function [appConfig subjectOrder]=nic_rsdd_dataReduction(appConfig,step,type,dataReductionOutDir)
% multiple-step data reduction using PCA
% Input: appConfig - structure containing all parameters 
% step: which step from 
% type: 1: change order in every runs
% dataReductionDir: output dir of each run
%-----------------------------------------------------------
%   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=20090726;
%-----------------------------------------------------------
%

reductionPara=appConfig.reductionPara;% parameter of each run
AllSubjectNum=length(appConfig.AppSubjectLabels);% sample number
numOfSubject=AllSubjectNum;
outPrefix=appConfig.outPrefix;% prefix
outDir=appConfig.outDir; % directory of output
numOfGroup=numOfSubject; % the number group of 1st step equals with the number of subject

if(~exist('step','var'))
    step=1;
end
if(~exist('type','var'))
    type=0;
end
if(~exist('dataReductionOutDir','var'))
    dataReductionOutDir='mat';
end


dataPCA=[];

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

subjectOrder=[1:AllSubjectNum];% default order

% permute random order 
if(type==1)
    rand('state', sum(100*clock));
    subjectOrder=randperm(AllSubjectNum);
end

lastDir='mat';

for numOfStep=step:reductionPara.numOfReductionSteps
   for iSubject=1:reductionPara.cfg(numOfStep).numOfGroupsAfterCAT
       if(numOfStep==1)
           % load data
           % mask data    
           disp(['--Extracting principal components for data reduction( time #',num2str(numOfStep),',subject #',num2str(iSubject),' )']);
           disp(sprintf(''));
           [data HInfo]=nic_rsdd_loadData(appConfig.AppFileList{iSubject,3}.name,appConfig.maskInd);
           appConfig.HInfo=HInfo;
%            load subject1;
           % PCA
           [dataPCA,vectorPCA,lamdaPCA,whitePCA,dewhitePCA]=nic_rsdd_pca(data,reductionPara.cfg(numOfStep).numOfPCs);
           % save pcasig_1PCA,V_1PCA,Lambda_1PCA,whiteM_1PCA,dewhiteM_1PCA
           save(fullfile(fullfile(outDir,'mat'),sprintf('%s_pca_r%d-%d',outPrefix,numOfStep,iSubject)),'dataPCA','vectorPCA','lamdaPCA','whitePCA','dewhitePCA','HInfo');
           clear vectorPCA lamdaPCA whitePCA dewhitePCA;
           if(numOfSubject==1)
               %copyfile(fullfile(outDir,lastDir,sprintf('%s_pca_r%d-%d.mat',outPrefix,1,1)),fullfile(outDir,dataReductionOutDir,sprintf('%s_pca_r%d-%d.mat',outPrefix,numOfStep,iSubject)));
           end % modified by Zhang Han @2012.7.24. to solve a bug that one subject can't be processed.
       elseif numOfStep>1 && numOfStep<reductionPara.numOfReductionSteps
           % merge data
           disp(['--Extracting principal components for data reduction( time #',num2str(numOfStep),',group #',num2str(iSubject),' )']);
           % random order
           
           numOfSubject=reductionPara.cfg(numOfStep).numSubjectInGroup;
           if(iSubject>1)
               s=sum(reductionPara.cfg(numOfStep).numStacked(1:iSubject-1));
           else
               s=0;
           end
           subjectNums=[s+1:s+reductionPara.cfg(numOfStep).numStacked(iSubject)];
           
%            if(type==1)
            subjectNums=subjectOrder(subjectNums);% get order
%            end
           
           if(length(subjectNums)==1)
               copyfile(fullfile(outDir,lastDir,sprintf('%s_pca_r%d-%d.mat',outPrefix,numOfStep-1,subjectNums)),fullfile(outDir,dataReductionOutDir,sprintf('%s_pca_r%d-%d.mat',outPrefix,numOfStep,iSubject)));
           else
               data=load_mat_pca(outDir,outPrefix,subjectNums,numOfStep-1,lastDir);
               % PCA
               % save PCA
               [dataPCA,vectorPCA,lamdaPCA,whitePCA,dewhitePCA]=nic_rsdd_pca(data,reductionPara.cfg(numOfStep).numOfPCs);
               % save pcasig_1PCA,V_1PCA,Lambda_1PCA,whiteM_1PCA,dewhiteM_1PCA
               save(fullfile(fullfile(outDir,dataReductionOutDir),sprintf('%s_pca_r%d-%d.mat',outPrefix,numOfStep,iSubject)),'dataPCA','vectorPCA','lamdaPCA','whitePCA','dewhitePCA');
                clear  vectorPCA lamdaPCA whitePCA dewhitePCA;
           end
           
           % save pcasig_2PCA,V_2PCA,Lambda_2PCA,whiteM_2PCA,dewhiteM_2PCA,strOfSu
    
       elseif numOfStep>1 && numOfStep==reductionPara.numOfReductionSteps
           % merge data
           disp(['--Extracting principal components for data reduction( time #',num2str(numOfStep),',group #',num2str(iSubject),' )']);
           numOfSubject=reductionPara.cfg(numOfStep).numSubjectInGroup;
           if(iSubject>1)
               s=sum(reductionPara.cfg(numOfStep).numStacked(1:iSubject-1));
           else
               s=0;
           end
           subjectNums=[s+1:s+reductionPara.cfg(numOfStep).numStacked(iSubject)];
           
           if(length(subjectNums)==1)
               copyfile(fullfile(outDir,lastDir,sprintf('%s_pca_r%d-%d.mat',outPrefix,numOfStep-1,subjectNums)),fullfile(outDir,dataReductionOutDir,sprintf('%s_pca_r%d-%d.mat',outPrefix,numOfStep,iSubject)));
           else
               data=load_mat_pca(outDir,outPrefix,subjectNums,numOfStep-1,lastDir);
               % PCA
               % save PCA
               [dataPCA,vectorPCA,lamdaPCA,whitePCA,dewhitePCA]=nic_rsdd_pca(data,reductionPara.cfg(numOfStep).numOfPCs);
               save(fullfile(fullfile(outDir,dataReductionOutDir),sprintf('%s_pca_r%d-%d.mat',outPrefix,numOfStep,iSubject)),'dataPCA','vectorPCA','lamdaPCA','whitePCA','dewhitePCA');
               clear vectorPCA lamdaPCA whitePCA dewhitePCA;
           end
           % eM_3PCA
       end   
   end
   if(numOfStep>1)
        lastDir=dataReductionOutDir;
    end
end

% save the order of subject
if(type==1)
    save(fullfile(outDir,dataReductionOutDir,'subjectOrder'),'subjectOrder');%save the subject order
end

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


appConfig.status.DataReductionStep=2;

saveFile=fullfile(appConfig.outDir,[appConfig.outPrefix,'_mica_parameter_info.mat']);
moicacfg=appConfig;
save(saveFile,'moicacfg');


function data=load_mat_pca(outDir,outPrefix,subjectNums,numOfStep,lastDir)
% load pca mat data

data=[];
if(~exist('lastDir','var'))
    lastDir='mat';
end

for i=1:length(subjectNums)
    iSubject=subjectNums(i);
    load(fullfile(outDir,lastDir,sprintf('%s_pca_r%d-%d.mat',outPrefix,numOfStep,iSubject)));
    data=[data,dataPCA];
end

return;
