function [sesInfo] = icatb_calculateICA(sesInfo, statusHandle)

% ICA is performed on the reduced data set
% Input: sesInfo - structure containing all parameters necessary for group
% ica analysis


if ~exist('sesInfo','var')
    %P=icatb_spm_get(1,'*.mat','Select Parameter File');
    [P] = icatb_selectEntry('typeEntity', 'file', 'title', 'Select Parameter File', 'filter', '*param*.mat');
    if isempty(P)
        error('Parameter file is not selected for analysis');
    end
    [pathstr,fileName]=fileparts(P);
    outputDir = pathstr;
    sesInfo.outputDir = outputDir;
    % Make sure parameter file exists
    load(P);
    if(~exist('sesInfo','var'))
        %         infoCell{1} = P;
        %         icatb_error('The selected file does not contain sesInfo variable', infoCell);
        error(['The selected file ', P, ' does not contain the sesInfo variable']);
    end
else
    outputDir = sesInfo.outputDir;
end

if ~exist('statusHandle', 'var')
    statusHandle = [];
end

if sesInfo.isInitialized == 0
    error('Parameter file has not been initialized');
end


icatb_defaults;
global INDIVIDUAL_ICA_INDEX;
global GROUP_ICA_INDEX;
global agg_FileName;
global WRITE_COMPLEX_IMAGES;
global NUM_RUNS_GICA;

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


% Modality type
[modalityType, dataTitle] = icatb_get_modality;

if isempty(NUM_RUNS_GICA)
    NUM_RUNS_GICA = 1;
end

if NUM_RUNS_GICA < 1
    numRuns = 1;
else
    numRuns = ceil(NUM_RUNS_GICA);
end

if ~isfield(sesInfo.userInput, 'dataType')
    dataType = 'real';
else
    dataType = sesInfo.userInput.dataType;
end

sesInfo.dataType = dataType;

pcain = [sesInfo.data_reduction_mat_file,num2str(sesInfo.numReductionSteps),'-',num2str(1), '.mat'];

load(fullfile(outputDir, pcain));

%contains concatenated PC's
data = pcasig;

%number of components to extract
numOfIC = sesInfo.numComp;

%get dimensions of data
xdim = sesInfo.HInfo.DIM(1);ydim = sesInfo.HInfo.DIM(2);zdim = sesInfo.HInfo.DIM(3);

% data = reshape(data,xdim*ydim*zdim,size(data,2));
mask_ind = sesInfo.mask_ind;

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

icaAlgo = icatb_icaAlgorithm; % available ICA algorithms

algoVal = sesInfo.algorithm; % algorithm index

% selected ICA algorithm
selected_ica_algorithm = deblank(icaAlgo(algoVal, :));

ICA_Options = {};
% get ica options
if isfield(sesInfo.userInput, 'ICA_Options')
    [ICA_Options] = sesInfo.userInput.ICA_Options;
end

% convert to cell
if isempty(ICA_Options)
    if ~iscell(ICA_Options)
        ICA_Options = {};
    end
end

if ~isempty(statusHandle)
    % get the status handles
    statusData = get(statusHandle, 'userdata');
    statusData.perCompleted = statusData.perCompleted + statusData.unitPerCompleted;
    set(statusHandle, 'userdata', statusData);
    set(statusHandle, 'name', [num2str(round(statusData.perCompleted*100)), '% analysis done']); waitbar(statusData.perCompleted, statusHandle);
end


% Changed this code to allow the users add their own ICA algorithm
% transpose data to equal components by volume
data = data';
if strcmpi(selected_ica_algorithm, 'semi-blind infomax')
    ICA_Options = icatb_sbica_options(ICA_Options, dewhiteM);
    sesInfo.userInput.ICA_Options = ICA_Options;
end

% calculate ICA
%[icaAlgo, W, A, icasig_tmp] = icatb_icaAlgorithm(selected_ica_algorithm, data, ICA_Options);


fprintf('\n');
disp(['Number of times ICA will run is ', num2str(numRuns)]);

% Loop over number of runs
for nRun = 1:numRuns
    
    fprintf('\n');
    
    disp(['Run ', num2str(nRun), ' / ', num2str(numRuns)]);
    
    fprintf('\n');
    
    % Run ICA
    [icaAlgo, W, A, icasig_tmp] = icatb_icaAlgorithm(selected_ica_algorithm, data, ICA_Options);
    
    %A = dewhiteM*pinv(W);
    
    if (nRun == 1)
        icasig2 = zeros(numRuns, size(icasig_tmp, 1), size(icasig_tmp, 2));
    end
    
    if (nRun > 1),
        
        rho = zeros(size(icasig_tmp, 1), size(icasig_tmp, 1));
        for k1 = 1:size(icasig_tmp, 1)
            for k2 = 1:size(icasig_tmp, 1)
                rho(k1, k2) = icatb_corr2(flatrow(icasig_tmp(k1, :)), flatrow(icasig2(1, k2, :)));
            end
        end
        % rho = (rho1+rho2)/2;
        
        Y = zeros(1, size(icasig_tmp, 1));
        I = Y;
        Ys = Y;
        
        for k = 1:size(icasig_tmp,1)
            [Y(k) I(k)] = max(abs(rho(:,k)));
            Ys(k) = rho(I(k),k);%get signed correlation
            rho(I(k), k) = 0;
        end;
        
        %reorder and force to be positively correlated
        icasig_tmp = sign(repmat(Ys', 1, size(icasig_tmp,2))).*icasig_tmp(I,:);
        A = sign(repmat(Ys, size(A,1),1)).*A(:,I);
        
    end
    
    % store icasig_tmp and A
    icasig2(nRun, :, :) = icasig_tmp;
    A2(nRun, :, :) = A;
    
end
% end loop over number of runs

if numRuns > 1
    icasig_tmp = squeeze(mean(icasig2));
    A = squeeze(mean(A2));
    clear W;
    W = pinv(A);
end

clear icasig2;
clear A2;

clear data;

fprintf('\n');

if ~isempty(statusHandle)
    
    % get the status handles
    statusData = get(statusHandle, 'userdata');
    statusData.perCompleted = statusData.perCompleted + statusData.unitPerCompleted;
    set(statusHandle, 'userdata', statusData);
    set(statusHandle, 'name', [num2str(round(statusData.perCompleted*100)), '% analysis done']); waitbar(statusData.perCompleted, statusHandle);
end

clear data;clear sphere;clear signs;clear bias;clear lrates;
%%icasig=zeros(numOfIC,xdim*ydim*zdim);
%%icasig(:,mask_ind)=icasig_tmp;
icasig = icasig_tmp;
clear icasig_tmp;


if strcmpi(modalityType, 'fmri')
    
    %force group images to be positive
    for compNum=1:numOfIC
        changeSign = 0;
        [Y I] =max(abs(icasig(compNum,:)));
        highestVal = icasig(compNum,I);
        sumICA = sum(icasig(compNum,find(abs(icasig(compNum,:))>mean(abs(icasig(compNum,:))))));
        moreNeg = sumICA < 0;
        % 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

% enforce dimensions to save in analyze format
if(size(icasig,1)>size(icasig,2))
    icasig = icasig';
end

% save 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
    
    %[pp, fileName] = fileparts(sesInfo.userInput.param_file);
    %save(fullfile(outputDir, [fileName, '.mat']), 'sesInfo');
    
    % 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');
disp('---------------------------------------------------------------------');
disp('');

if ~isempty(statusHandle)
    
    % get the status handles
    statusData = get(statusHandle, 'userdata');
    statusData.perCompleted = statusData.perCompleted + statusData.unitPerCompleted;
    set(statusHandle, 'userdata', statusData);
    set(statusHandle, 'name', [num2str(round(statusData.perCompleted*100)), '% analysis done']); waitbar(statusData.perCompleted, statusHandle);
end


function [data] = flatrow(data)

data = data(:);