%======================================================================
%
% This function performs ICA on input data
%
% the control display is
%       it  number of epoche
%       L   learning rate
%       wchange         change in w-w_old
%       angle change    of delta
%       KD  Kullback-Distance of w, how far is w away from orthogonality
%           after sphering (=0 for theoretically perfect separation)
%======================================================================

%------- important paramters to play with
B=20;            %--blocksize
L=0.001;          %--learning rate
LF=0.985;         %--learning factor
alpha=0.1;        %--momentum constant
iter=25;         %--number of iterations
step=5;           %--display after steps
ksize=1000;       %--blocksize for kurtosis estimation

% load all the training files
%
load hsbe1.dat;
load msbe1.dat;
load lsbe1.dat;

hsbe = hsbe1';
lsbe = lsbe1';
msbe = msbe1';


% load a file with all images
%
% load train1.dat;
% train = train1';

% load the test data
%
load test1.dat;
test = test1' ;


% perform required process for hsbe training
%
%------ permutate
[N,P]=size(hsbe);                    
permute=randperm(P);                
hsbe=hsbe(:,permute);                   
wh=1;
if (wh==1),
%------ shering
fprintf('\n sphering ... ');
hsbe=hsbe-mean(hsbe')'*ones(1,P);               % subtract means from mitraines
wz=2*inv(sqrtm(cov(hsbe')));              % get decorrelating matritrain
hsbe=wz*hsbe;                               % decorrelate mitraines
fprintf(' done \n');
end;

%------- initialize ----
w=eye(N);  wold=w;                  
dw_old=w;
olddelta=zeros(1,N*N);
oldchange=1;
degconst = 180./pi;

K=diag(ones(1,N));
if (P>ksize),
  ksize=ksize;
else
  kzise=P;
end;

Id=eye(N);
noblocks=fix(P/B);
BI=B*Id;

%------- start learning ----------------------------------------
ii=1;
sweep=1

for it=1:iter,          %--- begin iterations

t=1;

for t=t:B:t-1+noblocks*B,   %--- begin epoche
    u=w*hsbe(:,t:t+B-1);

    %------ Etraintention to sub and supergaussians
    tmp1=mean((sech(u').^2)).*mean(u'.^2);
    tmp2=mean(tanh(u').*u');
    K=diag(sign(tmp1-tmp2));

    %------
    dw=L*(BI-K*tanh(u)*u'-u*u')*w;
 

% Logistic ICA weight update

%   y = 1./(1+exp(-u));                                                
%   dw = L*(BI+(1-2*y)*u')*w;                  

   d_w=dw+alpha*dw_old;
    dw_old=d_w;
    w=w+d_w;
    
    sweep=sweep+1; 

    %------------------------------------CONTROL
    if (rem(sweep,step)==0),
    
      if (abs(w(1,1))>10000),
        L=L*.9;
        fprintf('\n weights blow up ... lowering learning rate %5f \n',L);
        w=0.1*eye(N);  wold=w;        
        dw_old=w; 
        olddelta=zeros(1,N*N); oldchange=1;
      end;
      
      d=w*w';

      tt=1;
      for k=1:N                                 %how far is w*w' from I
        tt=tt*d(k,k);
      end
      KD=log((tt)/det(d));                      %Kullback distance measure

      fprintf(1, '%f \n', KD); 

      delta=reshape(wold-w,1,N*N);
      change=delta*delta';
      angledelta=acos((delta*olddelta')/sqrt(change*oldchange));
      olddelta=delta;
      oldchange=change;
      
      fprintf('Epoche %d - lrate %5f, wchange %6.5f, angledelta %4.1f, deg K-distance %4.3f \n',it,L,change,degconst*angledelta,KD);
                
    end              
    %------------------------------------ end of control display

    wold=w;
    
  end;                
  %-------------------------------------- end of epoche

  L=LF*L

end;                  
%---------------------------------------- end of iterations

wi_hsbe = w;
[R, C] = size(test);   
u_hsbe=wi_hsbe*hsbe;
test=test-mean(test')'*ones(1, C);
u_test_hsbe = wi_hsbe*wz*test


% perform required process for msbe training
%
%------ permutate
[N,P]=size(msbe);                    
permute=randperm(P);                
msbe=msbe(:,permute);                   

wh=1;
if (wh==1),
%------ shering
fprintf('\n sphering ... ');
msbe=msbe-mean(msbe')'*ones(1,P);               % subtract means from mitraines
wz=2*inv(sqrtm(cov(msbe')));              % get decorrelating matritrain
msbe=wz*msbe;                               % decorrelate mitraines
fprintf(' done \n');
end;

%------- initialize ----
w=eye(N);  wold=w;                  
dw_old=w;
olddelta=zeros(1,N*N);
oldchange=1;
degconst = 180./pi;

K=diag(ones(1,N));
if (P>ksize),
  ksize=ksize;
else
  kzise=P;
end;

Id=eye(N);
noblocks=fix(P/B);
BI=B*Id;

%------- start learning ----------------------------------------
ii=1;
sweep=1

for it=1:iter,          %--- begin iterations

t=1;

for t=t:B:t-1+noblocks*B,   %--- begin epoche
    u=w*msbe(:,t:t+B-1);

    %------ Etraintention to sub and supergaussians
    tmp1=mean((sech(u').^2)).*mean(u'.^2);
    tmp2=mean(tanh(u').*u');
    K=diag(sign(tmp1-tmp2));

    %------
    dw=L*(BI-K*tanh(u)*u'-u*u')*w;
 

% Logistic ICA weight update

%   y = 1./(1+exp(-u));                                                
%   dw = L*(BI+(1-2*y)*u')*w;                  

   d_w=dw+alpha*dw_old;
    dw_old=d_w;
    w=w+d_w;
    
    sweep=sweep+1; 

    %------------------------------------CONTROL
    if (rem(sweep,step)==0),
    
      if (abs(w(1,1))>10000),
        L=L*.9;
        fprintf('\n weights blow up ... lowering learning rate %5f \n',L);
        w=0.1*eye(N);  wold=w;        
        dw_old=w; 
        olddelta=zeros(1,N*N); oldchange=1;
      end;
      
      d=w*w';

      tt=1;
      for k=1:N                                 %how far is w*w' from I
        tt=tt*d(k,k);
      end
      KD=log((tt)/det(d));                      %Kullback distance measure

      fprintf(1, '%f \n', KD); 

      delta=reshape(wold-w,1,N*N);
      change=delta*delta';
      angledelta=acos((delta*olddelta')/sqrt(change*oldchange));
      olddelta=delta;
      oldchange=change;
      
      fprintf('Epoche %d - lrate %5f, wchange %6.5f, angledelta %4.1f, deg K-distance %4.3f \n',it,L,change,degconst*angledelta,KD);
                
    end              
    %------------------------------------ end of control display

    wold=w;
    
  end;                
  %-------------------------------------- end of epoche

  L=LF*L

end;                  
%---------------------------------------- end of iterations

wi_msbe = w;
u_msbe=wi_msbe*msbe;
u_test_msbe = wi_msbe*wz*test;

% perform required process for lsbe training
%
%------ permutate
[N,P]=size(lsbe);                    
permute=randperm(P);                
lsbe=lsbe(:,permute);                   

wh=1;
if (wh==1),
%------ shering
fprintf('\n sphering ... ');
lsbe=lsbe-mean(lsbe')'*ones(1,P);               % subtract means from mitraines
wz=2*inv(sqrtm(cov(lsbe')));              % get decorrelating matritrain
lsbe=wz*lsbe;                               % decorrelate mitraines
fprintf(' done \n');
end;

%------- initialize ----
w=eye(N);  wold=w;                  
dw_old=w;
olddelta=zeros(1,N*N);
oldchange=1;
degconst = 180./pi;

K=diag(ones(1,N));
if (P>ksize),
  ksize=ksize;
else
  kzise=P;
end;

Id=eye(N);
noblocks=fix(P/B);
BI=B*Id;

%------- start learning ----------------------------------------
ii=1;
sweep=1

for it=1:iter,          %--- begin iterations

t=1;

for t=t:B:t-1+noblocks*B,   %--- begin epoche
    u=w*lsbe(:,t:t+B-1);

    %------ Etraintention to sub and supergaussians
    tmp1=mean((sech(u').^2)).*mean(u'.^2);
    tmp2=mean(tanh(u').*u');
    K=diag(sign(tmp1-tmp2));

    %------
    dw=L*(BI-K*tanh(u)*u'-u*u')*w;
 

% Logistic ICA weight update

%   y = 1./(1+exp(-u));                                                
%   dw = L*(BI+(1-2*y)*u')*w;                  

   d_w=dw+alpha*dw_old;
    dw_old=d_w;
    w=w+d_w;
    
    sweep=sweep+1; 

    %------------------------------------CONTROL
    if (rem(sweep,step)==0),
    
      if (abs(w(1,1))>10000),
        L=L*.9;
        fprintf('\n weights blow up ... lowering learning rate %5f \n',L);
        w=0.1*eye(N);  wold=w;        
        dw_old=w; 
        olddelta=zeros(1,N*N); oldchange=1;
      end;
      
      d=w*w';

      tt=1;
      for k=1:N                                 %how far is w*w' from I
        tt=tt*d(k,k);
      end
      KD=log((tt)/det(d));                      %Kullback distance measure

      fprintf(1, '%f \n', KD); 

      delta=reshape(wold-w,1,N*N);
      change=delta*delta';
      angledelta=acos((delta*olddelta')/sqrt(change*oldchange));
      olddelta=delta;
      oldchange=change;
      
      fprintf('Epoche %d - lrate %5f, wchange %6.5f, angledelta %4.1f, deg K-distance %4.3f \n',it,L,change,degconst*angledelta,KD);
                
    end              
    %------------------------------------ end of control display

    wold=w;
    
  end;                
  %-------------------------------------- end of epoche

  L=LF*L

end;                  
%---------------------------------------- end of iterations

wi_lsbe = w;
u_lsbe=wi_lsbe*lsbe;
u_test_lsbe = wi_lsbe*wz*test;


% compute the mean and covariance matrices
% for all the transformed data
%
mean_hsbe = mean(u_hsbe');
mean_msbe = mean(u_msbe');
mean_lsbe = mean(u_lsbe');

% the inverse matrices of covariance matrices
%
cov_hsbe = inv(cov(u_hsbe'));
cov_msbe = inv(cov(u_msbe'));
cov_lsbe = inv(cov(u_lsbe'));

% count how many images have to be tested
%
[no_images, features] = size(test');

% setup file for storing euclidean distance
%
FID = fopen('ica_output.dat', 'w');

% loop through each image and compute euclidean 
% distance
%

u_test1 = u_test_hsbe';
u_test2 = u_test_msbe';
u_test3 = u_test_lsbe';

for index=1:no_images


dist_1 = (u_test1(index,:)-mean_hsbe)*cov_hsbe*(u_test1(index,:)-mean_hsbe)'           
                
dist_2 = (u_test2(index,:)-mean_msbe)*cov_msbe*(u_test2(index,:)-mean_msbe)'
                
dist_3 = (u_test3(index,:)-mean_lsbe)*cov_lsbe*(u_test3(index,:)-mean_lsbe)'

 dist = dist_1;
 value = 3
 if dist_1>dist_2
   dist = dist_2;
   value = 2;
 end
 if  dist>dist_3
  value = 1;
 end;

% print the euclidean distances to the output file
%
fprintf(FID, 'image %d   %f  %f %f result = %d \n', index, dist_1, dist_2, dist_3, value);

end;

% close the output file
%
fclose(FID);

% end of file












