Added gitignore

Author: fabriciocravo

.DS_Store

@@ -0,0 +1,4 @@
+data/*
+power_calculator_results/*
+*.asv
+/test_power_calculator/*

.gitignore

@@ -80,19 +80,15 @@
     return;
 end
 
-
+%% CHECK IF NUMBER OF UNIQUE SUBJECTS
 if strcmp(GLM.test,'onesample') 
     % Calculate degrees of freedom for one-sample t-test
     df = GLM.n_observations - 1; 
     % Calculate uncorrected p-values from the t distribution
     p_uncorr = tcdf(-edge_stats__target, df);
 elseif strcmp(GLM.test,'ttest') 
-    %if strcmp(ttest_type,'paired') % TODO: create
-    % warning('Assuming paired sample and right-tailed t-test.');
-    df = GLM.n_observations/2-1; % observations are 2 per subject
-    %elseif strcmp(ttest_type,'unpaired')
-    %    df=GLM.n_observations-2; % assuming equal variances
-    %end
+    df = GLM.n_observations - 2; 
+    
     p_uncorr = tcdf(-edge_stats__target, df);
 elseif strcmp(GLM.test,'ftest') 
     error('Under development.');

NBS_addon/NBSedge_level_parametric_corr.m

@@ -47,31 +47,26 @@
 
 test_stat=zeros(1, GLM.n_GLMs);
 
-beta = zeros(GLM.n_predictors,GLM.n_GLMs);
 beta = GLM.X\GLM.y;
 
+% Round beta to the nearest multiple of 1e-14 to remove numerical issues
+beta = round(beta, 14);
+
 %Compute statistic of interest
 % TODO: consider moving to switch/case in glm_setup
-if strcmp(GLM.test,'onesample')
+if strcmp(GLM.test,'onesample') || strcmp(GLM.test,'ttest')
 
     % Compute the residuals
     resid = GLM.y - GLM.X * beta; 
     % Mean squared error
     mse = sum(resid.^2) / (GLM.n_observations - GLM.n_predictors); 
     % Compute the standard error using contrast
     se = sqrt(mse .* (GLM.contrast / (GLM.X' * GLM.X) * GLM.contrast'));
+    % Prevent division by zero - numerical fix
+    se(se < 1e-15) = 1e-15;
     % Compute the t-statistic
     test_stat(:) = (GLM.contrast * beta) ./ se;
 
-elseif strcmp(GLM.test,'ttest')
-
-    resid=zeros(GLM.n_observations,GLM.n_GLMs);
-    mse=zeros(GLM.n_observations,GLM.n_GLMs);
-    resid=GLM.y-GLM.X*beta;
-    mse=sum(resid.^2)/(GLM.n_observations-GLM.n_predictors);
-    se=sqrt(mse*(GLM.contrast*inv(GLM.X'*GLM.X)*GLM.contrast'));
-    test_stat(:)=(GLM.contrast*beta)./se;
-
 elseif strcmp(GLM.test,'ftest')
 
     sse=zeros(1,GLM.n_GLMs);

NBS_addon/NBSglm_smn.m

@@ -50,6 +50,25 @@
 create_test_contrast(RepParams.test_type, RepParams.n_subs_subset);
 
 
+%% Number of observations
+switch RepParams.test_type
+
+    case 't'
+        RepParams.observations = RepParams.n_subs_subset;
+
+    case 't2'
+        RepParams.observations = RepParams.n_subs_subset_c1 + RepParams.n_subs_subset_c2;
+
+    case 'r'
+        RepParams.observations = RepParams.n_subs_subset;
+
+    otherwise
+        error('Test type not supported')
+
+end
+
+
+
 %% Assign params to structures
 % Goal: should be able to run config file, load rep_params and UI from reference, and replicate reference results
 
@@ -90,7 +109,7 @@
 
 %% Set up DIMS
 UI.DIMS.nodes = RepParams.n_nodes;
-UI.DIMS.observations = RepParams.n_subs_subset;
+UI.DIMS.observations = RepParams.observations;
 UI.DIMS.predictors = size(RepParams.X_rep, 2);
 
 

NBS_benchmarking/setup_benchmarking.m

@@ -10,8 +10,8 @@
 
 % NBS toolbox
 Params.save_directory = './power_calculator_results/';
-%Params.data_dir = './data/s_abcd_fc_rosenblatt.mat';
-Params.data_dir = './data/s_hcp_fc_noble_tasks.mat';
+Params.data_dir = './data/s_abcd_fc_rosenblatt.mat';
+%Params.data_dir = './data/s_hcp_fc_noble_tasks.mat';
 Params.gt_data_dir = './power_calculator_results/ground_truth/';
 
 Params.gt_origin = 'power_calculator';
@@ -42,8 +42,8 @@
 Params.n_frames.REST2 = 1200;
 
 %%% Resampling parameters %%%
-Params.n_workers = 5; % num parallel workers for parfor, best if # workers = # cores
 Params.parallel = false; % run stuff sequentially or in parallel
+Params.n_workers = 8; % num parallel workers for parfor, best if # workers = # cores
 Params.n_repetitions = 500;  % 500 recommended
 
 %% List of subjects per subset

atlas_storage/test_hcp_fc_atlas.mat

@@ -0,0 +1,128 @@
+function plot_power_results(dataset_or_directory)
+    % If input is a dataset name, construct the expected directory path
+    base_dir = './power_calculator_results/power_calculation/';
+
+    if isfolder(fullfile(base_dir, dataset_or_directory))
+        data_dir = fullfile(base_dir, dataset_or_directory);
+    elseif isfolder(dataset_or_directory)
+        data_dir = dataset_or_directory; % Assume full directory path is given
+    else
+        error('Invalid dataset name or directory: %s', dataset_or_directory);
+    end
+
+    % Get list of all result files in directory
+    files = dir(fullfile(data_dir, '*.mat'));
+    if isempty(files)
+        error('No result files found in the specified directory: %s', data_dir);
+    end
+
+    % Initialize structure to store power results
+    power_results = struct();
+    unique_subject_numbers = [];   
+
+    % Process each file
+    for i = 1:numel(files)
+        file_path = fullfile(files(i).folder, files(i).name);
+        data = load(file_path);
+    
+        % Ensure meta-data and TPR exist
+        if ~isfield(data, 'power_data') || ~isfield(data.meta_data, 'subject_number') || ...
+           ~isfield(data, 'meta_data') || ~isfield(data.power_data, 'tpr')
+            warning('Skipping file (missing meta_data or brain_data.tpr): %s', files(i).name);
+            continue;
+        end
+
+        % Extract subject number from meta_data
+        n_subjects = data.meta_data.subject_number;
+
+        % Keep track of unique subject numbers
+        if ~ismember(n_subjects, unique_subject_numbers)
+            unique_subject_numbers = [unique_subject_numbers, n_subjects]; %#ok<AGROW>
+        end
+
+        % Extract power values (vector across tasks)
+        tpr_values = data.power_data.tpr(:);
+    
+        % Extract method name from meta_data
+        method_name = data.meta_data.test_type;
+        test_components = get_test_components_from_meta_data(data.meta_data.test_components);
+      
+        % Define field path
+        field_path = {method_name, sprintf('n_%d', n_subjects), test_components};
+
+        % Assign value safely using `setfield`
+        power_results = setfield(power_results, field_path{:}, mean(tpr_values));
+    
+    end
+    
+    % Check if we have results
+    if isempty(fieldnames(power_results))
+        error('No valid power results found.');
+    end
+
+    % Sort subject numbers and methods
+    unique_subject_numbers = sort(unique_subject_numbers);
+    subject_labels = arrayfun(@(x) sprintf('nsub %d', x), unique_subject_numbers, 'UniformOutput', false);
+    num_subjects = numel(unique_subject_numbers);
+    method_names = fieldnames(power_results);
+    plot_method_names = strrep(method_names, '_', ' ');
+    num_methods = numel(method_names);
+
+    % Generate a figure with subplots (one per subject number)
+    num_subplots = numel(unique_subject_numbers);
+    figure;
+    set(gcf, 'Position', [100, 100, 200 * num_subplots, 500]); % Adjust figure size
+
+
+    % Collect power data
+    for i = 1:num_subjects
+        n_subjects = unique_subject_numbers(i);
+        for j = 1:num_methods
+            method_name = method_names{j};
+            
+            if isfield(power_results.(method_name), sprintf('n_%d', n_subjects))
+                % Extract all task-specific power values
+                task_values = struct2cell(power_results.(method_name).(sprintf('n_%d', n_subjects)));
+    
+                % Convert to array
+                all_task_values = cell2mat(task_values);
+                
+                % Compute mean and standard error
+                mean_power(i, j) = mean(all_task_values);
+                error_power(i, j) = std(all_task_values) / sqrt(length(all_task_values));
+            end
+        end
+    end
+
+    % Create figure
+    figure;
+    hold on;
+    
+    % Define colors for each method
+    colors = lines(num_methods);
+    
+    % Create grouped bar plot
+    bar_handle = bar(subject_labels, mean_power, 'grouped');
+    
+    % Apply colors and add error bars
+    for j = 1:num_methods
+        bar_handle(j).FaceColor = colors(j, :); % Assign color per method
+        x_positions = bar_handle(j).XEndPoints; 
+        errorbar(x_positions, mean_power(:, j), error_power(:, j), 'k.', 'LineWidth', 1.5);
+    end
+    
+    % Formatting
+    xlabel('Number of Subjects');
+    ylabel('Average Power (%)');
+    ylim([0, 100]);
+    legend(plot_method_names, 'Location', 'northwest', 'Interpreter', 'none', 'FontSize', 12);
+    grid on;
+    hold off;
+    
+    dataset_or_directory = strrep(dataset_or_directory, '_', ' ');
+    % Main title
+    sgtitle(sprintf('Power Calculation Results for %s', dataset_or_directory), ...
+        'FontSize', 16, 'FontWeight', 'bold');
+
+end
+