function icatb_icasso(param_file, selMode, numRuns)
% ICASSO plugin for GIFT. We now use ICASSO centroid results and store in
% ica.mat which will be used later in back reconstruction, scaling
% components steps.
%
% Inputs:
% 1. param_file - ICA parameter file
% 2. selMode - Selection mode. Options are 'randinit', 'bootstrap' and
% 'both'.
% 3. numRuns - Number of ICA runs.
%

% Defaults
icatb_defaults;

global PARAMETER_INFO_MAT_FILE;
global NUM_RUNS_GICA;
global agg_FileName;
global WRITE_COMPLEX_IMAGES;

giftPath = fileparts(which('gift.m'));
if ~isempty(giftPath)
    addpath(fullfile(giftPath, 'toolbox', 'icasso122'));
end

% load parameters file
filterP = ['*', PARAMETER_INFO_MAT_FILE, '*.mat'];

if ~exist('param_file', 'var') || isempty(param_file)
    param_file = icatb_selectEntry('typeEntity', 'file', 'title', 'Select a valid parameter file', 'filter', filterP);
end

if isempty(param_file)
    error('Parameter file is not selected');
end

outputDir = fileparts(param_file);

if isempty(outputDir)
    outputDir = pwd;
end

cd(outputDir);

load(param_file);

if ~exist('sesInfo', 'var')
    error('Not a valid parameter file');
else
    disp('Parameters file succesfully loaded');
end

% Return modality type
modalityType = icatb_get_modality;

%% output directory
sesInfo.outputDir = outputDir; % set the output directory

% Check selection mode and number of ICA runs
if (nargin < 2)
    
    tempRuns = NUM_RUNS_GICA;
    if (tempRuns < 2)
        tempRuns = 2;
    end
    
    numParameters = 1;
    modes = {'RandInit', 'Bootstrap', 'Both'}; % All modes
    inputText(numParameters).promptString = 'Select Mode';
    inputText(numParameters).uiType = 'popup';
    inputText(numParameters).answerString = modes;
    inputText(numParameters).dataType = 'string';
    inputText(numParameters).tag = 'mode';
    inputText(numParameters).enable = 'on';
    
    numParameters = numParameters + 1;
    
    inputText(numParameters).promptString = 'Enter no. of times ( > 1) you want ICA to be run.';
    inputText(numParameters).uiType = 'edit';
    inputText(numParameters).answerString = num2str(tempRuns);
    inputText(numParameters).dataType = 'numeric';
    inputText(numParameters).tag = 'numRuns';
    inputText(numParameters).enable = 'on';
    
    % Input dialog box
    answer = icatb_inputDialog('inputtext', inputText, 'Title', 'Input parameters to determine ICA algorithm reliability', ...
        'handle_visibility', 'on', 'windowStyle', 'normal');
    
    drawnow;
    
    if isempty(answer)
        error('Input parameters to determine ICA algorithm reliability are not selected');
    end
    
    % Selected mode
    selMode = answer{1};
    
    if (isempty(answer{2}) && ~isnumeric(answer{2}))
        error('Number of ICA runs must be numeric');
    end
    
    % Number of ICA runs
    numRuns = answer{2};
    
    if (numRuns < 0)
        numRuns = 1;
    end
    
end
% End for checking selection mode and number of ICA runs

if (~exist('numRuns', 'var'))
    numRuns = NUM_RUNS_GICA;
end

if (numRuns < 2)
    error('Error:ICASSO', ['You need to run ICA algorithm atleast two times inorder to use ICASSO.', ...
        '\nPlease check variable NUM_RUNS_GICA in icatb_defaults.m']);
end

if ~sesInfo.isInitialized
    % Run parameter initialization
    sesInfo = icatb_runAnalysis(sesInfo, 2);
end

% Number of components
numOfIC = sesInfo.numComp;


% Load last PCA reduction step file
pcain = [sesInfo.data_reduction_mat_file, num2str(sesInfo.numReductionSteps), '-', num2str(1), '.mat'];

if ~exist(fullfile(outputDir, pcain), 'file')
    % Run data reduction
    sesInfo = icatb_runAnalysis(sesInfo, 3);
end

load(fullfile(outputDir, pcain));

if (size(pcasig, 2) ~= numOfIC)
    % Run data reduction
    sesInfo = icatb_runAnalysis(sesInfo, 3);
    load(fullfile(outputDir, pcain));
end

% data contains concatenated PC's
data = pcasig;

% Mask indices
mask_ind = sesInfo.mask_ind;

if (size(data, 1) == prod(sesInfo.HInfo.DIM(1:3)))
    data = data(mask_ind, :);
end

try
    icaOptions = sesInfo.userInput.ICA_Options;
catch
    icaOptions = {};
end

selMode = lower(selMode);

%%% Run ICA Algorithm multiple times using ICASSO %%%

algorithmName = icatb_icaAlgorithm;
algorithmName = deblank(algorithmName(sesInfo.algorithm, :));

% if strcmpi(algorithmName, 'semi-blind infomax')
%     icaOptions = icatb_sbica_options(icaOptions, dewhiteM);
% end


% Transpose data such that it has dimensions components by voxels
data = data';


%%%%%% Calculate PCA and Whitening matrix %%%%%
% PCA
[V, Lambda] = icatb_v_pca(data, 1, numOfIC, 0, 'transpose', 'yes');

% Whiten matrix
[w, White, deWhite] = icatb_v_whiten(data, V, Lambda, 'transpose');

clear V Lambda;


disp(' ');
disp('---------------------------------------------------------------------');
disp('STARTING GROUP ICA STEP USING ICASSO');
disp('---------------------------------------------------------------------');

sR = icatb_icassoEst(selMode, data, numRuns, 'numOfPC', numOfIC, 'algoIndex', sesInfo.algorithm, ...
    'dewhiteM', deWhite, 'whiteM', White, 'whitesig', w, 'icaOptions', icaOptions);

clear data w deWhite White;

%%%% Visualization %%%%%%

sR = icassoExp(sR);

%%% Visualization & returning results

disp(['Launch Icasso visualization supposing ', num2str(numOfIC), ' estimate-clusters.']);
disp('Show demixing matrix rows.');
icassoShow(sR, 'L', numOfIC, 'estimate', 'demixing');

disp(['Launch Icasso visualization supposing ', num2str(numOfIC), ' estimate-clusters.']);
disp('Show IC source estimates (default), reduce number of lines');
disp('Collect results.');

[iq, A, W, icasig] = icassoShow(sR, 'L', numOfIC, 'colorlimit', [.8 .9]);

icassoResultsFile = fullfile(outputDir, [sesInfo.userInput.prefix, '_icasso_results.mat']);

icatb_save(icassoResultsFile, 'iq', 'A', 'W', 'icasig', 'sR', 'algorithmName');

clear sR;

% End for running ICA several times using ICASSO %%%%

fprintf('\n');

%% Force group images to be positive
if strcmpi(modalityType, 'fmri')
    % Loop over components
    for compNum = 1:numOfIC
        [Y, I] = max(abs(icasig(compNum, :)));
        highestVal = icasig(compNum, I);
        % Use a better criteria to change the sign
        if (sign(highestVal) == -1)
            disp(['Changing sign of component ', num2str(compNum)]);
            icasig(compNum, :) = icasig(compNum, :)*-1;
        end
    end
    % End loop over components
end


%% Save ICA ouput in matlab format
icaout = [sesInfo.ica_mat_file, '.mat'];
icaout = fullfile(outputDir, icaout);

drawnow;

icatb_save(icaout, 'A', 'W', 'icasig', 'mask_ind');

% Naming of complex images
[sesInfo, complexInfo] = icatb_name_complex_images(sesInfo, 'write');

% global variable necessary for defining the real&imag or magnitude&phase
WRITE_COMPLEX_IMAGES = sesInfo.userInput.write_complex_images;

if strcmpi(modalityType, 'fmri')
    
    if isfield(sesInfo, 'zipContents')
        sesInfo = rmfield(sesInfo, 'zipContents');
    end
    
    sesInfo.zipContents.zipFiles = {};
    sesInfo.zipContents.files_in_zip.name = {};
    
    if(sesInfo.numOfSub ==1 && sesInfo.numOfSess==1)
        % complex data
        if ~isreal(icasig)
            % convert to complex data structure
            icasig = complex_data(icasig);
            A = complex_data(A);
        end
        
        outfile = sesInfo.icaOutputFiles(1).ses(1).name(1,:);
        [fileNames, zipfilename, files_in_zip] = icatb_saveICAData(outfile, icasig, A, mask_ind, numOfIC, sesInfo.HInfo, ...
            sesInfo.dataType, complexInfo, outputDir);
    else
        outfile = sesInfo.aggregate_components_an3_file;
        % complex data
        if ~isreal(icasig)
            % convert to complex data
            icasig = complex_data(icasig);
            A = complex_data(A);
        end
        [aggFileNames, zipfilename, files_in_zip] = icatb_saveICAData(outfile, icasig, A, mask_ind, numOfIC, sesInfo.HInfo, ...
            sesInfo.dataType, complexInfo, outputDir);
        drawnow;
        icatb_save(fullfile(outputDir, [sesInfo.userInput.prefix, agg_FileName, '.mat']), 'aggFileNames');
    end
    
    % store the zip filenames to a structure
    sesInfo.zipContents.zipFiles{1} = zipfilename;
    sesInfo.zipContents.files_in_zip(1).name = files_in_zip;
    
end

[pp, fileName] = fileparts(sesInfo.userInput.param_file);
drawnow;

icatb_save(fullfile(outputDir, [fileName, '.mat']), 'sesInfo');

disp('---------------------------------------------------------------------');
disp('DONE CALCULATING GROUP ICA USING ICASSO');
disp('---------------------------------------------------------------------');
disp('');

%% Run back reconstruction
sesInfo = icatb_runAnalysis(sesInfo, 5);

%% Run scaling components
sesInfo = icatb_runAnalysis(sesInfo, 6);

%% Run group stats
icatb_runAnalysis(sesInfo, 7);

