Changed output directory structure to be more organized

Author: fabriciocravo

NBS_benchmarking/summary_tools.m

@@ -6,15 +6,41 @@
 n = 268; % Full connectivity matrix for 268 ROIs (35,778 edges)
 img = zeros(n, n);
 
-% Create random functional connectivity values between -1 and 1
-% Fill the upper triangular part of the matrix with random FC values
+% Parameters for peak generation
+num_peaks = 8; % Number of functional network "hubs"
+peak_width = 15; % Width of Gaussian peaks (in ROI units)
+peak_strength = 0.99; % Maximum correlation strength
+background_noise = 0.15; % Background connectivity noise level
+
+% Generate random peak centers (representing functional network hubs)
+rng(20); % Set seed for reproducibility
+peak_centers = randi([1, n], num_peaks, 2);
+
+% Create the connectivity matrix
 for i = 1:n
     for j = (i+1):n
-        img(i,j) = (rand() * 2) - 1;  % Random value between -1 and 1
+        % Start with Gaussian background noise
+        value = randn() * background_noise;
+        
+        % Add Gaussian peaks centered at network hubs
+        for p = 1:num_peaks
+            % Distance from this edge to the peak center
+            dist_i = min(abs(i - peak_centers(p, 1)), abs(i - peak_centers(p, 2)));
+            dist_j = min(abs(j - peak_centers(p, 1)), abs(j - peak_centers(p, 2)));
+            dist = sqrt(dist_i^2 + dist_j^2);
+            
+            % Add Gaussian contribution from this peak
+            value = value + peak_strength * exp(-dist^2 / (2 * peak_width^2));
+        end
+        
+        % Add some distance-dependent decay (closer ROIs more correlated)
+        distance_factor = exp(-abs(i-j) / 100);
+        value = value + 0.2 * distance_factor * randn();
+        
+        img(i, j) = value;
     end
 end
 
-
 % Make it symmetric
 img = img + img';
 
@@ -87,6 +113,9 @@
 title('Absolute Difference', 'FontSize', 18, 'FontWeight', 'bold');
 colormap(gca, 'redblue');
 colorbar('FontSize', 12);
+tfce_max = max([max(tfced1(:)), max(tfced2(:))]);
+caxis([0, 1]);
 set(gca, 'FontSize', 14);
 xlabel('ROI', 'FontSize', 14);
-ylabel('ROI', 'FontSize', 14);
+ylabel('ROI', 'FontSize', 14);
+

TFCE_Benchmarkin/compare_full_tfce.m

@@ -66,15 +66,16 @@ function calculate_power(varargin)
     [Params.output, ~, ~] = get_data_set_name(Study_Info, Params);
 
     %% Process each repetition file one by one to reduce memory usage
-    rep_files = dir(fullfile(Params.save_directory, Params.output, '*.mat'));
+    rep_files = dir(fullfile(Params.save_directory, Params.output, "repetitions",'*.mat'));
 
     % If no files were found output an error
     if isempty(rep_files)
         error('No files found.')
     end
 
     %% Create output directory (only if it doesn't exist)
-    Params = create_power_output_directory(Params);
+    [Params.gt_data_dir, Params.gt_output] = setup_gt_directory(Params);
+    Params.save_directory = create_power_output_directory(Params);
 
     for i = 1:length(rep_files)
         % Load a single repetition data file
@@ -86,7 +87,6 @@ function calculate_power(varargin)
         method_list = rep_data.meta_data.method_list;
 
 
-
         [test_components, test_type, sub_number, testing_code] = get_data_for_file_naming(rep_data.meta_data);
         [~, file_name] = create_and_check_rep_file(Params.save_directory, Params.output, test_components, ...
                                                    test_type, ...
@@ -116,9 +116,9 @@ function calculate_power(varargin)
             method = method_list{j};
             method_data = rep_data.(method);
 
-            gt_filename = construct_gt_filename(rep_data.meta_data, Params.output);
-            gt_fullpath = fullfile(Params.gt_data_dir, Params.output, gt_filename);
-   
+            gt_filename = construct_gt_filename(rep_data.meta_data, Params.gt_output);
+            gt_fullpath = fullfile(Params.gt_data_dir, gt_filename);
+          
             if exist(gt_filename, 'file')
                 gt_data = load(gt_fullpath);
             else

calculate_power.m

@@ -31,11 +31,11 @@
 % 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 = 'tfce_power_comp';
+Params.output = 'test_new_output_struct';
 
 % Save specifications - if NaN output becomes dataset file name
 Params.save_directory = './power_calculator_results/';
-Params.gt_data_dir = './power_calculator_results/ground_truth/';
+% Params.gt_data_dir = './power_calculator_results/'; 
 
 
 
@@ -56,7 +56,7 @@
 Params.other_scripts_dir='./NBS_benchmarking/support_scripts/';
 
 %%% Resampling parameters %%%
-Params.parallel = true; % run stuff sequentially or in parallel
+Params.parallel = false; % run stuff sequentially or in parallel
 Params.n_workers = 5; % num parallel workers for parfor, best if # workers = # cores
 Params.n_repetitions = 500;  % 500 recommended
 Params.batch_size = 5;
@@ -92,8 +92,8 @@
 %%%%% DEVELOPERS ONLY %%%%%
 % Use a small subset of permutations for faster development -- inappropriate for inference
 
-Params.testing = false;
-Params.test_n_perms = 50;
+Params.testing = true;
+Params.test_n_perms = 10;
 Params.test_n_repetitions = 5;
 Params.test_n_workers = 1;
 Params.test_disable_save = false;

config_files/setparams.m

@@ -1,4 +1,12 @@
-function results = get_output_file_statistics(file_name)
+function results = get_output_file_statistics(file_name, varargin)
+
+    % Input Parser 
+    p = inputParser;
+
+    addOptional(p, 'attribute', 'tpr');
+    parse(p, varargin{:});
+    attribute = p.Results.attribute;
+    
 
     results = struct();
     results.mean = struct();
@@ -11,8 +19,8 @@
     for m = 1:length(method_list)
         method_name = method_list{m}; 
 
-        results.mean.(method_name) = mean(data.(method_name).tpr, 'omitnan');
-        results.deviation.(method_name) = std(data.(method_name).tpr, 'omitnan');
+        results.mean.(method_name) = mean(data.(method_name).(attribute), 'omitnan');
+        results.deviation.(method_name) = std(data.(method_name).(attribute), 'omitnan');
 
     end
 

data_fit_scripts/get_output_file_statistics.m

@@ -1,4 +1,12 @@
-function results_cell = multi_experiment_average(input_path)
+function results_cell = multi_experiment_average(input_path, varargin)
+    
+    % Input Parser 
+    p = inputParser;
+
+    addOptional(p, 'attribute_name_calculation', 'tpr');
+
+    parse(p, varargin{:});
+    attribute = p.Results.attribute_name_calculation;
 
     %% Check if input is multiple files or single
     multiple_files = false;
@@ -36,7 +44,7 @@
     for fi = 1:num_files
         full_file = files{fi};
 
-        file_results = get_output_file_statistics(full_file);
+        file_results = get_output_file_statistics(full_file, 'attribute', attribute);
 
         results_cell{fi} = file_results;
     end

data_fit_scripts/multi_experiment_average.m

@@ -11,18 +11,21 @@ function plot_aggregated_power_curve(varargin)
     default_dir = '/Users/f.cravogomes/Desktop/Cloned Repos/Power_Calculator/power_calculator_results/power_calculation/tfce_power_comp';
     default_undesired_sub_numbers = {};
     default_map = map_tfce_comp;
+    default_attribute_name = 'tpr';
 
      % Add optional parameter
     addParameter(p, 'dir', default_dir);
     addParameter(p, 'undesired_subject_numbers', default_undesired_sub_numbers);
     addParameter(p, 'map_function', default_map);
+    addParameter(p, 'attribute_name_calculation', default_attribute_name)
 
     % Parse input
     parse(p, varargin{:});
 
     directory = p.Results.dir;
     undesired_subject_numbers = p.Results.undesired_subject_numbers;
     map_function = p.Results.map_function;
+    attribute = p.Results.attribute_name_calculation;
 
     %% Check if input is a directory
     if ~isfolder(directory)
@@ -31,7 +34,7 @@ function plot_aggregated_power_curve(varargin)
     %%%%% Config end %%%%%%%%%%%%%%%%%%%%%%%%%%%
 
     %% Get averages 
-    multi_variable_average = multi_experiment_average(directory);
+    multi_variable_average = multi_experiment_average(directory, 'attribute_name_calculation', attribute);
 
     data_agregator = struct();
 

plot_scripts/plot_aggregated_power_curve.m

@@ -1,4 +1,4 @@
-function Params = create_output_directory(Params)
+function save_directory = create_output_directory(Params)
 %% create_output_directory
 % **Description**
 % Ensures the existence of the output directory for saving results. If the 
@@ -22,30 +22,26 @@
 % **Author**: Fabricio Cravo  
 % **Date**: March 2025
 
+    save_directory = [Params.save_directory, Params.output];
+
+    if ~exist(save_directory, 'dir') % Check if the directory does not exist
+        mkdir(save_directory);       % Create the directory
+    end
+
     if ~Params.ground_truth
 
-        if ~exist(Params.save_directory, 'dir') % Check if the directory does not exist
-            mkdir(Params.save_directory);       % Create the directory
-        end
-    
-        Params.save_directory = [Params.save_directory, Params.output, '/'];
+        save_directory = [save_directory, '/repetitions/'];
 
         if ~exist(Params.save_directory, 'dir') % Check if the directory does not exist
             mkdir(Params.save_directory);       % Create the directory
         end
 
     else
-
-        Params.save_directory = [Params.save_directory, 'ground_truth/'];
-
-        if ~exist(Params.save_directory, 'dir') % Check if the directory does not exist
-            mkdir(Params.save_directory);       % Create the directory
-        end
 
-        Params.save_directory = [Params.save_directory, Params.output, '/'];
+        save_directory = [save_directory, '/ground_truth/'];
 
-        if ~exist(Params.save_directory, 'dir') % Check if the directory does not exist
-            mkdir(Params.save_directory);       % Create the directory
+        if ~exist(save_directory, 'dir') % Check if the directory does not exist
+            mkdir(save_directory);       % Create the directory
         end
 
     end

power_calculator_tools/create_output_directory.m

@@ -1,15 +1,9 @@
-function Params = create_power_output_directory(Params)
-    
-    Params.save_directory = [Params.save_directory, 'power_calculation/'];
+function save_directory = create_power_output_directory(Params)
 
-    if ~exist(Params.save_directory, 'dir') % Check if the directory does not exist
-        mkdir(Params.save_directory);       % Create the directory
-    end
-
-    Params.save_directory = [Params.save_directory, Params.output, '/'];
+    save_directory = [Params.save_directory, Params.output, '/power_calculation/'];
 
-    if ~exist(Params.save_directory, 'dir') % Check if the directory does not exist
-        mkdir(Params.save_directory);       % Create the directory
+    if ~exist(save_directory, 'dir') % Check if the directory does not exist
+        mkdir(save_directory);       % Create the directory
     end
 
 end

power_calculator_tools/create_power_output_directory.m

@@ -0,0 +1,53 @@
+function positives = get_significance_vector(is_fp, ids_pos_vec, ids_neg_vec, PowerRes, ...
+        gt_size, stat_gt_level_str, pos_effect, neg_effect)
+    %%
+    % Calculates the absolute number of positives for both the tpr
+    % calculation and fpr
+    %
+    % is_fpr - bool - if we are calculating the fpr we require negation
+    % ids_pos_vec - bool array - varaibles with positive effect according to gt
+    % ids_neg_vec - bool array - variables with negative effect acordin to gt
+    % PowerRes - struct - current struct that stores the values (for gt accumulation)
+    % gt_size - array - size of gt variable
+    % stat_gt_level_str - level of statistics used
+    % pos_effect - for whole-brain only, if there is a pos_effect
+    % neg_effect - whole-brain only, if there is also a negative effect
+    
+    if is_fp
+        ids_pos = ~ids_pos_vec;
+        ids_neg = ~ids_neg_vec;
+    else
+        ids_pos = ids_pos_vec;
+        ids_neg = ids_neg_vec;
+    end
+
+    % calculate TPR
+    positives = nan(gt_size);
+    if contains(stat_gt_level_str, 'variable')
+
+        positives(ids_pos)=PowerRes.positives_total(ids_pos);
+        positives(ids_neg)=PowerRes.positives_total_neg(ids_neg);
+
+    elseif contains(stat_gt_level_str,'network')
+
+        positives(ids_pos)=PowerRes.positives_total(ids_pos);
+        positives(ids_neg)=PowerRes.positives_total_neg(ids_neg);
+
+    elseif contains(stat_gt_level_str,'whole_brain')
+
+        positives = 0;
+        if xor(pos_effect, is_fp)
+            positives = positives + PowerRes.positives_total;
+        end
+        if xor(neg_effect, is_fp)
+            positives = positives + PowerRes.positives_total_neg;
+        end
+        
+        %% Divide by two only if both positive and negative effects are present
+        if xor(pos_effect, is_fp) && xor(neg_effect, is_fp)
+           positives = floor(positives/2);
+        end
+
+    end
+
+end