%Notes for Jesse's Graph Analysis Workshop
% put this .m file inside NITP_network_lab_2012
clear all;
%Set binarized threshold for DTI graphs, fMRI graphs
THRESH_DTI=.2;
THRESH_FMRI=.2;

% Load in matrices
load('dti_matrix.txt');
load('fmri_matrix.txt');

%% Binarize matrices, thresholding to keep strongest % of connections
dti_matrix_bin=adjmatrix(dti_matrix,THRESH_DTI,false);
fmri_matrix_bin=adjmatrix(fmri_matrix,THRESH_FMRI,false);

%% Modules
% Calculate modules in DTI graph using Louvain algorithm
x_dti=0;
i_dti=0;
[Comm_index_dti	Q]=modularity_louvain_und(dti_matrix_bin);		
[x_dti i_dti]=sort(Comm_index_dti');	
dti_matrix_bin_sort=dti_matrix_bin(i_dti,i_dti);

% Calculate modules in fmri graph using Louvain algorithm
x_fmri=0;
i_fmri=0;
[Comm_index_fmri Q]=modularity_louvain_und(fmri_matrix_bin);		
[x_fmri i_fmri]=sort(Comm_index_fmri');	
fmri_matrix_bin_sort=fmri_matrix_bin(i_fmri,i_fmri);

% Get components
[comps_dti,comp_sizes_dti]=get_components(dti_matrix_bin);
[comps_fmri,comp_sizes_fmri]=get_components(fmri_matrix_bin);


%% Regional Measures
% Regional DTI Measures
d_dti=degrees_und(dti_matrix_bin);	
cc_dti=clustering_coef_bu(dti_matrix_bin);	
bc_dti=betweenness_bin(dti_matrix_bin);	
P_dti=participation_coef(dti_matrix_bin,Comm_index_dti);	

%region_sorted_ranks for DTI
sorted_ranks_dti=region_sorted_ranks(d_dti,cc_dti,bc_dti,P_dti);

% Regional fMRI Measures
d_fmri=degrees_und(fmri_matrix_bin);	
cc_fmri=clustering_coef_bu(fmri_matrix_bin);	
bc_fmri=betweenness_bin(fmri_matrix_bin);	
P_fmri=participation_coef(fmri_matrix_bin,Comm_index_fmri);	

%region_sorted_ranks for fMRI
sorted_ranks_fmri=region_sorted_ranks(d_fmri,cc_fmri,bc_fmri,P_fmri);

%% Global Measures
% DTI
% Global DTI Measures
dti_matrix_bin_distmat=distance_bin(dti_matrix_bin);
[cpl_dti,eglob_dti]=charpath(dti_matrix_bin_distmat);
mcc_dti=mean(clustering_coef_bu(dti_matrix_bin));

% Random DTI network
sim_dti_matrix=randmio_und(dti_matrix_bin,10);	
sim_dti_distmat=distance_bin(sim_dti_matrix);

[Comm_index_sim_dti	Q]=modularity_louvain_und(sim_dti_matrix);

% Random DTI network global measures
[sim_cpl_dti,sim_eglob_dti]=charpath(sim_dti_distmat);	
sim_mcc_dti=mean(clustering_coef_bu(sim_dti_matrix));

% Random DTI network Regional measures
d_Sim_dti=degrees_und(sim_dti_matrix);	
cc_sim_dti=clustering_coef_bu(sim_dti_matrix);	
bc_sim_dti=betweenness_bin(sim_dti_matrix);	
P_sim_dti=participation_coef(sim_dti_matrix,Comm_index_dti);

%Calculate Small worldness (sigma)
gamma_dti=mean(cc_dti./cc_sim_dti);
lambda_dti=mean(cpl_dti./sim_cpl_dti);
sigma_dti=gamma_dti/lambda_dti;

% fmri
% Global fmri Measures
fmri_matrix_bin_distmat=distance_bin(fmri_matrix_bin);
[cpl_fmri,eglob_fmri]=charpath(fmri_matrix_bin_distmat);
mcc_fmri=mean(clustering_coef_bu(fmri_matrix_bin));

% Random fmri network
sim_fmri_matrix=randmio_und(fmri_matrix_bin,10);	
sim_fmri_distmat=distance_bin(sim_fmri_matrix);

[Comm_index_sim_fmri Q]=modularity_louvain_und(sim_fmri_matrix);

% Random fmri network global measures
[sim_cpl_fmri,sim_eglob_fmri]=charpath(sim_fmri_distmat);	
sim_mcc_fmri=mean(clustering_coef_bu(sim_fmri_matrix));

% Random fmri network Regional measures
d_Sim_fmri=degrees_und(sim_fmri_matrix);	
cc_sim_fmri=clustering_coef_bu(sim_fmri_matrix);	
bc_sim_fmri=betweenness_bin(sim_fmri_matrix);	
P_sim_fmri=participation_coef(sim_fmri_matrix,Comm_index_fmri);

% Calculate Small worldness (sigma)
gamma_fmri=mean(cc_fmri./cc_sim_fmri); % cc_ratio
lambda_fmri=mean(cpl_fmri./sim_cpl_fmri); % cpl_ratio
sigma_fmri=gamma_fmri/lambda_fmri;

%% Display
% Display Graph for DTI
subplot(2,3,1), imagesc(dti_matrix), title('dti matrix weighted'); % original matrix
subplot(2,3,2), imagesc(dti_matrix_bin), title(['dti matrix ',num2str(THRESH_DTI),' binarized']); %binarized matrix
subplot(2,3,3), imagesc(dti_matrix_bin_sort), title('sorted by modularity');

% Display Graph for fmri
subplot(2,3,4), imagesc(fmri_matrix), title('fmri matrix weighted'); % original matrix
subplot(2,3,5), imagesc(fmri_matrix_bin), title(['fmri matrix ',num2str(THRESH_FMRI),' binarized']); %binarized matrix
subplot(2,3,6), imagesc(fmri_matrix_bin_sort), title('sorted by modularity');

dd=[mean(cpl_dti) mean(cpl_fmri) mean(sim_cpl_dti) mean(sim_cpl_fmri); mean(cc_dti) mean(cc_fmri) mean(cc_sim_dti) mean(cc_sim_fmri); sigma_dti sigma_fmri 0 0];

%figure;
f=figure;
set(f,'name','Global Measures','numbertitle','off')
%dd=rand(3,4); %# data
colnames = {'DTI' 'fMRI' 'Random DTI' 'Random fMRI'};
rownames = {'cpl' 'cc' 'sigma'};
t=uitable(f,'columnname',colnames,'rowname',rownames,'data',dd,'units','normalized','pos',[0 0 1 1]);
%t=uitable(f,'columnname',colnames,'rowname',rownames,'data',dd);

%% Random notes:
% print all components of a given module:
% com4_dti=find(ismember(Comm_index_dti,4));
% you can then paste this into a text file and cat it to run with
% fslview_regions.sh to visualize regions in module:
% sh ./fslview_regions.sh `cat com4.txt`

% If you wish to average the number of sigma values, you can modify the
% random network generation to iterate a number of times (1000)
% Save each sigma in a matrix then average over all.