function [aggCoefRanking,indCoefRanking]=nic_rsdd_aggComp(appConfig,labelResult)
% Calculate the mean component using multiple components from multiple runs
% to eventually get the final "aggregated" component, i.e., MICA component
% Input: 
%	appConfig - structure containing all parameters necessary for
%		    group
%	labelResult: a matrix specifying the component index corresponding in
% every runs
% Output:
%	aggCoefRanking: the ranking coefficient of aggragated components
%	indCoefRanking: the ranking coefficient of individual components
%-----------------------------------------------------------
%   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;
%   Last Revised by Zhang Han
%-----------------------------------------------------------
%

nic_rsdd_defaults;

global IsOutputWorstComp;

runTimes = appConfig.runTimes;
numOfIC = appConfig.compNum;

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

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


if(~exist('labelResult','var'))
    %load label Result
    load(fullfile(appConfig.outDir,'labelResult.mat'));
end


disp(' ');
disp('---------------------------------------------------------------------');
disp('STARTING AGGRAGATING COMPONENTS STEP');
disp('---------------------------------------------------------------------');

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

if 1
disp('Aggregating agg component');
% agg
icsCell={};
for j=1:runTimes
    icaDir=sprintf('mat%03d',j);
    icaIn=fullfile(outDir,icaDir,sprintf('%s_ica.mat',outPrefix));% load ICA of agg component
    ic = load(icaIn, 'icasig');
    A= load(icaIn, 'A');
    A=A.A;
    icasig = ic.icasig;
    icsCell{j}=icasig;
    clear icasig ic;
end

newicasig=[];
coefRanking=[];
dists={};
if(IsOutputWorstComp)
    worstComp=[];
end
for i=1:numOfIC
    % 
    disp(sprintf('Aggregating agg component #%d',i));
    dd=labelResult{i};
    icasigs=[];
    xxlabel=dd(2,:);
    for j=1:runTimes       
       % load data and average all agg and individual subject
        icasig= icsCell{j};
        icasigs=[icasigs,icasig(xxlabel(j),:)'];
    end
    meansig=mean(icasigs,2);
    
    coefDist=0;
    if(IsOutputWorstComp)
        theworst=1;
        theworstIndex=1;
    end
    for j=1:runTimes       
       % load data and average all agg and individual subject
        icasig= icsCell{j};
        coef=corrcoef(icasig(xxlabel(j),:),meansig);
        if(coef(1,2)<0)
               icasigs(:,j)=-icasigs(:,j);% reverse the sign of component               
        end
        if(IsOutputWorstComp)
            if(theworst>abs(coef(1,2)))
                theworst=abs(coef(1,2));
                theworstIndex=j;
            end
        end
        coefDist=coefDist+abs(coef(1,2));% to calc the average component
    end
    meansig=mean(icasigs,2);% the average component
    
%     worstComp{i}=icasigs(:,theworstIndex);
    
    corrdist=corrcoef([meansig,icasigs]); % distance
    corrdist=1-abs(corrdist(2:end,1));
       
    sigma=std(corrdist);
    mean1=mean(corrdist);
    
    % maybe this is not valid
    yindex=find(corrdist>=mean1-3*sigma);
    yindex1=find(corrdist<=mean1+3*sigma);
    yindex1=intersect(yindex,yindex1);
        
    if(appConfig.scatter.removeScatter==1)
        yindex2=find(corrdist<=appConfig.scatter.scatterThreshold);% the threshold
        yindex=intersect(yindex1,yindex2);
     else
         yindex=yindex1;
      end
      if(isempty(yindex))
         yindex=yindex1;
       end
%     length(yindex)
    
    % repeatedly calculating the average component after filtering the
    % index
    meansig=mean(icasigs(:,yindex),2);
    newicasig=[newicasig,meansig];
    
    if(IsOutputWorstComp)
        worstComp=[worstComp,icasigs(:,theworstIndex)];
    end
%     coefDist=sum(corrdist(yindex));
%     
%     coefValue=coefDist/length(yindex);
    coefDist=sum(corrdist);
    coefValue=coefDist/runTimes;
    
    coefRanking=[coefRanking coefValue];
    
    
    
    icasigs=[icasigs(:,yindex),meansig];
    D=pdist(icasigs','corr'); 
    D=1-abs(1-D);
    D(find(D<0))=0;
    % mdscale
    if(length(D)>2)            
        [Y stress]=mdscale(D,2,'Criterion', 'strain' );
        Dchanged=pdist(Y,'euclidean');
        dist.rawDist=D;
        dist.mds=Y;
        dist.scaleDist=Dchanged;
        dist.index=yindex;
        dists{i}=dist;
    else
        dists{i}=[];
    end
%     newicasig=newicasig;
    clear icasigs meansig;
end
clear icsCell;
icasig=newicasig';

icaOut=fullfile(outDir,'mat',sprintf('%s_ica.mat',outPrefix));% load ICA
save(icaOut,'icasig','coefRanking','dists');% save repeatedly
aggCoefRanking=coefRanking;
clear newicasig;
if(~isfield(appConfig,'HInfo'))
    [da pInfo]=nic_rsdd_loadData(appConfig.AppFileList{1,3}.name(1,:),[]);
    clear da;
else
    pInfo=appConfig.HInfo;
end
complexInfoRead = struct('complex_file_naming',{'R_'  'I_'},'complexType','real&imaginary');
    
outfile = [outPrefix,'_agg__component_ica_'];

% outDir=fullfile(outDir,'');

[aggFileNames, zipfilename, files_in_zip] = nic_icatb_saveICAData(outfile, icasig, A, appConfig.maskInd, numOfIC, pInfo, ...
         'real', complexInfoRead, outDir);% save the mean agg components
     
if(IsOutputWorstComp)
    outfile = [outPrefix,'_agg_component_ica_worst_'];
    [aggFileNames, zipfilename, files_in_zip] = nic_icatb_saveICAData(outfile, worstComp', A, appConfig.maskInd, numOfIC, pInfo, ...
             'real', complexInfoRead, outDir);% save the mean agg components
end
nic_waitbar(1/AllSubjectNum,sprintf('Aggregate components...'), 'Multiple-runs ICA (Post processing)','Child');
end

disp('Aggregating individual component');
indCoefRanking={};

for subj=1:AllSubjectNum
    icsCell={};% the spatial component
    tcsCell={};% the time courses
    disp(sprintf('Aggregating individual component (subject: #%d)',subj));
    
    % load files from every run
    for j=1:runTimes
        icaDir=sprintf('mat%03d',j);
        icaIn=fullfile(outDir,icaDir,sprintf('%s_ica_br%d.mat',outPrefix,subj));% load ICA
        compSet = load(icaIn, 'compSet');
        icasig = compSet.compSet(1).ic;
        icatc = compSet.compSet(1).tc;
        icsCell{j}=icasig;
        tcsCell{j}=icatc;
    end

    % calculate
    newicasig=[];
    newicatc=[];
    coefRanking=[];
    dists={};
    if(IsOutputWorstComp)
        worstComp=[];
        worstTC=[];
    end
    for i=1:numOfIC
        disp(sprintf('Aggregating individual component #%d',i));
        dd=labelResult{i};% the corresponding index of component
        icasigs=[];
        icatcs=[];
        xxlabel=dd(2,:);% a row vector containing the index of coponent
        
        for j=1:runTimes       
           % load data and average all agg and individual subject
            icasig= icsCell{j};
            icasigs=[icasigs,icasig(xxlabel(j),:)']; %corresponding component
            
            icatc= tcsCell{j};
            icatcs=[icatcs,icatc(:,xxlabel(j))];
            
        end
        
        meansig=mean(icasigs,2); % mean component
%         newicasig=[newicasig,meansig];
%         newicatc=[newicatc,mean(icatcs,2)];
        coefDist=0;
        
        theworst=1;
        theworstIndex=1;
        % adjust the direction of component
        for j=1:runTimes       
           % load data and average all agg and individual subject
            icasig= icsCell{j};
            coef=corrcoef(icasig(xxlabel(j),:),meansig);
            if(coef(1,2)<0)% reverse the direction
               icasigs(:,j)=-icasigs(:,j);
               icatcs(:,j)=-icatcs(:,j);
            end
            coefDist=coefDist+abs(coef(1,2));
            if(IsOutputWorstComp)
                if(theworst>abs(coef(1,2)))
                    theworst=abs(coef(1,2));
                    theworstIndex=j;
                end
            end
        end
        
        % calc the mean component again
        meansig=mean(icasigs,2);
        
        % valid? filter the component
        corrdist=corrcoef([meansig,icasigs]); % distance
        corrdist=1-abs(corrdist(2:end,1));
            
        sigma=std(corrdist);
        mean1=mean(corrdist);
        yindex=find(corrdist>=mean1-3*sigma);
        yindex1=find(corrdist<=mean1+3*sigma);
        yindex1=intersect(yindex,yindex1);
        
        if(appConfig.scatter.removeScatter==1)
            yindex2=find(corrdist<=appConfig.scatter.scatterThreshold);% the threshold
            yindex=intersect(yindex1,yindex2);
        else
            yindex=yindex1;
        end
        if(isempty(yindex))
            yindex=yindex1;
        end
             
        % select the component
        meansig=mean(icasigs(:,yindex),2);
        
        newicatc=[newicatc,mean(icatcs(:,yindex),2)];
        
        newicasig=[newicasig,meansig];% store every average component
        
        if(IsOutputWorstComp)
            worstComp=[worstComp,icasigs(:,theworstIndex)];

            worstTC=[worstTC,icatcs(:,theworstIndex)];
        end
        % calculate the correlation coefficient
        coefDist=sum(corrdist);
        coefValue=coefDist/runTimes;
        
        coefRanking=[coefRanking coefValue];
  
        
        % calculate representation map using mdscale
        icasigs=[icasigs(:,yindex),meansig];
        icatcs=[icatcs(:,yindex),mean(icatcs(:,yindex),2)];
        % we use the same parameter now to make it like simulated experiments
        D=pdist(icasigs','corr'); 
        D=1-abs(1-D);% calculate the abstract of distance
        D(find(D<0))=0;
        if(length(D)>2)            
            [Y stress]=mdscale(D,2,'Criterion', 'strain' );
            Dchanged=pdist(Y,'euclidean');
            dist.rawDist=D;
            dist.mds=Y;
            dist.scaleDist=Dchanged;
            dist.index=yindex;
            
            D1=pdist(icatcs','corr'); 
            D1=1-abs(1-D1);
            D1(find(D1<0))=0;
            [Y stress]=mdscale(D1,2,'Criterion', 'strain' );
            Dchanged=pdist(Y,'euclidean');
            dist.rawDistTC=D1;
            dist.mdsTC=Y;
            dist.scaleDistTC=Dchanged;
            
            dists{i}=dist;
        else
            dists{i}=[];
        end
    %     newicasig=newicasig;
        clear icasigs meansig;
    end
    clear icsCell tcsCell;
    
    compSet=struct();
     % use the component and raw data to obtain the time course
     % for every component
     [data]=nic_rsdd_loadData(appConfig.AppFileList{subj,3}.name,appConfig.maskInd);
     pinvTC=data'*pinv(newicasig');
%      compSet = setfield(compSet, compSetFields{3}, pinvTC);
     clear data;
                
    compSet.ic=newicasig';
    compSet.tc=newicatc;
    compSet.ptc=pinvTC;

    icaOut=fullfile(outDir,'mat',sprintf('%s_ica_br%d.mat',outPrefix,subj));% save path
    % compSet: spatial and temporal component
    % coefRanking: coefficient of every sample
    % dists: the mdscale structure
    save(icaOut,'compSet','coefRanking','dists');% save the file repeatedly
    
    if(IsOutputWorstComp)
        compSet=struct();
        icaOut=fullfile(outDir,'mat',sprintf('%s_ica_br%d_worst.mat',outPrefix,subj));% save path
        compSet.ic=worstComp';
        compSet.tc=worstTC;
        save(icaOut,'compSet');% save the file repeatedly
    end
    indCoefRanking{subj}=coefRanking;    
    nic_waitbar(subj/AllSubjectNum,sprintf('Aggregate components...'), 'Multiple-runs ICA (Post processing)','Child');
%     clear newicasig;
end



disp(' ');
disp('---------------------------------------------------------------------');
disp('DONE WITH AGGRAGATING COMPONENTS STEP');
disp('---------------------------------------------------------------------');

return;





