function W=kernel_ica(x,varargin)

% KERNEL_ICA - Runs the kernel ica algorithm on the given mixtures.
%              Returns a demixing matrix W such that s=W*x are the 
%              independent components
%
%	            It first centers and whitens the data and then minimizes
%              the given contrast function over orthogonal matrices.
%
%
%
% OPTIONAL PARAMETERS
%
%    Field name        Parameter                             Default
%
%    'contrast'        contrast function                     'kcca'
%                      values: 'kcca' 'kgv'
%
%    'opt'             Optimization procedure                'steepest'
%                      values:  'steepest' 'conjugate'
%
%    'polish'          1 if finishing with a half sigma      1
%                      value (better estimate)
%
%    'restarts'        number of restarts                    1
%    'sig'             kernel width                          1
%    'kap'             regularization parameter              0.01
%    'W0'              demixing matric initialization        rand_orth(m)
%    'disp'            verbose intermediate results          1







% default values

optimization='steepest';
contrast='kcca';
verbose=1;
[m,N]=size(x);
W0=rand_orth(m);

% currently fixed values but should depend on N and m
sigma=1;
kappa=.01;
eta=0.0001;
Nrestarts=1;
polish=1;


if (rem(length(varargin),2)==1)
   error('Optional parameters should always go by pairs');
else
   for i=1:2:(length(varargin)-1)
      switch varargin{i}
      case 'contrast'
         contrast= varargin{i+1};
      case 'sig'
         sigma=varargin{i+1};
      case 'kap'
         kappa=varargin{i+1};
      case 'W0'
         W0=varargin{i+1};
      case 'opt'
         optimization=varargin{i+1};
      case 'disp'
         verbose=varargin{i+1};
      case 'polish'
         polish=varargin{i+1};
      case 'restarts'
         Nrestarts=varargin{i+1};
         
      end
   end
end

% definition of parameters
kparam.sigmas=sigma*ones(1,m);
kparam.kappas=kappa*ones(1,m);
kparam.etas=kappa*1e-2*ones(1,m);
kparam.neigs=N*ones(1,m);
kparam.nchols=N*ones(1,m);


% first centers and scales data
if (verbose), fprintf('centering and scaling...'); end
xc=x-repmat(mean(x,2),1,N);  % centers data
covmat=xc*xc'/N;
sqcovmat=sqrtm(covmat);
invsqcovmat=inv(sqcovmat);
xc=invsqcovmat*xc;           % scales data
if (verbose), fprintf('done\n'); end

% making initial guess orthogonal
[U,S,V]=svd(W0*sqcovmat);
W0=U*V';

optparam.tolW=1e-2;
optparam.tolJ=1e-2;
optparam.maxit=20;
optparam.type=optimization;
optparam.verbose=verbose;

if (verbose), fprintf('Starting optimization, with %d restarts\n',Nrestarts); end
if (Nrestarts<=1)
   [J,W,details] = empder_search(contrast,xc,W0,kparam,optparam);
else
   optparam.Nrestarts=Nrestarts;
   optparam.Jaccept=0;
   [J,W,details]= global_mini(contrast,xc,W0,kparam,optparam);
   
end

% polishing
if (polish)
   if (verbose), fprintf('\nPolishing...\n'); end
   kparam.sigmas=kparam.sigmas/2;
   [J,W,details] = empder_search(contrast,xc,W,kparam,optparam);
end

W=W*invsqcovmat;