% Model
% 
% converted from quarto-docs/tutorial-fmri-spm/model.qmd

% select subject
sub = 'subj1';
% load information
load(fullfile(sub, 'info.mat'))
% BOLD data (contains realignment information)
fnBOLD = fullfile(sub, 'bold.nii');
% realigned BOLD data
fnrBOLD = fullfile(sub, 'rbold.nii');
% directory for model
modelDir = fullfile(sub, 'model');
% file containing model
fnSPM = fullfile(modelDir, 'SPM.mat');

%%

% ## Specify the model
% 
% Instead of the approach of the paper to use a single design matrix across all
% sessions, we here follow the standard approach of SPM, and use session-wise
% design matrices. Each of them contains eight HRF-convoluted regressors
% corresponding to the stimulus presentations of the eight conditions, six
% motion regressors derived from the realignment parameters, and a constant
% regressor.

% get realignment parameters for motion regressors
V = spm_vol(fnBOLD);
mat = cat(4, V.mat);

% create directory
mkdir(modelDir)

% SPM module 'fMRI model specification'
fmri_spec = struct;
fmri_spec.dir = {modelDir};         % model directory
fmri_spec.timing.units = 'secs';    % units for onsets and durations
fmri_spec.timing.RT = TR;           % repetition time (interscan interval)
fmri_spec.timing.fmri_t = 25;       % microtime resolution: 0.1 s
fmri_spec.timing.fmri_t0 = 12;      % microtime onset: middle of scan
% for each session
for si = 1 : nSessions
    % indices of scans in 4d BOLD image corresponding to session
    vi = (si - 1) * nVolsPerSession + (1 : nVolsPerSession);
    fmri_spec.sess(si).scans = arrayfun(@(i) ...
        sprintf('%s,%d', ...                % SPM's <filename>,<scan> syntax
            fullfile(pwd, fnrBOLD), i), ... % filename must include full path
            vi, ...
            'UniformOutput', false)';
    % for each condition
    for ci = 1 : nConds
        % name, onsets, and durations
        fmri_spec.sess(si).cond(ci).name = conditions{ci};
        fmri_spec.sess(si).cond(ci).onset = onsets{si, ci};
        fmri_spec.sess(si).cond(ci).duration = durations{si, ci};
        fmri_spec.sess(si).cond(ci).tmod = 0;   % no time modulation
        fmri_spec.sess(si).cond(ci).orth = 1;
    end
    fmri_spec.sess(si).multi = {''};        % no 'multiple conditions'
    % motion regressors from realignment parameters
    Q = nan(nVolsPerSession, 6);
    for i = 1 : nVolsPerSession
        qq = spm_imatrix(mat(:, :, vi(i)) / mat(:, :, vi(1)));
        Q(i, :) = qq(1 : 6);
    end
    for mi = 1 : 6
        fmri_spec.sess(si).regress(mi).name = sprintf('motion(%d)', mi);
        fmri_spec.sess(si).regress(mi).val = Q(:, mi);
    end
    fmri_spec.sess(si).multi_reg = {''};    % no 'multiple regressors'
    fmri_spec.sess(si).hpf = 128;           % high-pass cutoff time (def.)
end
fmri_spec.bases.hrf.derivs = [0 0]; % no HRF derivatives
fmri_spec.volt = 1;                 % do not model interactions
fmri_spec.global = 'None';          % no global intensity normalization
fmri_spec.mthresh = 0.8;            % masking threshold (def.)
fmri_spec.mask = {''};              % no explicit mask
fmri_spec.cvi = 'AR(1)';            % model serial correlations as AR(1) (def.)
% run module
matlabbatch = {};
matlabbatch{1}.spm.stats.fmri_spec = fmri_spec;
spm_jobman('run', matlabbatch(1))

%%

% ## Visualize the model

load(fullfile(modelDir, 'SPM.mat'), 'SPM')
filenames = SPM.xY.P(:, numel(pwd) + 2 : end);
filenames = reshape(cellstr(filenames), size(SPM.xY.VY));
spm_DesRep('DesMtx', SPM.xX, filenames, SPM.xsDes);

%%

% ## Estimate the model

% SPM module 'fMRI model estimation'
fmri_est = struct;
fmri_est.spmmat = {fnSPM};          % path to `SPM.mat`
fmri_est.write_residuals = 0;       % do not write residuals
fmri_est.method.Classical = 1;      % frequentist statistical estimation
% run module
matlabbatch = {};
matlabbatch{1}.spm.stats.fmri_est = fmri_est;
spm_jobman('run', matlabbatch(1))

%%

% The toolbox only needs the updates to the `SPM.mat` file made during
% estimation, not the parameter estimates themselves. We therefore delete the
% created `beta_####.nii` files.

load(fnSPM, 'SPM');
for i = 1 : numel(SPM.Vbeta)
    spm_unlink(fullfile(modelDir, SPM.Vbeta(i).fname))
end

%%

% <!-- Copyright © 2023 Carsten Allefeld
% SPDX-License-Identifier: GPL-3.0-or-later -->

