function [B, eta] = adaptica(X, k, varargin)
% ADAPTICA    Use the Bell and Sejnowski ICA algorithm to learn the weights of 
% the matrix B for ICA with Amari's dynamic learning rate.  
%
% B = adaptica(X, k, options)
%   X is the set of input vectors.  Each column of X is one sample.
%   k is the number of independent components to extract.
%   B is the seperating matrix, Y = BX are the independent signals
%
% options (specified by key/value pairs)
%   'rate' = 0.01       the learning rate (initial)
%   'B' = rand          the initial seperating matrix
%   'niter' = 1         the number of iterations to run
%   'alpha' = 0.002     Amari's alpha parameter for learning
%   'beta' = 20         Amari's beta parameter for learning
%   'delta' = 0.01      Amari's delta parameter for learning
%   'g' = tanh(t)       Bell and Senjowski's independence forcing function

%
% David Gleich
% CS 152 - Neural Networks
% 12 December 2003
%

% dimensionality assessment
[m n] = size(X);
if (k > m)
    error('Cannot extra more sources than sensors.');
end;

% subtract the mean
X = X - repmat(mean(X,2), 1, n);

% options
options = struct(...
    'rate', 0.01, ...
    'B', randn(k,m), ...
    'niter', 1, ...
    'alpha', 0.002, ...
    'beta', 20, ...
    'delta', 0.01, ...
    'g', inline('tanh(t)'));
options = getopt(options, varargin);

B = options.B;
u = options.rate;
g = options.g;
h = options.h;
d = options.delta;
a = options.alpha;
b = options.beta;

success = 0;
mdS = 0;
iter = 1;
eta = [u];

% set the initial r
S = pinv(B)*(eye-B);
r = S;
maxr = r;


% initialize the leaky averages


for niter=1:options.niter
    Bold = B;
    Sold = S;
    
    % update B, S
    for ii = 1:n
        y = B*X(:, ii);
        
        dB = (eye(k,k) - 2*g(y)*y')*B;
        B = B + u*dB;
        
        B = normr(B);
        
        if (any(~isfinite(B)))
            warning(sprintf('Lost convergence at iterator %i; lower learning rate?', iter));
            success = 11;
            break;
        end;
        
        % B = (I+S)^-1
        S = pinv(B)*(eye-B);
        dS = S - Sold;
       
        r = (1-d)*r + d*dS;
        if (norm(r, 'fro') > norm(maxr, 'fro'))
            maxr = r;
        end;
        
        u = u + a*u*(b*norm(r, 'fro')/norm(maxr, 'fro') - u);
        
        eta = [eta u];
       
        
        iter = iter + 1;
    end;
   
    
    if (success == 1)
        break;
    end;
end;