%% 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
n_size = size(S,1);
%% Properties of the sources
%Distribution
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=3; %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];
%alg_trick = [1,3,5,7,10];
num_alg = length(alg_trick);
[m,n] = size(normal_images);
ker_EMK = 3; % 3 is default number of kernels. Change it for different kernel number
enable_crossISI = 1 % 0- Single run (no crossISI); 1- multiple runs (crossISI)

%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;


% %% 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
%         A_MATRIX = MIX_MATRIX;
%% Algo_run
if enable_crossISI == 0
    %0- Single run (no crossISI); 1- multiple runs (crossISI)
    num_iter = 1;
end
for iter=1:num_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
        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);
                    W_store{1} = W_info; S_store{1} = S_info; t_store{1} = t_info;
                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);
                W_store{2} = W_fast_k; S_store{2} = S_fast_k; t_store{2} = t_fastICA_k;
                
            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);
                W_store{3} = W_fast_n; S_store{3} = S_fast_n; t_store{3} = t_fastICA_n;
                
                
            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);
                W_store{4} = W_JADE; S_store{4} = S_JADE; t_store{4} = t_JADE;
                
            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);
                W_store{5} = W_WASOBI; S_store{5} = S_WASOBI; t_store{5} = t_WASOBI;

                
            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);
                W_store{6} = W_EFICA; S_store{6} = S_EFICA; t_store{6} = t_EFICA;

            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);
                W_store{7} = W_EBM; S_store{7} = S_EBM; t_store{7} = t_EBM;


           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);
                W_store{8} = W_EMK3; S_store{8} = S_EMK3; t_store{8} = t_EMK;
                
            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);
                W_store{9} = W_MultiCOMBI; S_store{9} = S_MultiCOMBI; t_store{9} = t_MultiCOMBI;
                
            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);
                W_store{10} = W_ERBM; S_store{10} = S_ERBM; t_store{10} = t_ERBM;
            
                                
        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,crossISIA]= crossISIidx(A_MATRIX);
%% 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;];
%% alternative of above section
all_ISI = nan(10,num_iter);
for alg = 1:num_alg
    W_temp = W_store{alg_trick(alg)};
for isi_iter=1:num_iter
    all_ISI(alg_trick(alg),isi_iter) = calculate_ISI(W_temp(:,:,isi_iter), A_MATRIX(:,:,isi_iter));
end
clear W_temp
end
%%
alg_trick = [1,3,5,7,10];
xticksss = {'Infomax','FastICA(Kurtosis)','FastICA(NegEntropy)','JADE',...
        'WASOBI','EFICA','EBM','EMK','MultiCOMBI','ERBM'};
for i = 1:length(alg_trick)
x_temp{i} = xticksss{alg_trick(i)};
end
%% 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.2]);
%     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 correlation values
if num_iter > 1
    Data = squeeze(mean(Cor_max,2));
    figure('Name','Correlation Values')
    h = boxplot(Data(:,alg_trick));
    set(h,{'linew'},{3})
    %         Infomax --> FastICA (the two versions) --> JADE ---> WASOBI 
%         ---> EFICA ---> EBM ---> EMK ---> MULTI-COMBI ---> ERBM 

    xticklabels(x_temp)
    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(alg_trick));
    title('Average Correlation Values')
    ylim([0.2 1.2]);
    xticks([1:num_alg])

    xticklabels(x_temp)
    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 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(alg_trick,:)');
    set(h,{'linew'},{3})
%     xticklabels({'Infomax','FastICA(Kurtosis)','FastICA(NegEntropy)',...
%         'EBM','EMK','ERBM','JADE', 'WASOBI','MultiCOMBI','EFICA'})
    xticklabels(x_temp)
    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(alg_trick));
    title('Average execution time')
    xticks([1:num_alg])
    xticklabels(x_temp)
    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(all_ISI(alg_trick,:)');
  %  bar(1:num_alg, Data)
     set(h,{'linew'},{3})

    xticklabels(x_temp)
    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 1]);
    xticks([1:num_alg])
    xticklabels(x_temp)
    xtickangle(320)
    ylabel('ISI')
end
    %%
    S_hat_all(:,:,1) = S_info(:,:,consistentRunidx);
    S_hat_all(:,:,2) = S_fast_k(:,:,consistentRunidx); 
    S_hat_all(:,:,3) = S_fast_n(:,:,consistentRunidx);     
    S_hat_all(:,:,4) = S_JADE(:,:,consistentRunidx);
    S_hat_all(:,:,5) = S_WASOBI(:,:,consistentRunidx);    
    S_hat_all(:,:,6) = S_EFICA(:,:,consistentRunidx);
    S_hat_all(:,:,7) = S_EBM(:,:,consistentRunidx); 
    S_hat_all(:,:,8) = S_EMK3(:,:,consistentRunidx); 
    S_hat_all(:,:,9) = S_MultiCOMBI(:,:,consistentRunidx);
    S_hat_all(:,:,10) = S_ERBM(:,:,consistentRunidx);
%%
for i = 1:num_alg
S_hat_all(:,:,i) = S_store{i};
end
%%
 %Infomax --> FastICA (the two versions) --> JADE ---> WASOBI 
% %         ---> EFICA ---> EBM ---> EMK ---> MULTI-COMBI ---> ERBM 

    %A_name = ["FastICA(Kurtosis)", "FastICA(NegEntropy)","EBM","EMK","ERBM"];
    A_name = [{'Infomax'}, {'FastICA(Kurtosis)'}, {'FastICA(NegEntropy)'},{'JADE'},{'WASOBI'},{'EFICA'},{'EBM'},{'EMK'},{'MultiCOMBI'},{'ERBM'}];
    
    %%
    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 ..']);

    


    %% Not necessary. Only needed for sign correction of the estimated images
    zero_test = zeros(Num_images, Num_images,10);
    S_test = S_hat_all;
    for i = 1: size(S_hat_all,3)
        for j1  = 1:Num_images
            for j2 = 1:Num_images
                temp_corr(j1,j2,i) = corr2(S(j1,:),S_hat_all(j2,:,i));
            end
            [in,lo] = max(abs(temp_corr(j1,:,i)));
            if sign(temp_corr(j1,lo,i)) ~= 1
                zero_test(j1,lo,i) = 1;
                temp = -1* S_hat_all(lo,:,i);
                S_test(lo,:,i) = temp;
            end
            clear temp in lo
        end
    end
    %%
    figure
    imagesc(reshape(S_hat_all(1,:,10),168,168)); colormap(gray)
    figure
    imagesc(reshape(S_test(3,:,9),168,168)); colormap(gray)
    
    %% Plotting Algo results
    for i = 1:10
    new_S_hat(:,:,i) =  S_hat_all(id_max(consistentRunidx,:,alg_trick(i)),:,alg_trick(i));
    end
    %%
    %S_hat_all = S_test;
    for i = 1:length(alg_trick)
        figure('Name',A_name{i})
        
        for k = 1:Num_images
        subplot(2,Num_images,k)
        imagesc(reshape(S(k,:),168,168)); colormap(gray)
        xticklabels([])
        yticklabels([])
        title(['True-',num2str(k)],'FontSize',12)
        end
        %sgtitle('I am a super title')
        for j = 1:Num_images
            subplot(2,Num_images,j+4)
            imagesc(reshape(new_S_hat(j,:,alg_trick(i)),168,168)); colormap(gray)
            xticklabels([])
            yticklabels([])
            title(['Estimates-',num2str(j)],'FontSize',12)
        end
       
        %sgtitle('I am a super title')
    end
    
    %%
    %definput = {'0','0','0','0','0','0','0','0','0','0'};
%     fig = uifigure;
%     selection = uiconfirm(fig,'Do you want to do sign-correction?','Sign-correction',...
%                         'Options',{'Yes','No'},...
%            'DefaultOption',2,'CancelOption',2);

%%
for i = 1:length(alg_trick)
    A_new_Name{1,i} = A_name{alg_trick(i)};
end
%%
selection = questdlg('Do you want to do sign-correction?', ...
	'Sign-correction', ...
	'Yes','No','No');
% Handle response
switch selection
    case 'Yes'
        answer = inputdlg(A_new_Name,...
             'Sample', [1 50])
    case 'No'
        f = msgbox('Sign-correction canceled');
end
       
           
      
       
   for i = 1: length(answer)
       ope = str2num(answer{i});
       if length(ope) ~=0 & length(ope) <= Num_images
           ope
           new_S_hat(ope,:,i) = -1*new_S_hat(ope,:,i);
       end
       
   end
   
   %% delete later
    plot_count = [0:1:(length(alg_trick)+1)]*Num_images;    
    S_hat_all = S_test;
    for k = 1:Num_images
        subplot(length(alg_trick)+2,Num_images,k)
        imagesc(reshape(S(k,:),168,168)); colormap(gray)
        xticklabels([])
        yticklabels([])
    end
 
    for k = 5:Num_images*2
        subplot(length(alg_trick)+2,Num_images,k)
        imagesc(reshape(X(k,:),168,168)); colormap(gray)
        xticklabels([])
        yticklabels([])
    end
    for i = 1:length(alg_trick)
    %    figure('Name',A_name{i})
    
        for j = 1:Num_images
            subplot(length(alg_trick)+2,Num_images,plot_count(i+2)+j)
            imagesc(reshape(S_hat_all(id_max(consistentRunidx,j,alg_trick(i)),:,alg_trick(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

%%
