function [U, S, V] = lyngby_sopressvd(X, R, arg1, arg2, arg3, arg4, ...
    arg5, arg6, arg7, arg8, arg9, arg10, arg11, arg12, arg13, arg14, ...
    arg15, arg16)

% lyngby_sopressvd     - SOP with SVD on the residual
%
%	function [U,S,V] = lyngby_sopresvd(X, R,
%	    'PropertyName', 'PropertyValue')
%
%       Input:  X   The datamatrix
%               R   The run structure
%               PropertyName:
%                  'Paradigm'         The paradigm 
%                  'CanComponents'    { 5 } Number of canonical values
%                                     maintained
%                  'PrincComponents'  { (all) } Number of principal
%                                     components maintained 
%
%       Output: U   Canonical/Eigen sequences
%               S   Canonical/Eigen values
%               V   Canonical/Eigen images
%
%       Perform Strother Orthonormalized Partial Least Square (SOP) on
%       the data matrix X by using the run structure (group structure)
%       R. The analysis will only the return number of components
%       specified by 'CanComponents'.
%       The subspace spanned by the result of this analysis is then
%       returned as the first part of the result [U,S,V]. 
%       The second part of the result is found by an SVD of the
%       residual of the data matrix, ie, an SVD of the space
%       orthogonal to SOP space. 
%
%       See also: lyngby_sop_main, lyngby_opls, lyngby_cva, lyngby_svd

% cvs : $Id: lyngby_sopressvd.m,v 1.3 1999/07/09 14:03:31 fnielsen Exp $
%       $Revision: 1.3 $

    lyngby_global

    if nargin < 2 
      error('Wrong number of arguments');
    end
    
    if size(X,1) ~= size(R, 1)
      error(sprintf(['X and R do not have the same length. ',...
	    'X: %dx%d, R: %dx%d'], size(X), size(R)));
    end
    
    % Strother's design matrix
    Y = lyngby_sdesign(R);

    if (size(Y,2) < 2)
      [Runs, ScansPerRun] = lyngby_runinfo(R);
      error(sprintf(['The SOP model will not give any ', ...
	    'meaningfull results if there is only one run.\n', ...
	    'There are %d run with %d scans in.'], Runs, ...
	  ScansPerRun));
    end
    
    % Default parameters
    P = [];
    canComp = 5;
    princComp = size(X,1);
    
    % Parse Properties
    n = 1;
    while n <= nargin-2
      eval(sprintf('arg = lower(arg%d);', n));
      if strcmp(arg, 'paradigm')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if all(size(R) == size(arg))
	    P = arg;
	  else
	    error(sprintf(['The argument to ''Paradigm'' should ' ...
		  'be of the same size as R. Paradigm: %dx%d, ' ...
		  'R: %dx%d'], size(arg), size(R)));
	  end
	else
	  error('The argument with ''Paradigm'' should be a vector.'); 
	end
      elseif strcmp(arg, 'princcomponents')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if length(arg) == 1
	    if arg >= 1
	      princComp = arg;
	    else
	      error(['Argument to ''PrincComponents'' should ' ...
		    'be a larger than 1']);
	    end
	  else
	    error(['Argument to ''PrincComponents'' should ' ...
		  'be a single value']);
	  end
	else
	  error(['Argument to ''PrincComponents'' should ' ,...
		'be a value.']); 
	end
      elseif strcmp(arg, 'cancomponents')
	n = n + 1;
	eval(sprintf('arg = arg%d;', n));
	if isreal(arg)
	  if length(arg) == 1
	    if arg >= 1
	      canComp = arg;
	    else
	      error(['Argument to ''CanComponents'' should ' ...
		    'be a larger than 1']);
	    end
	  else
	    error(['Argument to ''CanComponents'' 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

    % OPLS analsysis
    [Usop,Ssop,Vsop] = lyngby_opls(X, Y, 'Components', canComp);
    
    % Projection to a subspace
    Ysop = X * Vsop;
    
    % SVD the residual from the OPLS analysis
    [Usvd, Ssvd, Vsvd] = lyngby_svd(X - Ysop * Vsop', ...
	'ReduceComponents', princComp);
    
    V = [ Vsop Vsvd ];
    lyngby_log('Projecting with the combined basis')
    U = X * V;
    S = [ Ssop ; Ssvd ];
    
    % Flip the sign of the first canonical sequence and image if we are
    % out of phase with the paradigm
    if ~isempty(P)
      if P'*Y*U(:,1) < 0
	U(:,1) = -U(:,1);
	V(:,1) = -V(:,1);
      end
    end



















