Homec4science
Diffusion iCAPs public version 2 (master) Browse

iCAPs_v2/packages/spm12/spm_ADEM_M_set.m
e5f5005bb2efmaster

Last Change
Raw File

spm_ADEM_M_set.m
View Options

function [M] = spm_ADEM_M_set(M)
% Set indices and perform checks on hierarchical action models
% FORMAT [M] = spm_ADEM_M_set(M)
%
% for each level (i); required fields
%
% M(i).g = y(t) = g(x,v,a,P) {inline function, string or m-file}
% M(i).f = dx/dt = f(x,v,a,P) {inline function, string or m-file}
%
% and
%
% M(i).m = number of inputs v(i + 1);
% M(i).n = number of states x(i);
% M(i).l = number of output v(i);
% M(i).k = number of action a(i);
%
% or
%
% M(i).x = hidden states;
% M(i).v = causal states;
% M(i).a = action states;
%
% for each level (i); optional fields
%
% M(i).pE = prior expectation of p model-parameters
% M(i).V = precision (input noise)
% M(i).W = precision (state noise)
% M(i).U = precision (action)
%
%
% sets fields, checks internal consistency of model specification and sets
% estimation parameters. If (V,W) are not specified infinite precision is
% assumed.
%--------------------------------------------------------------------------
%
% M(1).E.s; = smoothness (s.d. in time bins)
% M(1).E.d; = embedding order q(v) (i.e., number of derivatives)
% M(1).E.n; = embedding order q(x)
%
% If the highest level involves any dynamic or static transformation
% of its inputs a further level is added with flat priors
%__________________________________________________________________________
% Copyright (C) 2008-2014 Wellcome Trust Centre for Neuroimaging
% Karl Friston
% $Id: spm_ADEM_M_set.m 5962 2014-04-17 12:47:43Z spm $
% order
%--------------------------------------------------------------------------
g = length(M);
% set missing fields
%==========================================================================
% check for specification of hidden states
%--------------------------------------------------------------------------
if isfield(M,'f') && ~isfield(M,'n') && ~isfield(M,'x')
msgbox('please specify hidden states or their number')
error(' ')
end
% check supra-ordinate level and add one (with flat priors) if necessary
%--------------------------------------------------------------------------
try
fcnchk(M(g).g);
g = g + 1;
M(g).l = M(g - 1).m;
end
M(g).m = 0;
M(g).n = 0;
% default fields for static models (hidden states)
%--------------------------------------------------------------------------
if ~isfield(M,'f')
[M.f] = deal(inline('sparse(0,1)','x','v','a','P'));
[M.x] = deal(sparse(0,1));
[M.n] = deal(0);
end
for i = 1:g
try
fcnchk(M(i).f);
catch
M(i).f = inline('sparse(0,1)','x','v','a','P');
M(i).x = sparse(0,1);
end
end
% consistency and format check on states, parameters and functions
%==========================================================================
% prior expectation of parameters M.pE
%--------------------------------------------------------------------------
try
M.pE;
catch
% Assume fixed parameters
%----------------------------------------------------------------------
for i = 1:g
M(i).pE = sparse(0,0);
end
end
% get inputs
%--------------------------------------------------------------------------
try
v = M(g).v;
catch
v = sparse(0,1);
end
if isempty(v)
try
v = sparse(M(g - 1).m,1);
end
end
if isempty(v)
try
v = sparse(M(g).l,1);
end
end
M(g).l = length(spm_vec(v));
M(g).v = v;
% ensure action is specified
%--------------------------------------------------------------------------
for i = 1:g
try
a = M(i).a;
catch
a = sparse(0,1);
end
if isempty(a)
try
a = sparse(M(i).k,1);
end
end
M(i).k = length(spm_vec(a));
M(i).a = a;
end
% check functions
%--------------------------------------------------------------------------
for i = (g - 1):-1:1
try
x = M(i).x;
catch
x = sparse(M(i).n,1);
end
if isempty(x) && M(i).n
x = sparse(M(i).n,1);
end
% check f(x,v,P)
%----------------------------------------------------------------------
try
M(i).f = fcnchk(M(i).f);
if nargin(M(i).f) ~= 4
M(i).f = inline(char(M(i).f),'x','v','a','P');
end
end
try
f = feval(M(i).f,x,v,a,M(i).pE);
if length(spm_vec(x)) ~= length(spm_vec(f))
errordlg(sprintf('please check nargout: G(%i).f(x,v,a,P)',i));
end
catch
errordlg(sprintf('evaluation failure: G(%i).f(x,v,a,P)',i))
end
% check g(x,v,P)
%----------------------------------------------------------------------
try
M(i).g = fcnchk(M(i).g);
if nargin(M(i).g) ~= 4
M(i).g = inline(char(M(i).g),'x','v','a','P');
end
end
try
M(i).m = length(spm_vec(v));
v = feval(M(i).g,x,v,a,M(i).pE);
a = M(i).a;
M(i).k = length(spm_vec(a));
M(i).l = length(spm_vec(v));
M(i).n = length(spm_vec(x));
M(i).a = a;
M(i).v = v;
M(i).x = x;
catch
errordlg(sprintf('evaluation failure: G(%i).g(x,v,a,P)',i))
end
end
% remove empty levels
%--------------------------------------------------------------------------
try
g = find(~spm_vec(M.m),1);
M = M(1:g);
catch
errordlg('please specify number of variables')
end
% number of x (hidden states)
%--------------------------------------------------------------------------
nx = sum(spm_vec(M.n));
% check precisions
%==========================================================================
try, M.V; catch, M(1).V = []; end
try, M.W; catch, M(1).W = []; end
% check precisions
%--------------------------------------------------------------------------
for i = 1:g
% check V and assume unit precision if improperly specified
%----------------------------------------------------------------------
if isvector(M(i).V), M(i).V = diag(M(i).V); end
if length(M(i).V) ~= M(i).l
try
M(i).V = speye(M(i).l,M(i).l)*M(i).V(1);
catch
M(i).V = speye(M(i).l,M(i).l)*exp(32);
end
end
% check W and assume unit precision if improperly specified
%----------------------------------------------------------------------
if isvector(M(i).W), M(i).W = diag(M(i).W); end
if length(M(i).W) ~= M(i).n
try
M(i).W = speye(M(i).n,M(i).n)*M(i).W(1);
catch
M(i).W = speye(M(i).n,M(i).n)*exp(32);
end
end
end
% check restiction of precision for action (gain)
%==========================================================================
try, M(1).U; catch, M(1).U = []; end
if isvector(M(1).U), M(1).U = diag(M(1).U); end
if length(M(1).U) ~= M(1).l
try
M(1).U = speye(M(1).l,M(1).l)*M(1).U(1);
catch
M(1).U = speye(M(1).l,M(1).l)*exp(2);
end
end
% estimation parameters M(1).E.s, n,...
%==========================================================================
% E.s; % smoothness (seconds)
% E.dt; % time step
% E.d; % approximation order of q(x,v)
% E.n; % order of embedding (n >= d)
% temporal smoothness - s.d. of kernel
%--------------------------------------------------------------------------
try M(1).E.s; catch, if nx, M(1).E.s = 1/2; else M(1).E.s = 0; end, end
% time step
%--------------------------------------------------------------------------
try M(1).E.dt; catch, M(1).E.dt = 1; end
% embedding orders
%--------------------------------------------------------------------------
try M(1).E.d; catch, if nx, M(1).E.d = 2; else M(1).E.d = 0; end, end
try M(1).E.n; catch, if nx, M(1).E.n = 6; else M(1).E.n = 0; end, end
M(1).E.d = min(M(1).E.d,M(1).E.n);
% checks on smoothness hyperparameter
%==========================================================================
for i = 1:g
try, M(i).sv; catch, M(i).sv = M(1).E.s; end
try, M(i).sw; catch, M(i).sw = M(1).E.s; end
if ~isscalar(M(i).sv), M(i).sv = M(1).E.s; end
if ~isscalar(M(i).sw), M(i).sw = M(1).E.s; end
end
c4science · Help