function [U,S,V] = lyngby_opls(X,Y, arg1, arg2, arg3, arg4);

% lyngby_opls          - Orthonormalized PLS
%
%       function [U,S,V] = lyngby_opls(X,Y, 'PropertyName',
%           'PropertyValue');
%
%	Input:	X   First datamatrix (not orthonormalized)
%               Y   Second datamatrix (Orthonormalized)
%
%       Property: Components  [ {5} | Integer ] Number of eigenvalues
%                             maintained  
%
%       Output: U   Canonical sequence (canonical vectors - sequence)
%               S   Canonical values (canonical values)
%               V   Canonical image (canonical vectors - images)
%
%       Orthonormalized partial least square. Y is normalized with its
%       dispersionmatrix, and X is not normalized, - it has its full
%       variance maintained. (A design matrix is suitable as the Y,
%       and the data matrix is suitable as X.)
%
%       The data matrix X is not centered. 
%
%       Only canonical values larger than max(size(C)) * norm(C) * eps 
%       will be returned.
%
%       See also LYNGBY, LYNGBY_SVD, LYNGBY_CVA, LYNGBY_SOP_MAIN.
%
% $Id: lyngby_opls.m,v 1.6 2002/12/17 17:47:35 fnielsen Exp $


    lyngby_log('OPLS analysis');

    [Nxr, Nxc] = size(X);
    [Nyr, Nyc] = size(Y);
    
    % Default values
    samValues = min([Nxr Nxc Nyr Nyc 5]);
    
    % Parse Properties
    n = 1;
    while n <= nargin-2
      eval(sprintf('arg = lower(arg%d);', n));
      if strcmp(arg, 'components')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if length(arg) == 1
	    if arg >= 1
	      samValues = arg;
	    else
	      error(['Argument to ''Components'' should ' ...
		    'be a larger than 1']);
	    end
	  else
	    error(['Argument to ''Components'' should ' ...
		  'be a single value']);
	  end
	else
	  error('Argument to ''Components'' should be a value.'); 
	end
      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    samValues = min([Nxr Nxc Nyr Nyc samValues]);

    
    if Nxr >= Nxc
      % More scans than voxels
      
      % lyngby_log('OPLS: building data dispersionmatrix ...');
      % S22 = X'*X; 

      lyngby_log('OPLS: building design dispersionmatrix ...');
      S11 = Y'*Y; 

      lyngby_log('OPLS: building covariance matrix ...');
      S12 = Y'*X ; 
      
      % lyngby_log('Building OPLS normalization ...');
      C = S11^(-0.5) * S12;

      lyngby_log('Beginning OPLS svd ...');
      [U, S, V] = lyngby_svd(C, 'ReduceComponents', samValues);
      
      lyngby_log('OPLS completed!');

    else

      lyngby_log('OPLS: Computing first covariances ... 2-8');
      S22 = X*X';

      lyngby_log('OPLS: Computing eigenvalues ... 3-8');
      [Ux, Sx] = eig(S22);

      [Sx,I] = sort(-diag(Sx));
      Sx = -Sx;
      Ux = Ux(:,I); 

      % Zero all null-definite directions
      tol = length(Sx) * Sx(1) * eps;
      posdef = max(find(Sx>tol));
      if posdef < length(Sx)
	zeroIndices = (posdef+1):length(Sx);
	% Ux(:,zeroIndices) = 0;
	% Sx(zeroIndices) = 0;
	Ux = Ux(:,1:posdef);
	Sx = Sx(1:posdef);
      end
      

      % Diagonal pseudo inverse
      Sx = sqrt(Sx);
      iSx = 1 ./ Sx(find(Sx));
      
      lyngby_log('OPLS: Computing projection ... 4-8');
      Vx = X' * Ux * diag(iSx);    
    
      S11 = Y' * Y;
      S21 = Y' * Ux * diag(Sx);

      % lyngby_log('Building OPLS normalization ... 5-8');
      % C = S11^(-0.5) * S21;
      C = real(pinv(S11)^(0.5)) * S21;

      
      lyngby_log('Beginning OPLS svd ... 6-8');
      [U, S, V] = lyngby_svd(C, 'ReduceComponents', samValues);

      lyngby_log('Computing eigenimages ... 7-8');

      V = Vx * V;
      lyngby_log('OPLS completed! 8-8');
      
    end


















