function second_level = spm_cfg_dcm_peb
% SPM Configuration file for second-level DCM (PEB)
%__________________________________________________________________________

% Peter Zeidman
% Copyright (C) 2016-2022 Wellcome Centre for Human Neuroimaging


%==========================================================================
% Directory / filename selection
%==========================================================================

%--------------------------------------------------------------------------
% dir Directory
%--------------------------------------------------------------------------
dir         = cfg_files;
dir.tag     = 'dir';
dir.name    = 'Directory';
dir.help    = {'Select the directory where the output will be written.'};
dir.filter  = 'dir';
dir.ufilter = '.*';
dir.num     = [1 1];

%--------------------------------------------------------------------------
% name Model name
%--------------------------------------------------------------------------
name         = cfg_entry;
name.tag     = 'name';
name.name    = 'Name';
name.help    = {'Specify a name for the output.'};
name.strtype = 's';
name.num     = [0 Inf];

%==========================================================================
% DCM / PEB file selection
%==========================================================================

%--------------------------------------------------------------------------
% dcmmat Select DCM_*.mat
%--------------------------------------------------------------------------
dcmmat         = cfg_files;
dcmmat.tag     = 'dcmmat';
dcmmat.name    = 'Select DCM files';
dcmmat.help    = {'Select DCM_*.mat files.'};
dcmmat.filter  = 'mat';
dcmmat.ufilter = '^DCM_.*\.mat$';
dcmmat.num     = [1 Inf];

%--------------------------------------------------------------------------
% subj Create single subject
%--------------------------------------------------------------------------
subj      = cfg_branch;
subj.tag  = 'subj';
subj.name = 'Subject';
subj.val  = {dcmmat};
subj.help = {'Subject with one or more models.'};
    
%--------------------------------------------------------------------------
% peb_mat Select PEB_*.mat
%--------------------------------------------------------------------------
peb_mat         = cfg_files;
peb_mat.tag     = 'peb_mat';
peb_mat.name    = 'Select PEB file';
peb_mat.help    = {'Select PEB_*.mat file.'};
peb_mat.filter  = 'mat';
peb_mat.ufilter = '^PEB_.*\.mat$';
peb_mat.num     = [1 1];

%--------------------------------------------------------------------------
% peb_mat Select model_space_*.mat
%--------------------------------------------------------------------------
model_space_mat         = cfg_files;
model_space_mat.tag     = 'model_space_mat';
model_space_mat.name    = 'DCMs';
model_space_mat.help    = {['Select group DCM file (GCM_*.mat). This is a ' ...
                           'cell array with one row per subject and one ' ...
                           'column per DCM.']};
model_space_mat.filter  = 'mat';
model_space_mat.ufilter = '^GCM.*\.mat$';
model_space_mat.num     = [1 Inf];

%==========================================================================
% Covariates entry
%==========================================================================

%--------------------------------------------------------------------------
% design Design matrix
%--------------------------------------------------------------------------
cov_design         = cfg_entry;
cov_design.tag     = 'cov_design';
cov_design.name    = 'Design matrix';
cov_design.help    = {['Enter or paste the N x C design matrix for N ' ...
                      'subjects and C covariates. Note that a column of '...
                      'ones will automatically be added to the start of the '...
                      'matrix, to model the group mean, if one is not found.']};
cov_design.strtype = 'r';
cov_design.num     = [Inf Inf];

%--------------------------------------------------------------------------
% cov_name Name for a column in the design matrix
%--------------------------------------------------------------------------
cov_name         = cfg_entry;
cov_name.tag     = 'name';
cov_name.name    = 'Name';
cov_name.help    = {'Enter a name for a covariate.'};
cov_name.strtype = 's';
cov_name.num     = [0 Inf];

%--------------------------------------------------------------------------
% cov_names Contains the names for the covariates
%--------------------------------------------------------------------------
cov_names         = cfg_repeat;
cov_names.tag     = 'names';
cov_names.name    = 'Covariate names';
cov_names.values  = { cov_name };
cov_names.num     = [0 Inf];
cov_names.help    = {['Enter names for each covariate (excluding the mean ' ...
                     'regressor which is added automatically).']};

%--------------------------------------------------------------------------
% design_mtx Specify whole design matrix 
%--------------------------------------------------------------------------
design_mtx         = cfg_branch;
design_mtx.tag     = 'design_mtx';
design_mtx.name    = 'Specify design matrix';
design_mtx.val     = { cov_design cov_names };
design_mtx.help    = {'Specify the second-level design matrix.'};

%--------------------------------------------------------------------------
% cov_val Value
%--------------------------------------------------------------------------
cov_val         = cfg_entry;
cov_val.tag     = 'value';
cov_val.name    = 'Value';
cov_val.help    = {'Enter the vector of regressor values, one element ' ...
                   'per subject.'};
cov_val.strtype = 'r';
cov_val.num     = [Inf 1];

%--------------------------------------------------------------------------
% covariate A single covariate
%--------------------------------------------------------------------------
regressor         = cfg_branch;
regressor.tag     = 'regressor';
regressor.name    = 'Covariate';
regressor.val     = {cov_name cov_val };
regressor.help    = {'Regressor.'};

%--------------------------------------------------------------------------
% covariate Specify design matrix per covariate
%--------------------------------------------------------------------------
regressors         = cfg_repeat;
regressors.tag     = 'regressors';
regressors.name    = 'Specify covariates individually';
regressors.values  = { regressor };
regressors.help    = {'Specify the second-level design matrix one '...
                     'covariate (regressor) at a time. Note that a ' ...
                     'column of ones to model the mean across subjects ' ...
                     'is added automatically if one is not found.'};
regressors.num     = [1 Inf];

%--------------------------------------------------------------------------
% none No covariates
%--------------------------------------------------------------------------
cov_none         = cfg_branch;
cov_none.tag     = 'none';
cov_none.name    = 'None';
cov_none.val     = {};
cov_none.help    = {'Include no covariates (only the group mean for each '...
                'connection).'};

%--------------------------------------------------------------------------
% covariates Covariates branch
%--------------------------------------------------------------------------
covariates         = cfg_choice;
covariates.tag     = 'cov';
covariates.name    = 'Covariates';
covariates.values  = { cov_none design_mtx regressors };
covariates.help    = {['Specify between-subjects effects (covariates). The ' ...
                      'covariates may be entered all at once as a design ' ...
                      'matrix or individually. Note that if performing ' ...
                      'bayesian model comparison, only the first covariate ' ...
                      'will be treated as being of experimental interest.'] '' ...                      
                      ['Each parameter in the estimated PEB model will ' ...
                      'represent the influence of a covariate on a ' ...
                      'connection. If none is set, only the group mean will '...
                      'be estimated.']};
covariates.val    = {cov_none};
    
%==========================================================================
% PEB fields selection
%==========================================================================

%--------------------------------------------------------------------------
% field_default Select fields A,B
%--------------------------------------------------------------------------
field_default      = cfg_const;
field_default.tag  = 'default';
field_default.name = 'A- and B-matrix';
field_default.val  = {{'A','B'}};
field_default.help = {'A- and B-matrix.'};

%--------------------------------------------------------------------------
% field_all Select all fields
%--------------------------------------------------------------------------
field_all      = cfg_const;
field_all.tag  = 'all';
field_all.name = 'All';
field_all.val  = {'All fields'};
field_all.help = {'All fields.'};

%--------------------------------------------------------------------------
% field_entry Custom field entry
%--------------------------------------------------------------------------
field_entry  = cfg_entry;
field_entry.name = 'Enter manually';
field_entry.tag  = 'custom';
field_entry.help = {['Enter the fields as a cell array e.g. {''A''} ' ...
                     'or {''A'', ''C''} or {''B(:,:,1)'', ''B(:,:,3)''}']};
field_entry.strtype = 'e';
field_entry.num     = [0 Inf];

%--------------------------------------------------------------------------
% fields DCM fields to include
%--------------------------------------------------------------------------
fields        = cfg_choice;
fields.tag    = 'fields';
fields.name   = 'Fields';
fields.values = {field_default field_all field_entry};
fields.help   = {'Select the fields of the DCM to include in the model.' '' ...
                  'A- and B-matrix: Includes all A- and B- connections' ...
                  'All: Includes all fields' ...
                  'Enter manually: Enter a cell array e.g. {''A'',''C''}'};
fields.val    = {field_default};

%==========================================================================
% DCM model index selection
%==========================================================================

%--------------------------------------------------------------------------
% dcm_all Select all DCMs
%--------------------------------------------------------------------------
dcm_all      = cfg_const;
dcm_all.tag  = 'all';
dcm_all.name = 'All';
dcm_all.val  = {'All DCMs'};
dcm_all.help = {'All DCMs.'};

%--------------------------------------------------------------------------
% dcm_idx Single DCM index selection
%--------------------------------------------------------------------------
dcm_sel_idx  = cfg_entry;
dcm_sel_idx.name = 'Selected DCM index';
dcm_sel_idx.tag  = 'index';
dcm_sel_idx.help = {['Select index of the DCM (within subject) on which to ' ...
                 'to build the PEB - e.g. 1 to use each subject''s ' ...
                 'first DCM.']};
dcm_sel_idx.strtype = 'r';
dcm_sel_idx.num     = [1 Inf];

%--------------------------------------------------------------------------
% dcms Which DCMs to include
%--------------------------------------------------------------------------
dcm_idx_1     = dcm_sel_idx;
dcm_idx_1.val = {1};

dcm_idx        = cfg_choice;
dcm_idx.tag    = 'dcm';
dcm_idx.name   = 'DCM index';
dcm_idx.values = {dcm_idx_1 dcm_all};
dcm_idx.help   = {['If each subject has multiple DCMs, select which DCM to use ' ...
                  'for this analysis (or select all). For a standard PEB ' ... 
                  'analysis, it is recommended to specify the PEB over the ' ...
                  'first DCM only, then compare models at the second level.']};
dcm_idx.val    = {dcm_idx_1};
           
%==========================================================================
% Priors on log precision (between-subjects variability) entry
%==========================================================================

%--------------------------------------------------------------------------
% show_review Select whether to review results
%--------------------------------------------------------------------------
components = cfg_menu;
components.tag    = 'components';
components.name   = 'Precision components';
components.labels = {'Single','Fields','All','None'};
components.values = {'Single','Fields','All','None'};
components.val    = {'All'};
components.help   = {['The precision components to include. By default, ' ...
                     'the between-subjects variability is separately ' ...
                     'estimated for each DCM connectivity parameter.'], ...
                     'Single: a single component for all DCM parameters. ', ...
                     'Fields: one component per field (e.g. A,B,C)' ...
                     'All (default): one component per DCM parameter', ...
                     'None: all DCM parameters are treated as fixed effects.'};

%--------------------------------------------------------------------------
% priors_log_precision_mu Priors on log precision expectation
%--------------------------------------------------------------------------
priors_log_precision_mu      = cfg_entry;
priors_log_precision_mu.name = 'Expectation';
priors_log_precision_mu.tag  = 'expectation';
priors_log_precision_mu.help = {['Prior expectation of the log precision ' ...
      'parameters (M.hE), which scale each precision component. The default ' ...
      'is 0, which is translated internally to exp(0) = 1.']};
priors_log_precision_mu.strtype = 'r';
priors_log_precision_mu.num     = [1 1];
priors_log_precision_mu.val     = {0};

%--------------------------------------------------------------------------
% priors_log_precision_var Priors on log precision variance
%--------------------------------------------------------------------------
priors_log_precision_var      = cfg_entry;
priors_log_precision_var.name = 'Uncertainty';
priors_log_precision_var.tag  = 'var';
priors_log_precision_var.help = {['Uncertainty over the prior expectation ' ...
    'of the log precision parameters (M.hC), which scale each precision ' ...
    'component. The default is 1/16.']};
priors_log_precision_var.strtype = 'r';
priors_log_precision_var.num     = [1 1];
priors_log_precision_var.val     = {1/16};

%--------------------------------------------------------------------------
% group priors_parameters_ratio Priors on log precision variance
%--------------------------------------------------------------------------
group_priors_parameters_ratio      = cfg_entry;
group_priors_parameters_ratio.name = 'Group ratio';
group_priors_parameters_ratio.tag  = 'group_ratio';
group_priors_parameters_ratio.help = {['The group ratio (M.alpha) expresses ' ...
     'our prior uncertainty about second-level (group ' ...
     'level GLM) parameters. The default is 1, meaning our uncertainty about ' ...
     'the connection strengths at the second level is the same as our ' ...
     'uncertainty about the connection strengths at the first level. '] '' ...
     ['Internally, the prior covariance of second level parameters is set ' ...
     'to DCM{1}.M.pC / alpha, where DCM{1} is the first DCM provided.']};
group_priors_parameters_ratio.strtype = 'r';
group_priors_parameters_ratio.num     = [1 1];
group_priors_parameters_ratio.val     = {1};

%--------------------------------------------------------------------------
% priors_parameters_ratio Priors on log precision variance
%--------------------------------------------------------------------------
priors_parameters_ratio  = cfg_entry;
priors_parameters_ratio.name = 'Within:between ratio';
priors_parameters_ratio.tag  = 'ratio';
priors_parameters_ratio.help = {'Within:between variance ratio (M.beta). ' ...
    ['This ratio controls the expected between-subjects ' ...
    'variability for each connection (DCM parameter). The default is 16, ' ...
    'meaning we expect the variability in connection strengths across ' ...
    'subjects to be 1/16 of our uncertainty about connection strengths at ' ...
    'the first level.'] '' ...
    ['Internally, the prior second level covariance of DCM parameters (M.pC) ' ...
    'is set to DCM{1}.M.pC / M.beta, where DCM{1} is the first DCM provided.']};
priors_parameters_ratio.strtype = 'r';
priors_parameters_ratio.num     = [1 1];
priors_parameters_ratio.val     = {16};


%--------------------------------------------------------------------------
% priors_between Priors on log precision branch
%--------------------------------------------------------------------------
priors_between      = cfg_branch;
priors_between.tag  = 'priors_between';
priors_between.name = 'Between-subjects variability';
priors_between.val  = { components ...
                        priors_parameters_ratio ...
                        priors_log_precision_mu ...
                        priors_log_precision_var};
priors_between.help = {['Between-subjects variability over second-' ...
     'level parameters.'], '' ...
     ['A multi-component model is used. Each component is a ' ...
      '[p x p] precision matrix given p DCM parameters, where elements on ' ...
      'the diagonal represent the precision (inverse variance) across ' ...
      'subjects of each DCM connection. These precisions are set via the ' ...
      'Within:between ratio, below. Each precision component is scaled by ' ...
      'a hyper-parameter, which is estimated from the data. The prior ' ...
      'expectation and uncertainty of these hyper-parameters are also '...
      'set below.']};
  

%--------------------------------------------------------------------------
% priors_third Priors on log precision branch
%--------------------------------------------------------------------------
priors_glm         = cfg_branch;
priors_glm.tag     = 'priors_glm';
priors_glm.name    = 'Second level (GLM) priors';
priors_glm.val     = { group_priors_parameters_ratio};
priors_glm.help    = {['Priors on the expected values of the second ' ...
                       'level General Linear Model (GLM) parameters.']};

%--------------------------------------------------------------------------
% maxit Maximum VL iterations
%--------------------------------------------------------------------------
maxit      = cfg_entry;
maxit.name = 'Max iterations';
maxit.tag  = 'maxit';
maxit.help = {'Maximum iterations of the VL estimation procedure.'
              ['Only change this if you find the free energy is still ' ...
               'increasing non-trivially up to the default limit of 64.']};
maxit.strtype = 'r';
maxit.num     = [1 1];
maxit.val     = {64};

%--------------------------------------------------------------------------
% estimation Model estimation settings
%--------------------------------------------------------------------------
estimation         = cfg_branch;
estimation.tag     = 'estimation';
estimation.name    = 'Estimation';
estimation.val     = { maxit };
estimation.help    = {'Settings for estimating the PEB model'};

%--------------------------------------------------------------------------
% show_review Select whether to review results
%--------------------------------------------------------------------------
show_review        = cfg_menu;
show_review.tag    = 'show_review';
show_review.name   = 'Review PEB parameters';
show_review.labels = {'Yes','No'};
show_review.values = {1,0};
show_review.val    = {1};
show_review.help   = {'Review PEB parameters'};

%==========================================================================
% PEB specification batch
%==========================================================================
% Set show review default to off
sr      = show_review;
sr.val  = {0};

specify      = cfg_exbranch;
specify.tag  = 'specify';
specify.name = 'Specify / Estimate PEB';
specify.val  = { name model_space_mat dcm_idx covariates fields ...
                 priors_between priors_glm estimation sr };
specify.help = {['Specifies and estimates a second-level DCM (PEB) model. ' ...
                 'A PEB model will be created for each first level DCM.' ]};
            
specify.prog = @spm_run_create_peb;
specify.vout = @vout_peb;

%==========================================================================
% PEB reduce / average / compare batch
%==========================================================================
model_space_mat_op = model_space_mat;
model_space_mat_op.num = [0 Inf];
model_space_mat_op.val = {''};

null_prior_covariance  = cfg_entry;
null_prior_covariance.name = 'Null prior variance';
null_prior_covariance.tag  = 'nullpcov';
null_prior_covariance.help = {['This is the prior variance for any ''switched ' ...
    'off'' parameter. Setting this to 0 (zero) tests the point null ' ...
    'hypothesis that parameters are present vs absent. Setting this to a ' ...
    'small positive number tests the null hypothesis that the parameter is ' ...
    'trivially small. The smaller this setting, the fewer parameters will ' ...
    'survive the model search. (See setting gamma in spm_dcm_bmr_all.m)']};
null_prior_covariance.strtype = 'r';
null_prior_covariance.num     = [1 1];
null_prior_covariance.val     = {1/16};

peb_compare      = cfg_exbranch;
peb_compare.tag  = 'compare';
peb_compare.name = 'Compare / Average PEB models';
peb_compare.val  = { peb_mat model_space_mat_op show_review};
peb_compare.help = {['Addresses the question: which combination of ' ...
    'connections best explains the commonalities across subjects and ' ...
    'the group differences between subjects?'] '' ...
    ['If a group difference is to be investigated, this should be the ' ...
     'first manually entered between-subjects covariate in the PEB.'] ''  ...
    ['This analysis is performed by comparing a PEB model to nested ' ...
    'sub-models where ' ...
    'certain parameters have been disabled (fixed at their prior mean of ' ...
    'zero). Parameters are then averaged over reduced models to give an '...
    'averaged PEB (referred to as a Bayesian Model Average, BMA). Each ' ...
    'parameter in the PEB or BMA represents the effect of one between-' ...
    'subjects covariate on one connection.'] '' ...
    ['If only one first-level DCM is provided per subject, a search is made ' ...
    'over reduced PEB models to prune away any parameters not contributing ' ... 
    'to the model evidence. If multiple DCMs are provided per subject, ' ...
    'these are used to define the combinations of second level parameters ' ...
    'to be compared.']};
peb_compare.prog = @spm_run_bmc;
peb_compare.vout = @vout_bma;

reduce_all       = cfg_exbranch;
reduce_all.tag   = 'reduce_all';
reduce_all.name  = 'Search nested PEB models';
reduce_all.val   = { peb_mat model_space_mat_op null_prior_covariance show_review};
reduce_all.help  = {['Optimises a PEB model by trying different ' ...
                         'combinations of switching off parameters (fixing ' ...
                         'them at their prior value), where doing so does ' ...
                         'not reduce the model evidence. Any parameters ' ...
                         'not contributing will be set to zero.'] '' ...
                         ['This is equivalent to the function of '... 
                         'spm_dcm_post_hoc but on the second level (PEB) ' ...
                         'parameters.']};
reduce_all.prog = @spm_run_bmr_all;
reduce_all.vout = @vout_bma;

%==========================================================================
% PEB review batch
%==========================================================================
review      = cfg_exbranch;
review.tag  = 'peb_review';
review.name = 'Review PEB';
review.val  = { peb_mat model_space_mat_op };
review.help = {'Reviews PEB results'};
review.prog = @spm_run_dcm_peb_review;

%==========================================================================
% PREDICT leave-one-out cross validation
%==========================================================================

covariates_min1 = covariates;
covariates_min1.val{1} = covariates_min1.values{2};
covariates_min1.values(1) = [];

predict      = cfg_exbranch;
predict.tag  = 'predict';
predict.name = 'Predict (cross-validation)';
predict.val  = { name model_space_mat dcm_idx covariates_min1 fields ...
                 priors_between priors_glm estimation};
predict.help = {['Builds a PEB model on all but one subjects, and uses ' ...
                 'it to predict a between-subjects effect, such as ' ...
                 'group membership, in the remaining subject. This process ' ...
                 'is repeated, leaving out each subject in turn ' ...
                 '(leave-one-out cross-validation) to establish the cross- ' ...
                 'validation accuracy of the model.'] ...
                 ['The first covariate (after the automatically inserted ' ...
                  'mean regressor) is used as the predictor variable.']};            
predict.prog = @spm_run_dcm_loo;
predict.vout = @vout_loo;

%==========================================================================
% second_level Second level DCM batch
%==========================================================================
second_level         = cfg_choice; 
second_level.tag     = 'peb';
second_level.name    = 'Second level';
second_level.help    = {'Parametric Empirical Bayes for DCM.'};
second_level.values  = { specify peb_compare reduce_all review predict };


%==========================================================================
function dep = vout_bma(varargin)
%==========================================================================
dep(1)            = cfg_dep;
dep(1).sname      = 'BMA mat filename';
dep(1).src_output = substruct('.','bmamat');
dep(1).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});

%==========================================================================
function dep = vout_loo(varargin)
%==========================================================================
dep(1)            = cfg_dep;
dep(1).sname      = 'LOO mat File(s)';
dep(1).src_output = substruct('.','loo_mat');
dep(1).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});

%==========================================================================
function dep = vout_peb(varargin)
%==========================================================================
dep(1)            = cfg_dep;
dep(1).sname      = 'PEB mat File(s)';
dep(1).src_output = substruct('.','peb_mat');
dep(1).tgt_spec   = cfg_findspec({{'filter','mat','strtype','e'}});


%==========================================================================
function out = spm_run_dcm_peb_review(job)
%==========================================================================
% Run the PEB review batch
P   = job.peb_mat;
DCM = job.model_space_mat;
spm_dcm_peb_review(P{1},DCM);
out = job.peb_mat;

%==========================================================================
function out = spm_run_create_peb(job)
%==========================================================================
% Run the PEB specification / estimation batch

[GCM,M,field,dir_out] = prepare_peb_inputs(job);
    
% Specify / estimate PEB on full model only
name    = job.name;    
PEB     = spm_dcm_peb(GCM,M,field);

% Write PEB
peb_filename = fullfile(dir_out,['PEB_' name '.mat']);
save(peb_filename,'PEB', spm_get_defaults('mat.format'));

out.peb_mat = {peb_filename};

% Review PEB
if job.show_review == 1
    spm_dcm_peb_review(peb_filename,GCM);
end

%==========================================================================
function out = spm_run_dcm_loo(job)
%==========================================================================
% Run leave-one-out cross validation

[GCM,M,field,dir_out] = prepare_peb_inputs(job);

[qE,qC,Q] = spm_dcm_loo(GCM,M,field);

% Write output
name         = job.name;    
loo_filename = fullfile(dir_out,['LOO_' name '.mat']);

save(loo_filename,'qE','qC','Q', spm_get_defaults('mat.format'));

out.loo_mat = {loo_filename};

%==========================================================================
function [GCM,M,field,dir_out] = prepare_peb_inputs(job)
%==========================================================================
% Prepare the inputs needed to specify a PEB or run LOO

[GCM,dir_out] = load_dcm(job);

ns = size(GCM,1);

if ~isfield(GCM{1},'Ep')
    error('Please estimate DCMs before second-level analysis.');
end

% DCM field(s)
if isfield(job.fields,'default')
    field = {'A','B'};
elseif isfield(job.fields,'all')
    field = 'all';
else
    field = job.fields.custom;
end

% Covariates
if isfield(job.cov, 'none')
    % No covariates (except the mean, added later)
    Xnames = {};
    X = [];
           
elseif isfield(job.cov, 'design_mtx')
    % Whole design matrix entered
    x = job.cov.design_mtx.cov_design;
    
    if size(x,1) ~= ns
        error('Please ensure design matrix has one row per subject.');
    end
    
    X = x;
    Xnames = job.cov.design_mtx.name;    
                
elseif isfield(job.cov, 'regressor')
    % Design matrix entered per-regressor
    
    regressors = job.cov.regressor;
    
    X = [];
    Xnames = {};
    
    for r = 1:length(regressors)
        regressor = regressors(r).value;
        name      = regressors(r).name;
        
        if size(regressor,1) ~= ns
            error('Please ensure regressor %d has one row per subject.',r);
        end
                
        X = [X regressor];
        Xnames = [Xnames name];
    end
end

% Ensure there's a constant column
if isempty(X)
    has_constant = false;
else
    has_constant = ~any(diff(X(:,1)));
end
if ~has_constant
    Xnames = ['Commonalities' Xnames];
    X = [ones(ns,1) X];
end

% Ensure a column colinear with the mean wasn't entered accidentally
if size(X,2) > 1
    bad = find(~any(diff(X(:,2:end))));
    if ~isempty(bad)
        error('Please check regressor %d.', bad);
    end
end

if size(X,2) ~= length(Xnames)
    error('Please ensure there is one covariate name per covariate.');
end

% Option for precision components
Q = job.priors_between.components;
if isempty(Q)
    Q = 'all';
else
    Q = lower(Q);
end

% Priors / covariance components
M = struct();
M.alpha  = job.priors_glm.group_ratio;
M.beta   = job.priors_between.ratio;
M.hE     = job.priors_between.expectation;
M.hC     = job.priors_between.var;
M.Q      = Q;
M.X      = X;
M.Xnames = Xnames;
M.maxit  = job.estimation.maxit;

%==========================================================================
function out = spm_run_bmr_all(job)
%==========================================================================
% Run a search over all reduced PEB models and average

GCM = load_dcm(job);
nm  = size(GCM,2);

if nm > 1
    disp('Running search on the full DCM only.');        
    GCM = GCM(:,1);
end

out = run_peb_bmc_internal(job,GCM);

%==========================================================================
function out = spm_run_bmc(job)
%==========================================================================
% Run a model comparison between specified reduced PEB models

GCM = load_dcm(job);
nm  = size(GCM,2);

if nm < 2
    error('More than one DCM per subject is required for BMC.');
end

out = run_peb_bmc_internal(job,GCM);

%==========================================================================
function out = run_peb_bmc_internal(job, GCM)
%==========================================================================
% Search / model comparison both use spm_dcm_peb_bmc, which is called here

PEB = load(job.peb_mat{1});
PEB = PEB.PEB;

if isfield(job,'nullpcov')
    PEB.gamma = job.nullpcov;
end

nm  = size(GCM,2);

% Run BMA on defined reduced models or all submodels
if nm > 1
    BMA = spm_dcm_peb_bmc(PEB,GCM(1,:));
else
    BMA = spm_dcm_peb_bmc(PEB);
end

% Write BMA
[dir_out, name] = fileparts(job.peb_mat{1});
if nm == 1
    name = ['search_' name];
end
filename = fullfile(dir_out, ['BMA_' name '.mat']);
save(filename,'BMA', spm_get_defaults('mat.format'));

out.bmamat = {filename};

% Review BMA
if job.show_review == 1
    
    DCM = job.model_space_mat;
    if ~isempty(DCM)
        DCM = load(DCM{1});
        if isfield(DCM,'GCM')
            DCM = DCM.GCM;
        else
            DCM = DCM.DCM;
        end
    end    
    
    spm_dcm_peb_review(BMA,DCM);
end

%==========================================================================
function [GCM,dir_out] = load_dcm(job)
%==========================================================================
% Load and validate selected model space

if iscell(job.model_space_mat) && ischar(job.model_space_mat{1}) ...
        && numel(job.model_space_mat) > 1
    % An array of DCM filenames was provided (via the dependency system)    
    GCM     = job.model_space_mat;    
    dir_out = fileparts(GCM{1});
elseif isempty(job.model_space_mat)
    GCM = {};
    dir_out = pwd;
    return
else
    % The filename of a single GCM mat file was provided
    gcm_file = char(job.model_space_mat);
    dir_out  = fileparts(gcm_file);
    GCM      = load(gcm_file);
    if ~isfield(GCM,'GCM')
        error('Provided file is not a valid model space.');
    end
    GCM = GCM.GCM;
end

% If a GCM of filenames was provided, load the DCMs
if ischar(GCM{1})
    GCM = spm_dcm_load(GCM);
end

% Limit to specific model(s) if requested
if isfield(job,'dcm') && isfield(job.dcm,'index')
    index = job.dcm.index;
    GCM = GCM(:,index);
end
