function sesInfo = icatb_dataReduction(sesInfo, statusHandle)
% Reduces each subject's data using pca. If more than one reduction
% step is specified by parameters then a concatenation step is done along
% with another pca step.

% Input: sesInfo - structure containing all parameters nessisary for group
% ica analysis.

if ~exist('sesInfo', 'var')
    [P] = icatb_selectEntry('typeEntity', 'file', 'title', 'Select Parameter File', 'filter', '*param*.mat');
    if isempty(P)
        error('Parameter file is not selected for analysis');
    end
    outputDir = fileparts(P);
    % Make sure parameter file exists
    load(P);
    if ~exist('sesInfo', 'var')
        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)
    icatb_error('Parameter file has not been initialized');
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

%%%%%%%% Form mask as a logical array %%%%%%
mask_ind = sesInfo.mask_ind;

disp(' ');
disp('------------------------------------------------------------------------------------');
disp('STARTING DATA REDUCTION (PRINCIPAL COMPONENTS ANALYSIS)');
disp('------------------------------------------------------------------------------------');

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

[sesInfo, complexInfo] = icatb_name_complex_images(sesInfo, 'read');
numReductionSteps = sesInfo.numReductionSteps;

%% Loop through each reduction step
for j=1:numReductionSteps

    %#####################################################
    % LOAD AND CONCATENATE
    %#####################################################
    if(j > 1)

        disp(['--Starting to concatenate data from reduction( time #',num2str(j),')']);
        disp(['  The Data will be concatenated to ',num2str(sesInfo.reduction(j).numOfGroupsAfterCAT),' group/groups']);


        count=1;
        for i=1:sesInfo.reduction(j).numOfGroupsAfterCAT

            %% Load Group Data and Concatenate
            numInNewGroup = sesInfo.reduction(j).numOfPrevGroupsInEachNewGroupAfterCAT(i);
            %setup size of pca
            pca = zeros(length(mask_ind), numOfPC*numInNewGroup);
            for groupIndex=1:numInNewGroup
                %load pca data
                pcain = [sesInfo.data_reduction_mat_file,num2str(j-1),'-',num2str(count), '.mat'];
                load(fullfile(outputDir, pcain), 'HInfo', 'pcasig');
                numOfPC = sesInfo.reduction(j).numOfPCBeforeCAT;
                multi = mod(groupIndex-1,numInNewGroup);
                pca(:, (1:numOfPC) + (multi*numOfPC)) = pcasig;
                clear pcasig;
                count = count + 1;
            end

            pcaout = [sesInfo.data_reduction_mat_file,num2str(j),'-',num2str(i), '.mat'];

            pcasig = pca;
            clear pca;

            %% Save concatenated pca data
            pcaout = fullfile(outputDir, pcaout);
            drawnow;

            icatb_save(pcaout, 'pcasig', 'HInfo');
        end
        disp('Done with concatenation');
        disp(' ');
    end

    %###################################################
    % LOOP THROUGH EACH GROUP AND RUN PCA
    %#################################################
    disp(['--Extracting principal components for data reduction( time #',num2str(j),' )']);
    numOfGroupsAfter = sesInfo.reduction(j).numOfGroupsAfterCAT;
    for i=1:numOfGroupsAfter
        disp(' ');
        if(j==1)
            subNum = ceil(i/sesInfo.numOfSess);
            ses = mod(i-1,sesInfo.numOfSess)+1;
            string = ['--Doing pca on subject #',num2str(subNum),' Session #',num2str(ses)];
            clear subNum;clear ses;
        else
            string = ['--Doing pca time #',num2str(j),' on group #',num2str(i)];
        end
        disp(string);


        % Load Original Data if its the first reduction else load PC's
        % from previous reduction
        if(j == 1)

            % Load data
            [data, HInfo] = icatb_loadData(sesInfo.inputFiles(i).name, sesInfo.dataType, complexInfo, 'read', [], mask_ind);
            tdim = size(data, 2);
            HInfo.DIM(4) = tdim;

        else

            pcain = [sesInfo.data_reduction_mat_file,num2str(j),'-',num2str(i), '.mat'];
            load(fullfile(outputDir, pcain), 'pcasig');
            data = pcasig;
            clear pcasig;
        end


        %Calculate pca
        numOfPC = sesInfo.reduction(j).numOfPCAfterReduction;

        % Calculate pca step
        [pcasig, dewhiteM, Lambda, V, whiteM] = icatb_calculate_pca(data, numOfPC);

        drawnow;

        clear data;

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

        end

        %% Save PCA vars
        pcaout = [sesInfo.data_reduction_mat_file, num2str(j), '-', num2str(i), '.mat'];
        pcaout = fullfile(outputDir, pcaout);
        drawnow;
        icatb_save(pcaout, 'pcasig', 'dewhiteM', 'whiteM', 'Lambda', 'V', 'HInfo');
        clear pcasig dewhiteM whiteM Lambda V

    end


    disp(['Done with data reduction( time # ', num2str(j),')' ]);
    disp(' ');

end

%% Save parameter file
[pp, fileName] = fileparts(sesInfo.userInput.param_file);
drawnow;
icatb_save(fullfile(outputDir, [fileName, '.mat']), 'sesInfo');

disp('------------------------------------------------------------------------------------');
disp('ENDING DATA REDUCTION (PRINCIPAL COMPONENTS ANALYSIS)');
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