Merge pull request #2 from neuroprismlab/abcd_dataset_code

Abcd dataset code

Author: fabriciocravo

NBS_addon/NBSedge_level_parametric_corr.m

@@ -67,6 +67,20 @@
 GLM = NBSglm_setup_smn(GLM);
 edge_stats__target = NBSglm_smn(GLM);
 
+% If gt - we only need the result from the glm and the cluster stats
+if STATS.ground_truth
+    % Get the shape of edge_stats__target
+    shape = size(edge_stats__target);
+
+    % Create variables with the same shape
+    con_mat = {false(shape)};  % Empty logical array in a cell
+    pval = {NaN(shape)};  % NaN array inside a cell
+    any_significant = false; % Ensures downstream code doesn't break
+
+    return;
+end
+
+
 if strcmp(GLM.test,'onesample') 
     % Calculate degrees of freedom for one-sample t-test
     df = GLM.n_observations - 1; 

NBS_addon/NBSrun_smn.m

@@ -402,6 +402,12 @@
 % Number of nodes
 nbs.STATS.N = DIMS.nodes; 
 
+try
+   nbs.STATS.ground_truth = logical(UI.ground_truth);
+catch
+   error('Error setting up gt value')
+end
+
 % Do error checking on user inputs
 [msg,stop] = errorcheck(UI,DIMS,S);
 
@@ -467,6 +473,11 @@
 end
 
 
+% If gt change the method 
+% We only need the GLM result for the gt
+
+
+
 %Do NBS
 if strcmp(UI.method.ui, 'Run NBS')
 

NBS_addon/NBSstats_smn.m

@@ -1,4 +1,5 @@
-function [any_significant, con_mat, pval, edge_stats__target, cluster_stats__target] = NBSstats_smn(varargin)
+function [any_significant, con_mat, pval, edge_stats__target, cluster_stats__target] = ...
+NBSstats_smn(varargin)
 %NBSstats Computes network components among edges that survive a primary 
 %test statistic threshold. Assigns a corrected p-value to each component 
 %using permuted data that has been supplied as part of the structure STATS. 
@@ -187,6 +188,25 @@
 edge_stats__target = get_univariate_stats(STATS,GLM,precomputed,1);
 [cluster_stats__target, max_stat__target] = get_cluster_stats(edge_stats__target,STATS,ind_upper, N, Intensity,bgl);
 
+% If gt - we only need the result from the glm and the cluster stats
+if STATS.ground_truth
+     % Get the shape of cluster_stats__target
+    shape = size(cluster_stats__target);
+
+    % Create variables with the same shape
+    con_mat = false(shape);  % Logical array (not inside a cell)
+    
+    % Check if Omnibus method
+    if strcmp(STATS.statistic_type, 'Omnibus')
+        pval = 1;         % Scalar p-value
+    else
+        pval = NaN(shape);  % NaN array (not inside a cell)
+    end
+
+    any_significant = false; % Ensures downstream code doesn't break
+    return;
+end
+
 % 2, Permuted cluster statistics
 %First row of test_stat is the observed test statistics, so start at the second row
 for i=2:K+1

NBS_benchmarking/setup_benchmarking.m

@@ -85,10 +85,14 @@
 UI.mask.ui = RepParams.mask;
 % UI.do_Constrained_FWER_second_level.ui=do_Constrained_FWER_second_level;
 
+% Ground truth calculations do not require p-values
+UI.ground_truth = RepParams.ground_truth;
+
 %% Set up DIMS
 UI.DIMS.nodes = RepParams.n_nodes;
 UI.DIMS.observations = RepParams.n_subs_subset;
 UI.DIMS.predictors = size(RepParams.X_rep, 2);
 
+
 end
 

NBS_benchmarking/summarize_tprs.m

@@ -56,12 +56,16 @@
 % addOptional(p,'summary_type',summary_type_default);
 addRequired(p,'summary_type');
 addOptional(p,'tpr_dthresh', 0);
-addOptional(p,'save_directory', './power_calculation_results/')
+addOptional(p,'save_directory', NaN)
 parse(p, summary_type, varargin{:});
 
 % summary_type=p.Results.summary_type;
 tpr_dthresh = p.Results.tpr_dthresh;
-save_directory = [p.Results.save_directory, '/power_calculation/'];
+save_directory = p.Results.save_directory;
+
+if isnan(save_directory)
+    error('Save directory not specified')
+end
 
 % Make save directory if it does not exist
 
@@ -104,7 +108,7 @@
 
     %% Check if already calculated and get file_name
     data_set_name = strcat(rep_data.meta_data.dataset, '_', rep_data.meta_data.map);
-    test_components = strcat(rep_data.meta_data.test_components{1}, '_', rep_data.meta_data.test_components{2});
+    test_components = get_test_components_from_meta_data(rep_data.meta_data.test_components);
     [~, file_name] = create_and_check_rep_file(save_directory, data_set_name, test_components, ...
                                                rep_data.meta_data.test, rep_data.meta_data.test_type, ...
                                                rep_data.meta_data.omnibus, rep_data.meta_data.subject_number, ...

NBS_benchmarking/support_scripts/atlas_scripts/load_atlas_edge_groups.m

@@ -1,4 +1,4 @@
-function edge_groups=load_atlas_edge_groups(n_nodes, mapping_category, atlas_file)
+function edge_groups=load_atlas_edge_groups(atlas_file)
 
 % returns a n_nodes x n_nodes matrix where the value of each edge reflects
 % the network pair it belongs to

atlas_storage/test_hcp_fc_atlas.mat

@@ -22,34 +22,27 @@
 clear(vars{:});   % Clear all other variables
 clc;
 
-data_dir = '/Users/f.cravogomes/Desktop/Cloned Repos/NBS_Calculator/data/s_hcp_fc_noble_tasks.mat';
-if ~exist('Dataset', 'var')
-    Dataset = load(data_dir);
-else
-    % disp('Data already loaded')
-end
-% Extract dataset name
-[~, data_set_name, ~] = fileparts(data_dir);
-
 [current_path, ~, ~] = fileparts(mfilename('fullpath')); % assuming NBS_benchmarking is current folder
 addpath(genpath(current_path));
 
 %% Prepare parameters and dataset
 Params = setparams();
+Params.ground_truth = true;
 
+%% Load data
 if ~exist('Dataset', 'var')
     Dataset = load(Params.data_dir);
-else
-    % disp('Data already loaded')
 end
 
 %% Set n_nodes, n_var, n_repetitions 
 Params = setup_experiment_data(Params, Dataset);
 
-%% Create directory and get dataset name
-Params.save_directory = [Params.save_directory, 'ground_truth/'];
-Params = create_output_directory(Params);
-Params.data_set = get_data_set_name(Dataset);
+%% Create directory and get dataset name - get atlas
+[Params.data_set, Params.data_set_base, Params.data_set_map] = get_data_set_name(Dataset);
+Params.atlas_file = atlas_data_set_map(Params);
+
+%% Create gt output directory
+Params = create_gt_output_directory(Params);
 
 %% Paralle workers
 setup_parallel_workers(Params.parallel, Params.n_workers);
@@ -66,23 +59,30 @@
     RP = Params;
     RP.test_name = t;
 
-    RP = infer_test_from_data(RP, OutcomeData.(t), BrainData);
+    [RP, test_type_origin] = infer_test_from_data(RP, OutcomeData.(t), BrainData);
+        
+    % Important, X is calculated here for all test_types
+    % However, only the r test type will use this test
+    switch test_type_origin
+        
+        % Here: RP.test_name, RP.n_subs_1, RP.n_subs_2, RP.n_subs
+        case 'score_cond'
+            [X, Y , RP] = subs_data_from_score_condition(RP, OutcomeData.(t), BrainData, t);
+       
+        case 'contrast'
+            [X, Y , RP] = subs_data_from_contrast(RP, OutcomeData.(t).contrast, BrainData);
+        
+        otherwise
+            error('Test type origin not found')
+
+    end  
 
-    % y_and_x also extracts subject number and sub_ids
-    % Encapsulate this level of setup in a function if there is to much
-    % function calls here 
-    [~, Y , RP] = subs_data_from_contrast(RP, OutcomeData.(t).contrast, BrainData);
     [RP.triumask, RP.trilmask] = create_masks_from_nodes(size(RP.mask, 1));
 
     % Sets subject repetition to all subjects and others
     RP = setup_ground_truth_parameters(RP);
 
-    run_benchmarking(RP, Y)
-    
-    return;
-    if RP.testing == 1 && ti == 2
-        return;
-    end
+    run_benchmarking(RP, Y, X)
 
 end
 

calculate_gt.m

@@ -1,37 +1,49 @@
-%% Questions
-    % There appears to be something wrong in the network level power
-    % calculations - edge level and whole brain are fine
-    % I likely made an error in the atlas 
-    % Cluster has network-based stats not edge
-    % Found it - Likely problem in edge groups!
-    % 
-    %% Potential fix - More resonable results!
-    %   Used tril mask on extract_atlas_related_parameters for edge groups
-    %   Removed atlas reordering in setupbenchmarking 
-
-% Initial setup
+%% Initial setup
 scriptDir = fileparts(mfilename('fullpath'));
 addpath(genpath(scriptDir));
 cd(scriptDir);
-vars = who;       % Get a list of all variable names in the workspace
-vars(strcmp(vars, 'RepData')) = [];  % Remove the variable you want to keep from the list
-vars(strcmp(vars, 'GtData')) = [];
-clear(vars{:});   % Clear all other variables
+clearvars -except RepData GtData;
 clc;
 
-%% Directory to save and find rep data - TODO add them all
+%% Directory to save and find rep data
 Params = setparams();
 
-%% Create storage directory - only if it does not exist
-if ~exist(Params.save_directory, 'dir') % Check if the directory does not exist
-    mkdir(Params.save_directory);       % Create the directory
+% Load dataset information
+if ~exist('Dataset', 'var')
+    Study_Info = load(Params.data_dir, 'study_info');
 end
+[Params.data_set, ~, ~] = get_data_set_name(Study_Info);
+
+%% Process each repetition file one by one to reduce memory usage
+rep_files = dir(fullfile(Params.save_directory, Params.data_set, '*.mat'));
 
-if ~exist('RepData', 'var') || ~exist('GtData', 'var')
-    [GtData, RepData] = load_rep_and_gt_results(Params, 'gt_origin', Params.gt_origin);
-end 
+%% Create output directory (only if it doesn't exist)
+Params = create_power_output_directory(Params);
 
-power_calculation_tprs = @(x) summarize_tprs('calculate_tpr', x, GtData, ...
-                                             'save_directory', Params.save_directory);
-dfs_struct(power_calculation_tprs, RepData);
+for i = 1:length(rep_files)
+    % Load a single repetition data file
+    rep_file_path = fullfile(rep_files(i).folder, rep_files(i).name);
+    rep_data = load(rep_file_path);
+    
+    % Extract metadata to find corresponding GT file
+    if isfield(rep_data, 'meta_data') && isfield(rep_data.meta_data, 'test_components')
+        gt_filename = construct_gt_filename(rep_data.meta_data);
+        gt_file_path = [Params.gt_data_dir, Params.data_set, '/', gt_filename];
+        
+        % Load only the necessary GT data
+        if exist(gt_file_path, 'file')
+            gt_data = load(gt_file_path);
+        else
+            warning('GT file %s not found. Skipping...', gt_filename);
+            continue;
+        end
 
+        % Compute power using the extracted repetition and GT data
+        summarize_tprs('calculate_tpr', rep_data, gt_data, 'save_directory', Params.save_directory);
+    else
+        warning('Metadata missing in %s, skipping...', rep_files(i).name);
+    end
+    
+    % Free memory before next iteration
+    clear rep_data gt_data;
+end

calculate_power.m

@@ -0,0 +1,43 @@
+%% Questions
+    % There appears to be something wrong in the network level power
+    % calculations - edge level and whole brain are fine
+    % I likely made an error in the atlas 
+    % Cluster has network-based stats not edge
+    % Found it - Likely problem in edge groups!
+    % 
+    %% Potential fix - More resonable results!
+    %   Used tril mask on extract_atlas_related_parameters for edge groups
+    %   Removed atlas reordering in setupbenchmarking 
+
+% Initial setup
+scriptDir = fileparts(mfilename('fullpath'));
+addpath(genpath(scriptDir));
+cd(scriptDir);
+vars = who;       % Get a list of all variable names in the workspace
+vars(strcmp(vars, 'RepData')) = [];  % Remove the variable you want to keep from the list
+vars(strcmp(vars, 'GtData')) = [];
+clear(vars{:});   % Clear all other variables
+clc;
+
+%% Directory to save and find rep data - TODO add them all
+Params = setparams();
+
+% Load dataset
+if ~exist('Dataset', 'var')
+    Study_Info = load(Params.data_dir, 'study_info');
+end
+[Params.data_set, ~, ~] = get_data_set_name(Study_Info);
+
+if ~exist('RepData', 'var') || ~exist('GtData', 'var')
+    [GtData, RepData] = load_rep_and_gt_results(Params, 'gt_origin', Params.gt_origin, ...
+                                                'dataset', Params.data_set);
+end 
+
+%% Create storage directory - only if it does not exist
+% the previous directory is used to load the data
+Params = create_power_output_directory(Params);
+
+power_calculation_tprs = @(x) summarize_tprs('calculate_tpr', x, GtData, ...
+                                             'save_directory', Params.save_directory);
+dfs_struct(power_calculation_tprs, RepData);
+