Add glm_map object class: scikit-learn-style mass-univariate GLM estimator

CLAUDE.md

@@ -38,6 +38,7 @@ Other top-level classes (not subclasses of `image_vector`):
 - **`brainpathway` / `brainpathway_multisubject`** — connectivity / pathway-modeling objects.
 - **`canlab_dataset`** — generic subject × variable behavioral/clinical data container with its own `glm`, `mediation`, `scatterplot`, etc.
 - **`fmri_glm_design_matrix`**, **`fmri_timeseries`**, **`predictive_model`** — specialized containers for design matrices, raw timeseries, and ML model artifacts.
+- **`glm_map`** — scikit-learn-style estimator for mass-univariate GLM / multiple regression. Its property names mirror the fields of the results structure `fmri_data.regress` builds (the variable `out` in that method): the design spec (wrapping an `fmri_glm_design_matrix` in `.design`, or a direct `.X`), the fitted result maps (`.betas`/`.t`/`.contrast_estimates`/`.contrast_t` as `statistic_image`; `.df`/`.sigma`/`.residuals` as `fmri_data`), and three nested option/diagnostic structs `.input_parameters`, `.input_image_metadata`, `.diagnostics` (VIF, contrast VIF, leverage, condition number, collinearity). The historical out-struct field names are available as read/write aliases (`.b`→`.betas`, `.con_t`→`.contrast_t`, `.contrast_images`→`.contrast_estimates`, `.resid`→`.residuals`, `.variable_names`→`.regressor_names`, `.C`→`.contrasts`). Two workflows: (1) call `fmri_data.regress` directly — it now **returns a glm_map** (not a struct) — and query the object; (2) build an estimator: `g = glm_map(...)`, `add_contrasts`, `run_diagnostics`, `g = fit(g, fmri_data_obj)`, then `threshold`/`table`/`montage`. (The results live in the `.diagnostics` property; the method that computes them is `run_diagnostics` to avoid a property/method name clash.) `glm_map(out_struct)` re-casts a regress-style struct into an object. Supports 1st-level event designs (`build_design`, `import_SPM`) and 2nd-level/group designs; `fmri_data.regress` is the compute engine. Walkthrough: `docs/workflows/regression_with_glm_map.m`.
 
 ### MATLAB `@class` directories
 

CanlabCore/@fmri_data/regress.m

@@ -99,10 +99,13 @@
 % :Outputs:
 %
 %  **regression_results:**
-%        A structure containing stats_img and fmri_data objects.
-%        In addition to the main outputs below, the regression_results structure also has
-%        fields for input_parameters, the design matrix (X), variable
-%        names, and warnings.
+%        A glm_map object containing stats_img and fmri_data objects.
+%        (Historically this was a plain struct; it is now re-cast as a
+%        glm_map whose property names match the field names documented here,
+%        so out.b, out.t, out.contrast_images, out.con_t, etc. all still
+%        work. See help glm_map.) In addition to the main outputs below, it
+%        also carries input_parameters, input_image_metadata, the design
+%        matrix (X), variable names, diagnostics, and warnings.
 %
 %  **regression_results.b:**
 %        stats_img object of beta values estimated from regression
@@ -121,7 +124,8 @@
 %
 %  **regression_results.diagnostics:***
 %        A structure containing VIFs and leverage values for the design matrix
-%        regression_results.diagnostics.Variance_inflation_factors = VIFs
+%        regression_results.diagnostics.Variance_inflation_factors = VIFs (from VIF.m)
+%        regression_results.diagnostics.Contrast_variance_inflation_factors = contrast VIFs (cVIF.m), if C provided
 %        regression_results.diagnostics.Leverages = leverage values
 %
 % :Examples:
@@ -563,15 +567,28 @@
 
 % Variance inflation
 % ---------------------------------------------------------------------
-
-vifs = getvif(X);
+% Use canonical VIF.m (getvif is deprecated). When a contrast matrix C is
+% present, VIF also returns contrast VIFs (cVIF) for the contrasts. C is
+% stored column-wise [regressors x contrasts]; VIF/cVIF expect one contrast
+% per row, so pass C'.
+
+contrast_vifs = [];
+if ~isempty(C) && size(C, 1) == size(X, 2)
+    [vifs, contrast_vifs] = VIF(X, C');
+else
+    vifs = VIF(X);
+end
 
 if any(vifs > 4)
 
     mywarnings{end+1} = 'Warning!!!  Design multicolinearity. Some regressors have variance inflation factors > 4. Check regression_results.diagnostics';
 
 end
 
+if ~isempty(contrast_vifs) && any(contrast_vifs > 4)
+    mywarnings{end+1} = 'Warning!!!  Some contrasts have contrast variance inflation factors (cVIF) > 4. Check regression_results.diagnostics.Contrast_variance_inflation_factors';
+end
+
 % Leverages
 % ---------------------------------------------------------------------
 
@@ -1067,7 +1084,8 @@
 
 end
 
-out.diagnostics = struct('Variance_inflation_factors', vifs, 'Leverages', leverages);
+out.diagnostics = struct('Variance_inflation_factors', vifs, ...
+    'Contrast_variance_inflation_factors', contrast_vifs, 'Leverages', leverages);
 out.warnings = mywarnings;
 
 % Create objects
@@ -1199,6 +1217,25 @@
     warning on
 end
 
+% ---------------------------------------------------------------------
+% Re-cast the results structure as a glm_map object
+% ---------------------------------------------------------------------
+% The fields built above (b, t, df, sigma, contrast_images, con_t,
+% input_parameters, input_image_metadata, diagnostics, ...) map onto glm_map
+% properties of the same names. The object preserves struct-style field
+% access (out.b, out.t, out.con_t, ...) via Dependent aliases, so existing
+% callers continue to work. Guard the conversion so that any non-standard
+% output path (or a missing glm_map on the path) falls back to the struct.
+if isstruct(out) && isfield(out, 'b') && isfield(out, 't') && exist('glm_map', 'class') == 8
+    try
+        out = glm_map(out);
+        out = validate_object(out);   % ensure all nested-struct fields are present
+    catch ME
+        warning('fmri_data:regress:glm_mapCastFailed', ...
+            'Returning results as a struct; could not cast to glm_map: %s', ME.message);
+    end
+end
+
 % ---------------------------------------------------------------------
 % Subfunctions
 % ---------------------------------------------------------------------

CanlabCore/@fmri_glm_design_matrix/add.m

@@ -97,18 +97,31 @@
             isothermethod(i + 1) = false;
 
         case 'condition_names'
-
-            % special function for checking SPM-style format rules
-            [names, nsess, nconds] = check_req_and_length(obj, varargin{i+1});
-
-            indx = 1;
+
+            % Condition names are assumed to be the same for all sessions.
+            % Accept either nconds names (shared across sessions) or
+            % nsess*nconds names (session-major). Use the actual number of
+            % conditions from the onset structures when available.
+            names = varargin{i + 1};
+            if ~iscell(names), names = {names}; end
+
+            nsess = length(obj.nscan);
+            nconds = length(obj.Sess(1).U);
+            if nconds == 0
+                nconds = length(names) ./ nsess;   % legacy fallback (no onsets yet)
+            end
+
+            shared = (numel(names) == nconds);
             for ss = 1:nsess
                 for condname = 1:nconds
-                    obj.Sess(ss).U(condname).name = names{indx};
-                    indx = indx + 1;
+                    if shared
+                        obj.Sess(ss).U(condname).name = names{condname};
+                    else
+                        obj.Sess(ss).U(condname).name = names{(ss - 1) * nconds + condname};
+                    end
                 end
             end
-            
+
             isothermethod(i + 1) = false;
 
         case 'SPM'

CanlabCore/@fmri_glm_design_matrix/build.m

@@ -128,7 +128,10 @@
 end
 
 for i = 1:nconds
-    obj.xBF(i).name = sprintf('%s for Condition %3.0f', obj.xBF(i).name, i);
+    % Idempotent: strip any prior " for Condition N" suffix so repeated
+    % build() calls do not keep appending to the basis-set name.
+    basename = regexprep(obj.xBF(i).name, '\s*for Condition\s+\d+\s*$', '');
+    obj.xBF(i).name = sprintf('%s for Condition %3.0f', basename, i);
 end
 
 end

CanlabCore/@fmri_glm_design_matrix/build_single_trial.m

@@ -45,57 +45,57 @@
 [sess_delta,  sess_conditions,  sess_X, sess_C, sess_B] = deal(cell(1, nsess)); 
 
 
+res = 16;   % high-resolution samples per second (TR-independent)
+
 for s = 1:nsess
 
 [ons, name] = deal(cell(1, nconds));
 
 [ons{:}] = deal(obj.Sess(s).U(:).ons);
 %[name{:}] = deal(obj.Sess(s).U(:).name);
 
-% make sure onsets are in TRs, not secs
-switch obj.xBF(1).UNITS
-    case 'secs'
-        % onsets are in sec, convert
-        for i = 1:nconds
-            ons{i} = ons{i} ./ TR;
-        end
-
-    case {'tr', 'TR', 'trs'}
-        % ok, do nothing
+% Convert onsets to seconds (we build at high resolution, downsample to TR)
+to_sec = 1;
+if any(strcmpi(obj.xBF(1).UNITS, {'scans', 'tr', 'trs'})), to_sec = TR; end
+for i = 1:nconds
+    ons{i} = ons{i}(:) * to_sec;
 end
 
 % ----------------------------------------------
 % Build predictors for each session
 % ----------------------------------------------
 
-delta = onsets2delta(ons, obj.nscan(s));
-
-ns = size(delta, 1);
+nscan = obj.nscan(s);
+len_hires = round(nscan * TR * res);
 
 [trialdelta, cond_assignment] = deal(cell(1, nconds));
 
 % ----------------------------------------------
-% For each condition, parse into separate columns
-% for each onset (single-trial)
+% For each condition, parse into separate columns for each onset
+% (single-trial), at high resolution; convolve and downsample to TR.
 % ----------------------------------------------
 for j = 1:nconds
-
-    wh = find(delta(:, j));
-
-    trialdelta{j} = false(ns, length(wh));
-
-    for k = 1:length(wh)
-        trialdelta{j}(wh(k), k) = true;
+
+    os = ons{j};
+    ntrials = numel(os);
+
+    trialdelta{j} = zeros(len_hires, ntrials);   % one hi-res column per trial
+    for k = 1:ntrials
+        idx = round(os(k) * res) + 1;
+        if idx >= 1 && idx <= len_hires, trialdelta{j}(idx, k) = 1; end
     end
-
-    cond_assignment{j} = ones(length(wh), 1);
-
+
+    cond_assignment{j} = ones(ntrials, 1);
+
+    % custom per-condition HRF, sampled at high resolution (res samples/s)
     bf = obj.xBF(j);
     xvals = 0:TR:bf.length - TR;
-    hrf = interp1(xvals, inputhrf{j}(1:length(xvals)), 1:bf.dt:bf.length, 'spline', 'extrap');  
-
-    condX{j} = getPredictors(trialdelta{j}, hrf, 16);
-
+    hrf = interp1(xvals, inputhrf{j}(1:length(xvals)), 0:1/res:bf.length, 'spline', 'extrap');
+    hrf = hrf(:);
+
+    Xi = getPredictors(trialdelta{j}, hrf, 'dsrate', res * TR, 'force_delta');
+    condX{j} = Xi(1:nscan, :);
+
 end % conditions
 
 sess_delta{s} = cat(2, trialdelta{:}); % not needed?
@@ -116,13 +116,15 @@
 % Put the pieces together
 % ----------------------------------------------
 
-obj.xX.X = [blkdiag(sess_X{:}) blkdiag(sess_C{:}) blkdiag(sess_B{:})];
+% Index xX(1) explicitly: obj.xX starts as an empty (0x0) struct, and a
+% dot-assignment to an empty struct array is illegal.
+obj.xX(1).X = [blkdiag(sess_X{:}) blkdiag(sess_C{:}) blkdiag(sess_B{:})];
 
-obj.xX.cond_assignments = cat(1, sess_conditions{:});
+obj.xX(1).cond_assignments = cat(1, sess_conditions{:});
 
-obj.xX.iH = find(any(obj.xX.cond_assignments, 2));
-obj.xX.iC = find(~any(obj.xX.cond_assignments, 2));
-obj.xX.iB = size(obj.xX.X, 2) - nsess : size(obj.xX.X, 2);
+obj.xX(1).iH = find(any(obj.xX(1).cond_assignments, 2));
+obj.xX(1).iC = find(~any(obj.xX(1).cond_assignments, 2));
+obj.xX(1).iB = size(obj.xX(1).X, 2) - nsess + 1 : size(obj.xX(1).X, 2);
 
 
 end % function
@@ -174,7 +176,10 @@
 end
 
 for i = 1:nconds
-    obj.xBF(i).name = sprintf('%s for Condition %3.0f', obj.xBF(i).name, i);
+    % Idempotent: strip any prior " for Condition N" suffix so repeated
+    % build calls do not keep appending to the basis-set name.
+    basename = regexprep(obj.xBF(i).name, '\s*for Condition\s+\d+\s*$', '');
+    obj.xBF(i).name = sprintf('%s for Condition %3.0f', basename, i);
 end
 
 end

CanlabCore/@fmri_glm_design_matrix/create_orthogonal_contrast_set.m

@@ -0,0 +1,89 @@
+function obj = create_orthogonal_contrast_set(obj, varargin)
+% Assign an orthogonal set of contrasts spanning the regressors of interest.
+%
+% Uses create_orthogonal_contrast_set.m to build a set of mutually orthogonal,
+% zero-sum contrasts spanning the space of differences among the regressors of
+% interest (the H partition, obj.xX.iH), and stores them in obj.xX.contrasts
+% (one row per contrast, one column per design regressor) with placeholder
+% names in obj.xX.contrast_names. Nuisance covariates and session constants
+% get zero weight. The model must be built first (so obj.xX.iH is defined).
+%
+% :Usage:
+% ::
+%
+%     obj = create_orthogonal_contrast_set(obj)
+%     obj = create_orthogonal_contrast_set(obj, 'names', {'c1','c2',...})
+%
+% :Inputs:
+%
+%   **obj:**
+%        A built fmri_glm_design_matrix object (obj.xX.iH non-empty).
+%
+% :Optional Inputs:
+%
+%   **'names', {cellstr}:**  contrast names (default 'OrthC1', ...).
+%
+% :Outputs:
+%
+%   **obj:**
+%        with obj.xX.contrasts [n_contrasts x n_regressors] and
+%        obj.xX.contrast_names set.
+%
+% :See also:
+%   - create_orthogonal_contrast_set (function), build, glm_map.create_orthogonal_contrast_set
+%
+% ..
+%    2026 - Initial implementation.
+% ..
+
+names = {};
+wh = find(strcmpi(varargin, 'names'));
+if ~isempty(wh), names = varargin{wh(1) + 1}; end
+
+% Regressors of interest = H partition. Requires a built design.
+whI = [];
+if isstruct(obj.xX) && isscalar(obj.xX) && isfield(obj.xX, 'iH') && ~isempty(obj.xX.X)
+    whI = obj.xX.iH(:)';
+end
+
+if isempty(whI)
+    error('fmri_glm_design_matrix:NoInterestRegressors', '%s', local_interest_guidance());
+end
+if numel(whI) < 2
+    error('fmri_glm_design_matrix:TooFewInterest', ...
+        'Orthogonal contrasts require at least 2 regressors of interest (obj.xX.iH); found %d.', numel(whI));
+end
+
+nreg = size(obj.xX.X, 2);
+Cint = create_orthogonal_contrast_set(numel(whI));   % [(n-1) x n], rows = contrasts
+ncon = size(Cint, 1);
+
+C = zeros(ncon, nreg);
+C(:, whI) = Cint;
+
+if isempty(names)
+    names = arrayfun(@(i) sprintf('OrthC%d', i), 1:ncon, 'UniformOutput', false);
+elseif numel(names) ~= ncon
+    error('fmri_glm_design_matrix:ContrastNames', 'Number of names (%d) must match number of contrasts (%d).', numel(names), ncon);
+end
+
+obj.xX(1).contrasts      = C;
+obj.xX(1).contrast_names = names(:)';
+
+if ~iscell(obj.history), obj.history = {}; end
+obj.history{end + 1} = sprintf('create_orthogonal_contrast_set: %d contrast(s) over %d interest regressors', ncon, numel(whI));
+
+end % create_orthogonal_contrast_set
+
+
+% =========================================================================
+function s = local_interest_guidance()
+s = sprintf([ ...
+    'No regressors of interest are defined (obj.xX.iH is empty).\n' ...
+    'Regressors of interest are the event/task regressors (the H partition).\n' ...
+    'Add event onsets and build the design first, e.g.:\n' ...
+    '  d = fmri_glm_design_matrix(TR, ''nscan'', nscan, ''units'', ''secs'', ...\n' ...
+    '          ''onsets'', onsets, ''condition_names'', names);\n' ...
+    '  d = build(d);\n' ...
+    'After build, the event regressors are placed in obj.xX.iH.']);
+end