%% Demo code to compare different ICA algorithms
% Algorithms: Infomax, FastICA (kurtosis and negentropy nonlinearities), EBM, EMK, ERBM
% Dataset : Face images (Faces in the wild) http://vis-www.cs.umass.edu/lfw/
% Package to compare ICA codes are available in Machine Learning for Signal Processing (MLSP) Lab, UMBC website
% http://mlsp.umbc.edu/resources.html
% ICA codes:
% EBM,EMK,ERBM: MLSP website http://mlsp.umbc.edu/resources.html
% FastICA: https://www.mathworks.com/matlabcentral/fileexchange/38300-pca-and-ica-package
% Infmax: https://trendscenter.org/software/gift/
% Code written by Mohammad Abu Baker Siddique Akhonda (mo32@umbc.edu)
% ----------------------------
% UPDATE: 
% 
%   Based on the original package, JADE, WASOBI, MULTI-COMBI algorithms are
%   added. ISI is added as an extra performance metric.
%   The algorithms are sorted: 
%         Infomax --> FastICA (the two versions) --> JADE ---> WASOBI 
%         ---> EFICA ---> EBM ---> EMK ---> MULTI-COMBI ---> ERBM 
%   WASOBI, MULTI-COMBI: http://si.utia.cas.cz/downloadPT.htm
%   Last updated by Hanlu Yang (hyang3@umbc.edu) on 03/10/2020



clc
clear all
close all
set(groot, ...
'DefaultFigureColor', 'w', ...
'DefaultAxesLineWidth', 5.5, ...
'DefaultAxesXColor', 'k', ...
'DefaultAxesYColor', 'k', ...
'DefaultAxesFontUnits', 'points', ...
'DefaultAxesFontSize', 40, ...
'DefaultAxesFontName', 'Helvetica', ...
'DefaultLineLineWidth', 1, ...
'DefaultTextFontUnits', 'Points', ...
'DefaultTextFontSize', 45, ...
'DefaultTextFontName', 'Helvetica', ...
'DefaultAxesBox', 'off', ...
'DefaultAxesTickLength', [0.02 0.025]);
% %%
% startingFolder = pwd; % Whatever.
% addpath(genpath([startingFolder]));
% folder = uigetdir(startingFolder);
% filePattern = fullfile(folder, '*.jpg');
% DIR = dir(filePattern);
%%DIR = dir('/*.jpg');% change this one as needed
%% User input files
[file,path] = uigetfile('*.jpg','MultiSelect', 'on');
Num_images = length(file);
if isequal(file,0)
   disp('User selected Cancel');
else
   disp(['User selected ', fullfile(path,file)]);
end
%% Reading data and creating sources
for i = 1:Num_images
    normal_images = double(imread(fullfile(path,file{i})));
    S(i,:)= normal_images(:);
end
%% Properties of the sources
%Distribution
n_size = size(S,1);
temp_range = [1:(n_size*2)];
range = permute(reshape(temp_range,2,[]),[2,1]);
figure('Name','Original Signals')
for i = 1:n_size   
   subplot(n_size,2,range(i,1))
   imagesc(reshape(S(i,:),168,168)); colormap(gray)
   xticklabels([])
   yticklabels([])
   %title(['Source = ' num2str(i)])
   subplot(n_size,2,range(i,2))
   t = histfit(S(i,:),50,'kernel')
   t(1).FaceColor = [0.8 0.8 1]
   xlim([0 300]);
   xticklabels([])
   yticklabels([])
end
clear temp_range range
%% Dependence
for i = 1:n_size
    for j = 1:n_size
        Corr_S(i,j) = corr2(S(i,:),S(j,:));
    end
end
figure('Name','2D and 3D Correlation map')
subplot(1,2,1)
imagesc(abs(Corr_S));
subplot(1,2,2)
b = bar3(abs(Corr_S));
colorbar
for k = 1:length(b)
    zdata = b(k).ZData;
    b(k).CData = zdata;
    b(k).FaceColor = 'interp';
end
%% Necessary parameters
num_iter=1; %number of iterations to run test code, only thing that changes is the random mixing matrix
num_alg = 10; %number of algorithm used(Infomax, FastICA_kurtosis, FastICA_negentropy,EBM,EMK and ERBM)
alg_trick = [1:10];
[m,n] = size(normal_images);
ker_EMK = 3; % 3 is default number of kernels. Change it for different kernel number

%Remove overall mean of new source matrix:
tempVar = mean(S,2)*ones(1,size(S, 2));
S_zero_mean = S - tempVar;
mean_of_S=tempVar;
clear tempVar;

%%
for iter=1:num_iter
    iter
    %% Create random mixing matrix for every iteration if needed
    k=0;
    while (k ~= n_size)
        MIX_MATRIX=randn(n_size);
        c = cond(MIX_MATRIX);
        k = rank(MIX_MATRIX);
    end
%     MIX_MATRIX=randn(4,n_size);

    A_MATRIX(:,:,iter) = MIX_MATRIX;
    X=MIX_MATRIX* S; %X=AS
    %Remove mean of new mixture matrix (this step is not necessary in this example, since the sources were already zero mean)
    tempVar = mean(X,2)*ones(1,size(X, 2));
    X_zero_mean = X - tempVar;
    clear tempVar;
    
    for alg=1:num_alg
        
        switch alg_trick(alg)
            
            case 1
                % Infomax    
                if  exist('ica_fuse_runica') == 2
                    fprintf('\n')
                    disp(['Infomax now..']);
                    tic
                    [W_info(:,:,iter),~,~,~,~,~,~] = ica_fuse_runica(X_zero_mean);
                    t_info(:,iter) = toc;
                    S_info(:,:,iter) = W_info(:,:,iter)*X_zero_mean;
                    icasig=S_info(:,:,iter);
                else
                    disp(['Infomax code not available..']);
                end
            
            case 2
                % fastICA with kurtosis
                fprintf('\n')
                disp(['fastICA kurtosis now..']);
                tic
                [S_fast_k(:,:,iter), W_fast_k(:,:,iter),~,~] = fastICA(X_zero_mean,n_size,'kurtosis',0);
                t_fastICA_k(:,iter) = toc;
                icasig=S_fast_k(:,:,iter);
                
            case 3
                % fastICA with negentropy
                fprintf('\n')
                disp(['fastICA Neg Entropy now..']);
                tic
                [S_fast_n(:,:,iter), W_fast_n(:,:,iter),~,~] = fastICA(X_zero_mean,n_size,'negentropy',0);
                t_fastICA_n(:,iter) = toc;
                icasig=S_fast_n(:,:,iter);
                
            case 4
                % JADE
                fprintf('\n')
                disp(['JADE now..']);
                tic
                [W_JADE(:,:,iter),~]=jadem2(X_zero_mean);
                t_JADE(:,iter) = toc;
                S_JADE(:,:,iter) = W_JADE(:,:,iter)*X_zero_mean;
                icasig=S_JADE(:,:,iter);
                
                
            case 5
                % WASOBI
                fprintf('\n')
                disp(['WASOBI now..']);
                tic
                [W_WASOBI(:,:,iter),~,~]=iwasobi(X_zero_mean,5);
                t_WASOBI(:,iter) = toc;
                S_WASOBI(:,:,iter) = W_WASOBI(:,:,iter)*X_zero_mean;
                icasig=S_WASOBI(:,:,iter);
                

                
            case 6
                % EFICA
                fprintf('\n')
                disp(['EFICA now..']);
                tic
                [W_EFICA(:,:,iter),~,~,~,~,~]=efica(X_zero_mean);
                t_EFICA(:,iter) = toc;
                S_EFICA(:,:,iter) = W_EFICA(:,:,iter)*X_zero_mean;
                icasig=S_EFICA(:,:,iter);
                

            case 7
                % ICA-EBM
                fprintf('\n')
                disp(['ICA-EBM now..']);
                tic
                W_EBM(:,:,iter)=ICA_EBM(X_zero_mean);
                t_EBM(:,iter) = toc;
                S_EBM(:,:,iter)=W_EBM(:,:,iter)*X_zero_mean; %Estimated source
                icasig=S_EBM(:,:,iter);

           case 8
                % ICA-EMK
                fprintf('\n')
                disp(['ICA-EMK3 now..']);
                tic
                W_EMK3(:,:,iter)=ICA_EMK(X_zero_mean,ker_EMK);% change number of kernel by changing the value 'ker_EMK'
                t_EMK(:,iter) = toc;
                S_EMK3(:,:,iter)=W_EMK3(:,:,iter)*X_zero_mean; %Estimated source
                icasig=S_EMK3(:,:,iter);
            case 9
                % MultiCOMBI
                fprintf('\n')
                disp(['MultiCOMBI now..']);
                tic
                [W_MultiCOMBI(:,:,iter),~,~,~,~,~,~]=multicombi(X_zero_mean);
                t_MultiCOMBI(:,iter) = toc;
                S_MultiCOMBI(:,:,iter) = W_MultiCOMBI(:,:,iter)*X_zero_mean;
                icasig=S_MultiCOMBI(:,:,iter);
            case 10
                % ERBM
                fprintf('\n')
                disp(['ICA-ERBM now..']);
                tic
                W_ERBM(:,:,iter)=ERBM(X_zero_mean);
                t_ERBM(:,iter) = toc;
                S_ERBM(:,:,iter) = W_ERBM(:,:,iter)*X_zero_mean;
                icasig=S_ERBM(:,:,iter);
                
            
                                
        end
           
        %Re-shape estimated sources into 2D matrices (display them by un-commenting lines inside for loop):
        S_hat=zeros(n_size,m*n);
        for ij=1:n_size
            S_hat(ij,:)=icasig(ij,:);
        end
        %Find correspondence between estimated and original sources:
        zzz=cat(1,S,S_hat);
        R = corrcoef(zzz');
        R_imp=abs(R(1:n_size,n_size+1:end));
        
        max_val=zeros(n_size,1);
        max_idx=zeros(n_size,1);
        for ij=1:n_size
            [max_val(ij) max_idx(ij)]=max(R_imp(ij,:));
        end
        %In cases were source estimation was not one-to-one, we use auction algorithm
        %for assignment.
        pro_dec=0;
        if (length(unique(max_idx)) < n_size)
            %                 pause
            disp(['Using auction algorithm to determine correspondence between truth and estimated sources in alg ' num2str(alg) ' and iter ' num2str(iter)]);
            pro_dec=1;
            [colsol, rowsol] = auction((1-R_imp)');
            max_idx=colsol';
            for ij=1:n_size
                max_val(ij)=R_imp(ij,max_idx(ij));
            end
        end
        
        %Saving correlation values for each algorithm:
        Corr_ind(alg).max_idx=max_idx;
        Corr_val(alg).max_val=max_val;
    end
    
    for jj = 1:num_alg
        Cor_max(iter,:,jj) = Corr_val(jj).max_val';
        id_max(iter,:,jj) = Corr_ind(jj).max_idx';
    end
    i=1;
    
end
disp(['All algorithm completed..']);
%% Finding best run if multiple run used. For single run will return consistent run index as 1;
[consistentRunidx,crossISIW]= crossISIidx(W_fast_k);
[consistentRunidx0,crossISIW0]= crossISIidx(W_fast_n);
[consistentRunidx1,crossISIW1]= crossISIidx(W_EBM);
[consistentRunidx2,crossISIW2]= crossISIidx(W_EMK3);
[consistentRunidx3,crossISIW3]= crossISIidx(W_ERBM);
[consistentRunidx4,crossISIW4]= crossISIidx(W_info);
[consistentRunidx5,crossISIW5]= crossISIidx(W_JADE);
[consistentRunidx6,crossISIW6]= crossISIidx(W_WASOBI);
[consistentRunidx7,crossISIW7]= crossISIidx(W_MultiCOMBI);
[consistentRunidx8,crossISIW8]= crossISIidx(W_EFICA);


%% Calculate ISI for the best run 
% ISIW= calculate_ISI(W_fast_k(:,:,consistentRunidx), A_MATRIX(:,:,consistentRunidx));
% ISIW0= calculate_ISI(W_fast_n(:,:,consistentRunidx0), A_MATRIX(:,:,consistentRunidx0));
% ISIW1= calculate_ISI(W_EBM(:,:,consistentRunidx1), A_MATRIX(:,:,consistentRunidx1));
% ISIW2= calculate_ISI(W_EMK3(:,:,consistentRunidx2), A_MATRIX(:,:,consistentRunidx2));
% ISIW3= calculate_ISI(W_ERBM(:,:,consistentRunidx3), A_MATRIX(:,:,consistentRunidx3));
% ISIW4= calculate_ISI(W_info(:,:,consistentRunidx4), A_MATRIX(:,:,consistentRunidx4));
% ISIW5= calculate_ISI(W_JADE(:,:,consistentRunidx5), A_MATRIX(:,:,consistentRunidx5));
% ISIW6= calculate_ISI(W_WASOBI(:,:,consistentRunidx6), A_MATRIX(:,:,consistentRunidx6));
% ISIW7= calculate_ISI(W_MultiCOMBI(:,:,consistentRunidx7), A_MATRIX(:,:,consistentRunidx7));

%% Calculate ISI for all runs 
for isi_iter=1:num_iter
    all_ISIW(isi_iter)= calculate_ISI(W_fast_k(:,:,isi_iter), A_MATRIX(:,:,isi_iter));
    all_ISIW0(isi_iter)= calculate_ISI(W_fast_n(:,:,isi_iter), A_MATRIX(:,:,isi_iter));
    all_ISIW1(isi_iter)= calculate_ISI(W_EBM(:,:,isi_iter), A_MATRIX(:,:,isi_iter));
    all_ISIW2(isi_iter)= calculate_ISI(W_EMK3(:,:,isi_iter), A_MATRIX(:,:,isi_iter));
    all_ISIW3(isi_iter)= calculate_ISI(W_ERBM(:,:,isi_iter), A_MATRIX(:,:,isi_iter));
    all_ISIW4(isi_iter)= calculate_ISI(W_info(:,:,isi_iter), A_MATRIX(:,:,isi_iter));
    all_ISIW5(isi_iter)= calculate_ISI(W_JADE(:,:,isi_iter), A_MATRIX(:,:,isi_iter));
    all_ISIW6(isi_iter)= calculate_ISI(W_WASOBI(:,:,isi_iter), A_MATRIX(:,:,isi_iter));
    all_ISIW7(isi_iter)= calculate_ISI(W_MultiCOMBI(:,:,isi_iter), A_MATRIX(:,:,isi_iter));
    all_ISIW8(isi_iter)= calculate_ISI(W_EFICA(:,:,isi_iter), A_MATRIX(:,:,isi_iter));

end


all_ISI = [all_ISIW4;all_ISIW;all_ISIW0;all_ISIW5;all_ISIW6;all_ISIW8;...
    all_ISIW1;all_ISIW2;all_ISIW7;all_ISIW3;];
%% Plotting average correlation values
if num_iter > 1
    Data = squeeze(mean(Cor_max,2));
    figure('Name','Correlation Values')
    h = boxplot(Data);
    set(h,{'linew'},{3})
    %         Infomax --> FastICA (the two versions) --> JADE ---> WASOBI 
%         ---> EFICA ---> EBM ---> EMK ---> MULTI-COMBI ---> ERBM 

    xticklabels({'Infomax','FastICA(Kurtosis)','FastICA(NegEntropy)','JADE',...
        'WASOBI','EFICA','EBM','EMK','MultiCOMBI','ERBM'})
    xtickangle(320)
    ylabel('Average correlation')
else
    figure('Name','Average Correlation Values')
    y = [ mean(mean(Cor_max(:,:,1))) mean(mean(Cor_max(:,:,2))) ...
        mean(mean(Cor_max(:,:,3))) mean(mean(Cor_max(:,:,4))) ...
        mean(mean(Cor_max(:,:,5))) mean(mean(Cor_max(:,:,6))) ...
        mean(mean(Cor_max(:,:,7))) mean(mean(Cor_max(:,:,8))) ...
        mean(mean(Cor_max(:,:,9))) mean(mean(Cor_max(:,:,10)))];
    bar(y);
    title('Average Correlation Values')
    ylim([0.2 1]);
    xticks([1:num_alg])

    xticklabels({'Infomax','FastICA(Kurtosis)','FastICA(NegEntropy)','JADE',...
        'WASOBI','EFICA','EBM','EMK','MultiCOMBI','ERBM'})
    xtickangle(320)
    ylabel('Average correlation')
end
    %% Plotting average timing performance
if num_iter > 1
    t = [(t_info);(t_fastICA_k); (t_fastICA_n); (t_JADE); (t_WASOBI);...
        (t_EFICA); (t_EBM); (t_EMK); ...
        (t_MultiCOMBI); (t_ERBM)];
    figure('Name','Execution timing')
    h = boxplot(t');
    set(h,{'linew'},{3})
%     xticklabels({'Infomax','FastICA(Kurtosis)','FastICA(NegEntropy)',...
%         'EBM','EMK','ERBM','JADE', 'WASOBI','MultiCOMBI','EFICA'})
    xticklabels({'Infomax','FastICA(Kurtosis)','FastICA(NegEntropy)','JADE',...
        'WASOBI','EFICA','EBM','EMK','MultiCOMBI','ERBM'})
    xtickangle(320)
    ylabel('Execution time (s)')
else
    figure('Name','Average execution timing')
    t = [mean(t_info) mean(t_fastICA_k)  mean(t_fastICA_n) mean(t_JADE)...
        mean(t_WASOBI) mean(t_EFICA) mean(t_EBM) ...
        mean(t_EMK) mean(t_MultiCOMBI) mean(t_ERBM)];
    bar(t);
    title('Average execution time')
    xticks([1:num_alg])
    xticklabels({'Infomax','FastICA(Kurtosis)','FastICA(NegEntropy)','JADE',...
        'WASOBI','EFICA','EBM','EMK','MultiCOMBI','ERBM'})
    xtickangle(320)
    ylabel('Execution time (s)')
end

%% Plotting ISI values
if num_iter > 1
    Data = squeeze(mean(all_ISI,2));
    figure('Name','Average ISI Values')
%     h = boxplot(Data);
    bar(1:num_alg, Data)
%     set(h,{'linew'},{3})

    xticklabels({'Infomax','FastICA(Kurtosis)','FastICA(NegEntropy)', ...
        'JADE','WASOBI','EFICA','EBM','EMK','MultiCOMBI','ERBM'})
    xtickangle(320)
    ylabel('Average ISI')
else
    figure('Name','ISI Values')
    y = [ all_ISI(1,:) all_ISI(2,:) all_ISI(3,:) all_ISI(4,:) ...
        all_ISI(5,:) all_ISI(6,:) all_ISI(7,:) all_ISI(8,:) ...
        all_ISI(9,:) all_ISI(10,:)];
    bar(y);
    title('ISI Values')
    ylim([0.2 1]);
    xticks([1:num_alg])
    xticklabels({'Infomax','FastICA(Kurtosis)','FastICA(NegEntropy)', ...
        'JADE','WASOBI','EFICA','EBM','EMK','MultiCOMBI','ERBM'})
    xtickangle(320)
    ylabel('ISI')
end
    %%
    S_hat_all(:,:,1) = S_info(:,:,consistentRunidx4);
    S_hat_all(:,:,2) = S_fast_k(:,:,consistentRunidx); 
    S_hat_all(:,:,3) = S_fast_n(:,:,consistentRunidx0); 
    S_hat_all(:,:,4) = S_EBM(:,:,consistentRunidx1); 
    S_hat_all(:,:,5) = S_EMK3(:,:,consistentRunidx2); 
    S_hat_all(:,:,6) = S_ERBM(:,:,consistentRunidx3);
    S_hat_all(:,:,7) = S_JADE(:,:,consistentRunidx5);
    S_hat_all(:,:,8) = S_WASOBI(:,:,consistentRunidx6);
    S_hat_all(:,:,9) = S_MultiCOMBI(:,:,consistentRunidx7);
    S_hat_all(:,:,10) = S_EFICA(:,:,consistentRunidx8);


    %A_name = ["FastICA(Kurtosis)", "FastICA(NegEntropy)","EBM","EMK","ERBM"];
    A_name = [{'Infomax'}, {'FastICA(Kurtosis)'}, {'FastICA(NegEntropy)'},{'EBM'},{'EMK'},{'ERBM'},{'JADE'},{'WASOBI'},{'MultiCOMBI'},{'EFICA'}];
    
    %%
    disp(['Generating Figures ..']);
    disp(['Figure 1:Comparison ..']);
    disp(['Figure 2:Original Signals ..']);
    disp(['Figure 3:Mixture ..']);
    disp(['Figure 4:Infomax ..']);
    disp(['Figure 5:FastICA kurtosis ..']);
    disp(['Figure 6:FastICA Negentropy ..']);
    disp(['Figure 7:EBM ..']);
    disp(['Figure 8:EMK-3 ..']);
    disp(['Figure 9:ERBM ..']);
    disp(['Figure 10:JADE ..']);
    disp(['Figure 11:WASOBI ..']);
    disp(['Figure 12:MultiCOMBI ..']);

    


    %% No necessary. Only needed for sign correction of the estimated images
    for i = 1: size(S_hat_all,3)
        for t1  = 1:Num_images
            for t2 = 1:Num_images
                temp_corr(t1,t2,i) = corr2(S(t1,:),S_hat_all(t2,:,i));
                if sign(temp_corr(t1,t2,i)) ~= 1
                    S_hat_all(t2,:,i) = -1* S_hat_all(t2,:,i);
                end
            end
        end
    end  
    %% Plotting Algo results
    for i = 1:num_alg
        figure('Name',A_name{i})
        
        for j = 1:Num_images
            subplot(1,Num_images,j)
            imagesc(reshape(S_hat_all(id_max(1,j,i),:,i),168,168)); colormap(gray)
            xticklabels([])
            yticklabels([])
        end
    end

%% Plotting Signals
figure('Name','Original Signals')
for j = 1:Num_images
    subplot(1,Num_images,j)
    imagesc(reshape(S(j,:),168,168)); colormap(gray)
    xticklabels([])
    yticklabels([])
end

%% Plotting mixtures
figure('Name','Mixtures')
for j = 1:Num_images
    subplot(1,Num_images,j)
    imagesc(reshape(X(j,:),168,168)); colormap(gray)
    xticklabels([])
    yticklabels([])
end

%%
