%pls_train Partial Least Squares (training)
%
%  [B,XRes,YRes,Options] = pls_train(X,Y)
%  [B,XRes,YRes,Options] = pls_train(X,Y,Options)
%
% INPUT
%  X   [N -by- d_X]  the training (input)  data matrix, N samples, d_X variables
%  Y   [N -by- d_Y]  the training (output) data matrix, N samples, d_Y variables
%
%   Options.
%   maxLV        maximal number of latent variables (will be corrected
%                if > rank(X)); 
%                maxLV=inf means maxLV=min(N,d_X) -- theoretical maximum
%                number of LV; 
%                by default =inf
%   method       'NIPALS' or 'SIMPLS'; by default ='SIMPLS'
%
%   X_centering  do nothing (=[] or 0), do mean centering (=nan), center
%                around some vaue v (=v); 
%                by default  =[]
%   Y_centering  do nothing (=[] or 0), do mean centering (=nan), center
%                around some vaue v (=v); 
%                by default  =[]
%   X_scaling    do nothing (=[] or 1), divide each col by std (=nan),
%                divide by some v (=v); 
%                by default =[]
%   Y_scaling    do nothing (=[] or 1), divide each col by std (=nan),
%                divide by some v (=v); 
%                by default =[]
%
% OUTPUT
%  B      [d_X -by- d_Y -by- nLV]  collection of regression matrices:
%         Y_new = X_new*B(:,:,n) represents regression on the first n
%         latent variables (X_new here after preprocessing, Y_new before
%         un-preprocessing)
%  XRes.
%   ssq  [1 -by- nLV]    the part of explaind sum of squares of
%        (preprocessed) X matrix
%   T    [N -by- nLV]    scores   (transformed (preprocessed) X)
%   R    [d_X -by- nLV]  weights  (transformation matrix)
%   P    [d_X -by- nLV]  loadings (back-transformation matrix)
%                        P(:,k) are the regression coef of
%                        (preprocessed) X on T(:,k)
%   W    [d_X -by- nLV]  local weights (local transformation matrix);
%                        ONLY FOR NIPALS
%   V    [d_X -by- nLV]  first k columns of this matrix are the
%                        orthornormal basis of the space spanned by k
%                        first columns of P; ONLY FOR SIMPLS
%  YRes.
%   ssq  [1 -by- nLV]    the part of explained sum of squares of
%                        (preprocessed) Y matrix
%   U    [N -by- nLV]    scores   (transformed (preprocessed) Y)
%   Q    [d_Y -by- nLV]  weights  (transformation matrix)
%   C    [d_Y -by- nLV]  C(:,k) are the regression coeff of
%                        (preprocessed) Y on T(:,k)
%   bin  [1 -by- nLV]    bin(k) is the regression coeff of U(:,k) on T(:,k)
%
% Options. contains the same fields as in input, but the values can be changed 
% (e.g. after mean centering X_centering = mean(X,1))
%
% DESCRIPTION
% Trains PLS (Partial Least Squares) regression model
%
% Relations between matrices (X end Y are assumed to be preprocessed):
% NIPALS:
% T = X*R (columns of T are orthogonal)
% R = W*prinv(P'*W)
% P = X'*T*prinv(T'*T)
% U = Y*Q - T*(C'*Q-tril(C'*Q))
% C = Y'*T*prinv(T'*T) = Q*diag(bin)
% bin = sqrt(diag(T'*Y*Y'*T))'*prinv(T'*T))
% B = R*C' = W*prinv(P'*W)*C'
%
% SIMPLS:
% T = X*R (columns of T are orthonormal)
% P = X'*T
% U = Y*Q
% C = Y'*T = Q*diag(bin)
% bin = sqrt(diag(T'*Y*Y'*T))'
% B = R*C'
%
% BOTH:
% T_new = X_new*R
% Y_new = X_new*B
%  
% SEE ALSO (<a href="http://37steps.com/prtools">PRTools Guide</a>)
% PLS_APPLY, PLS_TRANSFORM

% Copyright: S.Verzakov, s.verzakov@ewi.tudelft.nl 
% Faculty EWI, Delft University of Technology
% P.O. Box 5031, 2600 GA Delft, The Netherlands

% $Id: pls_train.m,v 1.2 2010/02/08 15:29:48 duin Exp $

function [B, XRes, YRes, Options] = pls_train(X,Y,Options)

[N_X, d_X] = size(X);
[N_Y, d_Y] = size(Y);

if N_X ~= N_Y
  error('size(X,1) must be equal to size(Y,1)');
else
  N = N_X;
end

if nargin < 3
  Options  = [];
end

DefaultOptions.X_centering = [];
DefaultOptions.Y_centering = [];
DefaultOptions.X_scaling = [];
DefaultOptions.Y_scaling = [];
DefaultOptions.maxLV = inf;
DefaultOptions.method = 'SIMPLS';

Options = pls_updstruct(DefaultOptions, Options);

if isinf(Options.maxLV)
  Options.maxLV = min(N,d_X);
elseif Options.maxLV > min(N,d_X)
  error('PLS: The number of LV(s) cannot be greater then min(N,d_X)');
end

[X, Options.X_centering, Options.X_scaling] = pls_prepro(X, Options.X_centering, Options.X_scaling);
[Y, Options.Y_centering, Options.Y_scaling] = pls_prepro(Y, Options.Y_centering, Options.Y_scaling);

ssq_X = sum(X(:).^2);
ssq_Y = sum(Y(:).^2);

B = zeros(d_X,d_Y,Options.maxLV);
XRes.ssq = zeros(1,Options.maxLV);
XRes.T   = zeros(N,Options.maxLV);
XRes.R   = zeros(d_X,Options.maxLV);
XRes.P   = zeros(d_X,Options.maxLV);
XRes.W   = [];
XRes.V   = [];

YRes.ssq = zeros(1,Options.maxLV);  
YRes.U   = zeros(N,Options.maxLV);  
YRes.Q   = zeros(d_Y,Options.maxLV);
YRes.C   = zeros(d_Y,Options.maxLV);
YRes.bin = zeros(1,Options.maxLV);  

ev = zeros(1,Options.maxLV);
nLV = Options.maxLV;
opts.disp   = 0;
opts.issym  = 1;
opts.isreal = 1;

switch upper(Options.method)
case 'NIPALS'
  XRes.W = zeros(d_X,Options.maxLV);
  for LV = 1:Options.maxLV
    S = X'*Y;
    if d_X <= d_Y
      if d_X > 1
        [w, ev(LV)] = eigs(S*S',1,'LA',opts);
      else
        w = 1;
        ev(LV) = S*S';
      end
      
      t = X*w;
  	  t2 = (t'*t);
      proj = t/t2;
  	  p = X'*proj;

      c = Y'*proj;
      bin = norm(c);  %bin = sqrt(ev)/t2:
      q = c/bin;
      u = Y*q;

    else
      if d_Y > 1
        [q, ev(LV)] = eigs(S'*S,1,'LA',opts);
      else
        q = 1;
        ev(LV) = S'*S;
      end
      u = Y*q;

      w = X'*u;
  	  w = w/norm(w);
  	  t = X*w;
   	  t2 = (t'*t);
      proj = t/t2;
  	  p = X'*proj;
      
      bin = u'*proj;  %bin = sqrt(ev)/t2:
      c = q*bin;
  	end

	  if LV == 1 
      if ev(LV) == 0
        error('PLS: Rank of the covariation matrix X''*Y is zero.');
      end  
	  elseif ev(LV) <= 1e-16*ev(1)
      nLV = LV-1;
      WarnMsg = sprintf(['\nPLS: Rank of the covariation matrix X''*Y is exausted ' ...
                         'after removing %d Latent Variable(s).\n'...
                         'Results only for the %d LV(s) will be returned'],nLV,nLV);
      if exist('prwarning')
        prwarning(1, WarnMsg);
      else
        warning(WarnMsg);
      end
      break;
	  end

    %R = W*prinv(P'*W)
    %the next portion of the code makes use of the fact taht P'*W is the upper triangle matrix:
    % 
    %if LV == 1
    %  InvPTW(1,1) = 1/(p'*w);
    %  r = w*InvPTW(1,1);
    %else
    %  InvPTW(1:LV,LV) = [(-InvPTW(1:LV-1,1:LV-1)*(XRes.P(:,1:LV-1)'*w)); 1]/(p'*w);
    %  r = [XRes.W(:,1:LV-1), w] * InvPTW(1:LV,LV);
    %end  
    % 
    %some optimization of the above code (notice that ones(m,0)*ones(0,n) == zeros(m,n)):
    %
    r = (w - (XRes.R(:,1:LV-1)*(XRes.P(:,1:LV-1)'*w)))/(p'*w);

    B(:,:,LV) = r*c';

  	XRes.ssq(:,LV) = t2*(p'*p)/ssq_X;
  	XRes.T(:,LV)   = t;
	  XRes.R(:,LV)   = r;
  	XRes.P(:,LV)   = p;
  	XRes.W(:,LV)   = w;
         
	  YRes.ssq(:,LV) = t2*(bin.^2)/ssq_Y;
  	YRes.U(:,LV)   = u;
  	YRes.Q(:,LV)   = q;
  	YRes.C(:,LV)   = c;
  	YRes.bin(:,LV) = bin;
	  
  	X = X - t*p';
  	Y = Y - t*c';
  end

  if nLV < Options.maxLV
    XRes.W = XRes.W(:,1:nLV);
  end

case 'SIMPLS'
  XRes.V = zeros(d_X,Options.maxLV);
  S = X'*Y;
  for LV = 1:Options.maxLV
	  if d_X <= d_Y
		  if d_X > 1
        [r, ev(LV)] = eigs(S*S',1,'LA',opts);
      else
        r = 1;
        ev(LV) = S*S';
      end
		  
      t = X*r;
		  norm_t = norm(t);
		  r = r/norm_t;
		  t = t/norm_t;
		  p = X'*t;

		  c = Y'*t;      %c = S'*r;
		  bin = norm(c); %bin = sqrt(ev)/norm_t:
		  q = c/bin;
		  u = Y*q;

    else
      if d_Y > 1
        [q, ev(LV)] = eigs(S'*S,1,'LA',opts);
      else
        q = 1;
        ev(LV) = S'*S;
      end
      r = S*q;
		  t = X*r;
		  norm_t = norm(t);
		  r = r/norm_t;
		  t = t/norm_t;
		  p = X'*t;

			u = Y*q;
		  bin = u'*t;    %bin = ev/norm_t:
		  c = q*bin;
	  end

	  if LV == 1
		  if ev(LV) == 0
        error('PLS: Rank of the covariation matrix X''*Y is zero.');
      else
        v = p;
      end  
	  elseif ev(LV) <= 1e-16*ev(1)
      nLV = LV-1;
      WarnMsg = sprintf(['\nPLS: Rank of the covariation matrix X''*Y is exausted ' ...
                         'after removing %d Latent Variable(s).\n'...
                         'Results only for the %d LV(s) will be returned'],nLV,nLV);
      if exist('prwarning')
        prwarning(1, WarnMsg);
      else
        warning(WarnMsg);
      end
      break;
    else
  	  v = p - XRes.V(:,1:LV-1)*(XRes.V(:,1:LV-1)'*p);
	  end

	  v = v/norm(v);

    B(:,:,LV) = r*c';	  

    XRes.ssq(:,LV) = (p'*p)/ssq_X;
	  XRes.T(:,LV)   = t;
	  XRes.R(:,LV)   = r;
	  XRes.P(:,LV)   = p;
	  XRes.V(:,LV)   = v;

	  YRes.ssq(:,LV) = (bin.^2)/ssq_Y;
	  YRes.U(:,LV)   = u;
	  YRes.Q(:,LV)   = q;
	  YRes.C(:,LV)   = c;
	  YRes.bin(:,LV) = bin;
	  
	  S = S - v*(v'*S);
  end

  if nLV < Options.maxLV
    XRes.V = XRes.V(:,1:nLV);
  end
end

if nLV < Options.maxLV
  B = B(:,:,1:nLV);

  XRes.ssq = XRes.ssq(:,1:nLV);
  XRes.T   = XRes.T(:,1:nLV);  
  XRes.R   = XRes.R(:,1:nLV);    
  XRes.P   = XRes.P(:,1:nLV);  

  YRes.ssq = YRes.ssq(:,1:nLV);
  YRes.U   = YRes.U(:,1:nLV);  
  YRes.Q   = YRes.Q(:,1:nLV);  
  YRes.C   = YRes.C(:,1:nLV);  
  YRes.bin = YRes.bin(:,1:nLV);
 
  Options.maxLV = nLV;
end

B = cumsum(B,3);

return;



function A = MakeSym(A)
A = 0.5*(A+A');
%A = max(A,A');

return