function Ep1 = fraa(E,Ep,L,iter)
%Array estimate -- estimate missing values
%
%Usage: fraa(E,Ep,L,iter)
%E=matrix with missing values
%Ep must be initial solution
%L: parameter (number of significant singular values + 1)
%iter number of iterations to perform

%Modified 8 September 2005   Laura Chihara
%
%%%%%%%%%% THIS IS  THE SET-UP  
%Get size of E
    [N,M]=size(E);
     
      if (L> M)
        error('L must be less than or equal to the number of columns of E')
      end;

  %get index of missing values (vectorized)
   missing=find(isnan(E));

   %Number of missing values
    m=length(missing);
    m2=m*m;

%%%%%%%%%%% NOW WE WORK WITH THE ALGORITHM

  Xp1=zeros(N,M);

   track=iter;
  
 while(iter > 0)
   A=Ep'*Ep;
   
   %Find singular value decomposition of A
   [U,S,V]=svd(A);
   
   %Singular values of Ep
   sigma2=S(S~=0);
   singular=sqrt(sigma2);
   partial_sig2=sum(sigma2(L:M));
   total_sig2=sum(sigma2(1:M));

      
   fprintf('\n iteration %3.0f \n', track-iter+1)
   fraction=partial_sig2/total_sig2;
   fprintf(' partial sum/total sum of sq. singular values \n %1.8f', fraction)
   fprintf('\n')
   
%Construct B=Bp
  
    B=sparse(m,m);  %pre-allocate space
      [is,js]=ind2sub([N,M],missing(1:m));
          for s=1:m
         for t=s:m

            if (is(s)==is(t))
            B(s,t)=sum(U(js(s),L:M)*U(js(t),L:M)');%
            B(t,s)=B(s,t);    %B is symmetric
           end  %end if
             
        end %end For t
      end %end for s

%%%NOW CONSTRUCT THE VECTOR Wp

W=sparse(m,1); %pre-allocate space

    for t=1:m
 
     K=sparse(N,M);
     K(missing(t))=1;
     
     W(t)=sum(diag(U(:,L:M)'*Ep'*K*U(:,L:M))); 
     clear K;
   end %end for

   xp1=-B\W;

  %Create matrix B_{p+1}
   Xp1(missing)=xp1;

  %Update solution
   Ep=Ep+Xp1;
  %set counter
   iter=iter-1;

  end %End while

   fprintf('\n')
   fprintf(' singular values (final iteration):\n')
    % fprintf('%16.6f',sigma)
    fprintf('%16.6f',singular)

   Ep1=Ep;