From 8b62b1bcde22be437d5d5f882c6c8267f25e0de5 Mon Sep 17 00:00:00 2001 From: Bl0101 Date: Sat, 5 Oct 2019 09:33:15 -0400 Subject: [PATCH] add backbone of graph theory features --- {src => features}/dpli.m | 0 features/find_network_properties3.m | 111 ++++++++++++++++++ .../graph_theory/binary_small_worldness.m | 7 ++ .../graph_theory/clustering_coefficient.m | 7 ++ features/graph_theory/modularity.m | 7 ++ .../undirected_global_efficiency.m | 43 +++++++ {src => features}/hub_location.m | 0 {src => features}/permutation_entropy.m | 0 {src => features}/phase_amplitude_coupling.m | 0 {src => features}/spectral_power.m | 0 {src => features}/topographic_distribution.m | 0 {src => features}/wpli.m | 0 {features => src}/na_dpli.m | 0 {features => src}/na_hub_location.m | 0 {features => src}/na_permutation_entropy.m | 0 .../na_phase_amplitude_coupling.m | 0 {features => src}/na_spectral_power_ratio.m | 0 .../na_topographic_distribution.m | 0 {features => src}/na_wpli.m | 0 utils/binarize_matrix.m | 7 ++ utils/threshold_matrix.m | 25 ++++ 21 files changed, 207 insertions(+) rename {src => features}/dpli.m (100%) create mode 100644 features/find_network_properties3.m create mode 100644 features/graph_theory/binary_small_worldness.m create mode 100644 features/graph_theory/clustering_coefficient.m create mode 100644 features/graph_theory/modularity.m create mode 100644 features/graph_theory/undirected_global_efficiency.m rename {src => features}/hub_location.m (100%) rename {src => features}/permutation_entropy.m (100%) rename {src => features}/phase_amplitude_coupling.m (100%) rename {src => features}/spectral_power.m (100%) rename {src => features}/topographic_distribution.m (100%) rename {src => features}/wpli.m (100%) rename {features => src}/na_dpli.m (100%) rename {features => src}/na_hub_location.m (100%) rename {features => src}/na_permutation_entropy.m (100%) rename {features => src}/na_phase_amplitude_coupling.m (100%) rename {features => src}/na_spectral_power_ratio.m (100%) rename {features => src}/na_topographic_distribution.m (100%) rename {features => src}/na_wpli.m (100%) create mode 100644 utils/binarize_matrix.m create mode 100644 utils/threshold_matrix.m diff --git a/src/dpli.m b/features/dpli.m similarity index 100% rename from src/dpli.m rename to features/dpli.m diff --git a/features/find_network_properties3.m b/features/find_network_properties3.m new file mode 100644 index 00000000..84e03586 --- /dev/null +++ b/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) + diff --git a/features/graph_theory/binary_small_worldness.m b/features/graph_theory/binary_small_worldness.m new file mode 100644 index 00000000..2de4187c --- /dev/null +++ b/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 + diff --git a/features/graph_theory/clustering_coefficient.m b/features/graph_theory/clustering_coefficient.m new file mode 100644 index 00000000..c58ebe51 --- /dev/null +++ b/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 + diff --git a/features/graph_theory/modularity.m b/features/graph_theory/modularity.m new file mode 100644 index 00000000..57fe6df0 --- /dev/null +++ b/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 + diff --git a/features/graph_theory/undirected_global_efficiency.m b/features/graph_theory/undirected_global_efficiency.m new file mode 100644 index 00000000..92978a72 --- /dev/null +++ b/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 + diff --git a/src/hub_location.m b/features/hub_location.m similarity index 100% rename from src/hub_location.m rename to features/hub_location.m diff --git a/src/permutation_entropy.m b/features/permutation_entropy.m similarity index 100% rename from src/permutation_entropy.m rename to features/permutation_entropy.m diff --git a/src/phase_amplitude_coupling.m b/features/phase_amplitude_coupling.m similarity index 100% rename from src/phase_amplitude_coupling.m rename to features/phase_amplitude_coupling.m diff --git a/src/spectral_power.m b/features/spectral_power.m similarity index 100% rename from src/spectral_power.m rename to features/spectral_power.m diff --git a/src/topographic_distribution.m b/features/topographic_distribution.m similarity index 100% rename from src/topographic_distribution.m rename to features/topographic_distribution.m diff --git a/src/wpli.m b/features/wpli.m similarity index 100% rename from src/wpli.m rename to features/wpli.m diff --git a/features/na_dpli.m b/src/na_dpli.m similarity index 100% rename from features/na_dpli.m rename to src/na_dpli.m diff --git a/features/na_hub_location.m b/src/na_hub_location.m similarity index 100% rename from features/na_hub_location.m rename to src/na_hub_location.m diff --git a/features/na_permutation_entropy.m b/src/na_permutation_entropy.m similarity index 100% rename from features/na_permutation_entropy.m rename to src/na_permutation_entropy.m diff --git a/features/na_phase_amplitude_coupling.m b/src/na_phase_amplitude_coupling.m similarity index 100% rename from features/na_phase_amplitude_coupling.m rename to src/na_phase_amplitude_coupling.m diff --git a/features/na_spectral_power_ratio.m b/src/na_spectral_power_ratio.m similarity index 100% rename from features/na_spectral_power_ratio.m rename to src/na_spectral_power_ratio.m diff --git a/features/na_topographic_distribution.m b/src/na_topographic_distribution.m similarity index 100% rename from features/na_topographic_distribution.m rename to src/na_topographic_distribution.m diff --git a/features/na_wpli.m b/src/na_wpli.m similarity index 100% rename from features/na_wpli.m rename to src/na_wpli.m diff --git a/utils/binarize_matrix.m b/utils/binarize_matrix.m new file mode 100644 index 00000000..d368bb17 --- /dev/null +++ b/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 + diff --git a/utils/threshold_matrix.m b/utils/threshold_matrix.m new file mode 100644 index 00000000..2e76c2c6 --- /dev/null +++ b/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 +