function [MB] = spm_mb_ui(action,varargin)
% VOI extraction of adjusted data and Markov Blanket decomposition
% FORMAT [MB] = spm_mb_ui('specify',SPM)
% FORMAT [MB] = spm_mb_ui('blocking',MB)
% FORMAT [MB] = spm_mb_ui('results' ,MB)
%
% SPM    - structure containing generic analysis details
%
% MB.contrast -  contrast name
% MB.name     -  MB name
% MB.c        -  contrast weights
% MB.X        -  contrast subspace
% MB.Y        -  whitened and adjusted data
% MB.X0       -  null space of contrast
%
% MB.XYZ      -  locations of voxels (mm)
% MB.xyz      -  seed voxel location (mm)
% MB.VOX      -  dimension of voxels (mm)
%
% MB.V        -  canonical vectors  (data)
% MB.v        -  canonical variates (data)
% MB.W        -  canonical vectors  (design)
% MB.w        -  canonical variates (design)
% MB.C        -  canonical contrast (design)
%
% MB.chi      -  Chi-squared statistics testing D >= i
% MB.df       -  d.f.
% MB.p        -  p-values
%
% also saved in MB_*.mat in the SPM working directory
%
% FORMAT [MB] = spm_cva_ui('results',MB)
% Display the results of a MB analysis
%__________________________________________________________________________
%
% This routine 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.
%
% The following lists the various results options.  please see main body of
% this script for a description of the (graphical) output
% 
% 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,'connectogram');
% 
% characterise extrinsic coupling in terms of cross covariance functions
%--------------------------------------------------------------------------
% spm_mb_ui('results',MB,'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');
%__________________________________________________________________________

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


OPT.d    = 32;                         % maximum connection length (mm)
OPT.np   = 1024;                       % number of parcels (particles)
OPT.xyz  = [0;0;0];                    % centre of (spherical) VOI (mm)
OPT.spec = 128;                        % radius of (spherical) VOI (mm)
OPT.Ic   = 1;                          % contrast (for adjusted data)
OPT.T    = 1;                          % threshold for adiabatic reduction
OPT.N    = 8;                          % max modes for adiabatic reduction

switch lower(action)
    
    case 'specify'
        %==================================================================
        %                      M B  :  S P E C I F Y
        %==================================================================
        if isempty(varargin)
            [SPM,sts] = spm_select(1,'mat','Select SPM',[],[],'^SPM.*\.mat$');
            if ~sts, return; end
        else
            SPM = varargin{1};
        end
        if ischar(SPM)
            SPM = load(SPM);
            SPM = SPM.SPM;
        end
        
        %-Contrast specification
        %------------------------------------------------------------------
        con    = SPM.xCon(OPT.Ic).name;
        c      = SPM.xCon(OPT.Ic).c;
        c      = full(c);
               
        %-Extract required data from results files
        %==================================================================
        
        %-Get (normalised) explanatory variables (data)
        %------------------------------------------------------------------
        cd(SPM.swd)
        
        XYZ    = SPM.xVol.XYZ;
        F      = spm_get_data(SPM.xCon(1).Vspm,XYZ);
        XYZ    = XYZ(:,logical(F > 1));
        Y      = spm_get_data(SPM.xY.VY, XYZ);
        R      = spm_get_data(SPM.VResMS,XYZ);
        Y      = bsxfun(@times,Y,1./sqrt(R));
        
        XYZmm  = SPM.xVol.M(1:3,1:3)*XYZ;
        XYZmm  = bsxfun(@plus,XYZmm,SPM.xVol.M(1:3,4));

        %-Remove serial correlations and get design (note X := W*X)
        %------------------------------------------------------------------
        Y      = SPM.xX.W*Y;
        X      = SPM.xX.xKXs.X;
        M      = SPM.xVol.M(1:3,1:3);                 %-voxels to mm
        VOX    = diag(sqrt(diag(M'*M))');             %-voxel size
        
        %-Null-space
        %------------------------------------------------------------------
        X0      = [];
        try, X0 = [X0 blkdiag(SPM.xX.K.X0(:,1:min(16,end)))]; end  %-drift terms
        try, X0 = [X0 spm_detrend(SPM.xGX.gSF)];              end  %-global estimate
        X0      = full(spm_svd([X0, (X - X*c*pinv(c))]));          %-null space of c
        Y       = Y - X0*(pinv(X0)*Y);
                
        % exogenous inputs (decimated)
        %------------------------------------------------------------------
        for i = 1:numel(SPM.Sess.U)
            U(:,i)  = SPM.Sess.U(i).u(33:end,:);
            name{i} = SPM.Sess.U(i).name{1};
        end
        Dy      = spm_dctmtx(size(Y,1),size(Y,1));
        Du      = spm_dctmtx(size(U,1),size(Y,1));
        Dy      = Dy*sqrt(size(Y,1)/size(U,1));
        U       = Dy*(Du'*U);
        
        %-First level partition: compact support
        %==================================================================
        spm_figure('GetWin','Markov blanket');clf
        subplot(2,2,1), spm_mip(std(Y),XYZmm,VOX);
        axis image, axis off, title('Variance of voxels','FontSize',16)
                
        % local SVD
        %------------------------------------------------------------------
        V       = spm_mvb_U(Y,'compact',X0,XYZ,[],OPT.np);
        Y       = Y*V;
        nv      = size(Y,2);
        x       = cell(3,nv);      
        x(2,:)  = num2cell(1:nv);                 % indices of partition
        z       = num2cell(1:nv);                 % indices of states

        
        % distance matrix (mm)
        %------------------------------------------------------------------
        [R,S]   = spm_mb_distance(V,XYZmm);
        
        
        subplot(2,2,2), spm_mip(abs(V)*std(Y)',XYZmm,VOX);
        axis image, axis off, title('Variance of modes','FontSize',16)
        subplot(2,2,4), spm_mip(logical(V)*std(Y)',XYZmm,VOX);
        axis image, axis off, title('Variance of particles','FontSize',16)
        subplot(2,2,3), imagesc(R), axis image, axis off
        title('Disttance between particles','FontSize',16), drawnow
        
        % self inhibition based on spatial scale
        %------------------------------------------------------------------
        I       = -S^(-3/2)*16;
        
        % first order Jacobian
        %==================================================================
        [J,K]   = spm_dcm_J(Y,U,X0,SPM.xY.RT,R,I,OPT.d);
        
        % Populate first level of market blanket structure
        %==================================================================
        clear MB
        
        %-Assemble results
        %------------------------------------------------------------------
        MB{1}.swd  = SPM.swd;                    % results directory
        MB{1}.dt   = SPM.xY.RT;                  % sampling interval
        MB{1}.con  = con;                        % contrast name
        MB{1}.name = name;                       % input name
        MB{1}.XYZ  = XYZmm;                      % locations of voxels (mm)
        MB{1}.VOX  = VOX;                        % dimension of voxels (mm)
        
        MB{1}.J = J;                             % Jacobian (states)
        MB{1}.K = K;                             % Jacobian (inputs)
        MB{1}.z = z;                             % indices of states
        MB{1}.x = x;                             % indices of states
        MB{1}.s = num2cell(-abs(diag(J)));       % Lyapunov exponents
        MB{1}.U = U;                             % exogenous inputs
        MB{1}.V = V;                             % support in voxel space
        MB{1}.W = sparse(logical(abs(V)));       % support in voxel space
        MB{1}.Y = Y;                             % time-series
        
    case 'blocking'
        %==================================================================
        %                      M B  :  B L O C K I N G
        %==================================================================
        MB    = varargin{1};
        N     = 4;                               % maximum number of scales
        m     = ones(N,1);                       % # internal states
        
        % preclude long range coupling
        %------------------------------------------------------------------
        % R               = cell(N,1);
        % R{1}            = spm_mb_distance(MB{1}.V,MB{1}.XYZ);
        % R{1}(R{1} > 48) = 0;
        
        %% group renormalization (course graining or blocking)
        %------------------------------------------------------------------
        for i = 1:N
            
            % Markov blanket (particular) partition
            %--------------------------------------------------------------
            x = spm_Markov_blanket(MB{i}.J,MB{i}.z,m(i));
            
            % break if no states
            %--------------------------------------------------------------
            if size(x,2) < 1, break,  end
            
            % dimension reduction (eliminating internal states)
            %--------------------------------------------------------------
            [J,z,v,s]   = spm_A_reduce(MB{i}.J,x,OPT.T,OPT.N);
            MB{i + 1}.J = J;
            MB{i + 1}.z = z;
            MB{i + 1}.s = s;
            MB{i + 1}.x = x;
            MB{i + 1}.W = sparse(logical(abs(MB{i}.V)));
            
            % eigenmodes and variates (V and Y)
            %--------------------------------------------------------------
            for j = 1:numel(z)
                u                = [MB{i + 1}.x{1:2,j}];     % blanket states
                w                = MB{i + 1}.z{j};           % new states
                MB{i + 1}.V(:,w) = MB{i}.V(:,u)*v{j};
                MB{i + 1}.Y(:,w) = MB{i}.Y(:,u)*v{j};
            end
            
            % functional anatomy of states
            %--------------------------------------------------------------
            spm_figure('getwin',sprintf('Markov level %i',i)); clf;
            spm_mb_anatomy(MB{i},MB{1}.XYZ,MB{1}.VOX);
            
            % break if a single particle or parcel
            %--------------------------------------------------------------
            if size(MB{i}.x,2) < 2, break,  end
            
        end
        
        
        
        % Save results
        %==================================================================
        
        %-Save
        %------------------------------------------------------------------
        save(fullfile(MB{1}.swd,['MB_' MB{1}.con]),'MB');
        
        
        
    case 'results'
        %==================================================================
        %                      M B  :  R E S U L T S
        %==================================================================
        
        %-Get MB and ananlysis
        %------------------------------------------------------------------
        MB       = varargin{1};
        analysis = varargin{2};
        if ischar(MB)
            MB  = load(MB);
            MB  = MB.MB;
        end
        
        switch analysis
            
            case('anatomy')
                
%               This figure illustrates the partition of states at the
%               first and subsequent scales. The upper right panel shows
%               all the constituent particles as a maximum intensity
%               projection, where the spatial support of each particle has
%               been colour-coded according to the variance explained by
%               its eigenmode. The upper middle panel shows the
%               corresponding adjacency matrix or coupling among particles.
%               The coloured circles encode the identity of each particle.
%               In this instance, the particles have been arranged in order
%               of descending dissipation (i.e., the principal eigenvalue
%               of each particles Jacobian). The upper right panel shows
%               these eigenvalues above the corresponding particle (encoded
%               by coloured dots) in terms of rate constants (i.e., the
%               negative inverse of the eigenvalue of each particle). The
%               remaining panels show the first 12 particles as maximum
%               intensity projections. The colour of the background
%               corresponds to the colour that designates each particle. In
%               this first level, each particle has a single eigenstate.
%               The numbers in brackets above each maximum intensity
%               projection correspond to the number of internal, active and
%               sensory states, respectively, where the active and sensory
%               states comprise blanket states. At this lowest level, every
%               eigenstate is a sensory state because it can influence, and
%               be influenced by,the eigenstates of other particles.
                
                % functional anatomy (space)
                %----------------------------------------------------------
                VOX   = MB{1}.VOX;
                XYZ   = MB{1}.XYZ;
                for i = 1:numel(MB)
                    spm_figure('getwin',sprintf('Markov level %i',i)); clf;
                    spm_mb_anatomy(MB{i},XYZ,VOX)
                end
                
            case('dynamics')
                
%               Intrinsic timescales in the brain: This figure reports
%               intrinsic timescales at intermediate scales. The left
%               column shows the eigenmodes in terms of their principal
%               frequency, i.e., the largest complex eigenvalue (divided by
%               2*pi). The right column shows the corresponding eigenmodes
%               in terms of their principal time constants, i.e., the
%               reciprocal of the largest negative real part. These two
%               characterisations (principal frequency and time constant)
%               speak to different aspects of intrinsic timescales; both of
%               which contribute to the shape of an eigenstate's
%               correlation function of time. The first quantifies the
%               frequency of solenoidal flow, while the second reflects the
%               rate of decay associated with the dissipative flow.
                
                % functional anatomy (time)
                %----------------------------------------------------------
                VOX   = MB{1}.VOX;
                XYZ   = MB{1}.XYZ;
                
                spm_figure('getwin','Timing'); clf;
                N     = numel(MB) - 1;
                for i = 1:N
                    
                    % largest real and imaginary parts of each particle
                    %------------------------------------------------------
                    sr    = zeros(size(MB{i}.V,2),1);
                    si    = zeros(size(MB{i}.V,2),1);
                    for j = 1:numel(MB{i}.s)
                        try
                            z     = MB{i}.z{j};
                            sr(z) = max(real(MB{i}.s{j}));
                            si(z) = max(imag(MB{i}.s{j}));
                        end
                    end
                    
                    % imaginary part (frequency)
                    %------------------------------------------------------
                    si    = si/(2*pi);
                    subplot(N,2,2*(i - 1) + 1),cla,hold off
                    spm_mip(logical(abs(MB{i}.V))*si,XYZ,VOX);
                    axis image, axis off
                    str{1} = sprintf('Frequency - imaginary');
                    str{2} = sprintf('%2.1f to %1.3f Hz',min(si),max(si));
                    title(str,'FontSize',14,'FontWeight','bold')
                    
                    % real part (time constants)
                    %------------------------------------------------------
                    sr     = (-1./sr);
                    subplot(N,2,2*(i - 1) + 2),cla,hold off
                    spm_mip(logical(abs(MB{i}.V))*sr,XYZ,VOX);
                    axis image, axis off
                    str{1} = sprintf('Dissipation - real (scale %i)',i);
                    str{2} = sprintf('%2.1f to %2.1f sec',min(sr),max(sr));
                    title(str,'FontSize',14,'FontWeight','bold')
                    
                end
                 
            case('distance')
                
%               Local connectivity: This figure reports some of the
%               statistical characteristics of dynamical coupling among
%               particles at the first level. The upper left panel plots
%               each connection in terms of the real part of the
%               corresponding Jacobian in Hz, against the distance spanned
%               by the connection (i.e., Euclidean distance between the
%               centres of the two particles). The upper left panel plots
%               the log-coupling (real part) against log-distance. The
%               lower panel plots the strength of reciprocal connections
%               against each other, to illustrate the relative proportions
%               of recurrent excitatory and inhibitory coupling. The
%               rarefied region in the centre of this scatterplot reflects
%               the fact that connections wIf youith small coupling strengths
%               have been eliminated during Bayesian model reduction.
                
                % spatial distance and (first-level) coupling
                %----------------------------------------------------------
                spm_figure('getwin','Distance rules'); clf;
                col = {'r','b'};
                
                % get distance (from averge location, xyz)
                %----------------------------------------------------------
                J     = MB{1}.J;
                nv    = length(J);
                D     = abs(MB{1}.V);
                D     = bsxfun(@rdivide,D,sum(D));
                xyz   = MB{1}.XYZ*D;
                xyz(1,:) = abs(xyz(1,:));
                
                % characterise excitation and inhibitory coupling separately
                %----------------------------------------------------------
                u     = 1/1024;
                for e = 1:2
                    D = [];
                    E = [];
                    for i = 1:nv
                        for j = 1:nv
                            d = xyz(:,i) - xyz(:,j);
                            d = sqrt(d'*d);
                            if e > 1
                                if J(i,j) < -u && d
                                    D(end + 1,1) = sqrt(d'*d);
                                    E(end + 1,1) = abs(J(i,j));
                                end
                            else
                                if J(i,j) > +u && d
                                    D(end + 1,1) = sqrt(d'*d);
                                    E(end + 1,1) = abs(J(i,j));
                                end
                            end
                        end
                    end
                    
                    % coupling (cumulative)
                    %----------------------------------------------------------
                    n     = 16;
                    c     = zeros(n,1);
                    k     = zeros(n,1);
                    d     = linspace(8,max(D),n + 1);
                    for i = 1:n
                        j      = D > d(i) & D < d(i + 1);
                        c(i,1) = mean(E(j));                  % mean coupling
                        k(i,1) = sum(j);                      % numer of edges
                    end
                    d     = spm_vec(d(1:n));
                    j     = k > mean(k)/4;
                    d     = d(j);
                    c     = c(j);
                    k     = k(j);
                    
                    % regression
                    %----------------------------------------------------------
                    [F,df,B] = spm_ancova([ones(size(d)) log(d)],[],log(c));
                    
                    subplot(2,2,1)
                    plot(D,E,['.' col{e}],'MarkerSize',1), hold on
                    plot(d,c,['o' col{e}],'MarkerSize',8)
                    plot(d,exp(B(1))*d.^B(2),col{e},'LineWidth',2)
                    title('Distance scaling rule','FontSize',16)
                    xlabel('Distance (mm)'), ylabel('coupling strength (Hz}')
                    axis square, box off, set(gca,'YLim',[0, 1])
                    
                    subplot(2,2,2)
                    plot(log(d),log(c),['o' col{e}],'MarkerSize',8), hold on
                    plot(log(d),B(1) + B(2)*log(d),col{e},'LineWidth',1)
                    if e == 1
                        str = sprintf('Log scaling (%2.2f',B(2));
                    else
                        str = [str sprintf(', %2.2f)',B(2))];
                    end
                    title(str,'FontSize',16)
                    xlabel('log distance'), ylabel('log coupling')
                    axis square, box off
                    legend('excitatory',' ','inhibitory',' '),legend(gca,'boxoff')
                    
                    % subplot(2,2,3)
                    % p    = k./(d.^2);
                    % p    = p/sum(p);
                    % plot(d,p,col{e}), hold on
                    % title('Density','FontSize',16)
                    % xlabel('distance (mm)'), ylabel('Density')
                    % axis square, box off
                    
                end
                
                % asymmetry
                %----------------------------------------------------------
                D     = [];
                E     = [];
                for i = 1:nv
                    for j = 1:nv
                        if abs(J(i,j)) > exp(-8) && j ~= i
                            D(end + 1,1) = J(j,i);
                            E(end + 1,1) = J(i,j);
                        end
                    end
                end
                
                subplot(2,1,2)
                hold off, plot(E,D,'.r','MarkerSize',1)
                title('Reciprocal coupling','FontSize',16)
                xlabel('coupling (Hz)'), ylabel('coupling (Hz)')
                axis square, box off, axis([-1 1 -1 1])
                
                % percentages of different recurrent connections
                %----------------------------------------------------------
                R(1) = sum(D > 0 & E > 0);
                R(2) = sum(D < 0 & E > 0);
                R(3) = sum(D > 0 & E < 0);
                R(4) = sum(D < 0 & E < 0);
                R    = 100*R/sum(R);
                
                text( 1/4, 1/4,sprintf('%2.0f p.c.',R(1)),'FontWeight','bold');
                text( 1/4,-1/4,sprintf('%2.0f p.c.',R(2)),'FontWeight','bold');
                text(-1/4, 1/4,sprintf('%2.0f p.c.',R(3)),'FontWeight','bold');
                text(-1/4,-1/4,sprintf('%2.0f p.c.',R(4)),'FontWeight','bold');
                
                
            case('scaling')
                
%               Scale invariance: This figure illustrates scaling behaviour
%               across the scales of the particular decomposition. The
%               upper panel plots the real part of the eigenvalues of each
%               particle against its spatial scale; namely the Calliper
%               width of the particle's eigenmode. This is replicated for
%               each of the four scales, denoted by the different colours
%               (green, pink, cyan and puce, respectively). The expected
%               values are shown as encircled large dots. The lower left
%               panel plots the logarithms of these temporal and spatial
%               expectations against each other. The resulting regression
%               slope corresponds to the scaling exponent.  The light grey
%               circles correspond to what would have been seen at higher
%               and lower scales. The lower right plot shows the same
%               regression in terms of the implicit time constant, as a
%               function of spatial scale expressed in millimetres.
                
                % time constants over scales
                %----------------------------------------------------------
                spm_figure('getwin','Scaling'); clf;
                
                N     = numel(MB);
                SS    = cell(N,1);
                SR    = cell(N,1);
                for i = 1:N
                    for j = 1:numel(MB{i}.s)
                        
                        % get spatial scale (standard deviation)
                        %--------------------------------------------------
                        S     = sum(logical(abs(MB{i}.V(:,MB{i}.z{j})))).^(1/3);
                        SS{i} = [SS{i}; mean(S)];
                        
                        % get temporal scale (eigenvalue)
                        %--------------------------------------------------
                        S     = real(MB{i}.s{j});
                        SR{i} = [SR{i}; mean(S)];
                    end
                end
                
                
                % Expectations and graphics
                %----------------------------------------------------------
                for i = 1:N
                    j       = isfinite(SS{i}); 
                    SS{i}   = SS{i}(j);
                    ss(i,1) = mean(SS{i});
                    j       = isfinite(SR{i});
                    SR{i}   = SR{i}(j);
                    sr(i,1) = mean(SR{i});
                end
                
                subplot(2,1,1), hold off
                col   = spm_MB_col(N);
                for i = 1:N
                    plot(SS{i},SR{i},'.','MarkerSize', 8,'Color',col{i}), hold on
                    plot(ss(i),sr(i),'.','MarkerSize',32,'Color',col{i})
                    plot(ss(i),sr(i),'o','MarkerSize',32,'Color',col{i})
                end
                title('Scale invariance','FontSize',16)
                xlabel('spatial scale (mm)'), ylabel('eigenvalue (Hz)')
                axis square, box off, hold off
                
                % regression (Time)
                %----------------------------------------------------------
                ii        = (1:N)'; % scale
                [F,df,Br] = spm_ancova([ones(N,1) ii],     [],-log(-sr));
                [F,df,Bs] = spm_ancova([ones(N,1) ii],     [], log( ss));
                [F,df,A ] = spm_ancova([ones(N,1) log(ss)],[],-log(-sr));
                A         = [A(1), Br(2)/Bs(2)];
                
                fprintf('Beta (exp) time:  %2.2f (%2.2f)\n',Br(2),exp(Br(2)));
                fprintf('Beta (exp) space: %2.2f (%2.2f)\n',Bs(2),exp(Bs(2)));
                fprintf('Alpha: %2.2f\n',A(2));

                % extrapolate
                %----------------------------------------------------------
                bs  = exp(Bs(2)*(-4:4:8) + Bs(1))';
                br  = exp(A(1) + A(2)*log(bs));
                S   = linspace(min(bs),max(bs),64);
                
                subplot(2,2,3), cla, hold on
                plot(log(bs),log(br),'.','MarkerSize',32,'Color',[1,1,1]*.8)
                plot(log(bs),log(br),'o','MarkerSize',32,'Color',[1,1,1]*.8)           
                for i = 1:N
                    plot(log(SS{i}),-log(-SR{i}),'.','MarkerSize',4,'Color',col{i})
                    plot(log(ss(i)),-log(-sr(i)),'.','MarkerSize',32,'Color',col{i})
                    plot(log(ss(i)),-log(-sr(i)),'o','MarkerSize',32,'Color',col{i})
                    hold on
                end
                plot(log(S),A(1) + A(2)*log(S))
                
                title(sprintf('Log scaling - %2.2f',A(2)),'FontSize',16)
                xlabel('log distance (mm)'), ylabel('log time constant')
                axis square, box off, hold off
                
                subplot(2,2,4), hold off
                plot(S,exp(A(1) + A(2)*log(S))), hold on
                plot(bs,br,'.','MarkerSize',32,'Color',[1,1,1]*.8)
                plot(bs,br,'o','MarkerSize',32,'Color',[1,1,1]*.8)
                title('Scaling behaviour','FontSize',16)
                xlabel('distance (mm)'), ylabel('time constant (sec)')
                axis square, box off
                
                % some canonical scales
                %----------------------------------------------------------
                disp('mm')
                disp(num2str(bs,'%2.2e'))
                disp('Seconds')
                disp(num2str(br,'%2.2e'))
                for i = 1:numel(bs)
                    hold on, plot([bs(i),bs(i)],[0 br(i)],'r-.')
                end
                
            case('kernels')
                
%               Transfer functions: This figure characterises the dynamics
%               at successive scales in terms of transfer functions, as
%               quantified by the complex eigenvalues (c.f., a pole-zero
%               map). The left column shows the transfer functions over
%               frequency for all particles with a complex eigenvalue (at
%               successive scales). These eigenvalues are shown in the
%               right column in the complex number plane. As we ascend from
%               one scale to the next, the real part of the eigenvalue
%               approaches zero from the left and the number of eigenvalues
%               falls with the coarse-graining, i.e., number of particles.
%               The complex part of the eigenvalues corresponds to the peak
%               frequency of the associated transfer functions, while the
%               dispersion around this peak decreases as the real part
%               approaches zero. The emergence of spectral peaks in the
%               transfer functions inherit from the complex part of the
%               eigenvalues that emerge under asymmetric coupling with
%               solenoidal flow.
                
                % frequencies (Hz)
                %----------------------------------------------------------
                w    =  linspace(0,1/8,128);
                
                % transfer functions over scales
                %----------------------------------------------------------
                spm_figure('getwin','Transfer functions'); clf;
                for i = 2:min(4,numel(MB))
                    
                    % loop over particles
                    %------------------------------------------------------
                    nz    = numel(MB{i}.z);
                    col   = spm_MB_col(nz);
                    for p = 1:nz
                        
                        % loop over eigenmodes
                        %--------------------------------------------------
                        for j = 1:numel(MB{i}.s{p})
                            
                            % condition unstable eigenmodes
                            %----------------------------------------------
                            s   = MB{i}.s{p}(j);
                            if real(s) > 0
                                s  = 1j*imag(s) - 1/64;
                            end
                            
                            % Transfer functions (if imaginary)
                            %----------------------------------------------
                            if imag(s) ~= 0
                                
                                % Transfer function
                                %------------------------------------------
                                S      = 1./(1j*2*pi*w - s);
                                
                                subplot(3,2,(i - 2)*2 + 1), hold on
                                plot(w,abs(S),'color',col{p})
                                
                                % Pole-zero map
                                %------------------------------------------
                                subplot(3,2,(i - 2)*2 + 2), hold on
                                plot(MB{i}.s{p}(j),'.','MarkerSize',32,'color',col{p})
                                
                                
                            end
                        end
                    end
                    
                    subplot(3,2,(i - 2)*2 + 1)
                    title(sprintf('Scale %i',i - 1),'fontsize',16)
                    xlabel('frequency (Hz)'), ylabel('transfer tunctions')
                    axis square, box off, spm_axis tight, hold off
                    
                    subplot(3,2,(i - 2)*2 + 2)
                    plot([0,0],[-1,1],':')
                    title('Eigenvalues','fontsize',16)
                    xlabel('real (Hz)'), ylabel('imaginary (Hz)')
                    axis square, box off, axis([-1 1 -1/8 1/8]), hold off

                end
                
                
            case('connectivity')
                
%               Dynamic coupling: This figure characterises the coupling
%               between the two eigenstates, which (unless specified)
%               correspond to the strongest coupling at this highest scale.
%               This coupling is mediated by the corresponding element of
%               the (complex) Jacobian, circled in red in the upper middle
%               panel. The flanking panels on the left and right show the
%               corresponding eigenmodes in voxel space. The middle row
%               shows the auto-covariance functions of the two eigenstates,
%               illustrating serial correlations that can last for many
%               seconds. The lower panels report the cross-covariance
%               function between the two eigenstates, over 256 seconds
%               (lower left panel) and focusing on cross-covariances over
%               32 seconds (lower right panel). Asymmetrical
%               cross-covariance (and implicitly cross-correlation)
%               function reflects the solenoidal coupling and the implicit
%               breaking of detailed balance the particular decomposition
%               accommodates.
                
                
                % functional anatomy
                %----------------------------------------------------------
                VOX   = MB{1}.VOX;
                XYZ   = MB{1}.XYZ;
                spm_figure('getwin',' extrinsic connectivity'); clf;
                spm_mb_anatomy(MB{end},XYZ,VOX)
                
                
                % Jacobian and cross covariance function
                %----------------------------------------------------------
                dfdx      = full(MB{end}.J);
                [ccf,pst] = spm_ssm2ccf(dfdx);
                
                % get pair of eigenvariates
                %----------------------------------------------------------
                try
                    i     = varargin{3};
                    j     = varargin{4};
                catch
                    A     = dfdx;
                    A     = A - diag(diag(A));
                    [d,i] = max(real(A(:)));
                    [i,j] = ind2sub(size(A),i);
                end
                cf        = ccf(:,i,j);
                lim       = [min(cf) max(cf)];
                
                % if no correlations than plot conjugate modes
                %----------------------------------------------------------
                if max(abs(lim)) < 1e-6
                    [i,j] = sort(imag(spm_vec(MB{end}.s)),'descend');
                    i     = j(1);
                    j     = j(end);
                end
                
                % graphics
                %----------------------------------------------------------
                subplot(3,2,3), plot(pst,ccf(:,i,i))
                title(sprintf('covariance: state %i',i),'FontSize',16)
                xlabel('lag (sec)'),ylabel('covariance')
                axis square, box off, spm_axis tight

                subplot(3,2,4), plot(pst,ccf(:,j,j))
                title(sprintf('covariance: state %i',j),'FontSize',16)
                xlabel('lag (sec)'),ylabel('covariance')
                axis square, box off, spm_axis tight
                
                subplot(3,2,5), plot(pst,ccf(:,i,j),[0,0],lim,':b')
                title(sprintf('cross-covariance: states %i,%i',i,j),'FontSize',16)
                xlabel('lag (sec)'),ylabel('covariance')
                if max(abs(lim)) < 1e-6, set(gca,'YLim',[-1,1]), end
                axis square, box off, spm_axis tight

                % blow-up around zero lag
                %----------------------------------------------------------
                [d,m] = max(abs(ccf(:,i,j)));
                m     = pst(m);
                
                subplot(3,2,6), plot(pst,ccf(:,i,j),[0,0],lim,':b',[m,m],lim,'-.r')
                title(sprintf('cross-covariance: states %i,%i',i,j),'FontSize',16)
                xlabel('lag (sec)'),ylabel('covariance')
                if max(abs(lim)) < 1e-6, set(gca,'YLim',[-1,1]), end
                axis square, box off, spm_axis tight, set(gca,'XLim',[-1,1]*32)
                
                
                % modes
                %----------------------------------------------------------
                subplot(3,3,1),cla,hold off
                spm_mip(abs(MB{end}.V(:,i)),XYZ,VOX);
                axis image, axis off
                title(sprintf('Eigenmode %i',i),'FontSize',16)
                
                subplot(3,3,3),cla,hold off
                spm_mip(abs(MB{end}.V(:,j)),XYZ,VOX);
                axis image, axis off
                title(sprintf('Eigenmode %i',j),'FontSize',16)
                
                subplot(3,3,2), hold on, plot(j,i,'or','MarkerSize',16)
                
            case('eigenmodes')
                
%               Dissipative and solenoidal dynamics: This figure unpacks
%               the intrinsic coupling at the final (usually, single
%               particle) level. At this level, there are can be no
%               coupling between particles and (by construction) the
%               dynamics is completely characterised in terms of the
%               eigenstates that comprise the particle. In turn, these are
%               completely characterised by their complex eigenvalues;
%               namely, the intrinsic complex coupling. The upper panels
%               show the dissipative and solenoidal (kinetic) energy of the
%               eigenstates that comprise the particle. The corresponding
%               eigenmodes are shown in the subsequent panels as maximum
%               intensity projections (of their absolute values).
                
                % functional anatomy
                %----------------------------------------------------------
                spm_figure('getwin',' extrinsic coupling'); clf;
                VOX = MB{1}.VOX;
                XYZ = MB{1}.XYZ;
                

                % Jacobian
                %----------------------------------------------------------
                J = full(MB{end}.J);
                Y = full(MB{end}.Y);
                
                % complex form
                %----------------------------------------------------------
                J = diag(J);
                D = -real(J);                        % dissipative
                Q = (imag(J).^2)./D;                 % solenoidal
                
                % show dynamics
                %==========================================================
                try, dt = MB{1}.dt;    catch, dt    = 1;      end
                try,  i = varargin{3}; catch, [d,i] = max(Q); end
                
                
                subplot(4,2,1), bar(D)
                xlabel('eigenmode'),ylabel('kinetic energy'), axis square, box off
                title('Dissipative energy','FontSize',16)
                
                subplot(4,2,2), bar(Q)
                xlabel('eigenmode'),ylabel('kinetic energy'), axis square, box off
                title('Solenoidal energy','FontSize',16)
                
                % pst     = (1:size(Y,1))*dt;
                % subplot(3,2,4), plot(pst,abs(Y(:,i)),'-',pst,angle(Y(:,i)),':')
                % xlabel('time'),ylabel('state'), axis square, box off
                % legend({'amplitude','angle'}), legend(gca,'boxoff')
                % title('Dissipative dynamics','FontSize',16)
                
                % mode
                %----------------------------------------------------------
                for i = 1:min(9,numel(D))
                    subplot(4,3,i + 3),cla,hold off
                    spm_mip(abs(MB{end}.V(:,i)),XYZ,VOX);
                    axis image, axis off
                    title(sprintf('Eigenmode %i',i),'FontSize',12,'FontWeight','bold')
                end

                
            case('responses')
                
%               Induced responses: This figure illustrates the expression
%               of experimental or condition-specific effects at different
%               scales of the particular decomposition. The top panel is an
%               unusual form of statistical parametric mapping; namely, and
%               image of the F statistic, testing for the significance of
%               an effect of any of any exogenous inputs. Each row of the F
%               statistic map corresponds to a scale and comprises the F
%               statistic for each successive particle at that scale. This
%               map shows that the effect of (some linear mixture of)
%               exogenous inputs can be detected at several scales, as
%               evidenced by the dark bars.  unless specified, the most
%               significant eigenmode is shown on the lower left in voxel
%               space. It's expression over time (in terms of its real
%               value) is depicted in the middle panel (blue line), with
%               the best fitting prediction based upon exogenous input
%               (green line). This prediction is a contrast (i.e., linear
%               mixture) of the input functions shown in the design matrix
%               on the lower right. The coefficients of this contrast are
%               shown below the design matrix. The last column of the
%               design matrix is simply a column of ones.
                
                % functional anatomy
                %----------------------------------------------------------
                VOX   = MB{1}.VOX;
                XYZ   = MB{1}.XYZ;
                spm_figure('getwin','responses'); clf;
                                
                % ANOVA
                %----------------------------------------------------------
                U           = MB{1}.U;
                n           = size(U,2);
                U(:,n + 1)  = 1;
                c           = spm_speye(n + 1,n);
                
                % search over scles and eigenmodes
                %----------------------------------------------------------
                for j = 1:numel(MB)
                    for i = 1:size(MB{j}.Y,2)
                        Y     = real(MB{j}.Y(:,i));
                        F(j,i) = spm_ancova(U,[],Y,c);
                    end
                end
                
                % most significant mode
                %----------------------------------------------------------
                [Fj,j]       = max(F,[],2);
                [Fi,i]       = max(Fj);
                try, i       = varargin{3}; end
                j            = j(i);
                Y            = real(MB{i}.Y(:,j));
                [Fi,df,beta] = spm_ancova(U,[],Y,c);
                p            = 1 - spm_Fcdf(Fi,df);
                df           = round(df);
                
                % show dynamics
                %==========================================================
                try, dt = MB{1}.dt; catch, dt = 1; end
                pst     = (1:size(Y,1))*dt;
                str     = 'Eigenvariate and prediction';
                str     = sprintf('%s: F(%d,%d) = %2.1f, p = %2.3f',str,df(1),df(2),Fi,p);
                
                subplot(3,1,1), imagesc(1 - F)
                xlabel('eigenmode'),ylabel('scale'), spm_axis tight, box off
                title('F-statistic','FontSize',16), hold on
                plot(j,i,'or','MarkerSize',8)
                
                subplot(3,1,2), plot(pst,Y,':',pst,U*beta)
                xlabel('time'),ylabel('real part'), spm_axis tight, box off
                title(str,'FontSize',16)

                % mode
                %----------------------------------------------------------
                subplot(3,2,5),cla,hold off
                spm_mip(abs(MB{i}.V(:,j)),XYZ,VOX);
                axis image, axis off
                title(sprintf('Eigenmode %i',j),'FontSize',16)
                
                % inputs and contribution
                %----------------------------------------------------------
                subplot(6,2,10), imagesc(U)
                xlabel('input'),ylabel('time'), axis square, box off
                title('Inputs','FontSize',16)
                
                subplot(6,2,12), bar(beta)
                xlabel('input'),ylabel('contribution'), axis square, box off
                spm_axis tight
                
                
            case('connectogram')

%               Extrinsic connectivity: This figure illustrates asymmetric
%               extrinsic (between particle) coupling at the final scale
%               supplied in MB. The upper panels reproduce the results in
%               'anatomy', while the lower panel is a connectogram
%               illustrating the coupling among eigenstates that constitute
%               the  particles at this  scale. The width of each connector
%               reflects the strength of the coupling; after dividing the
%               strength into five bins and eliminating the lowest bin. The
%               colour of the dots corresponds to the colour of the
%               particle in the upper right panel. The colour of the
%               connectors corresponds to the source of the strongest
%               (reciprocal) connection. The coupling strength corresponds
%               to the real part of the Jacobian, in Hz.
                
                % functional anatomy of penultimate levels
                %----------------------------------------------------------
                VOX   = MB{1}.VOX;
                XYZ   = MB{1}.XYZ;
                i     = numel(MB);
                spm_figure('getwin',sprintf('connectogram level %i',i)); clf;
                spm_mb_anatomy(MB{i},XYZ,VOX,0)

                
                % assign colours to connections (afferent
                %----------------------------------------------------------
                MB   = MB{i};
                nz   = numel(MB.z);            % number of vector states
                col  = spm_MB_col(nz);         % colours
                COL  = [];
                for i = 1:nz
                    for j = 1:numel(MB.s{i})
                        COL(end + 1,:) = col{i}';
                    end
                end
                
                % connectivity of undirected connectivity
                %----------------------------------------------------------
                subplot(2,1,2)
                spm_circularGraph(MB.J,'colormap',COL)
                axis square off
                title('Connectogram','FontSize',16)
                
           case('inputs')
               
%              Induced responses over space and time: this figure
%              characterises induced responses in terms of first order
%              Volterra kernels (i.e., impulse response functions) of
%              particles at the first (finest) scale of course graining.
%              Each row corresponds to the inputs specified. The left
%              column shows the expression of these inputs over particles
%              (weighted by the absolute value of their eigenmodes). This
%              effect is the variance attributable to each input (i.e.,
%              square of the corresponding kernel, summed over time), shown
%              in the left row. These kernels are shown for the 32
%              particles with the greatest (absolute) magnitude.
                
               % functional anatomy (1st order kernels)
               %----------------------------------------------------------
               spm_figure('getwin','input effects'); clf;
               VOX   = MB{1}.VOX;
               XYZ   = MB{1}.XYZ;
               nu    = size(MB{1}.U,2);
               
               dfdx  = MB{1}.J;
               dfdu  = MB{1}.K;              
               pst   = (1:64)*4;
               ker   = spm_ssm2ker(dfdx,dfdu,[],pst);
               for i = 1:nu
                   
                   % spatial expression
                   %-------------------------------------------------------
                   subplot(nu,2,2*(i - 1) + 1)
                   K = sum(ker(:,:,i).^2,1);
                   
                   spm_mip(abs(MB{1}.V*K(:)),XYZ,VOX);
                   axis image, axis off
                   title(MB{1}.name{i},'FontSize',16)
                   
                   % temporal expression
                   %-------------------------------------------------------
                   [d,j] = sort(abs(ker(1,:,i)),'descend');
                   
                   subplot(nu,2,2*(i - 1) + 2)
                   plot(pst,ker(:,j(1:64),i))
                   xlabel('time (secs)'),ylabel('response'), axis square, box off
                   title('Neuronal responses','FontSize',16)
                   
               end
               
            otherwise
                error('Unknown action.');
                
        end
        
        
    otherwise
        error('Unknown action.');
        
end


% subroutines
%==========================================================================

% functional anatomy - distance (R)
%==========================================================================
function [R,S] = spm_mb_distance(V,XYZ)
% FORMAT [R,S] = spm_mb_distance(V,XYZ)
% V   - pattern vectors
% XYZ - locations
%
% R   - distance between patterns
% S   - spatial dispersion of patterns (sd; mm)
%__________________________________________________________________________

% distance matrix (mm)
%--------------------------------------------------------------------------
nv    = size(V,2);
R     = zeros(nv,nv);
D     = abs(V);
D     = bsxfun(@rdivide,D,sum(D));
xyz   = XYZ*D;
S     = (XYZ.^2)*D - (XYZ*D).^2;
S     = sqrt(mean(S(:)));

% interhemispheric coupling
%--------------------------------------------------------------------------
xyz(1,:) = abs(xyz(1,:));

for i = 1:nv
    for j = i:nv
        d      = xyz(:,i) - xyz(:,j);
        d      = sqrt(d'*d);
        R(i,j) = d;
        R(j,i) = d;
    end
end


% functional anatomy
%==========================================================================
function spm_mb_anatomy(MB,XYZ,VOX,N)
%__________________________________________________________________________


% set-up
%--------------------------------------------------------------------------
if nargin < 4, N = 12; end

nz   = numel(MB.z);                             % number of particles
ns   = size(MB.W,2);                            % number of states (n - 1)
nx   = size(MB.J,2);                            % number of states (n)
col  = spm_MB_col(nz);                          % colours

% functional anatomy of states
%--------------------------------------------------------------------------
subplot(3,3,1),cla,hold off
c    = var(abs(MB.Y));
spm_mip(logical(abs(MB.V))*c',XYZ,VOX);
axis image, axis off
str{1} = sprintf('%d states, %d particles',ns,nz);
str{2} = sprintf('(%d eigenstates)',nx);
title(str,'FontSize',16)

subplot(3,3,2),cla,hold off
if nx > 128
    imagesc(~abs(MB.J))
else
    J  = real(MB.J);
    imagesc(max(J,-1/8))
end
axis square
title({'Jacobian';''},'FontSize',16)
xlabel('States'), hold on

subplot(3,3,3),cla,hold off
bar(-1./real(spm_vec(MB.s)))
axis square, set(gca,'XLim',[0,nx + 1])
title('Time constants','FontSize',16)
xlabel('states'),ylabel('seconds')

% label vector states (states of particles previous level)
%--------------------------------------------------------------------------
xi    = 1;
for i = 1:nz
    for j = 1:numel(MB.s{i})
        hold on, subplot(3,3,2), plot(xi,0,'.','Color',col{i},'MarkerSize',24)
        hold on, subplot(3,3,3), plot(xi,0,'.','Color',col{i},'MarkerSize',24)
        xi = xi + 1;
    end
end


% functional anatomy of particles
%--------------------------------------------------------------------------
nn    = min(nz,N);
for j = 1:nn
    
    Z     = any(MB.W(:,MB.x{1,j}),2)*2;          % active states
    Z     = any(MB.W(:,MB.x{2,j}),2)*1 + Z;      % sensory states
    Z     = any(MB.W(:,MB.x{3,j}),2)*3 + Z;      % internal states
    
    
    % get subplot layout
    %----------------------------------------------------------------------
    if nn == 1
        ni = 3; nj = 1;
    elseif nn == 2
        ni = 3; nj = 2;
    elseif nn == 3
        ni = 3; nj = 3;
    elseif nn == 4
        ni = 3; nj = 2;
    elseif nn == 5
        ni = 4; nj = 2;
    elseif nn == 6
        ni = 4; nj = 3;
    else
        ni = 4; nj = 4;
    end
    
    subplot(ni,nj,j + nj)
    [z,i] = max(Z);
    if z < 3, Z(i)  = 3; end
    mip   = spm_mip(Z,XYZ,VOX);
    for k = 1:3
        c = col{j}(k)/2;
        MIP(:,:,k) = (1 - mip/64)*(1 - c) + c;
    end
    image(MIP), axis image, axis off
    str{1} = sprintf('Particle %d of %d (%d,%d,%d)',j,nz,...
        numel(MB.x{3,j}),numel(MB.x{1,j}),numel(MB.x{2,j}));
    str{2} = sprintf('%d eigenmodes',numel(MB.z{j}));
    title(str,'FontWeight','Bold')
end
