Homec4science
Diffusion iCAPs public (master) Browse

iCAPs/packages/spm12/spm_DEM_eval_diff.m
31022b3a13e6master

Last Change
Raw File

spm_DEM_eval_diff.m
View Options

function [D] = spm_DEM_eval_diff(x,v,qp,M,bilinear)
% evaluates derivatives for DEM schemes
% FORMAT [D] = spm_DEM_eval_diff(x,v,qp,M,bilinear)
% v{i} - casual states
% x(i) - hidden states
% qp - conditional density of parameters
% qp.p{i} - parameter deviates for i-th level
% qp.u(i) - basis set
% qp.x(i) - expansion point ( = prior expectation)
% M - model structure
% bilinear - optional flag to suppress second-order derivatives
%
% D - derivatives
% D.dgdv
% ...
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
% Karl Friston
% $Id: spm_DEM_eval_diff.m 6734 2016-03-02 12:02:46Z peter $
% check for evaluation of bilinear terms
%--------------------------------------------------------------------------
try
bilinear;
catch
bilinear = 1;
end
% get dimensions
%==========================================================================
nl = size(M,2); % number of levels
ne = sum(spm_vec(M.l)); % number of e (errors)
nx = sum(spm_vec(M.n)); % number of x (hidden states)
np = sum(spm_vec(M.p)); % number of p (parameters)
ny = M(1).l; % number of y (inputs)
nc = M(end).l; % number of c (prior causes)
% initialise cell arrays for hierarchical structure
%--------------------------------------------------------------------------
df.dv = cell(nl - 1,nl - 1);
df.dx = cell(nl - 1,nl - 1);
df.dp = cell(nl - 1,nl - 1);
dg.dv = cell(nl ,nl - 1);
dg.dx = cell(nl ,nl - 1);
dg.dp = cell(nl ,nl - 1);
for i = 1:(nl - 1)
dg.dv{i + 1,i} = sparse(M(i).m,M(i).m);
dg.dx{i + 1,i} = sparse(M(i).m,M(i).n);
dg.dp{i + 1,i} = sparse(M(i).m,M(i).p);
dg.dv{i ,i} = sparse(M(i).l,M(i).m);
dg.dx{i ,i} = sparse(M(i).l,M(i).n);
dg.dp{i ,i} = sparse(M(i).l,M(i).p);
df.dv{i ,i} = sparse(M(i).n,M(i).m);
df.dx{i ,i} = sparse(M(i).n,M(i).n);
df.dp{i ,i} = sparse(M(i).n,M(i).p);
end
if bilinear
for i = 1:(nl - 1)
dg.dvp{i} = cell(M(i).p,1);
dg.dxp{i} = cell(M(i).p,1);
df.dvp{i} = cell(M(i).p,1);
df.dxp{i} = cell(M(i).p,1);
[dg.dvp{i}{:}] = deal(dg.dv);
[dg.dxp{i}{:}] = deal(dg.dx);
[df.dvp{i}{:}] = deal(df.dv);
[df.dxp{i}{:}] = deal(df.dx);
end
end
% Derivatives at each hierarchical level
%==========================================================================
% inline function for evaluating projected parameters
%--------------------------------------------------------------------------
h = @(f,x,v,q,u,p) f(x,v,spm_unvec(spm_vec(p) + u*q,p));
for i = 1:(nl - 1)
% states level i
%----------------------------------------------------------------------
xvp = {x{i},v{i},qp.p{i},qp.u{i},M(i).pE};
% 1st and 2nd partial derivatives (states)
%----------------------------------------------------------------------
if bilinear && np
try
[dgdxp, dgdx] = spm_diff(h,M(i).gx,xvp{:},4,'q');
[dgdvp, dgdv] = spm_diff(h,M(i).gv,xvp{:},4,'q');
[dfdxp, dfdx] = spm_diff(h,M(i).fx,xvp{:},4,'q');
[dfdvp, dfdv] = spm_diff(h,M(i).fv,xvp{:},4,'q');
catch
[dgdxp, dgdx] = spm_diff(h,M(i).g,xvp{:},[2 4],'q');
[dgdvp, dgdv] = spm_diff(h,M(i).g,xvp{:},[3 4],'q');
[dfdxp, dfdx] = spm_diff(h,M(i).f,xvp{:},[2 4],'q');
[dfdvp, dfdv] = spm_diff(h,M(i).f,xvp{:},[3 4],'q');
end
else
try
dgdx = h(M(i).gx,xvp{:});
dgdv = h(M(i).gv,xvp{:});
dfdx = h(M(i).fx,xvp{:});
dfdv = h(M(i).fv,xvp{:});
catch
dgdx = spm_diff(h,M(i).g,xvp{:},2);
dgdv = spm_diff(h,M(i).g,xvp{:},3);
dfdx = spm_diff(h,M(i).f,xvp{:},2);
dfdv = spm_diff(h,M(i).f,xvp{:},3);
end
end
% 1st-order partial derivatives (parameters)
%----------------------------------------------------------------------
try
dfdp = h(M(i).fp,xvp{:});
dgdp = h(M(i).gp,xvp{:});
catch
dfdp = spm_diff(h,M(i).f,xvp{:},4);
dgdp = spm_diff(h,M(i).g,xvp{:},4);
end
% % check which dervatives need to be evaluated
% %====================================================================
% D(i).dgdv = nnz(dgdv) + nnz(spm_vec(dgdvp));
% D(i).dgdx = nnz(dgdx) + nnz(spm_vec(dgdxp));
% D(i).dfdv = nnz(dfdv) + nnz(spm_vec(dfdvp));
% D(i).dfdx = nnz(dfdx) + nnz(spm_vec(dfdxp));
% Constant terms (linking causes over levels)
%----------------------------------------------------------------------
dg.dv{i + 1,i} = -speye(M(i).m,M(i).m);
% place 1st derivatives in array
%----------------------------------------------------------------------
dg.dx{i,i} = dgdx;
dg.dv{i,i} = dgdv;
df.dx{i,i} = dfdx;
df.dv{i,i} = dfdv;
df.dp{i,i} = dfdp;
dg.dp{i,i} = dgdp;
% place 2nd derivatives in array
%----------------------------------------------------------------------
if bilinear && np
for j = 1:length(dgdxp)
dg.dxp{i}{j}{i,i} = dgdxp{j};
dg.dvp{i}{j}{i,i} = dgdvp{j};
df.dxp{i}{j}{i,i} = dfdxp{j};
df.dvp{i}{j}{i,i} = dfdvp{j};
end
end
end
% concatenate hierarchical forms
%==========================================================================
D.dgdv = spm_cat(dg.dv);
D.dgdx = spm_cat(dg.dx);
D.dfdv = spm_cat(df.dv);
D.dfdx = spm_cat(df.dx);
D.dfdp = spm_cat(df.dp);
D.dgdp = spm_cat(dg.dp);
% fixed derivatives w.r.t. prediction errors and states
%--------------------------------------------------------------------------
D.dfdy = sparse(nx,ny);
D.dfdc = sparse(nx,nc);
D.dedy = spm_speye(ne,ny);
D.dedc = -spm_speye(ne,nc,nc - ne);
% bilinear terms if required
%--------------------------------------------------------------------------
if bilinear
D.dgdvp = {};
D.dgdxp = {};
D.dfdvp = {};
D.dfdxp = {};
for i = 1:length(dg.dvp)
for j = 1:length(dg.dvp{i})
D.dgdvp{end + 1} = spm_cat(dg.dvp{i}{j});
D.dgdxp{end + 1} = spm_cat(dg.dxp{i}{j});
D.dfdvp{end + 1} = spm_cat(df.dvp{i}{j});
D.dfdxp{end + 1} = spm_cat(df.dxp{i}{j});
end
end
end
c4science · Help