function [Z,A] = lyngby_ica_ms_est(X, varargin)

% lyngby_ica_ms_est    - Molgedey-Schuster ICA
%
%       function [Z,A] = lyngby_ica_ms_est(X, ...
%           'PropertyName', 'PropertyValue')
% 
%       Input:    X           A data matrix (N x P)
%
%       Property: Components  [ integer | {10} ] 
%                 Tau         [ integer | {5} ]
%                 Paradigm    Paradigm function
%
%       Output:   Z           Sources (N x K)
%                 A           Mixing matrix (K x P)
%
%       This function perform Molgedey Schuster independent component
%       analysis on a singular  matrix:
%
%                       Z * A = X
%
%       Note that X, Z, and A are transposed compared to the "normal"
%       notation found in the literature: X = AS.
%
%       Example:
%         % load a data set
%         lyngby_ui_global
%         [Z,A] = lyngby_ica_ms_est(XN);
%         lyngby_ui_view([X_MEAN ; X_STD; A], 'Mean|Std|1|2|3|4|5|6|7|8|9|10')
%
%       Ref: Mardia A6.2, 
%            Hansen, Larsen, Kolenda. 
%
% $Id: lyngby_ica_ms_est.m,v 1.3 2002/03/11 18:42:25 fnielsen Exp $


    % Default values
    Components = 10;
    paradigm   = []; 
    tau        = 5;


    % Parse properties
    n = 1;
    while n < nargin
      arg = lower(varargin{n});

      if strcmp(arg, 'components')
	n = n + 1;
	arg = varargin{n};
	if isreal(arg) & arg > 0 & prod(size(arg))
	  Components = arg;
	else
	  error(sprintf(['Argument with ''Components'' should be '...
		'integers']));
	end

      elseif strcmp(arg, 'tau')
	n = n + 1;
	arg = lower(varargin{n});
	if isnumeric(arg) & isreal(arg) & arg > 0
	  tau = arg;
	else
	  error(sprintf(['Argument with ''Tau'' should be '...
		'an integer']));
	end
	
      elseif strcmp(arg, 'paradigm')
	n = n + 1;
	arg = lower(varargin{n});
	if isnumeric(arg) & isreal(arg) & length(arg) == size(X,1)
	  paradigm = arg(:);
	else
	  error(sprintf(['Argument with ''Paradigm'' should be '...
		'a vector with length  integer']));
	end

	
      else
	error(sprintf('Invalid property: %s', arg));
      end
      n = n + 1;
    end

    
    
    % SVD
    [U,L,V] = lyngby_svd(X, 'reducecomponents', Components);
    
    % Shifted eigensequences
    U0 = U(1:end-tau,:);
    Ut = U(tau+1:end,:);
    
    M = U0' * Ut;
    % [Atilde, Gammatilde] = reig(0.5(M + M')); No right eigenvalue decomposition
    [Abreve, Gammabreve] = eig(0.5*(M + M')');
    Atilde = Abreve';   % Atilde: right eigenvalues
    
    A = Atilde * diag(L) * V';  
    Z = U * pinv(Atilde);

    % Sort the components according to power
    [sortL, indexL] = sort(L);
    A = A(indexL,:);
    Z = Z(:,indexL);
      
    if ~isempty(paradigm)
      iNeg = find( sum(Z .* repmat(paradigm, 1, size(Z,2))) < 0);
      A(iNeg,:) = - A(iNeg,:);
      Z(:,iNeg) = - Z(:,iNeg);
      
    end
    













