Skip to content

Commit 23d7a7c

Browse files
fabriciocravofabriciocravo
committed
Maded TFCE method changes
1 parent e407c8d commit 23d7a7c

File tree

16 files changed

+320
-114
lines changed

16 files changed

+320
-114
lines changed

‎config_files/setparams.m‎

Lines changed: 9 additions & 10 deletions
Original file line numberDiff line numberDiff line change
@@ -29,9 +29,9 @@
2929

3030
% Datasets - Commented for easy use
3131
% Params.data_dir = './data/s_abcd_fc_rosenblatt.mat';
32-
% Params.data_dir = './data/s_hcp_fc_noble_tasks.mat';
33-
Params.data_dir = './data/s_hcp_act_noble_1.mat';
34-
Params.output = NaN;
32+
Params.data_dir = './data/s_hcp_fc_noble_tasks.mat';
33+
% Params.data_dir = './data/s_hcp_act_noble_1.mat';
34+
Params.output = 'test_tfce_comparions';
3535

3636
% Save specifications - if NaN output becomes dataset file name
3737
Params.save_directory = './power_calculator_results/';
@@ -57,9 +57,9 @@
5757

5858
%%% Resampling parameters %%%
5959
Params.parallel = true; % run stuff sequentially or in parallel
60-
Params.n_workers = 50; % num parallel workers for parfor, best if # workers = # cores
60+
Params.n_workers = 5; % num parallel workers for parfor, best if # workers = # cores
6161
Params.n_repetitions = 100; % 500 recommended
62-
Params.batch_size = 50;
62+
Params.batch_size = 25;
6363

6464

6565
%% Skip some tests - change ranges or the function
@@ -68,7 +68,7 @@
6868

6969

7070
%% List of subjects per subset
71-
Params.list_of_nsubset = {250, 500, 1000, 2000}; % To change this, add more when necessary
71+
Params.list_of_nsubset = {80}; % To change this, add more when necessary
7272
% size of subset is full group size (N=n*2 for two sample t-test or N=n for one-sample)
7373

7474
% Current model (see above design matrix) only designed for t-test
@@ -79,9 +79,8 @@
7979
Params.pthresh_second_level = 0.05; % FWER or FDR rate
8080
Params.tpr_dthresh = 0; % Threshold for true positives vs negatives
8181
Params.save_significance_thresh = 0.15;
82-
% Params.all_cluster_stat_types = {'Parametric', 'Size_cpp', 'Fast_TFCE_cpp', 'Constrained_cpp', 'Omnibus'};
83-
% Params.all_cluster_stat_types = {'Parametric', 'Size_cpp', 'Fast_TFCE_cpp', 'Constrained_cpp'};
84-
Params.all_cluster_stat_types = {'Parametric', 'Size_cpp', 'Fast_TFCE_cpp', 'Constrained_cpp', 'Omnibus_cNBS'};
82+
Params.all_cluster_stat_types = {'Exact_FC_TFCE', 'Fast_TFCE_dh1', 'Fast_TFCE_dh5', 'Fast_TFCE_dh10'};
83+
% Params.all_cluster_stat_types = {'Parametric', 'Size_cpp', 'Fast_TFCE_cpp', 'Constrained_cpp', 'Omnibus_cNBS'};
8584

8685
Params.all_submethods = {'FWER', 'FDR'};
8786

@@ -93,7 +92,7 @@
9392
%%%%% DEVELOPERS ONLY %%%%%
9493
% Use a small subset of permutations for faster development -- inappropriate for inference
9594

96-
Params.testing = false;
95+
Params.testing = true;
9796
Params.test_n_perms = 2;
9897
Params.test_n_repetitions = 5;
9998
Params.test_n_workers = 1;

‎create_flat_function.m‎

Lines changed: 14 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,14 @@
1+
function flat_function = create_flat_function(mask, varargin)
2+
% Parse variable_type like in unflatten_matrix
3+
p = inputParser;
4+
addParameter(p, 'variable_type', 'edge', @ischar);
5+
parse(p, varargin{:});
6+
variable_type = p.Results.variable_type;
7+
8+
switch variable_type
9+
case 'node'
10+
flat_function = @(x) x(:);
11+
otherwise
12+
flat_function = @(x) x(mask);
13+
end
14+
end

‎file_handlers/compact_file_append_sol_struct.m‎

Lines changed: 0 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -97,7 +97,6 @@ function compact_file_append_sol_struct(RP, output_file, all_pvals, all_pvals_ne
9797

9898
end
9999

100-
101100
% Update p-values for new repetitions
102101
for i_cell = 1:numel(reps_to_save)
103102
i = reps_to_save{i_cell}; % repetition index
@@ -121,18 +120,15 @@ function compact_file_append_sol_struct(RP, output_file, all_pvals, all_pvals_ne
121120

122121
method_struct.total_calculations = method_struct.total_calculations + 1;
123122
end
124-
125123

126124
has_field = cellfun(@(timing_struct) isfield(timing_struct, method_name), method_timing_all);
127-
128125

129126
if any(has_field)
130127
valid_indices = find(has_field);
131128
total_batch_time = sum(cellfun(@(timing_struct) timing_struct.(method_name), ...
132129
method_timing_all(valid_indices)));
133130
method_struct.total_time = method_struct.total_time + total_batch_time;
134131
end
135-
136132

137133
% Use eval to save method struct with dynamic name (unavoidable)
138134
eval([method_name ' = method_struct;']);

‎power_calculator_tools/process_repetition_batches.m‎

Lines changed: 1 addition & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -84,6 +84,7 @@ function process_repetition_batches(X, Y, RP, UI)
8484
STATS.edge_groups = RP.edge_groups;
8585
STATS.test_type = RP.test_type;
8686
STATS.unflatten_matrix = RP.unflat_matrix_fun;
87+
STATS.flatten_matrix = RP.flat_matrix_fun;
8788
STATS.mask = RP.mask;
8889
STATS.existing_repetitions = RP.existing_repetitions;
8990
STATS.all_submethods = RP.all_submethods;

‎power_calculator_tools/rep_cal_function.m‎

Lines changed: 2 additions & 1 deletion
Original file line numberDiff line numberDiff line change
@@ -95,9 +95,10 @@ function rep_cal_function(Params)
9595
[RP.flat_to_spatial, RP.spatial_to_flat] = create_spatial_flat_map(RP);
9696
[RP.triumask, RP.trilmask] = create_masks_from_nodes(size(RP.mask, 1));
9797

98-
% Create graph converter (flat to graph)
98+
% Create graph converter (flat to graph) and reverse
9999
RP.unflat_matrix_fun = unflatten_matrix(RP.mask, 'variable_type', ...
100100
RP.variable_type, 'flat_to_spatial', RP.flat_to_spatial, 'spatial_to_flat', RP.spatial_to_flat);
101+
RP.flat_matrix_fun = create_flat_function(RP.mask, 'variable_type', RP.variable_type);
101102

102103
[RP, test_type_origin] = infer_test_from_data(RP, OutcomeData.(t), BrainData);
103104

‎statistical_methods/.DS_Store‎

0 Bytes
Binary file not shown.

‎statistical_methods/Exact_FC_TFCE.m‎

Lines changed: 62 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,62 @@
1+
classdef Exact_FC_TFCE
2+
3+
properties
4+
level = "edge";
5+
permutation_based = true;
6+
permutations = 800; % Override permutation number
7+
method_params = Exact_FC_TFCE.get_tfce_params()
8+
end
9+
10+
methods (Static, Access = private)
11+
function method_params = get_tfce_params()
12+
method_params = struct();
13+
method_params.H = 3.0;
14+
method_params.E = 0.4;
15+
end
16+
end
17+
18+
methods
19+
20+
function pval = run_method(obj,varargin)
21+
22+
% Extract parameters
23+
params = struct(varargin{:});
24+
STATS = params.statistical_parameters;
25+
edge_stats = params.edge_stats;
26+
permuted_edge_stats = params.permuted_edge_data;
27+
28+
% Convert the edge statistics back into a matrix
29+
test_stat_mat = STATS.unflatten_matrix(edge_stats);
30+
31+
cluster_stats_target = apply_exact_tfce(test_stat_mat, ...
32+
'H', obj.method_params.H, 'E', obj.method_params.E);
33+
cluster_stats_target = STATS.flatten_matrix(cluster_stats_target);
34+
35+
% Ensure permutation data is provided
36+
if isempty(permuted_edge_stats)
37+
error('Permutation data is missing. Ensure precomputed permutations are provided.');
38+
end
39+
40+
% Number of permutations
41+
if size(permuted_edge_stats, 2) < obj.permutations
42+
K = size(permuted_edge_stats, 2);
43+
else
44+
K = obj.permutations;
45+
end
46+
null_dist = zeros(K, 1);
47+
48+
% Apply TFCE transformation to each permutation
49+
for i = 1:K
50+
perm_stat_mat = STATS.unflatten_matrix(permuted_edge_stats(:, i));
51+
tfce_null = apply_exact_tfce(perm_stat_mat, ...
52+
'H', obj.method_params.H, 'E', obj.method_params.E);
53+
null_dist(i) = max(tfce_null(:)); % Store max TFCE value for permutation
54+
end
55+
56+
% Compute p-values using permutation-based FWER correction
57+
pval = arrayfun(@(stat) (sum(stat <= null_dist)) /K, cluster_stats_target);
58+
59+
end
60+
61+
end
62+
end

‎statistical_methods/Fast_TFCE.m‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -40,12 +40,12 @@
4040
permuted_edge_stats = params.permuted_edge_data; % Explicitly using the new argument
4141

4242
% Convert the edge statistics back into a matrix
43-
test_stat_mat = unflatten_matrix(edge_stats, STATS.mask);
43+
test_stat_mat = STATS.unflatten_matrix(edge_stats);
4444

4545
% Apply TFCE transformation to the observed test statistics
4646
cluster_stats_target = apply_tfce(test_stat_mat, 'dh', obj.method_params.dh, ...
4747
'H', obj.method_params.H, 'E', obj.method_params.E);
48-
cluster_stats_target = flat_matrix(cluster_stats_target, STATS.mask);
48+
cluster_stats_target = STATS.flatten_matrix(cluster_stats_target);
4949

5050
% Ensure permutation data is provided
5151
if isempty(permuted_edge_stats)
@@ -62,14 +62,14 @@
6262

6363
% Apply TFCE transformation to each permutation
6464
for i = 1:K
65-
perm_stat_mat = unflatten_matrix(permuted_edge_stats(:, i), STATS.mask);
65+
perm_stat_mat = STATS.unflatten_matrix(permuted_edge_stats(:, i));
6666
tfce_null = apply_tfce(perm_stat_mat, 'dh', obj.method_params.dh, ...
6767
'H', obj.method_params.H, 'E', obj.method_params.E);
6868
null_dist(i) = max(tfce_null(:)); % Store max TFCE value for permutation
6969
end
7070

7171
% Compute p-values using permutation-based FWER correction
72-
pval = arrayfun(@(stat) (sum(stat <= null_dist)) /K, cluster_stats_target(:));
72+
pval = arrayfun(@(stat) (sum(stat <= null_dist)) /K, cluster_stats_target);
7373

7474
end
7575

‎statistical_methods/Fast_TFCE_cpp.m‎

Lines changed: 70 additions & 70 deletions
Original file line numberDiff line numberDiff line change
@@ -1,79 +1,79 @@
1-
classdef Fast_TFCE_cpp
2-
3-
properties
4-
level = "edge";
5-
permutation_based = true;
6-
permutations = 800; % Override permutation number
7-
method_params = Fast_TFCE_cpp.get_fast_tfce_params()
8-
end
9-
10-
methods (Static, Access = private)
11-
function method_params = get_fast_tfce_params()
12-
method_params = struct();
13-
method_params.dh = 0.1;
14-
method_params.H = 3.0;
15-
method_params.E = 0.4;
1+
classdef Fast_TFCE_cpp
2+
3+
properties
4+
level = "edge";
5+
permutation_based = true;
6+
permutations = 800; % Override permutation number
7+
method_params = Fast_TFCE_cpp.get_fast_tfce_params()
168
end
17-
end
189

19-
methods
20-
21-
function pval = run_method(obj,varargin)
22-
23-
% Applies Threshold-Free Cluster Enhancement (TFCE) and computes p-values
24-
% using a permutation-based approach.
25-
%
26-
% Inputs:
27-
% - STATS: Structure containing statistical parameters, including threshold.
28-
% - edge_stats: Raw test statistics for edges.
29-
% - permutation_edge_data: Precomputed permutation edge statistics.
30-
%
31-
% Outputs:
32-
% - pval: TFCE-corrected p-values.
33-
34-
params = struct(varargin{:});
35-
36-
% Extract relevant inputs
37-
STATS = params.statistical_parameters;
38-
edge_stats = params.edge_stats;
39-
permuted_edge_stats = params.permuted_edge_data; % Explicitly using the new argument
40-
41-
% Convert the edge statistics back into a matrix
42-
test_stat_mat = STATS.unflatten_matrix(edge_stats);
43-
44-
% Apply TFCE transformation to the observed test statistics
45-
cluster_stats_target = apply_tfce_cpp(test_stat_mat, obj.method_params.dh, ...
46-
obj.method_params.H, obj.method_params.E);
47-
cluster_stats_target = flat_matrix(cluster_stats_target, STATS.mask);
48-
49-
% Ensure permutation data is provided
50-
if isempty(permuted_edge_stats)
51-
error('Permutation data is missing. Ensure precomputed permutations are provided.');
10+
methods (Static, Access = private)
11+
function method_params = get_fast_tfce_params()
12+
method_params = struct();
13+
method_params.dh = 0.1;
14+
method_params.H = 3.0;
15+
method_params.E = 0.4;
5216
end
17+
end
5318

54-
% Number of permutations
55-
if size(permuted_edge_stats, 2) < obj.permutations
56-
K = size(permuted_edge_stats, 2);
57-
else
58-
K = obj.permutations;
59-
end
60-
null_dist = zeros(K, 1);
61-
62-
% Apply TFCE transformation to each permutation
63-
for i = 1:K
64-
perm_stat_mat = STATS.unflatten_matrix(permuted_edge_stats(:, i));
65-
tfce_null = apply_tfce_cpp(perm_stat_mat, obj.method_params.dh, ...
19+
methods
20+
21+
function pval = run_method(obj,varargin)
22+
23+
% Applies Threshold-Free Cluster Enhancement (TFCE) and computes p-values
24+
% using a permutation-based approach.
25+
%
26+
% Inputs:
27+
% - STATS: Structure containing statistical parameters, including threshold.
28+
% - edge_stats: Raw test statistics for edges.
29+
% - permutation_edge_data: Precomputed permutation edge statistics.
30+
%
31+
% Outputs:
32+
% - pval: TFCE-corrected p-values.
33+
34+
params = struct(varargin{:});
35+
36+
% Extract relevant inputs
37+
STATS = params.statistical_parameters;
38+
edge_stats = params.edge_stats;
39+
permuted_edge_stats = params.permuted_edge_data; % Explicitly using the new argument
40+
41+
% Convert the edge statistics back into a matrix
42+
test_stat_mat = STATS.unflatten_matrix(edge_stats);
43+
44+
% Apply TFCE transformation to the observed test statistics
45+
cluster_stats_target = apply_tfce_cpp(test_stat_mat, obj.method_params.dh, ...
6646
obj.method_params.H, obj.method_params.E);
67-
null_dist(i) = max(tfce_null(:)); % Store max TFCE value for permutation
47+
cluster_stats_target = STATS.flatten_matrix(cluster_stats_target);
48+
49+
% Ensure permutation data is provided
50+
if isempty(permuted_edge_stats)
51+
error('Permutation data is missing. Ensure precomputed permutations are provided.');
52+
end
53+
54+
% Number of permutations
55+
if size(permuted_edge_stats, 2) < obj.permutations
56+
K = size(permuted_edge_stats, 2);
57+
else
58+
K = obj.permutations;
59+
end
60+
null_dist = zeros(K, 1);
61+
62+
% Apply TFCE transformation to each permutation
63+
for i = 1:K
64+
perm_stat_mat = STATS.unflatten_matrix(permuted_edge_stats(:, i));
65+
tfce_null = apply_tfce_cpp(perm_stat_mat, obj.method_params.dh, ...
66+
obj.method_params.H, obj.method_params.E);
67+
null_dist(i) = max(tfce_null(:)); % Store max TFCE value for permutation
68+
end
69+
70+
% Compute p-values using permutation-based FWER correction
71+
pval = arrayfun(@(stat) (sum(stat <= null_dist)) /K, cluster_stats_target(:));
72+
6873
end
69-
70-
% Compute p-values using permutation-based FWER correction
71-
pval = arrayfun(@(stat) (sum(stat <= null_dist)) /K, cluster_stats_target(:));
72-
74+
7375
end
74-
76+
7577
end
76-
77-
end
78-
78+
7979

0 commit comments

Comments
 (0)