spm_LAP.m

function [DEM] = spm_LAP(DEM)
% Laplacian model inversion (see also spm_LAPS)
% FORMAT DEM   = spm_LAP(DEM)
%
% DEM.M  - hierarchical model
% DEM.Y  - response variable, output or data
% DEM.U  - explanatory variables, inputs or prior expectation of causes
%__________________________________________________________________________
%
% generative model
%--------------------------------------------------------------------------
%   M(i).g  = v     =  g(x,v,P)   {inline function, string or m-file}
%   M(i).f  = dx/dt =  f(x,v,P)   {inline function, string or m-file}
%
%   M(i).ph = pi(v) = ph(x,v,h,M) {inline function, string or m-file}
%   M(i).pg = pi(x) = pg(x,v,g,M) {inline function, string or m-file}
%                                  (assumed to be linear in v and x)
%
%   pi(v,x) = vectors of log-precisions; (h,g) = precision parameters
%
%   M(i).pE = prior expectation of p model-parameters
%   M(i).pC = prior covariances of p model-parameters
%   M(i).hE = prior expectation of h log-precision (cause noise)
%   M(i).hC = prior covariances of h log-precision (cause noise)
%   M(i).gE = prior expectation of g log-precision (state noise)
%   M(i).gC = prior covariances of g log-precision (state noise)
%   M(i).xP = precision (states)
%   M(i).Q  = precision components (input noise)
%   M(i).R  = precision components (state noise)
%   M(i).V  = fixed precision (input noise)
%   M(i).W  = fixed precision (state noise)
%
%   M(i).P  = optional initial value for parameters (defaults to M(i).pE)
%
%   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);
%
% conditional moments of model-states - q(u)
%--------------------------------------------------------------------------
%   qU.x    = Conditional expectation of hidden states
%   qU.v    = Conditional expectation of causal states
%   qU.w    = Conditional prediction error (states)
%   qU.z    = Conditional prediction error (causes)
%   qU.C    = Conditional covariance: cov(v)
%   qU.S    = Conditional covariance: cov(x)
%
% conditional moments of model-parameters - q(p)
%--------------------------------------------------------------------------
%   qP.P    = Conditional expectation
%   qP.C    = Conditional covariance
%
% conditional moments of hyper-parameters (log-transformed) - q(h)
%--------------------------------------------------------------------------
%   qH.h    = Conditional expectation (cause noise)
%   qH.g    = Conditional expectation (state noise)
%   qH.C    = Conditional covariance
%
% F         = log-evidence = log-marginal likelihood = negative free-energy
%
%__________________________________________________________________________
% Accelerated methods: To accelerate computations one can specify the 
% nature of the model equations using:
%
% M(1).E.linear = 0: full        - evaluates 1st and 2nd derivatives
% M(1).E.linear = 1: linear      - equations are linear in x and v
% M(1).E.linear = 2: bilinear    - equations are linear in x, v & x*v
% M(1).E.linear = 3: nonlinear   - equations are linear in x, v, x*v, & x*x
% M(1).E.linear = 4: full linear - evaluates 1st derivatives (for GF)
%
% similarly, for evaluating precisions:
%
% M(1).E.method.h = 0,1  switch for precision parameters (hidden causes)
% M(1).E.method.g = 0,1  switch for precision parameters (hidden states)
% M(1).E.method.v = 0,1  switch for precision (hidden causes)
% M(1).E.method.x = 0,1  switch for precision (hidden states)
%__________________________________________________________________________
%
% spm_LAP implements a variational scheme under the Laplace
% approximation to the conditional joint density q on states u, parameters 
% p and hyperparameters (h,g) of an analytic nonlinear hierarchical dynamic
% model, with additive Gaussian innovations.
%
%            q(u,p,h,g) = max <L(t)>q
%
% L is the ln p(y,u,p,h,g|M) under the model M. The conditional covariances
% obtain analytically from the curvature of L with respect to the unknowns.
%__________________________________________________________________________
% Copyright (C) 2010-2013 Wellcome Trust Centre for Neuroimaging
 
% Karl Friston
% $Id: spm_LAP.m 6508 2015-07-25 15:23:25Z karl $
 
 
% find or create a DEM figure
%--------------------------------------------------------------------------
try
    DEM.M(1).nograph;
catch
    DEM.M(1).nograph = 0;
end
if ~DEM.M(1).nograph
    Fdem = spm_figure('GetWin','DEM');
end
 
 
% check model, data and priors
%==========================================================================
[M,Y,U] = spm_DEM_set(DEM);


% set regularisation
%--------------------------------------------------------------------------
try
    dt = DEM.M(1).E.v;
catch
    dt = 0;
    DEM.M(1).E.v = dt;
end

 
% number of iterations
%--------------------------------------------------------------------------
try, nD = M(1).E.nD; catch, nD = 1;  end
try, nN = M(1).E.nN; catch, nN = 8;  end
 
 
% ensure integration scheme evaluates gradients at each time-step
%--------------------------------------------------------------------------
M(1).E.linear = 4;
 
% assume precisions are a function of, and only of, hyperparameters
%--------------------------------------------------------------------------
try
    method = M(1).E.method;
catch
    method.h = 1;
    method.g = 1;
    method.x = 0;
    method.v = 0;
end
try method.h; catch, method.h = 0; end
try method.g; catch, method.g = 0; end
try method.x; catch, method.x = 0; end
try method.v; catch, method.v = 0; end
 

% assume precisions are a function of, and only of, hyperparameters
%--------------------------------------------------------------------------
try
    form = M(1).E.form;
catch
    form = 'Gaussian';
end
 
% checks for Laplace models (precision functions; ph and pg)
%--------------------------------------------------------------------------
for i  = 1:length(M)
    try
        feval(M(i).ph,M(i).x,M(i + 1).v,M(i).hE,M(i)); method.v = 1;
    catch
        M(i).ph = inline('spm_LAP_ph(x,v,h,M)','x','v','h','M');
    end
    try
        feval(M(i).pg,M(i).x,M(i + 1).v,M(i).gE,M(i)); method.x = 1;
    catch
        M(i).pg = inline('spm_LAP_pg(x,v,h,M)','x','v','h','M');
    end
end

M(1).E.method = method;
 
 
% order parameters (d = n = 1 for static models) and checks
%==========================================================================
d   = M(1).E.d + 1;                          % embedding order of q(v)
n   = M(1).E.n + 1;                          % embedding order of q(x)
 
% number of states and parameters
%--------------------------------------------------------------------------
ns  = size(Y,2);                             % number of samples
nl  = size(M,2);                             % number of levels
nv  = sum(spm_vec(M.m));                     % number of v (casual states)
nx  = sum(spm_vec(M.n));                     % number of x (hidden states)
ny  = M(1).l;                                % number of y (inputs)
nc  = M(end).l;                              % number of c (prior causes)
nu  = nv*d + nx*n;                           % number of generalised states
ne  = nv*n + nx*n + ny*n;                    % number of generalised errors
 
 
% precision (R) of generalised errors and null matrices for concatenation
%==========================================================================
s     = M(1).E.s;
Rh    = spm_DEM_R(n,s,form);
Rg    = spm_DEM_R(n,s,form);

if ~nx, Rg = sparse(0,0); end
 
W     = sparse(nx*n,nx*n);
V     = sparse((ny + nv)*n,(ny + nv)*n);
 
 
% fixed priors on states (u)
%--------------------------------------------------------------------------
Px    = kron(spm_DEM_R(n,2),spm_cat(spm_diag({M(1:end).xP})));
Pv    = kron(spm_DEM_R(d,2),spm_cat(spm_diag({M(2:end).vP})));
Pu    = spm_cat(spm_diag({Px Pv}));
 
% hyperpriors
%--------------------------------------------------------------------------
ph.h  = spm_vec({M.hE M.gE});                % prior expectation of h,g
ph.c  = spm_cat(spm_diag({M.hC M.gC}));      % prior covariances of h,g
Ph    = spm_inv(ph.c);                       % prior precision of h,g
 
qh.h  = {M.hE};                              % conditional expectation h
qh.g  = {M.gE};                              % conditional expectation g
nh    = length(spm_vec(qh.h));               % number of hyperparameters h
ng    = length(spm_vec(qh.g));               % number of hyperparameters g
nb    = nh + ng;                             % number of hyperparameters
 
 
% priors on parameters (in reduced parameter space)
%==========================================================================
pp.c  = cell(nl,nl);
qp.p  = cell(nl,1);
for i = 1:(nl - 1)
 
    % eigenvector reduction: p <- pE + qp.u*qp.p
    %----------------------------------------------------------------------
    qp.u{i}   = spm_svd(M(i).pC,0);          % basis for parameters
    M(i).p    = size(qp.u{i},2);             % number of qp.p
    qp.p{i}   = sparse(M(i).p,1);            % initial deviates
    pp.c{i,i} = qp.u{i}'*M(i).pC*qp.u{i};    % prior covariance
 
end
Up    = spm_cat(spm_diag(qp.u));
 
% priors on parameters
%--------------------------------------------------------------------------
pp.p  = spm_vec(M.pE);
pp.c  = spm_cat(pp.c);
Pp    = spm_inv(pp.c);
 

% initialise conditional density q(p)
%--------------------------------------------------------------------------
for i = 1:(nl - 1)
    try
        qp.p{i} = qp.p{i} + qp.u{i}'*(spm_vec(M(i).P) - spm_vec(M(i).pE));
    end
end
np    = size(Up,2);
 
 
% initialise cell arrays for D-Step; e{i + 1} = (d/dt)^i[e] = e[i]
%==========================================================================
qu.x      = cell(n,1);
qu.v      = cell(n,1);
qu.y      = cell(n,1);
qu.u      = cell(n,1);
[qu.x{:}] = deal(sparse(nx,1));
[qu.v{:}] = deal(sparse(nv,1));
[qu.y{:}] = deal(sparse(ny,1));
[qu.u{:}] = deal(sparse(nc,1));
 
% initialise cell arrays for hierarchical structure of x[0] and v[0]
%--------------------------------------------------------------------------
x         = {M(1:end - 1).x};
v         = {M(1 + 1:end).v};
qu.x{1}   = spm_vec(x);
qu.v{1}   = spm_vec(v);
 
% derivatives for Jacobian of D-step
%--------------------------------------------------------------------------
Dx     = kron(spm_speye(n,n,1),spm_speye(nx,nx));
Dv     = kron(spm_speye(d,d,1),spm_speye(nv,nv));
Dy     = kron(spm_speye(n,n,1),spm_speye(ny,ny));
Dc     = kron(spm_speye(d,d,1),spm_speye(nc,nc));
Du     = spm_cat(spm_diag({Dx,Dv}));
Ib     = spm_speye(np + nb,np + nb);
dbdt   = sparse(np + nb,1);
 
 
% gradients of generalised weighted errors
%--------------------------------------------------------------------------
dedh   = sparse(nh,ne);
dedg   = sparse(ng,ne);
dedv   = sparse(nv,ne);
dedx   = sparse(nx,ne);
dedhh  = sparse(nh,nh);
dedgg  = sparse(ng,ng);
            
% curvatures of Gibb's energy w.r.t. hyperparameters
%--------------------------------------------------------------------------
dHdh   = sparse(nh,1);
dHdg   = sparse(ng,1);
dHdp   = sparse(np,1);
dHdu   = sparse(nu,1);

 
% preclude unnecessary iterations and set switches
%--------------------------------------------------------------------------
if ~np && ~nh && ~ng, nN = 1; end
mnx    = nx*~~method.x;
mnv    = nv*~~method.v;
 
% preclude very precise states from entering free-energy/action
%--------------------------------------------------------------------------
p      = spm_LAP_eval(M,qu,qh);
ih     = p.h < 8;
ig     = p.g < 8;
ie     = kron(ones(n,1),ih);
ix     = kron(ones(n,1),ig);
iv     = kron(ones(d,1),ih((1:nv) + ny));
je     = find([ie; ix]); ix(1:nx) = 1;
ju     = find([ix; iv]);
 
% and other useful indices
%--------------------------------------------------------------------------
ix     = (1:nx);
ih     = (1:nb);
iv     = (1:nv) + nx*n;
ip     = (1:np) + nu;


% number of iterations for convergence
%--------------------------------------------------------------------------
convergence = -4;

% Iterate Laplace scheme
%==========================================================================
F      = -Inf;
for iN = 1:nN
 
    % get time and clear persistent variables in evaluation routines
    %----------------------------------------------------------------------
    tic; clear spm_DEM_eval
 
    % [re-]set states & their derivatives
    %----------------------------------------------------------------------
    try, qu = Q(1).u; end
    
    
    % D-Step: (nD D-Steps for each sample)
    %======================================================================
    for is = 1:ns
 
        % D-Step: until convergence for static systems
        %==================================================================
        for iD = 1:nD
 
            % sampling time
            %--------------------------------------------------------------
            ts = is + (iD - 1)/nD;
 
            % derivatives of responses and inputs
            %--------------------------------------------------------------
            try
                qu.y(1:n) = spm_DEM_embed(Y,n,ts,1,M(1).delays);
                qu.u(1:d) = spm_DEM_embed(U,d,ts);
            catch
                qu.y(1:n) = spm_DEM_embed(Y,n,ts);
                qu.u(1:d) = spm_DEM_embed(U,d,ts);
            end
            
            
            % evaluate functions and derivatives
            %==============================================================
            
            % prediction errors (E) and precision vectors (p)
            %--------------------------------------------------------------
            [E,dE] = spm_DEM_eval(M,qu,qp);
            [p,dp] = spm_LAP_eval(M,qu,qh);
            
 
            % gradients of log(det(iS)) dDd...
            %==============================================================
            
            % get precision matrices
            %--------------------------------------------------------------
            iSh   = diag(exp(p.h));
            iSg   = diag(exp(p.g));
            iS    = blkdiag(kron(Rh,iSh),kron(Rg,iSg));
 
            
            % gradients of trace(diag(p)) = sum(p); p = precision vector
            %--------------------------------------------------------------
            dpdx  = n*sum(spm_cat({dp.h.dx; dp.g.dx}));
            dpdv  = n*sum(spm_cat({dp.h.dv; dp.g.dv}));
            dpdh  = n*sum(dp.h.dh);
            dpdg  = n*sum(dp.g.dg);
            dpdx  = kron(sparse(1,1,1,1,n),dpdx);
            dpdv  = kron(sparse(1,1,1,1,d),dpdv);
            dDdu  = [dpdx dpdv]';
            dDdh  = [dpdh dpdg]';
 
            
            % gradients precision-weighted generalised error dSd..
            %==============================================================
 
            % gradients w.r.t. hyperparameters
            %--------------------------------------------------------------
            for i = 1:nh
                diS       = diag(dp.h.dh(:,i).*exp(p.h));
                diSdh{i}  = blkdiag(kron(Rh,diS),W);
                dedh(i,:) = E'*diSdh{i};
            end
            for i = 1:ng
                diS       = diag(dp.g.dg(:,i).*exp(p.g));
                diSdg{i}  = blkdiag(V,kron(Rg,diS));
                dedg(i,:) = E'*diSdg{i};
            end
 
            % gradients w.r.t. hidden states
            %--------------------------------------------------------------
            for i = 1:mnx
                diV       = diag(dp.h.dx(:,i).*exp(p.h));
                diW       = diag(dp.g.dx(:,i).*exp(p.g));
                diSdx{i}  = blkdiag(kron(Rh,diV),kron(Rg,diW));
                dedx(i,:) = E'*diSdx{i};
            end
            
            % gradients w.r.t. causal states
            %--------------------------------------------------------------
            for i = 1:mnv
                diV       = diag(dp.h.dv(:,i).*exp(p.h));
                diW       = diag(dp.g.dv(:,i).*exp(p.g));
                diSdv{i}  = blkdiag(kron(Rh,diV),kron(Rg,diW));
                dedv(i,:) = E'*diSdv{i};
            end
 
            dSdx  = kron(sparse(1,1,1,n,1),dedx);
            dSdv  = kron(sparse(1,1,1,d,1),dedv);
            dSdu  = [dSdx; dSdv];
            dEdh  = [dedh; dedg];
            dEdp  = dE.dp'*iS;
            dEdu  = dE.du'*iS;
            
            % curvatures w.r.t. hyperparameters
            %--------------------------------------------------------------
            for i = 1:nh
                for j = i:nh
                    diS        = diag(dp.h.dh(:,i).*dp.h.dh(:,j).*exp(p.h));
                    diS        = blkdiag(kron(Rh,diS),W);
                    dedhh(i,j) = E'*diS*E;
                    dedhh(j,i) = dedhh(i,j);
                end
            end            
            for i = 1:ng
                for j = i:ng
                    diS        = diag(dp.g.dg(:,i).*dp.g.dg(:,j).*exp(p.g));
                    diS        = blkdiag(V,kron(Rg,diS));
                    dedgg(i,j) = E'*diS*E;
                    dedgg(j,i) = dedgg(i,j);
                end
            end
            
            % combined curvature
            %--------------------------------------------------------------
            dSdhh = spm_cat({dedhh  [] ;
                             [] dedgg});
                 
            
            % errors (from prior expectations) (NB pp.p = 0)
            %--------------------------------------------------------------
            Eu    = spm_vec(qu.x(1:n),qu.v(1:d));
            Ep    = spm_vec(qp.p);
            Eh    = spm_vec(qh.h,qh.g) - ph.h;
            
 
            % first-order derivatives of Gibb's Energy
            %==============================================================
            dLdu  = dEdu*E + dSdu*E/2 - dDdu/2 + Pu*Eu;
            dLdh  = dEdh*E/2          - dDdh/2 + Ph*Eh;
            dLdp  = dEdp*E                     + Pp*Ep;
            

            % and second-order derivatives of Gibb's Energy
            %--------------------------------------------------------------
            dLduu = dEdu*dE.du + Pu;
            dLdpp = dEdp*dE.dp + Pp;
            dLdhh = dSdhh/2    + Ph;
            dLduy = dEdu*dE.dy;
            dLduc = dEdu*dE.dc;
            dLdup = dEdu*dE.dp;
            dLdhp = dEdh*dE.dp;
            dLdpu = dLdup';
            dLdph = dLdhp';
 
 
            % precision and covariances for entropy
            %--------------------------------------------------------------                      
            dLdaa = spm_cat({dLduu dLdup  ;
                             dLdpu dLdpp});
            dLdbb = spm_cat({dLdpp dLdph  ;
                             dLdhp dLdhh});
            
            Cup   = spm_inv(dLdaa);
            Chh   = spm_inv(dLdhh);
            
            
            % first-order derivatives of Entropy term
            %==============================================================
            
            % log-precision
            %--------------------------------------------------------------
            for i = 1:nh
                Luub    = dE.du'*diSdh{i}*dE.du;
                Lpub    = dE.dp'*diSdh{i}*dE.du;
                Lppb    = dE.dp'*diSdh{i}*dE.dp;
                diCdh   = spm_cat({Luub Lpub';
                                   Lpub Lppb});
                dHdh(i) = spm_trace(diCdh,Cup)/2;
            end
            for i = 1:ng
                Luub    = dE.du'*diSdg{i}*dE.du;
                Lpub    = dE.dp'*diSdg{i}*dE.du;
                Lppb    = dE.dp'*diSdg{i}*dE.dp;
                diCdg   = spm_cat({Luub Lpub';
                                   Lpub Lppb});
                dHdg(i) = spm_trace(diCdg,Cup)/2;
            end
            
            % parameters
            %--------------------------------------------------------------
            for i = 1:np
                Luup    = dE.dup{i}'*dEdu';
                Lpup    = dEdp*dE.dup{i};
                Luup    = Luup + Luup';
                diCdp   = spm_cat({Luup Lpup';
                                   Lpup [] });
                dHdp(i) = spm_trace(diCdp,Cup)/2;
            end
            
            % hidden states and causes (disabled for stability)
            %--------------------------------------------------------------
            for i = 1:(nu - nu)
                Lppu    = dE.dpu{i}'*dEdp';
                Lupu    = dEdu*dE.dpu{i};
                Lppu    = Lppu + Lppu';
                diCdu   = spm_cat({[]    Lupu;
                                   Lupu' Lppu});
                dHdu(i) = spm_trace(diCdu,Cup)/2;
            end
 
            % and concatenate
            %--------------------------------------------------------------
            dHdb  = [dHdh; dHdg];
            dHdb  = [dHdp; dHdb];
            dLdb  = [dLdp; dLdh];

                        
            % save conditional moments (and prediction error) at Q{t}
            %==============================================================
            if iD == 1 || ns == 1
                
                % save means
                %----------------------------------------------------------
                Q(is).e   = E;
                Q(is).E   = diag(diag(iS))*E;
                Q(is).u   = qu;
                Q(is).p   = qp;
                Q(is).h   = qh;
                
                % and conditional covariances
                %----------------------------------------------------------
                Q(is).u.s = Cup(ix,ix);
                Q(is).u.c = Cup(iv,iv);
                Q(is).p.c = Cup(ip,ip);
                Q(is).h.c = Chh(ih,ih);              
                
                % Free-energy (components)
                %----------------------------------------------------------
                Fc(is,1)  = - E(je)'*iS(je,je)*E(je)/2;
                Fc(is,2)  = - Eu(ju)'*Pu(ju,ju)*Eu(ju)/2;
                Fc(is,3)  = - n*ny*log(2*pi)/2;
                Fc(is,4)  = spm_logdet(iS(je,je))/2;
                Fc(is,5)  = spm_logdet(Pu(ju,ju)*Cup(ju,ju))/2;

                                
                % Free-action (states and parameters)
                %----------------------------------------------------------
                AC(is)    = sum(Fc(is,:))    ...
                          - Ep'*Pp*Ep/2      ...
                          - Eh'*Ph*Eh/2      ...
                          + spm_logdet(Pp)/2 ...
                          + spm_logdet(Ph)/2 ...
                          - spm_logdet(dLdbb)/2;
  
            end
 
            % update conditional moments
            %==============================================================
            
            % prior precision of fluctuations on [hyper] parameters
            %--------------------------------------------------------------
            Kb    = ns*Ib;
            
            % accumulate curvatures of [hyper] parameters
            %--------------------------------------------------------------
            try
                dLdBB = dLdBB*(1 - 1/ns) + dLdbb/ns;
            catch
                dLdBB = dLdbb + Ib*32;
            end

            % whiten gradient (and curvatures) with regularised precision
            %--------------------------------------------------------------
            Cb    = spm_inv(dLdBB + diag(diag(dLdBB))*exp(dt));
            dLdb  = Cb*dLdb;
            dHdb  = Cb*dHdb;
            
            % assemble conditional means
            %--------------------------------------------------------------
            q.y  = qu.y(1:n);
            q.x  = qu.x(1:n);
            q.v  = qu.v(1:d);
            q.c  = qu.u(1:d);
            q.p  = qp.p;
            q.h  = qh.h;
            q.g  = qh.g;
            q.d  = dbdt;
            
                        
            % flow
            %--------------------------------------------------------------
            f.y  =  Dy*spm_vec(q.y)                  ;
            f.u  =  Du*spm_vec(q.x,q.v) - dLdu - dHdu;
            f.c  =  Dc*spm_vec(q.c)                  ;
            f.b  =     spm_vec(q.d)                  ;
            f.d  = -Kb*spm_vec(q.d)     - dLdb - dHdb;

            % and Jacobian
            %--------------------------------------------------------------
            dfdq = {Dy     []        []     []    [] ;
                   -dLduy  Du-dLduu -dLduc  []    [] ;
                    []     []        Dc     []    [] ;
                    []     []        []     []    Ib ;
                    []     []        []    -Ib   -Kb};
                
            dfdq = spm_cat(dfdq);
                
            % get eigenvalues of Jacobian (on hidden states)
            %--------------------------------------------------------------
            try DEM.options.eigenvalues;
                DEM.E(:,is) = eig(full(Du - dLduu));
            end
            
            
            % update conditional modes of states
            %==============================================================
            dq   = spm_dx(dfdq,spm_vec(f),1/nD);
            q    = spm_unvec(spm_vec(q) + dq,q);
            
            % unpack conditional means
            %--------------------------------------------------------------
            qu.x(1:n) = q.x;
            qu.v(1:d) = q.v;
            qp.p      = q.p;
            qh.h      = q.h;
            qh.g      = q.g;
            dbdt      = q.d;
 
        end % D-Step
 
    end % sequence (ns)
 
    
    % Bayesian parameter averaging
    %======================================================================
 
    % Conditional moments of time-averaged parameters
    %----------------------------------------------------------------------
    Ep    = 0;
    Qp    = 0;
    for i = 1:ns
        P   = spm_inv(Q(i).p.c);
        Ep  = Ep  + P*spm_vec(Q(i).p.p);
        Qp  = Qp + P;       
    end
    Ep    = spm_inv(Qp)*Ep;
    Cp    = spm_inv(Qp + (1 - ns)*Pp);
 
    % conditional moments of hyper-parameters
    %----------------------------------------------------------------------
    Eh    = 0;
    Qh    = 0;
    for i = 1:ns
        P   = spm_inv(Q(i).h.c);
        Eh  = Eh  + P*spm_vec(Q(i).h.h,Q(i).h.g);
        Qh  = Qh + P;
    end
    Eh    = spm_inv(Qh)*Eh - ph.h;
    Ch    = spm_inv(Qh + (1 - ns)*Ph);
 
    
    
    % Free-action of states plus free-energy of parameters
    %======================================================================
    FT    = sum(Fc,2);          %   instantaneous Free energy (of states)
    FC(1) = sum(Fc(:,1));       % - E'*iS*E/2;
    FC(2) = sum(Fc(:,2));       % - Eu'*Pu*Eu/2;
    FC(3) = sum(Fc(:,3));       % - n*ny*log(2*pi)/2;
    FC(4) = sum(Fc(:,4));       %   spm_logdet(iS)/2;
    FC(5) = sum(Fc(:,5));       %   spm_logdet(Pu*Cu)/2;
    FC(6) = -Ep'*Pp*Ep/2;
    FC(7) = -Eh'*Ph*Eh/2;
    FC(8) = spm_logdet(Pp*Cp)/2;
    FC(9) = spm_logdet(Ph*Ch)/2;
    
    CC(iN,:) = FC;
    S(iN)    = sum(AC);
    Fe       = sum(FC);
 
    % if F is decreasing, revert [hyper] parameters and slow down
    %----------------------------------------------------------------------
    if Fe < F(iN) && iN > 4
                
        % save free-energy
        %------------------------------------------------------------------
        F(iN + 1) = F(iN);
               
        % load current MAP estimates
        %------------------------------------------------------------------
        qp = PQ.qp;
        qh = PQ.qh;
        
        % decrease update time
        %------------------------------------------------------------------
        dt = max(dt + 2,2);
        
        % convergence
        %------------------------------------------------------------------
        if dt > 6; convergence = 1; end
        
    else
        
        % convergence
        %------------------------------------------------------------------
        if Fe - F(iN) < 1, convergence = convergence + 1; end
 
        % save free-energy
        %------------------------------------------------------------------
        F(iN)     = Fe;
        F(iN + 1) = Fe;
        
        % save current MAP estimates
        %------------------------------------------------------------------
        PQ.qp = qp;
        PQ.qh = qh;
        
        % increase update time
        %------------------------------------------------------------------
        dt    = max(dt - 1,-8);
        
    end
 
    % Convergence
    %======================================================================
    if convergence > 0; break, end
    
    % otherwise save conditional moments (for each time point)
    %======================================================================
    for t = 1:length(Q)
 
        % states and predictions
        %------------------------------------------------------------------
        v     = spm_unvec(Q(t).u.v{1},v);
        x     = spm_unvec(Q(t).u.x{1},x);
        z     = spm_unvec(Q(t).e(1:(ny + nv)),{M.v});
        Z     = spm_unvec(Q(t).E(1:(ny + nv)),{M.v});
        w     = spm_unvec(Q(t).e((1:nx) + (ny + nv)*n),{M.x});
        X     = spm_unvec(Q(t).E((1:nx) + (ny + nv)*n),{M.x});
        for i = 1:(nl - 1)
            if M(i).m, qU.v{i + 1}(:,t) = spm_vec(v{i});  end
            if M(i).n, qU.x{i}(:,t)     = spm_vec(x{i});  end
            if M(i).n, qU.w{i}(:,t)     = spm_vec(w{i});  end
            if M(i).l, qU.z{i}(:,t)     = spm_vec(z{i});  end
            if M(i).n, qU.W{i}(:,t)     = spm_vec(X{i});  end
            if M(i).l, qU.Z{i}(:,t)     = spm_vec(Z{i});  end
        end
        if    M(nl).l, qU.z{nl}(:,t)    = spm_vec(z{nl}); end
        if    M(nl).l, qU.Z{nl}(:,t)    = spm_vec(Z{nl}); end
 
        qU.v{1}(:,t)  = spm_vec(Q(t).u.y{1}) - spm_vec(z{1});
 
        % and conditional covariances
        %------------------------------------------------------------------
        qU.S{t} = Q(t).u.s;
        qU.C{t} = Q(t).u.c;
 
        % parameters
        %------------------------------------------------------------------
        qP.p{t} = spm_vec(Q(t).p.p);
        qP.c{t} = Q(t).p.c;
 
        % hyperparameters
        %------------------------------------------------------------------
        qH.p{t} = spm_vec(Q(t).h.h,Q(t).h.g);
        qH.c{t} = Q(t).h.c;
 
    end
 
    % graphics (states)
    %----------------------------------------------------------------------
    if exist('Fdem','var')
        spm_figure('Select', Fdem)
        spm_DEM_qU(qU)
        
        % graphics (parameters and log-precisions)
        %------------------------------------------------------------------
        if np && nb
            subplot(2*nl,2,4*nl - 2)
            plot(1:ns,full(spm_cat(qP.p)))
            set(gca,'XLim',[1 ns])
            title('parameters (modes)','FontSize',16)
            
            subplot(2*nl,2,4*nl)
            plot(1:ns,full(spm_cat(qH.p)))
            set(gca,'XLim',[1 ns])
            title('log-precision','FontSize',16)
            
        elseif nb
            subplot(nl,2,2*nl)
            plot(1:ns,full(spm_cat(qH.p)))
            set(gca,'XLim',[1 ns])
            title('log-precision','FontSize',16)
            
        elseif np
            subplot(nl,2,2*nl)
            plot(1:ns,full(spm_cat(qP.p)))
            set(gca,'XLim',[1 ns])
            title('parameters (modes)','FontSize',16)
            
        end
        drawnow
    end
    
    % report (EM-Steps)
    %----------------------------------------------------------------------
    try
        dF = F(iN) - F(iN - 1);
    catch
        dF = 0;
    end
    str{1} = sprintf('LAP: %i (%i)', iN,iD);
    if iN == 1
        str{2} = sprintf('  F0:%.4e', full(F(iN)));
    else
        str{2} = sprintf('F-F0:%.4e', full(F(iN) - F(1)));
    end
    str{3} = sprintf('dF:%.2e',      full(dF));
    str{4} = sprintf('(%.2e sec)',   full(toc));
    fprintf('%-16s%-20s%-14s%-16s\n',str{:})
 
end
fprintf('%-19sF:%.4e\n', 'LAP: Converged', F(end));
if exist('Fdem','var')
    spm_figure('Focus', Fdem)
end

% Place Bayesian parameter averages in output arguments
%==========================================================================
 
% Conditional moments of time-averaged parameters
%--------------------------------------------------------------------------
Qp = 0;
Ep = 0;
for i = 1:ns
    
    % weight in proportion to precisions
    %----------------------------------------------------------------------
    P  = spm_inv(qP.c{i});
    Ep = Ep + P*qP.p{i};
    Qp = Qp + P;
 
end
Ep     = spm_inv(Qp)*Ep;
Cp     = spm_inv(Qp + (1 - ns)*Pp);
qP.P   = spm_unvec(Up*Ep + pp.p,{M.pE});
qP.C   = Up*Cp*Up';
qP.V   = spm_unvec(diag(qP.C),{M.pE});
qP.U   = Up;
 
% conditional moments of hyper-parameters
%--------------------------------------------------------------------------
Qh = 0;
Eh = 0;
for i = 1:ns
    
    % weight in proportion to precisions
    %----------------------------------------------------------------------
    P  = spm_inv(qH.c{i});
    Eh = Eh + P*qH.p{i};
    Qh = Qh + P;
 
end
Eh     = spm_inv(Qh)*Eh;
Ch     = spm_inv(Qh + (1 - ns)*Ph);
P      = spm_unvec(Eh,{qh.h,qh.g});
qH.h   = P{1};
qH.g   = P{2};
qH.C   = Ch;
P      = spm_unvec(diag(qH.C),P);
qH.V   = P{1};
qH.W   = P{2};
 
 
 
% assign output variables
%--------------------------------------------------------------------------
DEM.M  = M;                   % model
DEM.U  = U;                   % causes
 
DEM.qU = qU;                  % conditional moments of model-states
DEM.qP = qP;                  % conditional moments of model-parameters
DEM.qH = qH;                  % conditional moments of hyper-parameters
 
DEM.F  = F(1:iN);             % [-ve] Free-energy
DEM.S  = S(1:iN);             % [-ve] Free-action
DEM.FC = FC;                  % Free-energy components
DEM.CC = CC;                  % over iterations
DEM.FT = FT;                  % over time