Addition of ADEM routines (active inference; i.e. DEM with active sam…

…pling)- alpha version SVN r1887
spm · Jul 4, 2008 · 41c3f3d · 41c3f3d
1 parent 20de94f
commit 41c3f3d
Show file tree

Hide file tree

Showing 9 changed files with 2,392 additions and 1 deletion.
diff --git a/spm_ADEM.m b/spm_ADEM.m
diff --git a/spm_ADEM_M_set.m b/spm_ADEM_M_set.m
diff --git a/spm_ADEM_diff.m b/spm_ADEM_diff.m
@@ -0,0  1,106 @@
 function [u dg df] = spm_ADEM_diff(M,u);
 % evaluates an active model given innovations z{i} and w{i}
 % FORMAT [u dgdv dgdx dfdv dfdx] = spm_ADEM_diff(M,u);
 %
 % M    - generative model
 %
 % u.a - active states
 % u.v - causal states - updated
 % u.x - hidden states - updated
 % u.z - innovation (causal state)
 % u.w - innovation (hidden states)
 %
 % dg.dv, ...  components of the Jacobian in generalised coordinates
 %
 % The system is evaluated at the prior expectation of the parameters
 %__________________________________________________________________________
 % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
 
 % Karl Friston
 % $Id: spm_ADEM_diff.m 1887 2008-07-04 17:48:42Z karl $
 
 % number of states and parameters
 %--------------------------------------------------------------------------
 nl    = size(M,2);                        % number of levels
 nv    = sum(spm_vec(M.l));                % number of v (causal states)
 na    = sum(spm_vec(M.k));                % number of a (active states)
 nx    = sum(spm_vec(M.n));                % number of x (hidden states)
 
 % order parameters (n = 1 for static models)
 %==========================================================================
 n     = M(1).E.n   1;                     % order of embedding
 
 % initialise arrays for hierarchical form
 %--------------------------------------------------------------------------
 dfdvi = cell(nl,nl);
 dfdxi = cell(nl,nl);
 dfdai = cell(nl,nl);
 dgdvi = cell(nl,nl);
 dgdxi = cell(nl,nl);
 dgdai = cell(nl,nl);
 for i = 1:nl
     dgdvi{i,i} = sparse(M(i).l,M(i).l);
     dgdxi{i,i} = sparse(M(i).l,M(i).n);
     dgdai{i,i} = sparse(M(i).l,M(i).k);
     dfdvi{i,i} = sparse(M(i).n,M(i).l);
     dfdxi{i,i} = sparse(M(i).n,M(i).n);
     dfdai{i,i} = sparse(M(i).n,M(i).k);
 end
 
 % partition states {x,v,z,w} into distinct vector arrays v{i}, ...
 %--------------------------------------------------------------------------
 vi    = spm_unvec(u.v{1},{M.v});
 xi    = spm_unvec(u.x{1},{M.x});
 ai    = spm_unvec(u.a{1},{M.a});
 zi    = spm_unvec(u.z{1},{M.v});
 wi    = spm_unvec(u.w{1},{M.x});
 
 % Derivatives for Jacobian
 %==========================================================================
 vi{nl} = zi{nl};
 for  i = (nl - 1):-1:1
 
     % evaluate
     %----------------------------------------------------------------------
     [dgdx g] = spm_diff(M(i).g,xi{i},vi{i   1},ai{i   1},M(i).pE,1);
     [dfdx f] = spm_diff(M(i).f,xi{i},vi{i   1},ai{i   1},M(i).pE,1);
     dgdv     = spm_diff(M(i).g,xi{i},vi{i   1},ai{i   1},M(i).pE,2);
     dfdv     = spm_diff(M(i).f,xi{i},vi{i   1},ai{i   1},M(i).pE,2);
     dgda     = spm_diff(M(i).g,xi{i},vi{i   1},ai{i   1},M(i).pE,3);
     dfda     = spm_diff(M(i).f,xi{i},vi{i   1},ai{i   1},M(i).pE,3);
 
     % g(x,v) & f(x,v)
     %----------------------------------------------------------------------
     gi{i}    = g;
     fi{i}    = f;
     vi{i}    = gi{i}   zi{i};
 
     % and partial derivatives
     %----------------------------------------------------------------------
     dgdxi{i,    i} = dgdx;
     dgdvi{i,i   1} = dgdv;
     dgdai{i,i   1} = dgda;
     dfdxi{i,    i} = dfdx;
     dfdvi{i,i   1} = dfdv;
     dfdai{i,i   1} = dfda;
 
 end
 
 % concatenate hierarchical arrays
 %--------------------------------------------------------------------------
 dg.da = spm_cat(dgdai);
 dg.dv = spm_cat(dgdvi);
 dg.dx = spm_cat(dgdxi);
 df.da = spm_cat(dfdai);
 df.dv = spm_cat(dfdvi);
 df.dx = spm_cat(dfdxi);
 
 % update generalised coordinates
 %--------------------------------------------------------------------------
 u.v{1}  = spm_vec(vi);
 u.x{2}  = spm_vec(fi)   u.w{1};
 for i = 2:(n - 1)
     u.v{i}     = dg.dv*u.v{i}   dg.dx*u.x{i}   dg.da*u.a{i}   u.z{i};
     u.x{i   1} = df.dv*u.v{i}   df.dx*u.x{i}   df.da*u.a{i}   u.w{i};
 end
 
diff --git a/spm_ADEM_set.m b/spm_ADEM_set.m
@@ -0,0  1,88 @@
 function [DEM] = spm_ADEM_set(DEM)
 % Performs checks on DEM structures for active inversion
 % FORMAT [DEM] = spm_DEM_set(DEM)
 %
 % DEM.G  - generative model
 % DEM.M  - recognition model
 % DEM.C  - causes
 % DEM.U  - explanatory variables, inputs or prior expectation of causes
 %__________________________________________________________________________
 % Copyright (C) 2005 Wellcome Department of Imaging Neuroscience
 
 % Karl Friston
 % $Id: spm_ADEM_set.m 1887 2008-07-04 17:48:42Z karl $
 
 % check recognition model
 % -------------------------------------------------------------------------
 try
     DEM.M = spm_DEM_M_set(DEM.M);
 catch
     errordlg('please check your inversion model')
 end
 try
     DEM.G = spm_ADEM_M_set(DEM.G);
 catch
     errordlg('please check your generative model')
 end
 
 % check data or generative model
 % -------------------------------------------------------------------------
 try
     N     = length(DEM.C);
 catch
     errordlg('please specify causes')
 end
 try
     DEM.class;
 catch
     DEM.class = 'active';
 end
 
 % ensure model and data dimensions check
 % -------------------------------------------------------------------------
 try
     if size(DEM.Y,1) ~= DEM.M(1).l
         errordlg('DCM and data are incompatible')
     end
 catch
     if size(DEM.C,1) ~= DEM.M(end).l
         errordlg('DCM and causes are incompatible')
     end
 end
 
 % Default priors and confounds
 % -------------------------------------------------------------------------
 n  = DEM.M(end).l;
 if ~isfield(DEM,'U')
     DEM.U = sparse(n,N);
 end
 
 % transpose causes and confounds, if specified in conventional fashion
 %--------------------------------------------------------------------------
 if size(DEM.U,2) < N, DEM.U = DEM.U';    end
 if size(DEM.C,2) < N, DEM.C = DEM.C';    end
 
 % check prior expectation of causes (at level n) and confounds
 %--------------------------------------------------------------------------
 if ~nnz(DEM.U), DEM.U = sparse(n,N); end
 if ~nnz(DEM.C), DEM.C = sparse(n,N); end
 
 % ensure inputs and cause dimensions check
 % -------------------------------------------------------------------------
 if size(DEM.U,1) ~= DEM.M(end).l
     errordlg('DCM inputs and priors are not compatible')
 end
 
 % ensure causes and data dimensions check
 % -------------------------------------------------------------------------
 if size(DEM.U,2) < N
     errordlg('priors and data have different lengths')
 end
 
 % check length of time-series
 %--------------------------------------------------------------------------
 if N < DEM.M(1).E.n
     errordlg('Please ensure time-series is longer than embedding order')
     return
 end
 
diff --git a/spm_DEM_qU.m b/spm_DEM_qU.m
@@ -13,7  13,7 @@ function spm_DEM_qU(qU,pU)
 % Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging
 
 % Karl Friston
-% $Id: spm_DEM_qU.m 1703 2008-05-21 13:59:23Z karl $
 % $Id: spm_DEM_qU.m 1887 2008-07-04 17:48:42Z karl $
 
 % unpack
 %--------------------------------------------------------------------------
@@ -175,4  175,23 @@ function spm_DEM_qU(qU,pU)
         end
     end
 end
 
 % plot action if specified
 %--------------------------------------------------------------------------
 try
     subplot(g,2,2*g)
     plot(qU.a{2});
     try
         hold on
         plot(pU.v{2},':'); hold off
     end
     xlabel('time','Fontsize',14)
     title('perturbation and action','Fontsize',16)
     grid on
     axis square
     set(gca,'XLim',[t(1) t(end)])
 catch
     delete(gca)
 end
 
 drawnow