add backbone of graph theory features

features/find_network_properties3.m

@@ -0,0 +1,111 @@
+function find_network_properties3
+
+%   This function defines a network like Joon's paper, but normalizes
+%   properties according to random networks
+
+samp_freq = 250;
+network_thresh = 0.05; %vary to see how stable our results are
+win = 10;   % number of seconds of EEG window
+% total_length = 300;    % total number of seconds of EEG epoch
+
+
+for subject = 5
+
+
+     Larray = zeros(1,floor(total_length/win)); %path length
+     Carray = zeros(1,floor(total_length/win)); %clustering coefficient
+     geffarray = zeros(1,floor(total_length/win)); %global efficiency
+     bswarray = zeros(1,floor(total_length/win)); %small worldness
+     Qarray = zeros(1,floor(total_length/win)); %modularity
+
+
+    for bp = 4
+
+
+        for state = 1:3
+
+
+            state
+%           EEG = pop_loadset('filename', [sname statename '.set'],'filepath',['F:\McDonnell Foundation study\University of Michigan\Anesthesia\' sname '\Resting state analysis']);
+            EEG = pop_loadset('filename', [sname statename '.set'],'filepath','C:\Users\Danielle\OneDrive - McGill University\Research\BIAPT Lab\DOC\Motif paper\WSAS09\DATA\5 min segments\');
+
+
+            [dataset, com, b] = pop_eegfiltnew(EEG, lp, hp);    
+            filt_data = dataset.data';
+
+            b_charpath = zeros(1,floor(total_length/win));
+            b_clustering = zeros(1,floor(total_length/win));
+            b_geff = zeros(1,floor(total_length/win));
+            bsw = zeros(1,floor(total_length/win));
+            Q = zeros(1,floor(total_length/win));
+
+            for i = 1:(floor((length(filt_data))/(win*samp_freq)))
+
+%                 EEG_seg = filt_data((i-1)*win*samp_freq + 1:i*win*samp_freq, EEG_chan);      % Only take win seconds length from channels that actually have EEG
+                EEG_seg = filt_data((i-1)*win*samp_freq + 1:i*win*samp_freq, :);
+
+                PLI = w_PhaseLagIndex(EEG_seg); %weighted PLI
+
+                A = sort(PLI);
+                B = sort(A(:));
+                C = B(1:length(B)-length(EEG_chan)); % Remove the 1.0 values from B (correlation of channels to themselves)
+
+                index = floor(length(C)*(1-network_thresh)); %top network_thresh% of data
+                thresh = C(index);  % Values below which the graph will be assigned 0, above which, graph will be assigned 1
+
+
+            % Create a (undirected, unweighted) network based on top network_thresh% of PLI connections    
+            for m = 1:length(PLI)
+                for n = 1:length(PLI)
+                    if (m == n)
+                        b_mat(m,n) = 0;
+                    else
+                        if (PLI(m,n) > thresh)
+                            b_mat(m,n) = 1;
+                        else
+                            b_mat(m,n) = 0;
+                        end
+                    end
+                end
+            end          
+
+                % Find average path length
+
+                D = distance_bin(b_mat);
+                [b_lambda,geff,~,~,~] = charpath(D,0,0);   % binary charpath
+                [W0,R] = null_model_und_sign(b_mat,10,0.1);    % generate random matrix
+
+                  % Find clustering coefficient
+
+                C = clustering_coef_bu(b_mat);  
+
+                % Find properties for random network
+
+                [rlambda,rgeff,~,~,~] = charpath(distance_bin(W0),0,0);   % charpath for random network
+                rC = clustering_coef_bu(W0); % cc for random network
+
+                b_clustering(i) = nanmean(C)/nanmean(rC); % binary clustering coefficient
+                b_charpath(i) = b_lambda/rlambda;  % charpath
+                b_geff(i) = geff/rgeff; % global efficiency
+
+                bsw(i) = b_clustering/b_charpath; % binary smallworldness
+
+                [M,modular] = community_louvain(b_mat,1); % community, modularity
+                Q(i) = modular;
+
+            end
+
+             Larray(state,:) = b_charpath(1,1:floor(total_length/win)); 
+             Carray(state,:) = b_clustering;
+             geffarray(state,:) = b_geff;
+             bswarray(state,:) = bsw;
+             Qarray(state,:) = Q;
+
+        end
+
+    end
+
+end
+
+ %figure; plot(Lnorm)
+

features/graph_theory/binary_small_worldness.m

@@ -0,0 +1,7 @@
+function [outputArg1,outputArg2] = binary_small_worldness(inputArg1,inputArg2)
+%BINARY_SMALL_WORLDNESS Summary of this function goes here
+%   Detailed explanation goes here
+outputArg1 = inputArg1;
+outputArg2 = inputArg2;
+end
+

features/graph_theory/clustering_coefficient.m

@@ -0,0 +1,7 @@
+function [outputArg1,outputArg2] = clustering_coefficient(inputArg1,inputArg2)
+%CLUSTERING_COEFFICIENT Summary of this function goes here
+%   Detailed explanation goes here
+outputArg1 = inputArg1;
+outputArg2 = inputArg2;
+end
+

features/graph_theory/modularity.m

@@ -0,0 +1,7 @@
+function [outputArg1,outputArg2] = modularity(inputArg1,inputArg2)
+%MODULARITY Summary of this function goes here
+%   Detailed explanation goes here
+outputArg1 = inputArg1;
+outputArg2 = inputArg2;
+end
+

features/graph_theory/undirected_global_efficiency.m

@@ -0,0 +1,43 @@
+function [g_efficiency,norm_g_efficiency,avg_path_length,norm_avg_path_length] = undirected_global_efficiency(matrix,t_level,number_null_network, bin_swaps, weight_frequency)
+%UNDIRECTED GLOBAL EFFICIENCY will calculate the undirected global
+%efficiency and path length for connectivity matrices like wpli
+%   matrix:
+%   t_level:
+%   number_null_network:
+%   bin_swaps:
+%   weight_frequency
+    %% Binarize and Threshold the matrix
+    b_matrix = binarize_matrix(threshold_matrix(matrix,t_level));
+
+    %% Create random null network
+    [num_row, num_col] = size(b_matrix); % These should be the same
+    null_matrices = zeros(number_null_network,num_row, num_col);
+    for i = 1:number_null_network
+        disp(i)
+        [null_matrix,~] = null_model_und_sign(b_matrix,bin_swaps,weight_frequency);    % generate random matrix
+        null_matrices(i,:,:) = null_matrix; % store all null matrix
+    end
+
+    %% Calculate the characteristic path length
+    input_distance = distance_bin(b_matrix);
+    [avg_path_length,g_efficiency,~,~,~] = charpath(input_distance,0,0);   % binary charpath
+
+    %% Calculate the characteristic path length for each null_matrix and average them
+    null_matrices_avg_path_length = zeros(1,number_null_network);
+    null_matrices_g_efficiency = zeros(1,number_null_network);
+    for i = 1:number_null_network
+        null_b_matrix = null_matrices(i,:,:);
+        null_input_distance = distance_bin(null_b_matrix);
+        [null_matrix_avg_path_length,null_matrix_g_efficiency,~,~,~] = charpath(null_input_distance,0,0);   % binary charpath    
+        null_matrices_avg_path_length(i) = null_matrix_avg_path_length;
+        null_matrices_g_efficiency(i) = null_matrix_g_efficiency;
+    end
+    % Calculate mean null path length and mean global efficiency
+    null_avg_path_length = mean(null_matrices_avg_path_length);
+    null_g_efficiency = mean(null_matrices_g_efficiency);
+
+    %% Normalizing average path length and global effiency
+    norm_g_efficiency = g_effiency/null_g_efficiency;
+    norm_avg_path_length = avg_path_length/null_avg_path_length;
+end
+

utils/binarize_matrix.m

@@ -0,0 +1,7 @@
+function [b_matrix] = binarize_matrix(matrix)
+%BINARIZE_MATRIX set the value of the matrix to 0 or 1
+%   matrix: a N*N matrix
+    b_matrix = matrix;
+    b_matrix(b_matrix > 0) = 1;
+end
+

utils/threshold_matrix.m

@@ -0,0 +1,25 @@
+function [t_matrix] = threshold_matrix(matrix,t_level)
+%THRESHOLD_MATRIX Threshold a matrix to have element below a significant
+%amount to 0
+%   Matrix: N*N matrix with value within any range
+%   t_level: value from 0 to 1 which set the ratio of highest value
+%   to keep. i.e. t_level = 0.05 -> keep only top 5% value
+
+    %% Flatten the matrix into an array
+    [num_row, num_col] = size(matrix);
+    num_element = num_row*num_col;
+    array = reshape(matrix, [1 num_element]);
+
+    %% Find the value to threshold on (the one at the limit of top t_level%
+    sorted_array = sort(array);
+    t_index = floor(num_element*(1 - t_level)) + 1;
+    t_element = sorted_array(t_index);
+
+    %% Threshold the matrix
+    t_matrix = matrix;
+    t_matrix(t_matrix < t_element) = 0;
+
+    %% Remove the diagonal elements
+    t_matrix = t_matrix - diag(diag(t_matrix));
+end
+