function [sesInfo] = icatb_backReconstruct(sesInfo, statusHandle)

% Back reconstruction is done to reconstruct individual subject image maps and time courses
% Input: sesInfo - structure containing all parameters for
% 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;
    % 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. You must Initialize the parameters file');
end


if ~isempty(statusHandle)

    %     sesInfo.perCompleted = sesInfo.perCompleted + sesInfo.unitPerCompleted;
    %
    %     set(statusHandle, 'name', [num2str(round(sesInfo.perCompleted*100)), '% analysis done']); waitbar(sesInfo.perCompleted, 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

disp(' ');
disp('---------------------------------------------------------------------');
disp('STARTING BACK RECONSTRUCTION STEP');
disp('---------------------------------------------------------------------');


%load ica data
icain =[sesInfo.ica_mat_file, '.mat'];
ic = load(fullfile(outputDir, icain), 'W');

%icasig = ic.icasig;
W=ic.W;

% Get modality, data title and compset fields
[modalityType, dataTitle, compSetFields] = icatb_get_modality;

if strcmpi(modalityType, 'fmri')
    compSetFields = {'ic', 'tc'};
else
    compSetFields = {'timecourse', 'topography'};
end

if ~isempty(statusHandle)

    %     sesInfo.perCompleted = sesInfo.perCompleted + sesInfo.unitPerCompleted;
    %
    %     set(statusHandle, 'name', [num2str(round(sesInfo.perCompleted*100)), '% analysis done']); waitbar(sesInfo.perCompleted, 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

% Get data from either 2 or 3 back reconstruct steps into the same format
% If only 1 step then back reconstruct and save.
subjectCounter=1;
fileCounter=1;
numReductionSteps = sesInfo.numReductionSteps;
if( numReductionSteps ==3)
    %3 reduction steps
    nPC=sesInfo.reduction(3).numOfPCAfterReduction;
    numOfGroups=sesInfo.reduction(3).numOfGroupsBeforeCAT;

    for groupIndex=1:sesInfo.reduction(3).numOfGroupsAfterCAT

        %load data
        pcain = [sesInfo.data_reduction_mat_file, num2str(3), '-', num2str(groupIndex), '.mat'];
        fileIn = fullfile(outputDir, pcain);
        [dewhiteM, whiteM] = load_vars(fileIn);
        stackPartition = nPC;

        %reconstruct using mixing matrix
        group(groupIndex).tc = dewhiteM * pinv(W(:,(1:stackPartition)+(groupIndex-1)*stackPartition,:));
        %group(groupIndex).ic =   dewhiteM * pinv(W(:,(1:stackPartition)+(groupIndex-1)*stackPartition,:));
        group(groupIndex).ic =   whiteM' * W(:,(1:stackPartition)+(groupIndex-1)*stackPartition,:)';
        clear dewhiteM whiteM;

    end
elseif(numReductionSteps ==2 )
    %2 reduction steps
    nPC=sesInfo.reduction(2).numOfPCAfterReduction;
    numOfGroups=sesInfo.reduction(2).numOfGroupsBeforeCAT;
    for groupIndex=1:sesInfo.reduction(2).numOfGroupsAfterCAT

        %load data
        pcain = [sesInfo.data_reduction_mat_file, num2str(2), '-', num2str(groupIndex), '.mat'];
        %pc = load(fullfile(outputDir, pcain),'dewhiteM');
        fileIn = fullfile(outputDir, pcain);
        [dewhiteM, whiteM] = load_vars(fileIn);
        %dewhiteM = pc.dewhiteM;
        stackPartition = nPC;

        %reconstruct using mixing matrix
        group(groupIndex).tc =   dewhiteM * pinv(W(:,(1:stackPartition)+(groupIndex-1)*stackPartition,:));
        %group(groupIndex).ic =   dewhiteM * pinv(W(:,(1:stackPartition)+(groupIndex-1)*stackPartition,:));
        % using unming matrix to reconstruct spatial maps
        group(groupIndex).ic =   whiteM' * W(:,(1:stackPartition)+(groupIndex-1)*stackPartition,:)';
        group(groupIndex).tc = group(groupIndex).tc';
        clear dewhiteM whiteM;

    end
elseif(numReductionSteps ==1)
    %1 reduction step

    nPC=sesInfo.reduction(1).numOfPCAfterReduction;
    numOfGroups=sesInfo.reduction(1).numOfGroupsBeforeCAT;



    %load data
    pcain = [sesInfo.data_reduction_mat_file, num2str(1), '-', num2str(1), '.mat'];
    %     pc = load(fullfile(outputDir, pcain),'dewhiteM','pcasig');
    %     dewhiteM = pc.dewhiteM;
    %     pcasig = pc.pcasig;

    fileIn = fullfile(outputDir, pcain);
    [dewhiteM, whiteM, pcasig] = load_vars(fileIn);

    stackPartition = nPC;

    %reconstruct using mixing matrix
    group(1).tc =   dewhiteM * pinv(W(:,(1:stackPartition)+(1-1)*stackPartition,:));
    %group(1).ic =   pcasig * pinv(W(:,(1:stackPartition)+(1-1)*stackPartition,:));
    % using unmixing matrix to reconstruct spatial maps
    group(1).ic = W*pcasig';
    group(1).tc=group(1).tc';
    clear dewhiteM whiteM pcasig;

    %save
    compSet = struct(compSetFields{1}, [], compSetFields{2}, []);
    compSet = setfield(compSet, compSetFields{1}, group(1).ic);
    compSet = setfield(compSet, compSetFields{2}, group(1).tc);
    %compSet(1).tc = group(1).tc;
    %compSet(1).ic = group(1).ic;
    %compSet(1).ic = group(1).ic';

    subFile = [sesInfo.back_reconstruction_mat_file, num2str(1), '.mat'];
    drawnow;
    icatb_save(fullfile(outputDir, subFile), 'compSet');

end

% If 3 back reconstruction steps
if(numReductionSteps ==3)

    %unstack group into temp
    numStacked = sesInfo.reduction(2).numOfGroupsAfterCAT;
    for i=1:numStacked

        %load pca info
        stackPartition = sesInfo.reduction(2).numOfPCAfterReduction(1);
        temp(i).tc = group(1).tc((1:stackPartition)+(i-1)*stackPartition,:);
        temp(i).ic = group(1).ic((1:stackPartition)+(i-1)*stackPartition,:);
    end
    clear group;

    %reconstruct
    for i=1:numStacked
        pcain = [sesInfo.data_reduction_mat_file, num2str(2), '-', num2str(i), '.mat'];
        %pc = load(fullfile(outputDir, pcain), 'dewhiteM', 'pcasig');
        %dewhiteM = pc.dewhiteM;

        fileIn = fullfile(outputDir, pcain);
        [dewhiteM, whiteM, pcasig] = load_vars(fileIn);

        %group(i).tc = detrend( temp(i).tc' * dewhiteM',0);
        group(i).tc = temp(i).tc' * dewhiteM';
        %group(i).ic = pinv(  temp(i).ic'  *  dewhiteM' );
        group(i).ic = whiteM'*temp(i).ic;

        clear dewhiteM whiteM pcasig;
    end
    clear temp;
end

% For both 2 and 3 backreconstruct steps
if(numReductionSteps >=2)

    % unstack
    numOfGroups = sesInfo.reduction(2).numOfGroupsAfterCAT;
    for j=1:numOfGroups
        numStacked = sesInfo.reduction(2).numOfPrevGroupsInEachNewGroupAfterCAT(j);
        stackPartition = sesInfo.reduction(2).numOfPCBeforeCAT;
        for i=1:numStacked
            gr(j).temp(i).tc = group(j).tc(:,(1:stackPartition)+(i-1)*stackPartition);
            gr(j).temp(i).ic = group(j).ic((1:stackPartition)+(i-1)*stackPartition,: );

        end
    end
    clear group;

    % reconstruct
    numOfGroups = sesInfo.reduction(2).numOfGroupsAfterCAT;
    totalSubjects = sesInfo.reduction(1).numOfGroupsBeforeCAT;
    subjectCount = 1;
    fileCount = 1;
    for j=1:numOfGroups
        numStacked = sesInfo.reduction(2).numOfPrevGroupsInEachNewGroupAfterCAT(j);
        for i=1:numStacked


            %------back reconstruct
            %load pca info
            pcain = [sesInfo.data_reduction_mat_file, num2str(1), '-', num2str(subjectCount), '.mat'];
            %pc = load(fullfile(outputDir, pcain), 'dewhiteM', 'pcasig');
            %dewhiteM = pc.dewhiteM;
            %pcasig = pc.pcasig;

            fileIn = fullfile(outputDir, pcain);
            [dewhiteM, whiteM, pcasig] = load_vars(fileIn);

            %finish back reconstruction
            %compSet(1).tc = dewhiteM * gr(j).temp(i).tc';
            %compSet(1).ic = (pcasig * gr(j).temp(i).ic)';

            tempIC = (pcasig * gr(j).temp(i).ic)';
            compSet = struct(compSetFields{1}, [], compSetFields{2}, []);
            compSet = setfield(compSet, compSetFields{1}, tempIC);
            compSet = setfield(compSet, compSetFields{2}, dewhiteM * gr(j).temp(i).tc');

            % enforce dimensions to be components by voxels
            if(size(tempIC, 1) > size(tempIC, 2))
                tempIC = tempIC';
                compSet = setfield(compSet, compSetFields{1}, tempIC);
            end

            disp(['-done back reconstructing  set ',num2str(subjectCount)]);
            drawnow;


            %---------save subject's data
            %save and give feedback
            subFile = [sesInfo.back_reconstruction_mat_file,num2str(fileCount), '.mat'];
            string = ['-saving back reconstructed ica data for set ', num2str(subjectCount),' -> ',subFile];
            disp(string);
            drawnow;
            icatb_save(fullfile(outputDir, subFile), 'compSet');

            %keep track of counting variables
            fileCount = fileCount+1;
            clear compSet;
            clear dewhiteM whiteM pcasig;

            subjectCount = subjectCount+1;
        end
    end
end

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

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


disp('---------------------------------------------------------------------');
disp('DONE WITH BACK RECONSTRUCTION STEP');
disp('---------------------------------------------------------------------');
disp(' ');

if ~isempty(statusHandle)

    %     sesInfo.perCompleted = sesInfo.perCompleted + sesInfo.unitPerCompleted;
    %
    %     set(statusHandle, 'name', [num2str(round(sesInfo.perCompleted*100)), '% analysis done']); waitbar(sesInfo.perCompleted, 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 [dewhiteM, whiteM, pcasig] = load_vars(fileIn)
% Load dewhiteM and whiteM variables

varsIn = whos('-file', fileIn);

matchIndex = strmatch('whiteM', str2mat(varsIn.name));

status = 0;
if ~isempty(matchIndex)
    status = 1;
end

load(fileIn);

if ~status
    whiteM = pinv(dewhiteM);
end
