function nic_rsdd_backReconstruct(appConfig,datareductionOutDir,ICAOutDir,subjectOrder)
% Back reconstruction to individual component
% Input: appConfig - structure containing all parameters 
% datareductionOutDir: the folder of pca 
% subjectOrder: subject order
% ICAOutDir: output dir of ICA
%-----------------------------------------------------------
%   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=20090730;
%-----------------------------------------------------------

% Please edit this file to support no order mode

% the number of subject
AllSubjectNum=length(appConfig.AppSubjectLabels);

% the dir of pca reduction
if(~exist('datareductionOutDir','var'))
    datareductionOutDir='mat';
end

if(~exist('subjectOrder','var'))
    subjectOrder=[1:AllSubjectNum];
end

if(~exist('ICAOutDir','var'))
    ICAOutDir='mat';
end

compSetFields = {'ic', 'tc'};
disp(' ');
disp('---------------------------------------------------------------------');
disp('STARTING BACK RECONSTRUCTION STEP');
disp('---------------------------------------------------------------------');

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

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

%load the final ICA data
icaout=fullfile(outDir,ICAOutDir,sprintf('%s_ica.mat',outPrefix));% load ICA result
ic = load(icaout, 'W');

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


% Get data from either 2 or 3 back reconstruct steps into the same format
% If there is only 1 step then back reconstruct and save.
subjectCounter=1;
fileCounter=1;
reductionPara=appConfig.reductionPara;
numReductionSteps = reductionPara.numOfReductionSteps;
group=[];
tempGroupCount=1;
ICSigs={};
%%%% Back reconstruction
for k=numReductionSteps:-1:1
%     if(k==1 && numReductionSteps>1)
%         break;
%     end
    nPC=reductionPara.cfg(k).numOfPCs;
    numOfGroups=reductionPara.cfg(k).numOfGroupsAfterCAT;
    if(k>1)
        stackPartition = reductionPara.cfg(k-1).numOfPCs;
    end
    
    if(k==numReductionSteps)
        for groupIndex=1:numOfGroups
            %load data
            if(k>1)
                fileIn=fullfile(outDir,datareductionOutDir,sprintf('%s_pca_r%d-%d.mat',outPrefix,k,groupIndex));
            else
                fileIn=fullfile(outDir,'mat',sprintf('%s_pca_r%d-%d.mat',outPrefix,k,groupIndex));
            end

            if(k==1)
                [dewhiteM, whiteM, pcasig] = load_vars(fileIn);

                stackPartition = nPC;% load the pcasig

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

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

                group(groupIndex,1).ic =   whiteM' * W(:,(1:stackPartition)+(groupIndex-1)*stackPartition,:)';

                clear dewhiteM whiteM;
            end
            tempGroupCount=1;     
        end    
    end
    
    %Split the data
    for groupIndex=1:numOfGroups
       numStacked=reductionPara.cfg(k).numStacked(groupIndex);
       [m n]=getGroupIndex(groupIndex,k,reductionPara);
        for i=1:numStacked
            if(k>1)
                gr(groupIndex).temp(i).tc = group(m,n).tc((1:stackPartition)+(i-1)*stackPartition,:);
                gr(groupIndex).temp(i).ic = group(m,n).ic((1:stackPartition)+(i-1)*stackPartition,: );
            else
                gr(groupIndex).temp(i).tc = group(m,n).tc;
                gr(groupIndex).temp(i).ic = group(m,n).ic;
            end
        end
    end
    clear group;
    tempGroupCount=numOfGroups;
     
    %reconstruct next
    if(k>=2)
        subjectCount=1;
        for groupIndex=1:numOfGroups
            %reconstruct
            numStacked=reductionPara.cfg(k).numStacked(groupIndex);

            for i=1:numStacked
                if(k>2)
                    fileIn=fullfile(outDir,datareductionOutDir,sprintf('%s_pca_r%d-%d.mat',outPrefix,k-1,subjectCount));%load 2 step result
                else
                    index=subjectOrder(subjectCount);
                    fileIn=fullfile(outDir,'mat',sprintf('%s_pca_r%d-%d.mat',outPrefix,k-1,index));%load 2 step result
                end

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


                group(groupIndex,i).tc = dewhiteM * gr(groupIndex).temp(i).tc;
            
                if(k==2)
                    group(groupIndex,i).ic = (pcasig*gr(groupIndex).temp(i).ic)';
                else
                     group(groupIndex,i).ic = whiteM'*gr(groupIndex).temp(i).ic;
                end
                
                subjectCount=subjectCount+1;

                clear dewhiteM whiteM pcasig;
            end % end for numStacked
        end % end for groupIndex
       clear gr;
    else %if k<2
        subjectCount=1;
        for groupIndex=1:numOfGroups
            %reconstruct
            numStacked=reductionPara.cfg(k).numStacked(groupIndex);
            
            for i=1:numStacked
                 
                  tempIC = gr(groupIndex).temp(i).ic;%(gr(groupIndex).temp(i).pcasig * gr(groupIndex).temp(i).ic)';
                  tempTC=gr(groupIndex).temp(i).tc;%gr(groupIndex).temp(i).dewhiteM * gr(groupIndex).temp(i).tc;
                  
                 
                 compSet = struct(compSetFields{1}, [], compSetFields{2}, []);
                 compSet = setfield(compSet, compSetFields{1}, tempIC);
                 compSet = setfield(compSet, compSetFields{2}, tempTC);
                 
                 % tempIC
                 % read the raw individual data and caculate the W
                 % add 2009-12-31
                 ICSigs={ICSigs tempIC}; % every subject
                 
                    % enforce dimensions to be components by voxels
                 if(size(tempIC, 1) > size(tempIC, 2))
                     tempIC = tempIC';
                     compSet = setfield(compSet, compSetFields{1}, tempIC);
                 end
                       
                 index=subjectOrder(subjectCount);
                 fileOut=fullfile(outDir,ICAOutDir,sprintf('%s_ica_br%d.mat',outPrefix,index));               
                
                save(fileOut, 'compSet');
                
                
                disp(['-done back reconstructing  set ',num2str(groupIndex)]);
                subjectCount=subjectCount+1;
            end
            
        end
        
    end
end
clear group;  

% for i=1:numOfIC
%     nsample=size(data,2);% 样本数目
%     meand=mean(data,2);%均值
%     stdd=std(data')';%标准差
%     tind=find(stdd>eps);%取非零的序号
%     ttest=zeros(size(meand));
%     
%     meand=meand*sqrt(nsample);
%     
%     ttest(tind)=meand(tind)./stdd(tind);%t-map
% end
% calculate tmap

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





function [dewhitePCA, whitePCA, dataPCA] = load_vars(fileIn)
% Load dewhiteM and whiteM and pcasig variables

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

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

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

load(fileIn);

if ~status
    whitePCA = pinv(dewhitePCA);
end


function [m n]=getGroupIndex(groupIndex,step,reductionPara)
% get the group index in each step

if(step==reductionPara.numOfReductionSteps)
    m=groupIndex;
    n=1;
    return;
end

s=0;
for i=1:reductionPara.cfg(step+1).numOfGroupsAfterCAT
    s=s+reductionPara.cfg(step+1).numStacked(i);
    if(s>=groupIndex)
        m=i;
        break;
    end
end
if(i>1)
    n=groupIndex-sum(reductionPara.cfg(step+1).numStacked(1:i-1));
else
    n=groupIndex;
end



