%
%
%   flexica.m
%
%		seungjin@postech.ac.kr
%
%		Last updated:
%		February 6, 2003
%
%		function for batch version of flexible ICA
%
%		function [W]=flexica(x,n,nb,ga_W,maxits)
%
%		W: demixing matrix
%		x: multivariate input data matrix
%			m by N where m is the number of sensors and
%			N is the number of data points
%		n: number of sources (default=m)
%     nb: number of blocks (default=5)
%		ga_W: learning rate (default=.07)
%		maxits: number of maximal iterations (default=100)
%
%		You can find independent components by y=Wx.
%
%     Tips: 1. Better to choose the number of blocks such that
%              the number of data points in the block is about 3000.
%           2. Sometimes you might want to increase maxits.
%  
%


function [W]=flexica(x,n,nb,ga_W,maxits);

if nargin==0; 
   fprintf('function [W]=flexica(x,n,nb,ga_W,maxits) \n');
   break; 
end;
if nargin<5; maxits=100; end;
if nargin<4; ga_W=.07; end;
if nargin<3; nb=5; end;
if nargin<2; n=length(x(:,1)); end;


% figure out the size of data matrix x
m=length(x(:,1));
N=length(x(1,:));


% compute the block size
bsize=floor((2*N)/(nb+1));
if mod(bsize,2) ~= 0
   bsize=bsize-1;
end


% preprocessing: data sphering
for i=1:m
   x(i,:)=x(i,:)-mean(x(i,:));
end
Rx=cov(x');
[u d v]=svd(Rx);
Q=sqrt(inv(d(1:n,1:n)))*u(:,1:n)';
tempx=Q*x;
clear x;
x=tempx;
clear tempx;


% initial conditions
W=eye(n,n);
I=eye(n,n);
nu=2*ones(n,1);
loop=1;
dellik=1;
tol=0.0005;

% initial likelihood
yy=W*x;
templik=0;
for i=1:n
   templik=templik+log(nu(i))-log(2*gamma(1/nu(i)))-mean(abs(yy(i,:)).^nu(i));
end
lik=log(det(W))+templik;
oldlik=lik;
oldoldlik=lik; 


% main algorithm
% big loop (sweep)
while (dellik > tol) & (loop <= maxits),

   oldoldlik=oldlik;
   oldlik=lik;

   % small loop depending on the number of blocks
   for k=1:nb
    
      itsnum=(loop-1)*(nb-1)+k;
      y=W*x(:,(k-1)*bsize/2+1:(k+1)*bsize/2);

      % nonlinear function
      for i=1:n
         fy(i,:)=sign(y(i,:)).*abs(y(i,:)).^(nu(i)-1);
      end

      % update W
      scalingW=1/(n*bsize)*abs(trace(fy*y'));
      change=1/N*(y*y'+1/(1+ga_W*scalingW)*(fy*y'-y*fy'));
      W=W+ga_W*( I-change )*W;

   end

   % update nonlinear functions
   yy=W*x;
   % estimate moments
   m2=mean(yy.^2');
   m4=mean(yy.^4');
   kurt=m4./(m2.^2)-3*ones(1,n);
   % choose Gaussian exponent
   for i=1:n
      if kurt(i) > 20
	      nu(i)=.8;
	   elseif kurt(i) > 5
         nu(i)=1;
	   elseif kurt(i) > 0
         nu(i)=1.3;
	   else
	      nu(i)=4;
      end
   end
   % likelihood calculation
   templik=0;
   for i=1:n
      templik=templik+log(nu(i))-log(2*gamma(1/nu(i)))-mean(abs(yy(i,:)).^nu(i));
   end
   loop=loop+1;
   lik=log(det(W))+templik;
   dellik=(abs(lik-oldlik)+abs(oldlik-oldoldlik))/2;

   % anealing learning rate
   if dellik<.01
      ga_W=.05;
   elseif dellik<.005
      ga_W=.03
   else
      ga_W=.07;
   end

   % display
   fprintf('sweep= %d, difference of likelihood= %f\n', loop-1, dellik);

end % end of big loop


% demixing matrix
W=W*Q;