%Simulation study given in
% Chung, M.K., Lee, H. Ombao. H., Solo, V. 2019 Exact topological inference 
% of the resting-state brain networks in twins, Network Neuroscience, in
% press.  http://www.stat.wisc.edu/~mchung/papers/chung.2018.NN.pdf
%
% This code is downloaded from 
% http://www.cs.wisc.edu/~mchung/twins or 
% http://www.cs.wisc.edu/~mchung/TDA 
%
% If you are using the code, please reference the above paper. If we iterate 
% the below procedure, we obtain the results given in the paper. Due to randomness, 
% results may change. 
%
%
% (C) 2018 Moo K. Chung, Universtiy of Wisconsin-Madison   
%     Hyekyoung Lee, Seoul National University
%     mkchung@wisc.edu
%
%
% Tested in iMAC 2012 (macOS 10.12.5) with MATLAB R2016a
%
% Update history
% 2018 June 6 Error related to SLM1=PH_SLM(1-C1) fixed
% 2018 August 8
% 2019 Febuary 15 Codes were run with 5000 nodes. 


%----------------------------------------------
%PERMUTATION TEST 
%Generage exactly 252 permutations

n=5; m=5; %There are two groups and each group has 3 networks.
p=120; % 20 to 500 nodes in a network
sigma =0.1; %sigma controls amount of noise. 0.01 and 0.1 are used.
%
pL1     = [];
pL2     = [];
pLinf   = [];
pGH     = [];
pKS0    = [];
pKS1    = [];
pQ      = [];

for nsim=1:20
    
    nsim
    %C1 and C2 are randomly simulated connectivity matrices that will be feed
    %into network distance functions.
    
    tic
    %-----------------------
    %4 vs 4 modules
    A = normrnd(0, 1, n,p);
    [X C1] = RG_module(A,10,sigma);
    [Y C2] = RG_module(A,10,sigma);
    
    %------------------------
    %Connectivity matrix = sqrt(1-corr), which is a metric.
    %Note 0<= sqrt(1-corr)<=sqrt(2)
    %
    %Numerically it is possible to have imaginary number like 0 - i 0.000000001.
    %Taking abs(real( )) fix this numerical problem related to sqrt.
    
    C1= abs(real(sqrt(1-C1)));  %figure; imagesc(C1); colorbar
    C2 = abs(real(sqrt(1-C2))); %figure; imagesc(C2); colorbar
    
    %---------------------------
    %pvalue computation for KS-distace
    [beta0C1,beta1C1,biggest0C1,biggest1C1] = PH_betti3(C1, [0:0.01:sqrt(2)]);
    [beta0C2,beta1C2,biggest0C2,biggest1C2] = PH_betti3(C2, [0:0.01:sqrt(2)]);
    
    %Betti-0 plot
    d0 =max(abs(beta0C1-beta0C2));
    d0 = ceil(d0);
    q0=size(C1,1);
    pvalue0  = PH_Aasym(q0,d0);
    %pvalue0  = PH_A(q0,d0);
    pKS0=[pKS0 pvalue0];
    
    %Betti-1 plot
    %q (number of features) is given as follows:
    %beta1 = beta0 + # edges - # of nodes
    %6555 = 1 + 6670 - 116 for 116 node network
    d1 =max(abs(beta1C1-beta1C2));
    d1 = ceil(d1);
    q1= 1 + (q0^2-q0)/2 - q0;
    pvalue1  = PH_Aasym(q1,d1);
    %pvalue1  = PH_A(q1,d1);
    pKS1=[pKS1 pvalue1];
    
    %KS-distance on the size of the largest cluster: gamma
    %This is shown in Chung et al. (2017) but not in this study.
    %d0= max(abs(C1biggest-C2biggest));
    %d0 = ceil(d0); q=length(C1biggest);
    %pvalue = PH_Aasym(q,d0)
    
    %-----------------------------
    %pvalue computation for GH-distance and L1, L2 and L_infity distances
    %based on the permutation test.
    %WARNING: hklee edited from SLM1=PH_SLM(1-C1); SLM2=PH_SLM(1-C2)
    
    
    SLM1=PH_SLM2(C1); %Single linkage matrix for C1
    SLM2=PH_SLM2(C2); %Signle linkage matrix for C2
    GH_observed= max(max(abs(SLM1-SLM2))); %observed values for GH-distance
    
    %observed values for L1, L2 and Linfty distances
    L1_observed = sum(sum(abs(C1-C2)));
    L2_observed = sqrt(sum(sum(abs (C1-C2).^2)));
    Linf_observed = max(max(abs(C1-C2)));
    
    %observed value for bottleneck distance
    %bottle_observed = PH_bottleneck2(C1,C2);
    
    %-----------------------------
    %Observed Q-modularity
    [Ci Q1] = modularity_und(C1);
    [Ci Q2] = modularity_und(C2);
    Qobs = abs(Q1-Q2);
    
    toc
    %-----------------------------
    %Setting for the permutation test. The permutation test is applied to
    %GH, L1, L2 and L_infity distances.
    %For each random simulation, we genearte every
    % possible (n+m choose n) permuations.
    % For n=m=5, there are exactly 252 permutations.
    
    toc
    
    totalperm = nchoosek(n+m,n);  %Total number of permutation for sample sizes n and m
    
    a=1:(n+m); %node index
    perm1 = nchoosek(a,n);
    
    perm2=[];
    for i=1:totalperm
        perm2(i,:)=setdiff(a,perm1(i,:));  %[perm1 perm2] gives the every possible permutations
    end
    
    XYcombined=[X;Y]; %combined data matrix
    
    GH=[];L1=[]; L2=[]; Linf=[]; bottle=[]; Qdiff=[];

    
    for i=1:totalperm
        
        Xperm = XYcombined(perm1(i,:),:);  %permuted data matrix for the first group
        Yperm = XYcombined(perm2(i,:),:);  %permuted data matrix for the second group
        
        C1=corr(Xperm);
        C1= abs(real(sqrt(1-C1)));
        
        C2=corr(Yperm);
        C2 = abs(real(sqrt(1-C2)));
        
        %-----------
        %Matrix norms
        L1(i) = sum(sum(abs(C1-C2))); %L1 distance
        L2(i) = sqrt(sum(sum(abs (C1-C2).^2))); %L2 distance
        Linf(i) = max(max(abs(C1-C2))); %L_infinity distance
        
        %----------
        %GH-distance
        SLM1=PH_SLM2(C1); %Sinle linkage matrices
        SLM2=PH_SLM2(C2);
        GH(i)= max(max(abs(SLM1-SLM2))); %GH-distance
        
        %-----------
        %Q-modularity
        [Ci Q1] = modularity_und(C1);
        [Ci Q2] = modularity_und(C2);
        Qdiff(i) = abs(Q1-Q2);
        
    end
    
    % the fraction of GH that is bigger than the observation is the pvalue
    pvalue_L1 = length(find(L1>L1_observed))/252; %pvalue for L1 distance
    pvalue_L2 = length(find(L2>L2_observed))/252; %pvalue for L2 distance
    pvalue_Linf = length(find(Linf>Linf_observed))/252; %pvalue for L_infity distance
    pvalue_GH = length(find(GH>GH_observed))/252; %pvalue for bottle neck distance
    pvalue_Q = length(find(Qdiff>Qobs))/252; %pvalue for Q-modularity
         
    pL1     = [pL1 pvalue_L1];
    pL2     = [pL2 pvalue_L2];
    pLinf   = [pLinf pvalue_Linf];
    pGH     = [pGH pvalue_GH];
    pQ =[pQ pvalue_Q];
    
    
end

%average p-values
[mean(pL1) std(pL1)]   
[mean(pL2) std(pL2)]
[mean(pLinf) std(pLinf)]  
[mean(pGH) std(pGH)]
[mean(pKS0) std(pKS0)]
[mean(pKS1) std(pKS1)] 
[mean(pQ) std(pQ)]

%error rate: false positives. Use it testing no network difference
[length(find(pL1<0.05))/100]   
[length(find(pL2<0.05))/100]
[length(find(pLinf<0.05))/100]  
[length(find(pGH<0.05))/100 ]
[length(find(pKS0<0.05))/100]
[length(find(pKS1<0.05))/100]
[length(find(pQ<0.05))/100]
    
%error rate: false negatives. Use it in testing network differece
[length(find(pL1>0.05))/100]   
[length(find(pL2>0.05))/100]
[length(find(pLinf>0.05))/100]  
[length(find(pGH>0.05))/100 ]
[length(find(pKS0>0.05))/100]
[length(find(pKS1>0.05))/100]
[length(find(pQ>0.05))/100]

%---------------
% VISUALIZATION

% Displaying modular structures. 

n=5; m=5; %There are two groups and each group has 5 networks.
p=20; % 100 nodes in a network
sigma =0.1; %sigma controls amount of noise. 

A = normrnd(0, 1, n,p);
[X C1] = RG_module(A,2,sigma);
[Y C2] = RG_module(A,10,sigma);

figure; imagesc(C1); colorbar; colormap('hot'); axis square; caxis([-1 1])
title('2 moduels'); axis off;  
figure; imagesc(C2); colorbar; colormap('hot'); axis square; caxis([-1 1])
title('10 modules'); axis off;

 
%Betti plots 
C1= abs(real(sqrt(1-C1))); C2 = abs(real(sqrt(1-C2)));
[beta0C1,beta1C1,biggest0C1,biggest1C1] = PH_betti2(C1, [0:0.01:sqrt(2)]);
[beta0C2,beta1C2,biggest0C2,biggest1C2] = PH_betti2(C2, [0:0.01:sqrt(2)]);

%Displaying monotonic feature functions (0-th Betti number)
figure; plot(0:0.01:sqrt(2), beta0C1,'--r','LineWidth',2);
hold on; plot(0:0.01:sqrt(2), beta0C2,'-k','LineWidth',2);
xlabel('\lambda'); ylabel('B(\lambda)')
legend('boxoff')

%Displaying monotonic feature functions (1-st Betti number)
figure; plot(0:0.01:sqrt(2), beta1C1,'--r','LineWidth',2);
hold on; plot(0:0.01:sqrt(2), beta1C2,'-k','LineWidth',2);
xlabel('\lambda'); ylabel('B(\lambda)')
legend('boxoff')





