function [S, A] = lyngby_ica_bs_est(X, varargin)

% lyngby_ica_bs_est     - Bell & Sejnowski ICA estimation
%
%       function [S, A] = lyngby_ica_bs_est(X, varargin)
%
%       Input:  X    Matrix with data, objects x sensors
%
%       Output: S    Source matrix, objects x sources
%               A    Mixing matrix, sources x sensors
%
%       Independent component analysis as Bell and
%       Sejnowski. Optimization of the mixing matrix is done by
%       conjugate gradient optimization. The estimated mixing matrix
%       is square and is applied on the right side of the source
%       matrix: 
%                          X = S * A
%
%       This function is based on the icaML.m function in the 
%       ICA:DTU Toolbox by Thomas Kolenda.
%
%       Example:
%         Strue = randn(500, 2).^3;
%         Atrue = [ 3 4 ; 1 4 ];
%         X = Strue * Atrue;
%         [S, A] = lyngby_ica_bs_est(X);
%         figure, plot(X(:,1), X(:,2), '.', ...
%               5*[-A(1,1) A(1,1)], 5*[-A(1,2) A(1,2)], 'r-', ...
%               5*[-A(2,1) A(2,1)], 5*[-A(2,2) A(2,2)], 'g-')
% 
%       Ref: Kolenda T, Winther O, LK Hansen, ICA:DTU Toolbox,
%               http://isp.imm.dtu.dk/toolbox/ 
%
%       See also LYNGBY, LYNGBY_SVD.
%
% $Id: lyngby_ica_bs_est.m,v 1.1 2003/02/20 16:28:46 fnielsen Exp $


    % Check arguments
    if nargin < 1
      error('Too few input arguments')
    end
    [N, P] = size(X);
    
    % Default values
    maxIterations = 1000;
    sortType = 'energy';
    W = eye(P);
    ucminf_opt= [1  1e-4  1e-8  maxIterations];    

    % Scale X to avoid numerical problems
    scaleX = max(abs(max(X(:))), abs(min(X(:))));
    X = X / scaleX;

    
    % Estimate inverse mixing matrix
    % W = ucminf( 'ica_MLf' , par , W(:) , ucminf_opt );
    func = @cost;
    W = lyngby_opt_cg(func, W(:), 'parameters', { X });
    W = reshape(W, P, P); 

    % Estimates
    A = pinv(W);
    S = X * W;

    % Sort components according to energy
    if strcmp(sortType, 'energy')
      Aenergy = sum(A.^2,2) / P;
      Senergy = sum(S.^2) / N;
      energy = Aenergy' .* Senergy;
      [dummy, i] = sort(-energy);
      S = S(:,i);
      A = A(i,:);
    elseif strcmp(sortType, 'none')
      % Nothing
    else
      error('Internal error')
    end

    % Scale S back
    S = S .* scaleX;
    
    % variance one for sources
    stdS = std(S);
    A = A .* repmat(stdS', 1, size(A,1));
    S = S ./ repmat(stdS, size(S,1), 1);



function [f, df] = cost(W, X)
    [N, P] = size(X);
    W = reshape(W, P, P);
    S = X * W;

    % Negative log likelihood function
    f = -N*log(abs(det(W))) + sum(sum(log( cosh(S) ))) + N*P*log(pi);
  
    if nargout > 1
      df = -(N*inv(W') - X'*tanh(S));
      df = df(:);
    end














