function [MB] = DEMO_DCM_MB
% This demo uses the notion of Markov blankets and the renormalisation
% group to evaluate the coupling among neuronal systems at increasing
% spatial scales. The underlying generative model is based upon the
% renormalisation group: a working definition of renormalization involves
% three elements: vectors of random variables, a course-graining operation
% and a requirement that the operation does not change the functional form
% of the Lagrangian. In our case, the random variables are neuronal states;
% the course graining operation corresponds to the grouping (G) into a
% particular partition and adiabatic reduction (R) - that leaves the
% functional form of the dynamics unchanged.
%
% Here, the grouping operator (G) is based upon coupling among states as
% measured by the Jacobian. In brief, the sparsity structure of the
% Jacobian is used to recursively identify Markov blankets around internal
% states to create a partition of states - at any level - into particles;
% where each particle comprises external and blanket states. The ensuing
% reduction operator (R) eliminates the internal states and retains the slow
% eigenmodes of the blanket states. These then constitute the (vector)
% states at the next level and the process begins again.
%
% This routine starts using a simple form of dynamic causal modelling
% applied to the principal eigenvariate of local parcels (i.e., particles)
% of voxels with compact support. The Jacobian is estimated using a
% linearised dynamic causal (state space) model, where observations are
% generated by applying a (e.g., haemodynamic) convolution operator to
% hidden (e.g., neuronal) states. This estimation uses parametric empirical
% Bayes (PEB: spm_PEB). The ensuing estimates of the Jacobian (i.e.,
% effective connectivity) are reduced using Bayesian model reduction (BMR:
% spm_dcm_BMR_all) within a bespoke routine (spm_dcm_J).
%
% The Jacobian is then partitioned using the course graining operator into
% particles or parcels (using spm_markov blanket). The resulting partition
% is then reduced by eliminating internal states and retaining slow
% eigenmodes with the largest (real) eigenvalues (spm_A_reduce). The
% Jacobian of the reduced states is then used to repeat the process -
% recording the locations of recursively coarse-grained particles - until
% there is a single particle.
%
% The result of this recursive decomposition (i.e., renormalisation)
% affords a characterisation of directed coupling, as characterised by a
% complex Jacobian; namely, a multivariate coupling matrix, describing the
% coupling between eigenmodes of Markov blankets at successive scales. This
% can be regarded as a recursive parcellation scheme based upon effective
% connectivity and a generative (dynamic causal) model of multivariate
% (neuronal) timeseries.
% 
%__________________________________________________________________________

% Karl Friston
% Copyright (C) 2008-2022 Wellcome Centre for Human Neuroimaging


SPM = load('C:\home\spm\SPM\analysis_set9\SPM.mat','SPM');
MB  = spm_mb_ui('specify',SPM.SPM);
MB  = spm_mb_ui('blocking',MB(1));

% display the results in terms of particular partitions and eigenmodes
%--------------------------------------------------------------------------
spm_mb_ui('results',MB,'anatomy');

% characterise connectivity at the smallest scale
%--------------------------------------------------------------------------
spm_mb_ui('results',MB,'distance');

% characterise scaling behaviour in terms of scaling exponent
%--------------------------------------------------------------------------
spm_mb_ui('results',MB,'scaling');

% characterise intrinsic coupling in terms of transfer functions
%--------------------------------------------------------------------------
spm_mb_ui('results',MB,'kernels');

% display the results in terms of particular partitions and eigenmodes
%--------------------------------------------------------------------------
spm_mb_ui('results',MB,'dynamics');

% characterise extrinsic coupling with a connectogram
%--------------------------------------------------------------------------
spm_mb_ui('results',MB(1:end - 1),'connectogram');

% characterise extrinsic coupling in terms of cross covariance functions
%--------------------------------------------------------------------------
spm_mb_ui('results',MB(1:end - 1),'connectivity');

% characterise intrinsic coupling in terms of dissipative flow
%--------------------------------------------------------------------------
spm_mb_ui('results',MB,'eigenmodes');

% characterise eigenmodes in terms of design or inputs
%--------------------------------------------------------------------------
spm_mb_ui('results',MB,'responses');

% input effects as active states at base level
%--------------------------------------------------------------------------
spm_mb_ui('results',MB,'inputs');


return


% Anti-symmetry (i.e., the emergence of anti-correlations)
%==========================================================================
for i = 1:numel(MB)
    J    = full(MB{i}.J);
    A(i) = norm(J - J')/norm(J + J');
end

% Anti-correlations
%==========================================================================
for i = 1:numel(MB)
    
    % Jacobian and cross covariance function
    %----------------------------------------------------------------------
    dfdx      = full(MB{i}.J);
    [ccf,pst] = spm_ssm2ccf(dfdx);
    cor       = spm_ccf2cor(ccf);
    
    %  predicted and recorded correlations (functional connectivity)
    %----------------------------------------------------------------------
    C{i} = spm_cov2corr(real(cor));
    R{i} = corr(real(MB{i}.Y));
    
end




