function [S,A,U,ll,info]=icaML(X,K,par,debug_draw)

% icaML     : ICA by ML (Infomax) with square mixing matrix and no noise.

%

% function [S,A,U,ll,info]=icaML(X,[K],[par]) Independent component analysis (ICA) using

%                                       maximum likelihood, square mixing matrix and 

%                                       no noise [1] (Infomax). Source prior is assumed 

%                                       to be p(s)=1/pi*exp(-ln(cosh(s))). For optimization

%                                       the BFGS algorithm is used [2]. See code for 

%                                       references.

%                                       

%                                       X  : Mixed signals

%                                       K  : Number of source components.

%                                           For K=0 (default) number of sources are equal to number of

%                                           observations.

%                                           For K < number of observations, SVD is used to reduce the

%                                           dimension.

%                                       par:  Vector with 4 elements:

%                                           (1)  :  Expected length of initial step

%                                           Stopping criteria:

%                                           (2)  :  Gradient  ||g||_inf <= par(2) 

%                                           (3)  :  Parameter changes ||dW||_2  <= par(3)*(par(3) + ||W||_2)

%                                           (4)  :  Maximum number of iterations 

%                                           Any illegal element in  opts  is replaced by its

%                                           default value,  [1  1e-4*||g(x0)||_inf  1e-8  100]

%                                       debug_draw : Draw debug information

%                                       

%                                       S  : Estimated source signals with variance

%                                            scaled to one.

%                                       A  : Estimated mixing matrix

%                                       U  : Principal directions of preprocessing PCA. %                                            If K (the number of sources) is equal to the number%                                            of observations then no PCA is performed and U=eye(K).  %                                       ll : Log likelihood for estimated sources

%                                       info :  Performance information, vector with 6 elements:

%                                           (1:3)  : final values of [ll  ||g||_inf  ||dx||_2] 

%                                           (4:5)  : no. of iteration steps and evaluations of (ll,g)

%                                           (6)    : 1 means stopped by small gradient

%                                                    2 means stopped by small x-step

%                                                    3 means stopped by max number of iterations.

%                                                    4 means stopped by zero step.

%

%

%                                       Eks. Separate with number of sources equal number of

%                                       observations.

%                                       

%                                               [S,A] = icaML(X);

%                                       

%                                       Eks. Separate with number of sources equal to k, using SVD

%                                       as pre-processing.

%                                       

%                                               [S,A,U] = icaML(X,k);%                                       % - version 1.5 (Revised 9/9-2003)% - IMM, Technical University of Denmark% - version 1.4

% - by Thomas Kolenda 2002 - IMM, Technical University of Denmark

% Revised: 9/9-2003, Mads, mad@imm.dtu.dk%   * Fixed help message to inform the user about U.%   * Removed the automatic use of PCA in the quadratic case.% Bibtex references:

% [1]  

%   @article{Bell95,

%       author =       "A. Bell and T.J. Sejnowski",

%       title =        "An Information-Maximization Approach to Blind Separation and Blind Deconvolution",

%       journal =      "Neural Computation",

%       year =         "1995",

%       volume =       "7",

%       pages =        "1129-1159",

%   }

%

% [2]

%   @techreport{Nielsen01:unopt,

%       author =        "H.B. Nielsen",

%       title =         "UCMINF - an Algorithm for Unconstrained, Nonlinear Optimization ",

%       institution =   "IMM, Technical University of Denmark",

%       number =        "IMM-TEC-0019",

%       year =          "2001",

%       url =           "http://www.imm.dtu.dk/pubdb/views/edoc_download.php/642/ps/imm642.ps",

%   }

% algorithm parameter settings

MaxNrIte = 1000;



try

  debug = debug_draw;

catch

  debug = 0;    

end



if debug==1,fprintf('\n** Start icaML ***************************************\n');tic;end



% Scale X to avoid numerical problems

Xorg = X;

scaleX = max([abs(max(X(:))),abs(min(X(:)))]);

X = X./scaleX;



% set number of source parametersif nargin<2, K=size(X,1); endif ((K>0) & (K<size(X,1))),  [U,X] = callSVD(X,K,debug);else  U=eye(size(X,1));  if debug==1,disp('Don''t use SVD');endend% initialize optimize parameters

try 

  ucminf_opt(1) = par(1);

  ucminf_opt(2) = par(2);

  ucminf_opt(3) = par(3);

  ucminf_opt(4) = par(4);

catch

  if debug==1,disp('Use default parameters to optimize');end

  ucminf_opt = [1  1e-4  1e-8  MaxNrIte];

end



% initialize variables

[M,N] = size(X);

W = eye(M);



if debug==1,fprintf('Number of samples %d - Number of sources %d\n',N,M);end



% optimize

if debug==1,fprintf('Optimize ICA ... ');end

par.X=X; par.M=M; par.N=N;

[W,info] = ucminf( 'ica_MLf' , par , W(:) , ucminf_opt );

W = reshape(W,M,M); 



% estimates

A=pinv(W);

S=W*X;



if debug==1,disp('done optimize ICA!');end



% sort components according to energy

Avar=diag(A'*A)/M;

Svar=diag(S*S')/N;

vS=var(S');

sig=Avar.*Svar;

[a,indx]=sort(sig);

S=S(indx(M:-1:1),:);

A=A(:,indx(M:-1:1));



% scale back

A=A.*scaleX;



% log likelihood

if nargout>3

  ll= N*log(abs(det(inv(A)))) - sum(sum(log( cosh(S) ))) - N*M*log(pi);

end



if debug==1,fprintf('** End of icaML ****** time %0.2f sec******************\n\n',toc/60);end







function [f,dW]=ica_MLf(W,par)

% returns the negative log likelihood and its gradient w.r.t. W



X=par.X; M=par.M; N=par.N;

W=reshape(W,M,M);



S=W*X;



% negative log likelihood function

f=-( N*log(abs(det(W))) - sum(sum(log( cosh(S) ))) - N*M*log(pi) );



if nargout>1

  % gradient w.r.t. W

  dW=-(N*inv(W') - tanh(S)*X');

  dW=dW(:);

end







function [U,DV] = callSVD(X,K,draw)

% Reduce dimension with SVD



[M,N] = size(X);



if N>M % Transpose if matrix is flat, to speed up the svd and later transpose back again

    if draw==1,disp('Do Transpose SVD');end

    [V,D,U] = svd(X',0);

else;  

    if draw==1,disp('Do SVD');end

    [U,D,V] = svd(X,0);

end;

  

DV = D(1:K,1:K)*V(:,1:K)';







function  [X, info, perf, D] = ucminf(fun,par, x0, opts, D0)

%UCMINF  BFGS method for unconstrained nonlinear optimization:

% Find  xm = argmin{f(x)} , where  x  is an n-vector and the scalar

% function  F  with gradient  g  (with elements  g(i) = DF/Dx_i )

% must be given by a MATLAB function with with declaration

%           function  [F, g] = fun(x, par)

% par  holds parameters of the function.  It may be dummy.