function [U,S,V] = lyngby_cva(X, Y, arg1, arg2, arg3, arg4, arg5, ...
    arg6, arg7, arg8, arg9, arg10, arg11);

% lyngby_cva           - Canonical variate analysis (Canonical Ridge)
%
%       function [U,S,V] = lyngby_cva(X, Y, ...
%           'Propertyname', 'PropertyValue');
%
%	Input:	X   First datamatrix
%               Y   Second datamatrix
%               Property:
%                  'RidgeX'     { 0 } Canonical ridge parameter for X
%                  'RidgeY'     { 0 } Canonical ridge parameter for Y
%                  'Components' { 5 } Max number of canonical
%                               components returned
%                  'SubsetComp' { size(X,2) } Number of subset
%                               component, ie, SVD components piped to
%                               CVA 
%                  'NormRidgeX' [ {1} (true) | 0 ] Normalization of the
%                               ridge parameter with the trace of X's
%                               covariance: trace(X'X)/length(X'X)
%                  'NormRidgeY' [ {1} (true) | 0 ] Normalization of the
%                               ridge parameter with the trace of Y's
%                               covariance: trace(Y'Y)/length(Y'Y)
%                  'InitSVD'    [ {1} (true) | 0 ] Perform initial SVD
%                               on a rank-deficient X matrix. 
%
%       Output: U   Canonical sequence (canonical correlation vectors -
%                   sequence)
%               S   Canonical values (squared canonical correlation
%                   coefficient)
%               V   Canonical image (canonical correlation vectors -
%                   images)
%
%       Performs canonical variate analysis, or rather canonic
%       correlation analysis, on full ranked or rank-deficient
%       matrices. Ridge terms can be applied forming canonical ridge
%       analysis.
%       If 'InitSVD' is true a initial SVD is performed on the first
%       matrix (X) if the matrix is rank-deficient and then the ridge
%       parameter is only applied in the SVD subspace of X. 
%
%       Only canonic values larger than max(size(K)) * norm(K) * eps
%       will be returned.
%
%       Ref: Mardia, Multivariate Analysis
%
%       See also: lyngby_svd, lyngby_opls

% cvs : $Id: lyngby_cva.m,v 1.11 1998/08/15 14:46:06 fnielsen Exp $
%       $Revision: 1.11 $

    lyngby_log('CVA analysis');

    [Nxr, Nxc] = size(X);
    [Nyr, Nyc] = size(Y);
    
    % Default values
    ridgeX = 0;
    ridgeY = 0;
    canValues = min([Nxr Nxc Nyr Nyc 5]);
    subsetComp = Nxc;       % Subset components
    bNormRidgeX = 1;
    bNormRidgeY = 1;
    bInitSVD = 1;
    type = 1;     % ridge = 1, svd size subset = 2,
    

    % Parse Properties
    n = 1;
    while n <= nargin-2
      eval(sprintf('arg = lower(arg%d);', n));

      if strcmp(arg, 'type')
	n = n + 1;
	eval(sprintf('arg = lower(arg%d);', n));
	if strcmp(arg, 'ridge')
	  type = 1;
	elseif strcmp(arg, 'subset')
	  type = 2;
	else
	  error(['Argument to ''Type'' should a string ']);
	end
	
      elseif strcmp(arg, 'ridgex')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if length(arg) == 1
	    if arg <= 1 & arg >= 0 
	      ridgeX =  arg;
	    else
	      error(['Argument to ''RidgeX'' should ' ...
		    'be in the interval [0;1]']);
	    end
	  else
	    error(['Argument to ''RidgeX'' should ' ...
		  'be a single value']);
	  end
	else
	  error('Argument to ''RidgeX'' should be a value.'); 
	end
      elseif strcmp(arg, 'ridgey')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if length(arg) == 1
	    if arg <= 1 & arg >= 0 
	      ridgeY =  arg;
	    else
	      error(['Argument to ''RidgeY'' should ' ...
		    'be in the interval [0;1]']);
	    end
	  else
	    error(['Argument to ''RidgeY'' should ' ...
		  'be a single value']);
	  end
	else
	  error('Argument to ''RidgeY'' should be a value.'); 
	end

      elseif strcmp(arg, 'components')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if length(arg) == 1
	    if arg >= 1
	      canValues = 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

      elseif strcmp(arg, 'subsetcomp')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if length(arg) == 1
	    if arg >= 1
	      subsetComp = arg;
	    else
	      error(['Argument to ''SubsetComp'' should ' ...
		    'be a larger than 1']);
	    end
	  else
	    error(['Argument to ''SubsetComp'' should ' ...
		  'be a single value']);
	  end
	else
	  error('Argument to ''SubsetComp'' should be a value.'); 
	end

	
      elseif strcmp(arg, 'normridgex')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if length(arg) == 1
	    bNormRidgeX = (arg ~= 0);
	  else
	    error(['Argument to ''NormRidgeX'' should ' ...
		  'be a single value']);
	  end
	else
	  error('Argument to ''NormRidgeX'' should be a value.'); 
	end
      
      elseif strcmp(arg, 'initsvd')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if length(arg) == 1
	    bInitSVD = (arg ~= 0);
	  else
	    error(['Argument to ''InitSVD'' should ' ...
		  'be a single value']);
	  end
	else
	  error('Argument to ''InitSVD'' should be a value.'); 
	end
      
      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    % Largest dimension of canonical space
    canValues = min([Nxr Nxc Nyr Nyc canValues]);

    lyngby_log('CVA: centering of Y ...');
    Y = Y - ones(size(Y,1),1) * mean(Y);
    
    if Nxr >= Nxc | bInitSVD == 0
      % More scans than voxels

      lyngby_log('CVA: Building design covariance matrix ...');
      if ridgeY == 0
	S11 = Y'*Y; 
      elseif ridgeY == 1
	S11 = eye(Nyc);
      else
	if bNormRidgeY
	  S11 = Y'*Y;
	  normRidgeY = trace(S11) / length(S11);
	  S11 = (1-ridgeY) * S11 + normRidgeY * ridgeY * eye(Nyc);
	else
	  S11 = (1-ridgeY) * Y'*Y + ridgeY * eye(Nyc);
	end
      end

      lyngby_log('CVA: Inverting design covariance matrix ...');
      if ridgeY == 1
	nS11 = S11;
      else
	nS11 = mpower(S11, -0.5);
      end

      lyngby_log('CVA: Building data covariance matrix ...');
      S22 = X'*X;
      if ridgeX ~= 0
	if bNormRidgeX
	  normRidgeX = trace(S22) / length(S22);
	  S22 = (1-ridgeX) * S22 + normRidgeX * ridgeX * eye(Nxc);
	else
	  S22 = (1-ridgeX) * S22 + ridgeX * eye(Nxc);
	end
      end
      
      lyngby_log('CVA: Inverting data convariance matrix ...');
      [Ux,Sx] = eig(S22);
      [Sx,I] = sort(-diag(Sx));
      Sx = -Sx;

      % Zero all null-definite directions
      tol = length(Sx) * Sx(1) * eps;
      posdef = max(find(Sx>tol));
      Sx = Sx(1:posdef);
      Ux = Ux(:,I(1:posdef));

      % Invert and take square root
      iSx = 1 ./ sqrt(Sx);
      nS22 = Ux * diag(iSx) * Ux'; 

      lyngby_log('CVA: Building design/data covariance matrix ...');
      S12 = Y'*X ; 
      
      lyngby_log('CVA: Building normalized covariance matrix ...');
      K = nS11 * S12 * nS22;
      
      lyngby_log('CVA: Computing CVA SVD ...');
      [U, S, V] = lyngby_svd(K, 'ReduceComponents', canValues);

      lyngby_log('CVA: Canonical correlation vectors ...');
      U = nS11 * U;
      V = nS22 * V;
      
    else

      % Rank deficient: more voxels than scans
      % initial singular value decomposition before canonical ridge
      % analysis
      
      lyngby_log('CVA: Building design covariance matrix ...');
      if ridgeY == 0
	S11 = Y'*Y; 
      else
	if bNormRidgeY 
	  S11 = Y'*Y;
	  normRidgeY = trace(S11) / length(S11);
	  S11 = (1-ridgeY) * S11 + normRidgeY * ridgeY * eye(Nyc);
	else
	  S11 = (1-ridgeY) * Y'*Y + ridgeY * eye(Nyc);
	end
      end

      lyngby_log('CVA: Inverting design covariance matrix ...');
      if ridgeY == 1
	nS11 = S11;
      else	
	nS11 = mpower(S11, -0.5);
      end
      
      lyngby_log('CVA: Initial SVD on data matrix ...');
      S22 = X*X';

      lyngby_log('CVA: Computing eigenvalues ...');
      [Ux, Sx2] = eig(S22);

      [Sx2,I] = sort(-diag(Sx2));
      Sx2 = -Sx2;        % Eigenvalues of data matrix

      % Zero all null-definite directions
      tol = length(Sx2) * Sx2(1) * eps;
      posdef = max(find(Sx2 > tol));
      comp = min([ posdef subsetComp ]);
      
      Ux = Ux(:,I(1:comp));
      Sx2 = Sx2(1:comp);
      Sx = sqrt(Sx2);
      iSx = 1 ./ Sx;
      iSx2 = 1 ./ Sx2;
      
      lyngby_log('CVA: Computing projection ...');
      Vx = X' * Ux * diag(iSx);    

      if ridgeX == 0
	nS22 = diag(iSx2);
	S21 = Y' * Ux;
      else
	if bNormRidgeX
	  normRidgeX = sum(Sx2) / length(Sx);
	  nS22 = ( (1-ridgeX)*(Sx2) + normRidgeX * ridgeX ).^(-0.5); 
	else
	  nS22 = ( (1-ridgeX)*(Sx2) + ridgeX ).^(-0.5); 
	end
	S21 = Y' * Ux * diag(Sx .* nS22);
	nS22 = diag(nS22);
      end

      lyngby_log('CVA: Building design/data covariance matrix ...');
      K = nS11 * S21;
      
      lyngby_log('CVA: Computer CVA SVD ...');
      [U, S, V] = lyngby_svd(K, 'ReduceComponents', canValues);

      lyngby_log('CVA: Canonical correlation vectors ...');
      U = nS11 * U;
      V = nS22 * V;
      V = Vx * V;
      
    end

    lyngby_log('CVA completed!');

















