% Using Graphical Lasso to get the matrix with wanted sparsity
% Perform Graphical Analysis on these matrices.
%__________________________________________________________________________
% Written by YAN Chao-Gan 130101.
% The Nathan Kline Institute for Psychiatric Research, 140 Old Orangeburg Road, Orangeburg, NY 10962, USA
% The Phyllis Green and Randolph Cowen Institute for Pediatric Neuroscience, New York University Child Study Center, New York, NY 10016, USA
% ycg.yan@gmail.com



%Get the corresponding Lambda  -- Using Negative sign, correspondes to using positive correlation

FromUpDir = '/AtlasROI_SeedSignal_ExtractedInMNISpace_Detrend_Filter01_1';
ToUpDir = '/SparsityThreshold';

ConditionList={'NoGSR','GSR'};

NodeNumber=90;
SparsityRange = [0.02:0.02:0.5]';
%LassoLambdaTestSet = [0.005:0.001:1];
LassoLambdaTestSet = [0.0001:0.0001:0.001, 0.002:0.001:1];

for iCondition=1:length(ConditionList)
    
    FisherZDir = [FromUpDir,filesep,ConditionList{iCondition}];
    OutputDir = [ToUpDir,filesep,ConditionList{iCondition},'_GLasso_NegativeTheta'];
    mkdir(OutputDir);

    nSparsity = length(SparsityRange);
    
    LassoSparsitySet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity);
    LassoLambdaSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity);
    
    ThetaMatrixSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity,NodeNumber,NodeNumber);
    BinarizedMatrixSet = zeros(AutoDataProcessParameter.SubjectNum,NodeNumber,NodeNumber);

    parfor i=1:AutoDataProcessParameter.SubjectNum
        
        fprintf('Subject %s ',AutoDataProcessParameter.SubjectID{i});
            
        MAT = load ([FisherZDir,filesep,'ROISignals_ROITC',AutoDataProcessParameter.SubjectID{i},'.mat']);
        TC = MAT.ROISignals(:,493:608);
        TC = TC(:,1:NodeNumber);
        
        
        TC = (TC-repmat(mean(TC),size(TC,1),1))./repmat(std(TC),size(TC,1),1);

        
        
        [wList, thetaList, lambdaList, errors] = GraphicalLassoPath(TC,LassoLambdaTestSet);
        
        fprintf('Lasso OK.\n');
        
        Temp = thetaList<0; %Negative theta!!!
        %Temp = thetaList~=0; %Not Zero!!!
        EstimatedSparsity = squeeze((sum(sum(Temp,1),2)-NodeNumber)/NodeNumber/(NodeNumber-1));
        

        LassoLambdaSetTemp=zeros(1,length(SparsityRange));
        LassoSparsitySetTemp=zeros(1,length(SparsityRange));
        ThetaMatrixSetTemp=zeros(1,length(SparsityRange),NodeNumber,NodeNumber);
                
        for iSparsity=1:length(SparsityRange)
            

            
            DeltaSparsity = EstimatedSparsity - SparsityRange(iSparsity);
            [MinDelta,IndexSparsity] = min(abs(DeltaSparsity));
            
            LassoLambdaSetTemp(1,iSparsity) = LassoLambdaTestSet(IndexSparsity(1));
            LassoSparsitySetTemp(1,iSparsity) = EstimatedSparsity(IndexSparsity(1));
            
            
            Matrix = thetaList(:,:,IndexSparsity(1));
            
            ThetaMatrixSetTemp(1,iSparsity,:,:) = Matrix;
            
            
        end
        
        LassoLambdaSet(i,:) = LassoLambdaSetTemp;
        LassoSparsitySet(i,:) = LassoSparsitySetTemp;
        ThetaMatrixSet(i,:,:,:) = ThetaMatrixSetTemp;
        
    end
    
    save([OutputDir,filesep,'GTA_BinarizedUndirected_LassoLambdaSet.mat'],'LassoLambdaSet','LassoSparsitySet','ThetaMatrixSet');
    
end




%Perform graph theoretical analysis  -- Using Negative sign, correspondes to using positive correlation

FromUpDir = '/AtlasROI_SeedSignal_ExtractedInMNISpace_Detrend_Filter01_1';
ToUpDir = '/SparsityThreshold';

ConditionList={'NoGSR','GSR'};

NodeNumber=90;
SparsityRange = [0.02:0.02:0.5]';
%LassoLambdaTestSet = [0.005:0.001:1];
LassoLambdaTestSet = [0.0001:0.0001:0.001, 0.002:0.001:1];


for iCondition=1:length(ConditionList)
    
    %FisherZDir = [FromUpDir,filesep,ConditionList{iCondition}];
    OutputDir = [ToUpDir,filesep,ConditionList{iCondition},'_GLasso_NegativeTheta'];
    mkdir(OutputDir);
   

    
    %Perform graph theoretical analysis (binarized undirected)
   
    nSparsity = length(SparsityRange);
    
    MaxSparsitySet = ones(AutoDataProcessParameter.SubjectNum,1);
    
    CpSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity);
    LpSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity);
    GammaSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity);
    LambdaSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity);
    SigmaSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity);
    ElocSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity);
    EglobSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity);
    AssortativitySet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity);
    ModularitySet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity);
    
    DegreeSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity,NodeNumber);
    NodalEfficiencySet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity,NodeNumber);
    BetweennessSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity,NodeNumber);
    ClusteringCoefficientSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity,NodeNumber);
    ParticipantCoefficientSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity,NodeNumber);
    SubgraphCentralitySet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity,NodeNumber);
    EigenvectorCentralitySet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity,NodeNumber);
    PageRankCentralitySet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity,NodeNumber);

    LassoSparsitySet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity);
    LassoLambdaSet = zeros(AutoDataProcessParameter.SubjectNum,nSparsity);

    %[Cp, Lp, Gamma, Lambda, Sigma, Eloc, Eglob, Assortativity, NO!!Hierarchy, NO!!Synchronization, Modularity]
    %[Degree, NodalEfficiency, Betweenness, ClusteringCoefficient, ParticipantCoefficient, NO!!NormalizedParticipantCoefficient, SubgraphCentrality, EigenvectorCentrality, PageRankCentrality]
    
    for i=1:AutoDataProcessParameter.SubjectNum
         %MAT = load ([FisherZDir,filesep,'ROICorrelation_FisherZ_ROITC',AutoDataProcessParameter.SubjectID{i},'.mat']);
        load ([OutputDir,filesep,'GTA_BinarizedUndirected_LassoLambdaSet.mat']);
        
  
                
        parfor iSparsity=1:length(SparsityRange)
            
            fprintf('Sparsity: %d for Subject %s',SparsityRange(iSparsity),AutoDataProcessParameter.SubjectID{i});
            
            
            
            Matrix = squeeze(ThetaMatrixSet(i,iSparsity,:,:));
            

            Matrix(find(eye(size(Matrix))))=0;
                        
            Matrix=double(Matrix<0); % Binarized %USE negative theta!!!
            %Matrix=double(Matrix~=0); % Binarized %Not Zero!!!
            %Matrix=Matrix.*(Matrix>MatrixValueThreshold);
            
            [GTA] = y_GraphTheoreticalAnalysis_bu(Matrix,100);  %[GTA] = y_GraphTheoreticalAnalysis_bu(Matrix,100);
            
            CpSet(i,iSparsity) = GTA.Cp;
            LpSet(i,iSparsity) = GTA.Lp;

            GammaSet(i,iSparsity) = GTA.Gamma;
            LambdaSet(i,iSparsity) = GTA.Lambda;
            SigmaSet(i,iSparsity) = GTA.Sigma;
            
            
            ElocSet(i,iSparsity) = GTA.Eloc;
            EglobSet(i,iSparsity) = GTA.Eglob;
            AssortativitySet(i,iSparsity) = GTA.Assortativity;
            ModularitySet(i,iSparsity) = GTA.Modularity;
            
            DegreeSet(i,iSparsity,:) = GTA.Degree;
            NodalEfficiencySet(i,iSparsity,:) = GTA.NodalEfficiency;
            BetweennessSet(i,iSparsity,:) = GTA.Betweenness;
            ClusteringCoefficientSet(i,iSparsity,:) = GTA.ClusteringCoefficient;
            ParticipantCoefficientSet(i,iSparsity,:) = GTA.ParticipantCoefficient;
            SubgraphCentralitySet(i,iSparsity,:) = GTA.SubgraphCentrality;
            EigenvectorCentralitySet(i,iSparsity,:) = GTA.EigenvectorCentrality;
            PageRankCentralitySet(i,iSparsity,:) = GTA.PageRankCentrality;    
        end
    end
    save([OutputDir,filesep,'GTA_BinarizedUndirected.mat'],'CpSet','LpSet','GammaSet','LambdaSet','SigmaSet','ElocSet','EglobSet','AssortativitySet','ModularitySet','DegreeSet','NodalEfficiencySet','BetweennessSet','ClusteringCoefficientSet','ParticipantCoefficientSet','SubgraphCentralitySet','EigenvectorCentralitySet','PageRankCentralitySet','MaxSparsitySet','SparsityRange','LassoSparsitySet','LassoLambdaSet'); 
end



