Homec4science
Diffusion iCAPs public (master) Browse

iCAPs/packages/spm12/spm_DEM_int.m
31022b3a13e6master

Last Change
Raw File

spm_DEM_int.m
View Options

function [V,X,Z,W] = spm_DEM_int(M,z,w,c)
% Integrates/evaluates a hierarchical model given innovations z{i} and w{i}
% FORMAT [V,X,Z,W] = spm_DEM_int(M,z,w,c);
%
% M{i} - model structure
% z{i} - innovations (causes)
% w{i} - innovations (states)
% c{i} - exogenous causes
%
% V{i} - causal states (V{1} = y = response)
% X{i} - hidden states
% Z{i} - fluctuations (causes)
% W{i} - fluctuations (states)
%
% The system is evaluated at the prior expectation of the parameters
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
% Karl Friston
% $Id: spm_DEM_int.m 6132 2014-08-06 19:59:46Z karl $
% set model indices and missing fields
%--------------------------------------------------------------------------
M = spm_DEM_M_set(M);
% innovations
%--------------------------------------------------------------------------
z = spm_cat(z(:)) + spm_cat(c(:));
w = spm_cat(w(:));
% number of states and parameters
%--------------------------------------------------------------------------
nt = size(z,2); % number of time steps
nl = size(M,2); % number of levels
nv = sum(spm_vec(M.l)); % number of v (casual states)
nx = sum(spm_vec(M.n)); % number of x (hidden states)
% order parameters (n= 1 for static models)
%==========================================================================
dt = M(1).E.dt; % time step
n = M(1).E.n + 1; % order of embedding
nD = M(1).E.nD; % number of iterations per sample
td = dt/nD; % integration time for D-Step
% initialize cell arrays for derivatives z{i} = (d/dt)^i[z], ...
%--------------------------------------------------------------------------
u.v = cell(n,1);
u.x = cell(n,1);
u.z = cell(n,1);
u.w = cell(n,1);
[u.v{:}] = deal(sparse(nv,1));
[u.x{:}] = deal(sparse(nx,1));
[u.z{:}] = deal(sparse(nv,1));
[u.w{:}] = deal(sparse(nx,1));
% hyperparameters
%--------------------------------------------------------------------------
ph.h = {M.hE};
ph.g = {M.gE};
% initialize with starting conditions
%--------------------------------------------------------------------------
vi = {M.v};
xi = {M.x};
u.v{1} = spm_vec(vi);
u.x{1} = spm_vec(xi);
% derivatives for Jacobian of D-step
%--------------------------------------------------------------------------
Dx = kron(spm_speye(n,n,1),spm_speye(nx,nx,0));
Dv = kron(spm_speye(n,n,1),spm_speye(nv,nv,0));
D = spm_cat(spm_diag({Dv,Dx,Dv,Dx}));
dfdw = kron(eye(n,n),eye(nx,nx));
% initialize conditional estimators of states to be saved (V and X)
%--------------------------------------------------------------------------
mnx = 0;
mnv = 0;
for i = 1:nl
V{i} = sparse(M(i).l,nt);
X{i} = sparse(M(i).n,nt);
Z{i} = sparse(M(i).l,nt);
W{i} = sparse(M(i).n,nt);
% check for state-dependent precision
%----------------------------------------------------------------------
mnx = mnx | length(M(i).pg);
mnv = mnv | length(M(i).ph);
end
% defaults for state-dependent precision
%--------------------------------------------------------------------------
Sz = 1;
Sw = 1;
% iterate over sequence (t) and within for static models
%==========================================================================
for t = 1:nt
for iD = 1:nD
% Get generalised motion of random fluctuations
%==================================================================
% sampling time
%------------------------------------------------------------------
ts = (t + (iD - 1)/nD)*dt;
% evaluate state-dependent precision
%------------------------------------------------------------------
if mnx || mnv
vi{nl} = vi{nl} + c{nl}(:,t);
pu.x = {spm_vec(xi(1:end - 1))};
pu.v = {spm_vec(vi(1 + 1:end))};
p = spm_LAP_eval(M,pu,ph);
if mnv, Sz = sparse(diag(exp(-p.h/2))); end
if mnx, Sw = sparse(diag(exp(-p.g/2))); end
end
% derivatives of innovations (and exogenous input)
%------------------------------------------------------------------
u.z = spm_DEM_embed(Sz*z,n,ts,dt);
u.w = spm_DEM_embed(Sw*w,n,ts,dt);
% Evaluate and update states
%==================================================================
% evaluate functions
%------------------------------------------------------------------
[u,dg,df] = spm_DEM_diff(M,u);
% tensor products for Jacobian
%------------------------------------------------------------------
dgdv = kron(spm_speye(n,n,1),dg.dv);
dgdx = kron(spm_speye(n,n,1),dg.dx);
dfdv = kron(spm_speye(n,n,0),df.dv);
dfdx = kron(spm_speye(n,n,0),df.dx);
% Save realization
%==================================================================
vi = spm_unvec(u.v{1},{M.v});
xi = spm_unvec(u.x{1},{M.x});
zi = spm_unvec(u.z{1},{M.v});
wi = spm_unvec(u.w{1},{M.x});
if iD == 1
for i = 1:nl
if M(i).l, V{i}(:,t) = spm_vec(vi{i}); end
if M(i).n, X{i}(:,t) = spm_vec(xi{i}); end
if M(i).l, Z{i}(:,t) = spm_vec(zi{i}); end
if M(i).n, W{i}(:,t) = spm_vec(wi{i}); end
end
end
% break for static models
%------------------------------------------------------------------
if nt == 1, break, end
% Jacobian for update
%------------------------------------------------------------------
J = spm_cat({dgdv dgdx Dv [] ;
dfdv dfdx [] dfdw;
[] [] Dv [] ;
[] [] [] Dx});
% update states u = {x,v,z,w}
%------------------------------------------------------------------
du = spm_dx(J,D*spm_vec(u),td);
% and unpack
%------------------------------------------------------------------
u = spm_unvec(spm_vec(u) + du,u);
end % iterations over iD
end % iterations over t
c4science · Help