Code prepared for cluster TFCE power test

Author: fabriciocravo

config_files/setparams.m

@@ -31,7 +31,7 @@
 % Params.data_dir = './data/s_abcd_fc_rosenblatt.mat';
 Params.data_dir = './data/s_hcp_fc_noble_tasks.mat';
 % Params.data_dir = './data/s_hcp_act_noble_1.mat';
-Params.output = 'test_tfce_comparison';
+Params.output = 'tfce_power_comp';
 
 % Save specifications - if NaN output becomes dataset file name
 Params.save_directory = './power_calculator_results/';
@@ -57,9 +57,9 @@
 
 %%% Resampling parameters %%%
 Params.parallel = false; % run stuff sequentially or in parallel
-Params.n_workers = 25; % num parallel workers for parfor, best if # workers = # cores
-Params.n_repetitions = 100;  % 500 recommended
-Params.batch_size = 25;
+Params.n_workers = 10; % num parallel workers for parfor, best if # workers = # cores
+Params.n_repetitions = 500;  % 500 recommended
+Params.batch_size = 10;
 
 
 %% Skip some tests - change ranges or the function
@@ -68,7 +68,7 @@
 
 
 %% List of subjects per subset
-Params.list_of_nsubset = {80}; % To change this, add more when necessary
+Params.list_of_nsubset = {20, 80, 200}; % To change this, add more when necessary
                     % size of subset is full group size (N=n*2 for two sample t-test or N=n for one-sample)
 
                             % Current model (see above design matrix) only designed for t-test
@@ -80,8 +80,8 @@
 Params.tpr_dthresh = 0; % Threshold for true positives vs negatives
 Params.save_significance_thresh = 0.15;
 % Params.all_cluster_stat_types = {'Parametric', 'Size_cpp', 'Fast_TFCE_cpp', 'Constrained_cpp', 'Omnibus_cNBS'};
-Params.all_cluster_stat_types = {'Exact_FC_TFCE', 'TFCE_dh1', 'TFCE_dh5', 'TFCE_dh10', 'TFCE_dh25', ...
-    'Fast_TFCE_dh1', 'Fast_TFCE_dh5', 'Fast_TFCE_dh10', 'Fast_TFCE_dh25'};
+Params.all_cluster_stat_types = {'IC_TFCE_FC_cpp_dh1','IC_TFCE_FC_cpp_dh5','IC_TFCE_FC_cpp_dh10', ...
+    'IC_TFCE_FC_cpp_dh25'};
 
 Params.all_submethods = {'FWER', 'FDR'};
 
@@ -94,8 +94,8 @@
 % Use a small subset of permutations for faster development -- inappropriate for inference
 
 Params.testing = true;
-Params.test_n_perms = 1;
-Params.test_n_repetitions = 1;
+Params.test_n_perms = 2;
+Params.test_n_repetitions = 2;
 Params.test_n_workers = 1;
 Params.test_disable_save = false;
 

statistical_methods/tfce_dh_methods/IC_TFCE_FC_cpp_dh25.m

@@ -0,0 +1,79 @@
+    classdef IC_TFCE_FC_cpp_dh25
+        
+        properties 
+            level = "edge";
+            permutation_based = true;
+            permutations = 800; % Override permutation number
+            method_params = IC_TFCE_FC_cpp_dh25.get_fast_tfce_params()
+        end
+    
+        methods (Static, Access = private)
+            function method_params = get_fast_tfce_params()
+                method_params = struct();
+                method_params.dh = 0.25;
+                method_params.H = 3.0;
+                method_params.E = 0.4;
+            end
+        end
+        
+        methods
+    
+            function pval = run_method(obj,varargin)
+    
+                % Applies Threshold-Free Cluster Enhancement (TFCE) and computes p-values
+                % using a permutation-based approach.
+                %
+                % Inputs:
+                %   - STATS: Structure containing statistical parameters, including threshold.
+                %   - edge_stats: Raw test statistics for edges.
+                %   - permutation_edge_data: Precomputed permutation edge statistics.
+                %
+                % Outputs:
+                %   - pval: TFCE-corrected p-values.
+            
+                params = struct(varargin{:});
+            
+                % Extract relevant inputs
+                STATS = params.statistical_parameters;
+                edge_stats = params.edge_stats;
+                permuted_edge_stats = params.permuted_edge_data; % Explicitly using the new argument
+            
+                % Convert the edge statistics back into a matrix
+                test_stat_mat = STATS.unflatten_matrix(edge_stats);
+            
+                % Apply TFCE transformation to the observed test statistics
+                cluster_stats_target = apply_tfce_cpp(test_stat_mat, obj.method_params.dh, ...
+                    obj.method_params.H, obj.method_params.E);
+                cluster_stats_target = STATS.flatten_matrix(cluster_stats_target);
+            
+                % Ensure permutation data is provided
+                if isempty(permuted_edge_stats)
+                    error('Permutation data is missing. Ensure precomputed permutations are provided.');
+                end
+            
+                % Number of permutations
+                if size(permuted_edge_stats, 2) < obj.permutations
+                    K = size(permuted_edge_stats, 2);
+                else
+                    K = obj.permutations;
+                end
+                null_dist = zeros(K, 1);
+            
+                % Apply TFCE transformation to each permutation
+                for i = 1:K
+                    perm_stat_mat = STATS.unflatten_matrix(permuted_edge_stats(:, i));
+                    tfce_null = apply_tfce_cpp(perm_stat_mat, obj.method_params.dh, ...
+                        obj.method_params.H, obj.method_params.E);
+                    null_dist(i) = max(tfce_null(:)); % Store max TFCE value for permutation
+                end
+            
+                % Compute p-values using permutation-based FWER correction
+                pval = arrayfun(@(stat) (sum(stat <= null_dist)) /K, cluster_stats_target(:));
+    
+            end
+            
+        end
+    
+    end
+    
+