function [Xfactors,Yfactors,B,ypred,ssx,ssy] = npls(X,y,DimX,DimY,Fac,show);

% $ Version 1.02 $ Date 28. July 1998 $ Not compiled $
% $ Version 1.03 $ Date 4. December 1998 $ Not compiled $ Cosmetic changes
%
% See also:
% 'parafac' 'tucker'
%
% Copyright, 1998 - 
% This M-file and the code in it belongs to the holder of the
% copyrights and is made public under the following constraints:
% It must not be changed or modified and code cannot be added.
% The file must be regarded as read-only. Furthermore, the
% code can not be made part of anything but the 'N-way Toolbox'.
% In case of doubt, contact the holder of the copyrights.
%
% Rasmus Bro
% Chemometrics Group, Food Technology
% Department of Food and Dairy Science
% Royal Veterinary and Agricultutal University
% Rolighedsvej 30, DK-1958 Frederiksberg, Denmark
% Phone  +45 35283296
% Fax    +45 35283245
% E-mail rb@kvl.dk
%
%
% MULTILINEAR PLS  -  N-PLS
%
% INPUT
% X        Array of independent variables
% y        Array of dependent variables
% DimX     Dimensions of X, 
%          e.g., if X is a 10 x 20 matrix, DimX = [10 20]
%          if X is a 10 x 20 x 5 array, DimX = [10 20 5]
% DimY     Dimensions of Y
% Fac      Number of factors to compute
% 
% OPTIONAL
% show	   If show = NaN, no outputs are given: npls(X,y,DimX,DimY,Fac,XCent,XScal,YCent,YScal,show);
%
%
% OUTPUT
% Xfactors Holds the parameters of the model of X in one vector.
%          Use fac2let to convert the parameters to scores and
%          weight matrices. I.e., for a three-way array do
%          [T,Wj,Wk]=fac2let(Xfactors,DimX,Fac);
%
% Yfactors Holds the parameters of the model of Y in one vector.
%          Use fac2let to convert the parameters to scores and
%          weight matrices. I.e., for a two-way array do
%          [U,Q]=fac2let(Yfactors,DimY,Fac);
% B        The regression coefficients from which the scores in
%          the Y-space are estimated from the scores in the X-
%          space (U = TB);
% ypred    The predicted values of Y
% ssx      Variation explained in the X-space.
%          ssx(f+1,1) is the sum-squared residual after first f factors.
%          ssx(f+1,2) is the percentage explained by first f factors.
% ssy      As above for the Y-space
%
%
% AUXILIARY
%
% If missing elements occur these must be represented by NaN.
%
%
% [Xfactors,Yfactors,B,ypred,ssx,ssy] = npls(X,y,DimX,DimY,Fac);
% or short
% [Xfactors,Yfactors,B] = npls(X,y,DimX,DimY,Fac);
%	Copyright
%	Rasmus Bro 1995
%	Denmark
%	E-mail rb@kvl.dk


 if nargin==0
   disp(' ')
   disp(' ')
   disp(' THE N-PLS REGRESSION MODEL')
   disp(' ')
   disp(' Type <<help npls>> for more info')
   disp('  ')
   disp(' [Xfactors,Yfactors,B,ypred,ssx,ssy] = npls(X,y,DimX,DimY,Fac);')
   disp(' or short')
   disp(' [Xfactors,Yfactors,B] = npls(X,y,DimX,DimY,Fac);')
   disp(' ')
   break
 elseif nargin<5
   error(' The inputs X, y, DimX, DimY, and Fac must be given')
 end



if ~exist('show')==1
  show=1;
end

maxit=20;

chkpfdim(X,DimX,NaN); % Check DimX
chkpfdim(y,DimY,NaN); % Check DimY 

 
[I,Jx]=size(X);
[I,Jy]=size(y);
if any(isnan(X(:)))
  Missing=1;
  if show~=0
    disp(' ')
    disp(' Don''t worry, missing values will be taken care of')
    disp(' ')
  end
  missX=abs(1-isnan(X));
  missy=abs(1-isnan(y));
else
  Missing=0;
end
[order]=length(DimX);
crit=1e-8;
B=zeros(Fac,Fac);
T=[];
U=[];
W=[];
Q=[];
if Missing
  SSX=sum(sum(X(find(missX)).^2));
  SSy=sum(sum(y(find(missy)).^2));
else
  SSX=sum(sum(X.^2));
  SSy=sum(sum(y.^2));
end
ssx=[];
ssy=[];
Xres=X;
yres=y;
xmodel=zeros(size(X));
q=zeros(Jy,1);
Q=[];
W=[];

for num_lv=1:Fac

  %init
  u=rand(DimX(1),1);t=rand(DimX(1),1);tgl=t+2;it=0;
    while (norm(t-tgl)/norm(t))>crit&it<maxit
    tgl=t;
    it=it+1;

      % w=X'u
      if Missing
         for i=1:Jx
           ww=Xres(find(missX(:,i)),i)'*u(find(missX(:,i)))/ ...
           (u(find(missX(:,i)))'*u(find(missX(:,i))));
           if length(ww)==0
             w(i)=0;
           else
             w(i)=ww;
          end
        end
      else
        w=Xres'*u;
      end
      wW=reshape(w,DimX(2),prod(DimX(3:length(DimX))));
%      w=npca(wW,DimX(2:length(DimX)));
%      [www,g]=tucker(wW,DimX(2:length(DimX)),ones(1,length(DimX)-1));
      if length(DimX)==3
        [w1,s,w2]=svd(wW);
        www=[w1(:,1);w2(:,1)];
      else
        www=parafac(wW,DimX(2:length(DimX)),1,[0 0 0 0 -1]');
     end
     w=fac2let(www,DimX(2:length(DimX)),1,1);

      % t=Xw
      if Missing
        for i=1:I,
         t(i)=Xres(i,find(missX(i,:)))*w(find(missX(i,:)))/(w(find(missX(i,:)))'*w(find(missX(i,:))));
        end
      else
        t=Xres*w;
      end

      % q=y't
      if Missing
        for i=1:Jy
          q(i)=yres(find(missy(:,i)),i)'*t(find(missy(:,i)))/(t(find(missy(:,i)))'*t(find(missy(:,i))));
        end
      else
        q=yres'*t;
      end
      if length(DimY)>2
        qQ=reshape(q,DimY(2),prod(DimY(3:length(DimY))));
        %        q=npca(qQ,DimY(2:length(DimY)));
        %        [qqq,g]=tucker(qQ,DimY(2:length(DimY)),ones(1,length(DimY)-1));
        if length(DimX)==3
          [q1,s,q2]=svd(qQ);
          qqq=[q1(:,1);q2(:,1)];
        else
          qqq=parafac(qQ,DimY(2:length(DimY)),ones(1,length(DimY)-1),1,[0 0 0 0 -1]');
        end
        qqq=qqq';
        q=fac2let(qqq,DimY(2:length(DimY)),1,1);
      else
        q=q/norm(q);
        qqq=q;
      end

      % u=yq
      if Missing
        for i=1:I
          u(i)=yres(i,find(missy(i,:)))*q(find(missy(i,:)))/(q(find(missy(i,:)))'*q(find(missy(i,:))));
        end
      else
        u=yres*q;
      end
    end
    T=[T;t];
%    W=[W;www];
    if length(DimX)>2&num_lv>1
      Z=[];
      for i=2:length(DimX);
        start=sum(DimX(2:i-1)*(num_lv-1))+1;
        endd=sum(DimX(2:i)*(num_lv-1));
        Z=[Z;W(start:endd)];
        start=sum(DimX(2:i-1))+1;
        endd=sum(DimX(2:i));
        Z=[Z;www(start:endd)];
      end
      W=Z;
    else
      W=[W;www];
    end
    U=[U;u];
    if length(DimY)>2&num_lv>1
      Z=[];
      for i=2:length(DimY);
        start=sum(DimY(2:i-1)*(num_lv-1))+1;
        endd=sum(DimY(2:i)*(num_lv-1));
        Z=[Z;Q(start:endd)];
        start=sum(DimY(2:i-1))+1;
        endd=sum(DimY(2:i));
        Z=[Z;qqq(start:endd)];
      end
      Q=Z;
    else
      Q=[Q;qqq];
    end

  tT=reshape(T,I,num_lv);
  uU=reshape(U,I,num_lv);
  B(1:num_lv,num_lv)=inv(tT'*tT)*tT'*uU(:,num_lv);
  if length(DimY)>2
    qQ=fac2let(Q,DimY(2:length(DimY)),num_lv,1);
  else
    qQ=reshape(Q,Jy,num_lv);
  end
  ypred=tT*B(1:num_lv,1:num_lv)*qQ';
  if Jy > 1
    if show~=0
      disp(' ') 
      fprintf('number of iterations: %g',it);
      disp(' ')
    end
  end
  xmodel=xmodel+tT(:,num_lv)*w';
  Xres=X-xmodel; 
  yres=y-ypred;
  if Missing
    ssx=[ssx;sum(sum(Xres(find(missX)).^2))];
    ssy=[ssy;sum(sum((y(find(missy))-ypred(find(missy))).^2))];
  else
    ssx=[ssx;sum(sum(Xres.^2))];
    ssy=[ssy;sum(sum((y-ypred).^2))];
  end
end

ssx= [ [SSX(1);ssx] [0;100*(1-ssx/SSX(1))]];
ssy= [ [SSy(1);ssy] [0;100*(1-ssy/SSy(1))]];

if show~=0
  disp('  ')
  disp('   Percent Variation Captured by N-PLS Model   ')
  disp('  ')
  disp('   LV      X-Block    Y-Block')
  disp('   ----    -------    -------')
  ssq = [(1:Fac)' ssx(2:Fac+1,2) ssy(2:Fac+1,2)];
  format = '   %3.0f     %6.2f     %6.2f';
  for i = 1:Fac
    tab = sprintf(format,ssq(i,:)); disp(tab)
  end
end

Xfactors=[T;W];
Yfactors=[U;Q];