function [T,P]=missnipals(X,w,T,Mode)
%[T]=missnipals(X,w,T,Mode)
%
% 'missnipals.m'
% $ Version 0.01 $ Date 10. July 1997 $ Not compiled $
%
% This algorithm requires the presence of:
% 'gsm.m' 'missmean.m'
%
% Copyright
% Claus A. Andersson 1995-
% E-mail claus@andersson.dk
%
% ----------------------------------------------------
%        Find eigenvectors according to NIPALS
% ----------------------------------------------------
%
% [T]=missnipals(X,w,T);
% [T]=missnipals(X,w);
%
% T is found so that X = T*P', s.t ||T||=1 and T'T=I
%
% X        : The matrix to be decomposed.
% w        : Number of factors to extract.
%            If w is high (perhaps>20) consider using SVD.
% T        : Initial guess of the solution, optional.
%            If T is not specified, a little time will
%            be used on finding orthogonal random 
%            starting values.
% Mode     : 'EM', uses an expectation maximization scheme.
%            This is often slower, but may be able to handle
%            systems with more missing values than 'WM'.
%            'WM', uses a weighing method. This approach may
%            be faster than 'EM' but it is also more 
%            vulnerable to the number and dispersion of
%            the missing values. Also, if w is close to, or greater,
%            than the rank of X, we must expect errorneous results.
%            The 'Mode' is only valid when there are missing data.
%           'EM' is the default choice
%
% It handles missing values NaNs (very dispersed, say, less than 15%)

% scalar ConvLim WarnLim ItMax a b i

ConvLim=1e-12;
WarnLim=1e-4;
ConvLimMiss=100*ConvLim;
ItMax=800;

P=[];
[a b]=size(X);

if (w>a | w>b) | w<1,
   T=NaN;
   disp('Error in MISSNIPALS.M: Number of factors to extract is invalid!')
   break;
end;

MissingExist=any(any(isnan(X)));

if ~exist('T') | isempty(T),
   T=orth(randn(a,w));
end;

if ~exist('Mode') | isempty(Mode),
   Mode='EM';
end;

if ~MissingExist, %no missing
   if all(size(T)==[a w]),
      if a>b,
         P=X'*T;
         l2=Inf;
         Z=X'*X;
         for i=1:w,
            p=P(:,i);
            d=1;
            it=0;
            while (d>ConvLim) & (it<ItMax),
               it=it+1;
               p=Z*p;
               l1=sqrt(sum(p.^2));
               p=p/l1;
               d=(l1-l2)^2;
               l2=l1;
            end;
            P(:,i)=sqrt(l1)*p;
            Z=Z-P(:,i)*P(:,i)';
            WarnLim=sqrt(l1)/1000;
            if it>=ItMax,
               disp('Warning from MISSNIPALS, High-X: Iterated up to the ItMax limit!')
            end;
            if d>WarnLim,
               disp('Warning from MISSNIPALS, High-X: The solution has not converged!')
            end;
         end;
         T=X*P;
      else
         l2=Inf;
         Z=X*X';
         for i=1:w,
            t=T(:,i); 
            d=1;
            it=0;
            while (d>ConvLim) & (it<ItMax),
               it=it+1;
               t=Z*t;
               l1=sqrt(sum(t.^2));
               t=t/l1;
               d=(l1-l2).^2;
               l2=l1;
            end;
            T(:,i)=sqrt(l1)*t;
            Z=Z-T(:,i)*T(:,i)';
            WarnLim=sqrt(l1)/1000;
            if it>=ItMax,
               disp('Warning from MISSNIPALS, Wide-X: Iterated up to the ItMax limit!')
            end;
            if d>WarnLim,
               disp('Warning from MISSNIPALS, Wide-X: The solution has not converged!')
            end;
         end;
      end;
   else
      disp('Error in MISSNIPALS.M (real): Inconsistency between input X and T');
      return;
   end;
elseif strcmp('WM',upper(Mode)) %Handles missing data by weighing
   if all(size(T) == [a w]),
      if a>b,
         P=missmult(X',T);
         l2=Inf;
         Z=missmult(X',X);
         for i=1:w,
            p=P(:,i);
            d=1;
            it=0;
            while (d>ConvLim) & (it<ItMax),
               it=it+1;
               p=Z*p;
               l1=sqrt(sum(p.^2));
               p=p/l1;
               d=(l1-l2)^2;
               l2=l1;
            end;
            P(:,i)=sqrt(l1)*p;
            Z=Z-P(:,i)*P(:,i)';
            WarnLim=sqrt(l1)/1000;
            if it>=ItMax,
               disp('Warning from MISSNIPALS (W), High-X: Iterated up to the ItMax limit!')
            end;
            if d>WarnLim,
               disp('Warning from MISSNIPALS (W), High-X: The solution has not converged!')
            end;
         end;
         T=missmult(X,P);
      else
         l2=Inf;
         Z=missmult(X,X');
         for i=1:w,
            t=T(:,i); 
            d=1;
            it=0;
            while (d>ConvLim) & (it<ItMax),
               it=it+1;
               t=Z*t;
               l1=sqrt(sum(t.^2));
               t=t/l1;
               d=(l1-l2).^2;
               l2=l1;
            end;
            T(:,i)=sqrt(l1)*t;
            Z=Z-T(:,i)*T(:,i)';
            WarnLim=sqrt(l1)/1000;
            if it>=ItMax,
               disp('Warning from MISSNIPALS (W), Wide-X: Iterated up to the ItMax limit!')
            end;
            if d>WarnLim,
               disp('Warning from MISSNIPALS (W), Wide-X: The solution has not converged!')
            end;
         end;
      end;
   else
      disp('Error in MISSNIPALS.M (W): Inconsistency between input X and T')
   end;
elseif strcmp('EM',upper(Mode)) %Handles missing data by EM
   if all(size(T) == [a w]),
      MissIdx=find(isnan(X));
      [i j]=find(isnan(X));
      mnx=missmean(X)/2;
      mny=missmean(X')/2;
      n=size(i,1);
      for k=1:n,
         i_=i(k);
         j_=j(k);
         X(i_,j_) = mny(i_) + mnx(j_);
      end;
      mnz=(missmean(mnx)+missmean(mny))/2;
      p=find(isnan(X));
      X(p)=mnz;
      
      T=missnipals(X,w,T);
      P=X'*T;
      Xm=T*P';
      X(MissIdx)=Xm(MissIdx);
      ssmisold=sum(sum( Xm(MissIdx).^2 ));
      sstotold=sum(sum( X.^2 ));
      ssrealold=sstotold-ssmisold;
      iterate=1;
      while iterate
         T=missnipals(X,w,T);
         P=X'*T;
         Xm=T*P';
         X(MissIdx)=Xm(MissIdx);
         ssmis=sum(sum( Xm(MissIdx).^2 ));
         sstot=sum(sum( X.^2 ));
         ssreal=sstot-ssmis;
         if ( (abs(ssreal-ssrealold)<ConvLim*ssrealold) & (abs(ssmis-ssmisold)<ConvLimMiss*ssmisold) ),
            iterate=0;
         end;
         ssrealold=ssreal;
         ssmisold=ssmis;
      end;
   else
      disp('Error in MISSNIPALS.M (E): Inconsistency between input X and T')
   end;
end;

T=gsm(T);
