function nic_rsdd_groupStats(appConfig)
%  Group stats is performed to calculate the mean, t-map, std
% Input: appConfig - structure containing all parameters necessary for
% group
% filein: the components after calibration
% fileout: the statistical components: t_map, mean, std_map
%-----------------------------------------------------------
%   Copyright (C) 2009 Neuroimage Computing Group, State Key Laboratory of
%   Cognitive Neuroscience and Learning, Beijing normal university
%	Written by Gwen Chen
%	Version=1.0;
%	Release=2009803;
%  Last Revised by Guiwen Chen
%-----------------------------------------------------------


numOfIC = appConfig.compNum;
mask_ind=appConfig.maskInd;


outPrefix=appConfig.outPrefix;
outDir=appConfig.outDir;
outputDir=outDir;

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

%%%%%%%% Test if the data sets have different number of time points %%%
% Count the number of files in each of the datasets

%--get mean of each subject's back reconstructed ica component
%--------------------------------------------------------------

disp('--calculating mean ica component and timecourse');
numSub =  length(appConfig.AppSubjectLabels);
numSess = 1;
meanICASig = zeros(length(mask_ind), numOfIC, numSess);
numTimePoints=appConfig.AppFileList{1,2};
meanA = zeros(numTimePoints, numOfIC, numSess);
compSetFields={'ic','tc'};

nic_waitbar(0,sprintf('Statistics...'), 'Multiple-runs ICA (Post processing)','Child');

for i=1:numSub
    for ses=1:numSess

        %load next file
        subFile = [outPrefix,'_ica_c',num2str(i),'-',num2str(ses), '.mat'];
        load(fullfile(outputDir, 'mat',subFile));
        eval(['ic = ', compSetFields{1}, ';']);
        eval(['tc = ', compSetFields{2}, ';']);
        if size(ic, 2) ~= length(mask_ind)
            ic = ic(:, mask_ind);
        end
        %get mean

        meanICASig(:, :, ses) = squeeze(meanICASig(:, :, ses)) + ic';
        meanA(:, :, ses) = squeeze(meanA(:, :, ses)) + tc(1:numTimePoints, :);
        eval(['clear ', compSetFields{1}, ' ', compSetFields{2}]);
    end
end

for ses=1:numSess
    meanICASig(:, :, ses) = meanICASig(:, :, ses)./numSub;
    meanA(:, :, ses) = meanA(:, :, ses)./numSub;
    disp(['done calculating mean for session ', num2str(ses)]);
end

meanICASig_all = zeros(length(mask_ind), numOfIC, 1);
meanA_all = zeros(numTimePoints, numOfIC, 1);

% sum the mean of each session
for ses = 1:numSess
    % mean for different sessions
    meanICASig_all = meanICASig(:, :, ses) + meanICASig_all;
    meanA_all = meanA(:, :, ses) + meanA_all;
end

% mean for different sessions
meanICASig_all = meanICASig_all ./ numSess;
meanA_all = meanA_all ./ numSess;
disp('done calculating mean for different sessions');
nic_waitbar(0.2,sprintf('Statistics...'), 'Multiple-runs ICA (Post processing)','Child');

%--get variance and standard deviation
%-----------------------------------------------------
disp('--calculating variance and standard deviation of components');
varICASig = zeros(size(meanICASig));
varA = zeros(size(meanA));

% numSub = sesInfo.numOfSub;


for i=1:numSub

    for ses=1:numSess
        %load next file
        subFile = [outPrefix,'_ica_c',num2str(i),'-',num2str(ses), '.mat'];
        load(fullfile(outputDir,'mat', subFile));
        eval(['ic = ', compSetFields{1}, ';']);
        eval(['tc = ', compSetFields{2}, ';']);
        if size(ic, 2) ~= length(mask_ind)
            ic = ic(:, mask_ind);
        end
        varICASig(:, :, ses) = squeeze(varICASig(:, :, ses)) + (ic' - squeeze(meanICASig(:, :, ses))).^2;
        varA(:, :, ses) = squeeze(varA(:, :, ses)) + (tc(1:numTimePoints, :) - squeeze(meanA(:, :, ses))).^2;
        eval(['clear ', compSetFields{1}, ' ', compSetFields{2}]);
    end


end

% implement std for different timepoints

for ses =1:numSess
    %variance and standard deviation for ica components
    varICASig(:, :, ses) = squeeze(varICASig(:, :, ses))./(numSub-1);
    stdICASig(:, :, ses) = sqrt(squeeze(varICASig(:, :, ses)));
    %variance and standard deviation for timecourse
    varA(:, :, ses) = squeeze(varA(:, :, ses))./(numSub-1);
    stdA(:, :, ses) = sqrt(squeeze(varA(:, :, ses)));
end
disp('done calculating variance and standard deviation');
nic_waitbar(0.6,sprintf('Statistics...'), 'Multiple-runs ICA (Post processing)','Child');

% implement tmaps for different timepoints

%--calculate tmaps
%-----------------------------------------------------
disp('--calculating tmaps');
tmapICASig = zeros(size(meanICASig));
for ses=1:numSess
    for comp=1:numOfIC
        tmap_ind = find(stdICASig(:, comp, ses) > eps);
        divisor = squeeze(stdICASig(tmap_ind, comp, ses))./sqrt(numSub-1);
        tmapICASig(tmap_ind, comp, ses) = squeeze(meanICASig(tmap_ind, comp, ses)) ./ divisor;
    end
end
disp('done calculating tmaps');
nic_waitbar(0.8,sprintf('Statistics...'), 'Multiple-runs ICA (Post processing)','Child');

% Transpose all of these so that the dimensions are components
% by voxels
meanICASig = permute(meanICASig, [3, 2, 1]);
meanA = permute(meanA, [3, 1, 2]);

meanICASig_all = permute(meanICASig_all, [3, 2, 1]);
meanA_all = permute(meanA_all, [3, 1, 2]);

stdICASig = permute(stdICASig, [3, 2, 1]);
stdA = permute(stdA, [3, 1, 2]);

varICASig  = permute(varICASig , [3, 2, 1]);
varA = permute(varA, [3, 1, 2]);

tmapICASig = permute(tmapICASig, [3, 2, 1]);


%-SAVE RESULTS
%--------------------------------------------------
disp('...saving group stats data...');

complexInfo = struct('complex_file_naming',{'R_'  'I_'},'complexType','real&imaginary');

if(~isfield(appConfig,'HInfo'))
    [da pInfo]=nic_rsdd_loadData(appConfig.AppFileList{1,3}.name(1,:),[]);
    clear da;
else
    pInfo=appConfig.HInfo;
end

for ses=1:numSess
    %save mean images
    
    outfile = [outPrefix,'_mean__component_ica_'];
    [fileNames, zipFileName, files_in_zip] = nic_icatb_saveICAData(outfile, squeeze(meanICASig(ses,:,:)), ...
        squeeze(meanA(ses,:,:)), mask_ind, numOfIC,pInfo, 'real', complexInfo, outputDir);

    %save standard deviation images
    outfile = [outPrefix,'_std__component_ica_'];
    [fileNames, zipFileName, files_in_zip] = nic_icatb_saveICAData(outfile, ...
        squeeze(stdICASig(ses,:,:)), squeeze(stdA(ses,:,:)), mask_ind, numOfIC,pInfo, ...
        'real', complexInfo, outputDir);

    %save tmap
    outfile = [outPrefix,'_tmap__component_ica_'];
    [fileNames, zipFileName, files_in_zip] = nic_icatb_saveICAData(outfile, ...
        squeeze(tmapICASig(ses,:,:)), squeeze(meanA(ses,:,:)), mask_ind, numOfIC,pInfo, ...
        'real', complexInfo, outputDir);

end


%%% add mean for different sessions
%save mean for different sessions images
outfile = [outPrefix,'_mean_all_component_ica_'];
[meanFiles, zipFileName, files_in_zip] = nic_icatb_saveICAData(outfile, squeeze(meanICASig_all(1,:,:)), ...
    squeeze(meanA_all(1,:,:)), mask_ind, numOfIC,pInfo,'real', ...
    complexInfo, outputDir);

nic_waitbar(1,sprintf('Statistics...'), 'Multiple-runs ICA (Post processing)','Child');

disp('---------------------------------------------------------------------');
disp('ENDING GROUP STATS STEP');
disp('---------------------------------------------------------------------');
disp(' ');



