Fix bug that made IC-TFCE a little slower

Author: fabriciocravo

.DS_Store

@@ -1,4 +1,4 @@
-%% This script compares TFCE speeds between different implementations 
+ %% This script compares TFCE speeds between different implementations 
 
 % Guarantee there is no behind the scenes parallelization
 maxNumCompThreads(1);

TFCE_Benchmarkin/compare_TFCE_speeds.m

@@ -1,139 +1,127 @@
-%% TFCE vs Fast_TFCE Speed Comparison Script
-% This script compares the speed of TFCE and Fast_TFCE methods across different dh values
-% to assess the performance gains.
-
-% Clear workspace but keep important variables
-vars = who;
-vars(strcmp(vars, 'data_matrix')) = [];
-vars(strcmp(vars, 'testing_yml_workflow')) = [];
-clear(vars{:});
-clc;
-
-% Initialize methods to compare
-fast_tfce_methods = {'Fast_TFCE_dh1', 'Fast_TFCE_dh5', 'Fast_TFCE_dh50', 'Fast_TFCE_dh100', 'Fast_TFCE_dh250'};
-tfce_methods = {'TFCE_dh1', 'TFCE_dh5', 'TFCE_dh50', 'TFCE_dh100', 'TFCE_dh250'};
-dh_values = [1, 5, 50, 100, 250]; % dh values for each method
-
-% Combine all methods
-all_methods = [fast_tfce_methods, tfce_methods];
-
-% Get the directory of the current script
-current_script_dir = fileparts(mfilename('fullpath'));
-parent_dir = fileparts(current_script_dir);
-
-% Change the current working directory to the parent directory
-cd(parent_dir);
-
-% Create the full path for the results directory - using a different folder
-speed_results_dir = './power_calculator_results/tfce_speed_dh_test/';
-if ~exist(speed_results_dir, 'dir')
-    mkdir(speed_results_dir);
-end
-
-% Define path for speed comparison results file
-speed_comparison_filepath = fullfile(speed_results_dir, 'tfce_speed_test_.mat');
-
-% Check if results already exist
-if exist(speed_comparison_filepath, 'file')
-    error('Delete or remove old file to continue')
-else
-    fprintf('No existing speed comparison results found. Starting fresh.\n');
-    speed_results = struct();
-end
-
-% Define path to clear
-clear_dir = fullfile(speed_results_dir, 'hcp_fc');
-if ~exist(clear_dir, 'dir')
-    mkdir(clear_dir);
-end
-
-% Set up parameters
-fprintf('\n========================================\n');
-fprintf('Starting TFCE vs Fast_TFCE dh Speed Comparison\n');
-fprintf('========================================\n\n');
-
-% Set parameters
-Params = setparams();
-Params.testing = false;
-Params.save_directory = speed_results_dir;
-Params.data_dir = './data/s_hcp_fc_noble_tasks.mat';
-Params.all_cluster_stat_types = all_methods; % All methods at once
-Params.n_perms = 100; % Reduced number of permutations since we're just comparing speed
-Params.n_repetitions = 10; % Reduced number of repetitions
-Params.list_of_nsubset = {20}; % Using 20 subjects
-Params.save_significance_thresh = 0.15;
-
-% ADD HERE: Skip tests on parameters if needed
-% Params.skip_param_tests = true; 
-
-% Run the analysis
-rep_cal_function(Params);
-
-% Find and load the results files
-result_files = dir(fullfile(speed_results_dir, 'hcp_fc', 'hcp_fc*.mat'));
-if isempty(result_files)
-    error('No result files found for %d repetitions', Params.n_repetitions);
+%% TFCE Speed Comparison - LaTeX Table Generator
+% Generates LaTeX tables for IC-TFCE vs TFCE speed comparison
+
+clear; clc;
+
+% Directory with result files
+results_dir = '/Users/f.cravogomes/Desktop/Cloned Repos/PRISME-Brain-Power-Calculator/power_calculator_results/re_tfce_speed_comp';
+
+% Method names
+ic_tfce_methods = {'IC_TFCE_FC_cpp_dh1', 'IC_TFCE_FC_cpp_dh5', 'IC_TFCE_FC_cpp_dh10', 'IC_TFCE_FC_cpp_dh25'};
+tfce_methods = {'TFCE_cpp_dh1', 'TFCE_cpp_dh5', 'TFCE_cpp_dh10', 'TFCE_cpp_dh25'};
+exact_method = 'Exact_FC_TFCE_cpp';
+dh_values = [1, 5, 10, 25];
+n_subs_list = [20, 80, 200];
+
+% Load all result files
+result_files = dir(fullfile(results_dir, '**/*.mat'));
+fprintf('Found %d result files\n', length(result_files));
+
+% Initialize storage
+ic_tfce_times = struct();
+tfce_times = struct();
+exact_tfce_times = struct();
+
+for ns = n_subs_list
+    key = sprintf('n%d', ns);
+    ic_tfce_times.(key) = cell(1, length(dh_values));
+    tfce_times.(key) = cell(1, length(dh_values));
+    exact_tfce_times.(key) = [];
+    for i = 1:length(dh_values)
+        ic_tfce_times.(key){i} = [];
+        tfce_times.(key){i} = [];
+    end
 end
 
-fprintf('Found %d result files. Loading all tasks...\n', length(result_files));
-
-% Initialize the results structure
-speed_results = struct();
-speed_results.tasks = {};
-speed_results.dh_values = dh_values;
-
-% Create structures to store timing data across all tasks
-all_fast_tfce_times = zeros(length(fast_tfce_methods), length(result_files));
-all_tfce_times = zeros(length(tfce_methods), length(result_files));
-task_names = cell(1, length(result_files));
-
-% Process all result files
-for file_idx = 1:length(result_files)
-    temp_file_path = fullfile(result_files(file_idx).folder, result_files(file_idx).name);
-    temp_results = load(temp_file_path);
+% Process each file
+for f = 1:length(result_files)
+    filepath = fullfile(result_files(f).folder, result_files(f).name);
+    data = load(filepath);
 
-    % Get the task name
-    task_name = strcat(temp_results.meta_data.test_components{1}, '_', ...
-        temp_results.meta_data.test_components{2});
-    task_names{file_idx} = task_name;
+    if ~isfield(data, 'meta_data') || ~isfield(data.meta_data, 'n_subs')
+        continue;
+    end
 
-    % Extract Fast_TFCE times for this task
-    for m = 1:length(fast_tfce_methods)
-        method = fast_tfce_methods{m};
-        if isfield(temp_results, method) && isfield(temp_results.(method), 'total_time')
-            all_fast_tfce_times(m, file_idx) = temp_results.(method).total_time;
-        else
-            all_fast_tfce_times(m, file_idx) = NaN;
-        end
+    n_subs = data.meta_data.n_subs;
+    key = sprintf('n%d', n_subs);
+    
+    if ~ismember(n_subs, n_subs_list)
+        continue;
     end
 
-    % Extract TFCE times for this task
-    for m = 1:length(tfce_methods)
-        method = tfce_methods{m};
-        if isfield(temp_results, method) && isfield(temp_results.(method), 'total_time')
-            all_tfce_times(m, file_idx) = temp_results.(method).total_time;
-        else
-            all_tfce_times(m, file_idx) = NaN;
+    % Extract times for each dh value
+    for i = 1:length(dh_values)
+        if isfield(data, ic_tfce_methods{i}) && isfield(data.(ic_tfce_methods{i}), 'total_time')
+            ic_tfce_times.(key){i}(end+1) = data.(ic_tfce_methods{i}).total_time;
+        end
+        if isfield(data, tfce_methods{i}) && isfield(data.(tfce_methods{i}), 'total_time')
+            tfce_times.(key){i}(end+1) = data.(tfce_methods{i}).total_time;
         end
     end
-
-    speed_results.all_fast_tfce_times = all_fast_tfce_times;
-    speed_results.all_tfce_times = all_tfce_times;
+    
+    % Extract Exact TFCE time
+    if isfield(data, exact_method) && isfield(data.(exact_method), 'total_time')
+        exact_tfce_times.(key)(end+1) = data.(exact_method).total_time;
+    end
 end
 
-% Store all task names
-speed_results.tasks = task_names;
+% Generate LaTeX output
+fprintf('\n\n%% ========== COPY BELOW ==========\n\n');
 
-% Calculate average times across all tasks
-speed_results.fast_tfce_times = nanmean(all_fast_tfce_times, 2)';
-speed_results.tfce_times = nanmean(all_tfce_times, 2)';
+fprintf('\\begin{table}[h!]\n');
+fprintf('\\centering\n');
+fprintf('\\caption{Runtime comparison of TFCE implementations across varying sample sizes.}\n');
+fprintf('\\label{tab:tfce_comparison}\n\n');
 
-% Calculate speed gain (time difference and percentage)
-speed_results.speed_gain = zeros(1, length(dh_values));
+labels = {'(a)', '(b)', '(c)'};
 
-for i = 1:length(dh_values)
-    speed_results.speed_gain(i) = speed_results.tfce_times(i)/speed_results.fast_tfce_times(i);     
+for ns_idx = 1:length(n_subs_list)
+    ns = n_subs_list(ns_idx);
+    key = sprintf('n%d', ns);
+    
+    fprintf('\\vspace{0.5em}\n');
+    fprintf('\\textbf{%s N = %d}\n\n', labels{ns_idx}, ns);
+    fprintf('\\begin{tabular}{cccc}\n');
+    fprintf('\\hline\n');
+    fprintf('dh & IC\\_TFCE (s) & TFCE (s) & Speedup \\\\\n');
+    fprintf('\\hline\n');
+    
+    for i = 1:length(dh_values)
+        ic_mean = mean(ic_tfce_times.(key){i});
+        tfce_mean = mean(tfce_times.(key){i});
+        speedup = tfce_mean / ic_mean;
+        fprintf('%d & %.4f & %.4f & %.2fx \\\\\n', dh_values(i), ic_mean, tfce_mean, speedup);
+    end
+    
+    fprintf('\\hline\n');
+    fprintf('\\end{tabular}\n\n');
 end
 
-% Save the speed comparison results
-save(speed_comparison_filepath, 'speed_results');
+% Exact TFCE table - by subject size
+fprintf('\\vspace{0.5em}\n');
+fprintf('\\textbf{(d) Exact TFCE}\n\n');
+fprintf('\\begin{tabular}{cccc}\n');
+fprintf('\\hline\n');
+fprintf('N & Exact TFCE (s) & \\begin{tabular}[c]{@{}c@{}}Speedup\\\\IC-TFCE\\end{tabular} & \\begin{tabular}[c]{@{}c@{}}Speedup\\\\TFCE\\end{tabular} \\\\\n');
+fprintf('\\hline\n');
+
+for ns_idx = 1:length(n_subs_list)
+    ns = n_subs_list(ns_idx);
+    key = sprintf('n%d', ns);
+    
+    exact_mean = mean(exact_tfce_times.(key));
+    % Compare to dh=1 (finest discretization)
+    ic_mean_dh1 = mean(ic_tfce_times.(key){1});
+    tfce_mean_dh1 = mean(tfce_times.(key){1});
+    
+    speedup_ic = exact_mean / ic_mean_dh1;
+    speedup_tfce = exact_mean / tfce_mean_dh1;
+    
+    fprintf('%d & %.4f & %.2fx & %.2fx \\\\\n', ns, exact_mean, speedup_ic, speedup_tfce);
+end
+
+fprintf('\\hline\n');
+fprintf('\\end{tabular}\n\n');
+fprintf('\\end{table}\n');
+
+fprintf('\n%% ========== END COPY ==========\n');

TFCE_Benchmarkin/compare_different_dhs.m

@@ -19,40 +19,41 @@
     lb = p.Results.lower_bounds;
     ub = p.Results.upper_bounds;
 
-    % Initial parameter estimates
-    % a: roughly where we expect decent power (median of x range)
-    % b: start with exponential shape (b=1)
     a_init = median(results_x);
     b_init = 1;
     c_init = 50;
 
-    % Initial guess: [scale, shape]
+    % Initial parameter estimates
+    % a: roughly where we expect decent power (median of x range)
+    % b: start with exponential shape (b=1)
     initial_params = [a_init, b_init, c_init];
 
     % Define cost function (sum of squared residuals)
     cost_func = @(params) sum((results_y - power_func(params, results_x)).^2);
 
     % Set optimization options with bounds to ensure positive parameters
-    options = optimset('Display', 'off', 'MaxFunEvals', 10000);
-
-    if exist('fmincon', 'file')
-        [fitted_params, ~] = fmincon(cost_func, initial_params, [], [], [], [], lb, ub, [], options);
-    else
-        % Fallback to fminsearch if Optimization Toolbox not available
-        [fitted_params, ~] = fminsearch(cost_func, initial_params, options);
-    end
-
+    options = optimoptions('fmincon', ...
+        'Display', 'off', ...
+        'MaxFunctionEvaluations', 10000, ...
+        'MaxIterations', 5000, ...
+        'OptimalityTolerance', 1e-8, ...
+        'StepTolerance', 1e-10, ...
+        'Algorithm', 'interior-point');
+
+    [fitted_params, ~] = fmincon(cost_func, initial_params, [], [], [], [], lb, ub, [], options);
+    
 
     % Generate smooth curve for plotting
     x = linspace(0, max(results_x), 200);
     y = power_func(fitted_params, x);
-    
+
+    % Clip power values to avoid weird fits
+    y = max(0, min(100, y));
+
     % Calculate R-squared
     y_pred = power_func(fitted_params, results_x);
     ss_res = sum((results_y - y_pred).^2);
     ss_tot = sum((results_y - mean(results_y)).^2);
     r_squared = 1 - (ss_res / ss_tot);
 
-
-
 end

data_fit_scripts/fit_power_curve.m

@@ -0,0 +1,8 @@
+function variable_type = get_variable_type_from_meta_data(meta_data)
+
+    if isfield(meta_data, 'variable_type')
+        variable_type = meta_data.variable_type;
+        return;
+    end
+
+end

file_handlers/get_variable_type_from_meta_data.m

@@ -2,8 +2,8 @@
     average_heat_map_gen(method, n_subs, directory, varargin)   
 %
 % For the paper: Shitf - down to select for paste
-% - average_heat_map_gen('Constrained_FDR',80,'/Users/f.cravogomes/Desktop/Pc_Res_Updated/SHOCK Paper Results/power_calculation/hcp_fc')
-% - average_heat_map_gen('Size_cpp',500,'/Users/f.cravogomes/Desktop/Pc_Res_Updated/SHOCK Paper Results/power_calculation/abcd_100_sex/')
+% - average_heat_map_gen('Constrained_FDR',80,'/Users/f.cravogomes/Desktop/Pc_Res_Updated/PRISME Paper Results/power_calculation/hcp_fc')
+% - average_heat_map_gen('Size_cpp',500,'/Users/f.cravogomes/Desktop/Pc_Res_Updated/PRISME Paper Results/power_calculation/abcd_100_sex')
 
     % Create parser object
     p = inputParser;
@@ -123,9 +123,12 @@
     end
 
     % Figure config
-    figure('Position', [100, 100, 800, 700]);
+    figure('Position', [100, 100, 800, 700], 'Color', 'white');
 
     imagesc(power_data);
+
+    set(gca, 'XTickLabel', []);
+    set(gca, 'YTickLabel', []);
 
     c = colorbar;
     c.FontSize = 12;