clear; clc;

cd /home/k1514724/subhadip/BBC/DTI_rsfMRI_ConnectivityMatrices/connectivity_matrices/
% Suppy the following: N, n_perm, files_rmfmri and files_DTI
%Parcellation resolution of rsfMRI and DTI network/connectivity matrices. It is assumes that both of them have same parcellation resolutions 
N=90;

% Total number of permutations
n_reps=10000;

% load rsfMRI and DTI network matrices
files_rsfmri=dir('Netmat*.txt');
files_DTI=dir('Connection*.txt');

% Computes significant correlation between DTI and rsFMRI connectivity
% networks within a population using spearmann's correlation 

Netmatrix_rsfMRI(1:N,1:N,1:length(files_rsfmri))=zeros;
Netmatrix_DTI(1:N,1:N,1:length(files_DTI))=zeros;

for l=1:length(files_rsfmri)
Netmatrix_rsfMRI(:,:,l)=load(files_rsfmri(l).name);
Netmatrix_DTI(:,:,l)=load(files_DTI(l).name);
end

%% 
c_observed(1:N,1:N)=zeros; % observed correlation coefficient
pval_permuted(1:N,1:N)=zeros; % p-values obtained after permutation tests

for i=1:N
    for j=i+1:N
        [a,b,rnk_x]=unique(Netmatrix_rsfMRI(i,j,:)); % coumputing the rank of rsfMRI connectivity values
        [a,b,rnk_y]=unique(Netmatrix_DTI(i,j,:));    % computing the rank of DTI connectivity values
        
        % converting the rnk_x and rnk_y into column vectors if they are row vectors
        size_x=size(rnk_x);size_y=size(rnk_y);
        if size_x(1,2) & size_y(1,2) >1
            rnk_x=rnk_x'; rnk_y=rnk_y';
        else
        end
        X=[rnk_x rnk_y];
        [temp_c,temp_p]=corr(X,'type','Spearman','tail','both');
        c_observed(i,j)=temp_c(1,2);
        clear temp_c; clear temp_p; clear X;
        
        %% randomly permute the values of rnk_x
        corrcoef_permuted(1:n_reps,1)=zeros;
        for m=1:n_reps
            perm_index_x=transpose(randperm(length(files_rsfmri))); % randomly permuted array indices
            Y=[rnk_x(perm_index_x) rnk_y]; % just permute/reshuffle the values of one array
            
            [c_temp,p_temp]=corr(Y,'type','Spearman','tail','both');
            corrcoef_permuted(m,1)=c_temp(1,2);
            clear p_temp; clear c_temp;
        end
        clear Y; clear clear rnk_x; clear rnk_y;
        % determine how many permuted correlation coefficients exceeds the
        % value of the observed correlation coefficient
        pval_permuted(i,j)=sum(corrcoef_permuted>c_observed(i,j))/n_reps;
        clear corrcoef_permuted;
    end
end

%% Now perform FDR testing on those pval_permuted