In [35]:
# TO DO WITH THIS FUNC: 

# Try it working up to spm_diff
# ...checking vs. spm func using spm debugger breakpoints

# then come up with alternative for spm_diff

# then can do full simulations

In [187]:
import scipy.sparse as sps

SPM fx erp port

code


In [405]:
def fx_erp(x,u,P,M):
    """
    % State equations for a neural mass model of erps
  
    Port of spm_fx_erp.m
  
    Usage: 
  
      f,J,D = spm_fx_erp(x,u,P,M,returnJ=True,returnD=True)
   
    x      - state vector
    x(:,1) - voltage (spiny stellate cells)
    x(:,2) - voltage (pyramidal cells) +ve
    x(:,3) - voltage (pyramidal cells) -ve
    x(:,4) - current (spiny stellate cells)    depolarizing
    x(:,5) - current (pyramidal cells)         depolarizing
    x(:,6) - current (pyramidal cells)         hyperpolarizing
    x(:,7) - voltage (inhibitory interneurons)
    x(:,8) - current (inhibitory interneurons) depolarizing
    x(:,9) - voltage (pyramidal cells)

    f        - dx(t)/dt  = f(x(t))
    J        - df(t)/dx(t)
    D        - delay operator dx(t)/dt = f(x(t - d)) = D(d)*f(x(t))

    Prior fixed parameter scaling [Defaults]

    M.pF.E = [32 16 4];           % extrinsic rates (forward, backward, lateral)
    M.pF.H = [1 4/5 1/4 1/4]*128; % intrinsic rates (g1, g2 g3, g4)
    M.pF.D = [2 16];              % propogation delays (intrinsic, extrinsic)
    M.pF.G = [4 32];              % receptor densities (excitatory, inhibitory)
    M.pF.T = [8 16];              % synaptic constants (excitatory, inhibitory)
    M.pF.S = [1 1/2];             % parameters of activation function


    JG 2017
    """

    
    
    """
    % get dimensions and configure state variables
    %--------------------------------------------------------------------------   
    n = length(P.A{1});         % number of sources
    x = spm_unvec(x,M.x);       % neuronal states
    """;

    n = len(P['A'][0]);        # % number of sources      # TO DO
    
    
    x = spm_unvec(x,M['x']);       # % neuronal states

    
    """
    #% [default] fixed parameters
    #%--------------------------------------------------------------------------
    E = [1 1/2 1/8]*32;         % extrinsic rates (forward, backward, lateral)
    G = [1 4/5 1/4 1/4]*128;    % intrinsic rates (g1 g2 g3 g4)
    D = [2 16];                 % propogation delays (intrinsic, extrinsic)
    H = [4 32];                 % receptor densities (excitatory, inhibitory)
    T = [8 16];                 % synaptic constants (excitatory, inhibitory)
    R = [2 1]/3;                % parameters of static nonlinearity
    """;
    
    E = np.array([1., 1/2., 1/8.,])*32;         # extrinsic rates (forward, backward, lateral)
    G = np.array([1, 4/5., 1/4., 1/4*128.]);   # % intrinsic rates (g1 g2 g3 g4)    
    D = np.array([2, 16]);  #               % propogation delays (intrinsic, extrinsic)
    H = np.array([4, 32]);  #               % receptor densities (excitatory, inhibitory)
    T = np.array([8, 16]);  #               % synaptic constants (excitatory, inhibitory)
    R = np.array([2, 1])/3.;                #% parameters of static nonlinearity

    
    """
    % [specified] fixed parameters
    %--------------------------------------------------------------------------
    if isfield(M,'pF')
     try, E = M.pF.E; end
     try, G = M.pF.H; end
     try, D = M.pF.D; end
     try, H = M.pF.G; end
     try, T = M.pF.T; end
     try, R = M.pF.R; end
    end   
    """   
    
    
    if 'pF' in M:             # if isfield(M,'pF')
        try: E = M['pF']['E']; 
        except: _
        try: G = M['pF']['H']; 
        except: _
        try: D = M['pF']['D']; 
        except: _
        try: H = M['pF']['G']; 
        except: _
        try: T = M['pF']['T']; 
        except: _  
        try: R = M['pF']['R']; 
        except: _  

    """
    #% test for free parameters on intrinsic connections
    #%--------------------------------------------------------------------------
    try
      G = G.*exp(P.H);
    end
    G     = ones(n,1)*G;
    """
    try:
        G = G*np.exp(P['H']);
    except: _
    G     = np.ones([n,1])*G;
    
    
    """
    % exponential transform to ensure positivity constraints
    %--------------------------------------------------------------------------
    A{1}  = exp(P.A{1})*E(1);
    A{2}  = exp(P.A{2})*E(2);
    A{3}  = exp(P.A{3})*E(3);
    C     = exp(P.C);
    """
    A[0] = np.exp(P['A'][0])*E[0]
    A[1] = np.exp(P['A'][1])*E[1]    
    A[2] = np.exp(P['A'][2])*E[2]
    C = np.exp(P['C'])
    
    """
    % intrinsic connectivity and parameters
    %--------------------------------------------------------------------------
    Te    = T(1)/1000*exp(P.T(:,1));         % excitatory time constants
    Ti    = T(2)/1000*exp(P.T(:,2));         % inhibitory time constants
    He    = H(1)*exp(P.G(:,1));              % excitatory receptor density
    Hi    = H(2)*exp(P.G(:,2));              % inhibitory receptor density
    """
    Te = T[0]/1000*np.exp(P['T'][:,0])
    Ti = T[1]/1000*np.exp(P['T'][:,1])
    He = H[0]*np.exp(P['G'][:,0])
    Hi = H[1]*np.exp(P['G'][:,1])    

    """
    % pre-synaptic inputs: s(V)
    %--------------------------------------------------------------------------
    R     = R.*exp(P.S);
    S     = 1./(1 + exp(-R(1)*(x - R(2)))) - 1./(1 + exp(R(1)*R(2)));
    """
    R = R*np.exp(P['S'])
    S = 1./(1 + np.exp(-R[0]*(x - R[1]))) - 1./(1 + np.exp(R[0]*R[1]));
    
    """
    % input
    %==========================================================================
    if isfield(M,'u')
    
      % endogenous input
      %----------------------------------------------------------------------
      U = u(:)*64;
    
    else
      % exogenous input
      %----------------------------------------------------------------------
      U = C*u(:)*2;
    end
    """
    if 'u' in M:
        #% endogenous input
        #%----------------------------------------------------------------------
        #U = u(:)*64;
        U = u[:]*64.
        
    else:
        #% exogenous input
        #%----------------------------------------------------------------------
        #U = C*u(:)*2;
        U = C*u[:]*2
    
    
    #% State: f(x)
    #%==========================================================================

    """
    % Supragranular layer (inhibitory interneurons): Voltage & depolarizing current
    %--------------------------------------------------------------------------
    f(:,7) = x(:,8);
    f(:,8) = (He.*((A{2} + A{3})*S(:,9) + G(:,3).*S(:,9)) - 2*x(:,8) - x(:,7)./Te)./Te;
    """
    f[:,6] = x[:,7]
    f[:,7] = (He*( (A[1] + A[2]) * S[:,8] + G[:,2]*S[:,8]) - 2*x[:,7] - x[:,6]/Te)/Te
    
    """
    % Granular layer (spiny stellate cells): Voltage & depolarizing current
    %--------------------------------------------------------------------------
    f(:,1) = x(:,4);
    f(:,4) = (He.*((A{1} + A{3})*S(:,9) + G(:,1).*S(:,9) + U) - 2*x(:,4) - x(:,1)./Te)./Te;
    """
    f[:,0] = x[:,3];
    f[:,3] = (He*((A[0] + A[2])*S[:,8] + G[:,0]*S[:,8] + U) - 2*x[:,3] - x[:,0]/Te)/Te;

    """
    % Infra-granular layer (pyramidal cells): depolarizing current
    %--------------------------------------------------------------------------
    f(:,2) = x(:,5);
    f(:,5) = (He.*((A{2} + A{3})*S(:,9) + G(:,2).*S(:,1)) - 2*x(:,5) - x(:,2)./Te)./Te;
    """
    f[:,1] = x[:,4];
    f[:,4] = (He*((A[1] + A[2])*S[:,8] + G[:,1]*S[:,0]) - 2*x[:,4] - x[:,1]/Te)/Te;

    
    """
    % Infra-granular layer (pyramidal cells): hyperpolarizing current
    %--------------------------------------------------------------------------
    f(:,3) = x(:,6);
    f(:,6) = (Hi.*G(:,4).*S(:,7) - 2*x(:,6) - x(:,3)./Ti)./Ti;
    """
    f[:,2] = x[:,5];
    f[:,5] = (Hi*G[:,3]*S[:,6] - 2*x[:,5] - x[:,2]/Ti)/Ti;

    """
    % Infra-granular layer (pyramidal cells): Voltage
    %--------------------------------------------------------------------------
    f(:,9) = x(:,5) - x(:,6);
    f      = spm_vec(f);
    """
    f[:,8] = x[:,4] - x[:,5];
    
    
    #f = spm_vec(f);
    f = spm_vec(f)


    #if nargout < 2; return, end
    if returnJ == False and returnD == False:
        return f
    else: 
    
        """
        % Jacobian
        %==========================================================================
        J  = spm_diff(M.f,x,u,P,M,1);
        """

        # TO DO
        #J  = spm_diff(M.f,x,u,P,M,1);
        J = diff_dfdx(M['f'],x,u,P,M,1)
    
        """
        % delays
        %==========================================================================
        % Delay differential equations can be integrated efficiently (but
        % approximately) by absorbing the delay operator into the Jacobian 
        %
        %    dx(t)/dt     = f(x(t - d))
        %                 = Q(d)f(x(t)) 
        %
        %    J(d)         = Q(d)df/dx
        %--------------------------------------------------------------------------
        De = D(2).*exp(P.D)/1000;
        Di = D(1)/1000;
        De = (1 - speye(n,n)).*De;
        Di = (1 - speye(9,9)).*Di;
        De = kron(ones(9,9),De);
        Di = kron(Di,speye(n,n));
        D  = Di + De;
        """
    
        De = D[1] * np.exp(P['D'])/1000.
        Di = D[0] / 1000.
    
        De = np.array*([1-sps.eye(n,n)])*De
        Di = np.array*([1-sps.eye(9,9)])*Di

        De = np.kron(np.ones([9,9]),De)
        Di = np.kron(Di,sps.eye([n,n]))
                  
        D = Di + De
    
        

        """
        % Implement: dx(t)/dt = f(x(t - d)) = inv(1 + D.*dfdx)*f(x(t))
        %                     = Q*f = Q*J*x(t)
        %--------------------------------------------------------------------------
        Q  = inv(speye(length(J)) + D.*J);
        """;
        
        Q = np.linalg.inv(sps.eye(len(J)) + D*J)

        return f,j,D

tests:

spm_diff port


In [409]:
%load /home/jgriffiths/Code/libraries_of_others/misc/SPM/spm12/spm_dfdx.m


---------------------------------------------------------------------------
ValueError                                Traceback (most recent call last)
<ipython-input-409-80e170245257> in <module>()
----> 1 get_ipython().magic(u'load /home/jgriffiths/Code/libraries_of_others/misc/SPM/spm12/toolbox/DEM/spm_dfdx.m')

/home/jgriffiths/Software/anaconda2/envs/jupyter_release/lib/python2.7/site-packages/IPython/core/interactiveshell.pyc in magic(self, arg_s)
   2156         magic_name, _, magic_arg_s = arg_s.partition(' ')
   2157         magic_name = magic_name.lstrip(prefilter.ESC_MAGIC)
-> 2158         return self.run_line_magic(magic_name, magic_arg_s)
   2159 
   2160     #-------------------------------------------------------------------------

/home/jgriffiths/Software/anaconda2/envs/jupyter_release/lib/python2.7/site-packages/IPython/core/interactiveshell.pyc in run_line_magic(self, magic_name, line)
   2077                 kwargs['local_ns'] = sys._getframe(stack_depth).f_locals
   2078             with self.builtin_trap:
-> 2079                 result = fn(*args,**kwargs)
   2080             return result
   2081 

<decorator-gen-45> in load(self, arg_s)

/home/jgriffiths/Software/anaconda2/envs/jupyter_release/lib/python2.7/site-packages/IPython/core/magic.pyc in <lambda>(f, *a, **k)
    186     # but it's overkill for just that one bit of state.
    187     def magic_deco(arg):
--> 188         call = lambda f, *a, **k: f(*a, **k)
    189 
    190         if callable(arg):

/home/jgriffiths/Software/anaconda2/envs/jupyter_release/lib/python2.7/site-packages/IPython/core/magics/code.pyc in load(self, arg_s)
    348         search_ns = 'n' in opts
    349 
--> 350         contents = self.shell.find_user_code(args, search_ns=search_ns)
    351 
    352         if 's' in opts:

/home/jgriffiths/Software/anaconda2/envs/jupyter_release/lib/python2.7/site-packages/IPython/core/interactiveshell.pyc in find_user_code(self, target, raw, py_only, skip_encoding_cookie, search_ns)
   3183         except Exception:
   3184             raise ValueError(("'%s' was not found in history, as a file, url, "
-> 3185                                 "nor in the user namespace.") % target)
   3186 
   3187         if isinstance(codeobj, string_types):

ValueError: '/home/jgriffiths/Code/libraries_of_others/misc/SPM/spm12/toolbox/DEM/spm_dfdx.m' was not found in history, as a file, url, nor in the user namespace.

In [ ]:
# %load /home/jgriffiths/Code/libraries_of_others/misc/SPM/spm12/spm_diff.m

In [ ]:
def spm_diff(f,x,n,V=None,q=None):
  """
  
  matrix high-order numerical differnetiation
  
  port of spm_diff.m
  
  
  JG June 2017
  
  """    
% matrix high-order numerical differentiation
% FORMAT [dfdx] = spm_diff(f,x,...,n)
% FORMAT [dfdx] = spm_diff(f,x,...,n,V)
% FORMAT [dfdx] = spm_diff(f,x,...,n,'q')
%
% f      - [inline] function f(x{1},...)
% x      - input argument[s]
% n      - arguments to differentiate w.r.t.
%
% V      - cell array of matrices that allow for differentiation w.r.t.
% to a linear transformation of the parameters: i.e., returns
%
% df/dy{i};    x = V{i}y{i};    V = dx(i)/dy(i)
%
% q      - (char) flag to preclude default concatenation of dfdx
%
% dfdx          - df/dx{i}                     ; n =  i
% dfdx{p}...{q} - df/dx{i}dx{j}(q)...dx{k}(p)  ; n = [i j ... k]
%
%
% - a cunning recursive routine
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_diff.m 6110 2014-07-21 09:36:13Z karl $

% create inline object
%--------------------------------------------------------------------------
f     = spm_funcheck(varargin{1});

% parse input arguments
%--------------------------------------------------------------------------
if iscell(varargin{end})
    x = varargin(2:(end - 2));
    n = varargin{end - 1};
    V = varargin{end};
    q = 1;
elseif isnumeric(varargin{end})
    x = varargin(2:(end - 1));
    n = varargin{end};
    V = cell(1,length(x));
    q = 1;
elseif ischar(varargin{end})
    x = varargin(2:(end - 2));
    n = varargin{end - 1};
    V = cell(1,length(x));
    q = 0;
else
    error('improper call')
end

% check transform matrices V = dxdy
%--------------------------------------------------------------------------
for i = 1:length(x)
    try
        V{i};
    catch
        V{i} = [];
    end
    if isempty(V{i}) && any(n == i);
        V{i} = speye(spm_length(x{i}));
    end
end

% initialise
%--------------------------------------------------------------------------
m     = n(end);
xm    = spm_vec(x{m});
dx    = exp(-8);
J     = cell(1,size(V{m},2));

% proceed to derivatives
%==========================================================================
if length(n) == 1
    
    % dfdx
    %----------------------------------------------------------------------
    f0    = f(x{:});
    for i = 1:length(J)
        xi    = x;
        xi{m} = spm_unvec(xm + V{m}(:,i)*dx,x{m});
        J{i}  = spm_dfdx(f(xi{:}),f0,dx);
    end

    
    % return numeric array for first-order derivatives
    %======================================================================
    
    % vectorise f
    %----------------------------------------------------------------------
    f  = spm_vec(f0);
    
    % if there are no arguments to differentiate w.r.t. ...
    %----------------------------------------------------------------------
    if isempty(xm)
        J = sparse(length(f),0);
        
        % or there are no arguments to differentiate
        %----------------------------------------------------------------------
    elseif isempty(f)
        J = sparse(0,length(xm));
    end
    
    % or differentiation of a vector
    %----------------------------------------------------------------------
    if isvector(f0) && isnumeric(f0) && q
        
        % concatenate into a matrix
        %------------------------------------------------------------------
        if size(f0,2) == 1
            J = spm_cat(J);
        else
            J = spm_cat(J')';
        end
    end
    
    % assign output argument and return
    %----------------------------------------------------------------------
    varargout{1} = J;
    varargout{2} = f0;
    
else
    
    % dfdxdxdx....
    %----------------------------------------------------------------------
    f0        = cell(1,length(n));
    [f0{:}]   = spm_diff(f,x{:},n(1:end - 1),V);
    
    for i = 1:length(J)
        xi    = x;
        xmi   = xm + V{m}(:,i)*dx;
        xi{m} = spm_unvec(xmi,x{m});
        fi    = spm_diff(f,xi{:},n(1:end - 1),V);
        J{i}  = spm_dfdx(fi,f0{1},dx);
    end
    varargout = [{J} f0];
end

In [ ]:
def spm_dfdx(f,f0,dx):

    """
    port of spm_dfdx
    ...which is defined at the bottom of spm_diff.m
    
    """
    
    """
    function dfdx = spm_dfdx(f,f0,dx)
    % cell subtraction
    %--------------------------------------------------------------------------
    """
    
    """
    if iscell(f)
      dfdx  = f;
      for i = 1:length(f(:))
        dfdx{i} = spm_dfdx(f{i},f0{i},dx);
      end  
    elseif isstruct(f)
      dfdx  = (spm_vec(f) - spm_vec(f0))/dx;
    else
      dfdx  = (f - f0)/dx;
    end
    """;
    
    # Check if input is a dict(cell)
    # ...if it is, loop through elements and 
    # re-call this function
    # ...if it isn't, 
    # 

    dfdx = []
    if type(f) == list:
        dfdx  = f;
        for i = in f[:]: 
            dfdx{i} = spm_dfdx(f{i},f0{i},dx);
    elif isstruct(f)
      dfdx  = (spm_vec(f) - spm_vec(f0))/dx;
    else
      dfdx  = (f - f0)/dx;
    end
    """;

In [ ]:
# %load /home/jgriffiths/Code/libraries_of_others/misc/SPM/spm12/spm_diff.m
function [varargout] = spm_diff(varargin)
% matrix high-order numerical differentiation
% FORMAT [dfdx] = spm_diff(f,x,...,n)
% FORMAT [dfdx] = spm_diff(f,x,...,n,V)
% FORMAT [dfdx] = spm_diff(f,x,...,n,'q')
%
% f      - [inline] function f(x{1},...)
% x      - input argument[s]
% n      - arguments to differentiate w.r.t.
%
% V      - cell array of matrices that allow for differentiation w.r.t.
% to a linear transformation of the parameters: i.e., returns
%
% df/dy{i};    x = V{i}y{i};    V = dx(i)/dy(i)
%
% q      - (char) flag to preclude default concatenation of dfdx
%
% dfdx          - df/dx{i}                     ; n =  i
% dfdx{p}...{q} - df/dx{i}dx{j}(q)...dx{k}(p)  ; n = [i j ... k]
%
%
% - a cunning recursive routine
%__________________________________________________________________________
% Copyright (C) 2008 Wellcome Trust Centre for Neuroimaging

% Karl Friston
% $Id: spm_diff.m 6110 2014-07-21 09:36:13Z karl $

% create inline object
%--------------------------------------------------------------------------
f     = spm_funcheck(varargin{1});

% parse input arguments
%--------------------------------------------------------------------------
if iscell(varargin{end})
    x = varargin(2:(end - 2));
    n = varargin{end - 1};
    V = varargin{end};
    q = 1;
elseif isnumeric(varargin{end})
    x = varargin(2:(end - 1));
    n = varargin{end};
    V = cell(1,length(x));
    q = 1;
elseif ischar(varargin{end})
    x = varargin(2:(end - 2));
    n = varargin{end - 1};
    V = cell(1,length(x));
    q = 0;
else
    error('improper call')
end

% check transform matrices V = dxdy
%--------------------------------------------------------------------------
for i = 1:length(x)
    try
        V{i};
    catch
        V{i} = [];
    end
    if isempty(V{i}) && any(n == i);
        V{i} = speye(spm_length(x{i}));
    end
end

% initialise
%--------------------------------------------------------------------------
m     = n(end);
xm    = spm_vec(x{m});
dx    = exp(-8);
J     = cell(1,size(V{m},2));

% proceed to derivatives
%==========================================================================
if length(n) == 1
    
    % dfdx
    %----------------------------------------------------------------------
    f0    = f(x{:});
    for i = 1:length(J)
        xi    = x;
        xi{m} = spm_unvec(xm + V{m}(:,i)*dx,x{m});
        J{i}  = spm_dfdx(f(xi{:}),f0,dx);
    end

    
    % return numeric array for first-order derivatives
    %======================================================================
    
    % vectorise f
    %----------------------------------------------------------------------
    f  = spm_vec(f0);
    
    % if there are no arguments to differentiate w.r.t. ...
    %----------------------------------------------------------------------
    if isempty(xm)
        J = sparse(length(f),0);
        
        % or there are no arguments to differentiate
        %----------------------------------------------------------------------
    elseif isempty(f)
        J = sparse(0,length(xm));
    end
    
    % or differentiation of a vector
    %----------------------------------------------------------------------
    if isvector(f0) && isnumeric(f0) && q
        
        % concatenate into a matrix
        %------------------------------------------------------------------
        if size(f0,2) == 1
            J = spm_cat(J);
        else
            J = spm_cat(J')';
        end
    end
    
    % assign output argument and return
    %----------------------------------------------------------------------
    varargout{1} = J;
    varargout{2} = f0;
    
else
    
    % dfdxdxdx....
    %----------------------------------------------------------------------
    f0        = cell(1,length(n));
    [f0{:}]   = spm_diff(f,x{:},n(1:end - 1),V);
    
    for i = 1:length(J)
        xi    = x;
        xmi   = xm + V{m}(:,i)*dx;
        xi{m} = spm_unvec(xmi,x{m});
        fi    = spm_diff(f,xi{:},n(1:end - 1),V);
        J{i}  = spm_dfdx(fi,f0{1},dx);
    end
    varargout = [{J} f0];
end


function dfdx = spm_dfdx(f,f0,dx)
% cell subtraction
%--------------------------------------------------------------------------
if iscell(f)
    dfdx  = f;
    for i = 1:length(f(:))
        dfdx{i} = spm_dfdx(f{i},f0{i},dx);
    end
elseif isstruct(f)
    dfdx  = (spm_vec(f) - spm_vec(f0))/dx;
else
    dfdx  = (f - f0)/dx;
end

spm_vec and spm_unvec ports

code:


In [381]:
def spm_vec(X):
    # Vectorise a numeric, cell or structure array - a compiled routine
    # FORMAT [vX] = spm_vec(X)
    #X  - numeric, cell or stucture array[s]
    #vX - vec(X)
    #X = [np.eye(2),3]
    
    #  Usage:
    #vX = vec(X)
    
    for v_it,v in enumerate(X):
        if v_it == 0:
            try:
                nr,nc = v.shape
                vX = np.reshape(v,[nr*nc,1])
            except:
                vX = np.array([v][:,np.newaxis])
        else:
            try:
                nr,nc = v.shape
                vX = np.concatenate([vX,np.reshape(v,[nr*nc,1])])
            except:
                vX = np.concatenate([vX,np.array([v])[:,np.newaxis]])            
                
    return np.squeeze(vX)[:,np.newaxis]



def unvec(Xflat,X):

    #nr,nc = X.shape
    #vX = np.reshape(X,[nr,nc])    
    
    cnt = 0
    vXlist = []
    for v_it,v in enumerate(X):
        try:
            nr,nc = v.shape
            rng = np.arange(cnt,nr*nc)
            vX = np.reshape(Xflat[rng],[nr,nc])
            vXlist.append(vX)
            cnt+=1
        except: 
            vXlist.append(v)
    
    return vXlist

tests:

matlab spm_vec


In [385]:
%%matlab
X = {eye(2) 3};
vX = spm_vec(X);
disp(['shape: ' num2str(size(vX))]);
vX


shape: 5  1

vX =

     1
     0
     0
     1
     3

python spm_vec


In [386]:
X = [np.eye(2),3]
vX = spm_vec(X)
print 'shape: %s,%s' %vX.shape
print 'value: \n%s' %vX


shape: 5,1
value: 
[[ 1.]
 [ 0.]
 [ 0.]
 [ 1.]
 [ 3.]]

matlab spm_unvec


In [387]:
%%matlab
X = {eye(2) 3};
vX = spm_unvec(spm_vec(X),X);
vX{:}


ans =

     1     0
     0     1


ans =

     3

python spm_unvec


In [388]:
X = [np.eye(2),3]
vX = spm_unvec(vec(X),X)
#print 'shape: %s,%s' %vX.shape
#print 'value: \n%s' %vX
vX


Out[388]:
[array([[ 1.,  0.],
        [ 0.,  1.]]), 3]

spm_cat port

code

tests

spm_length port

code

tests



In [ ]:
_,_,D = fx_erp()

Test matlab vs. python functions


In [2]:
%load_ext pymatbridge


Starting MATLAB on ZMQ socket ipc:///tmp/pymatbridge-587dd83b-4fb5-4066-a12b-7f566714207e
Send 'exit' command to kill the server
...................MATLAB started and connected!

In [3]:
%%matlab

addpath(genpath('/home/jgriffiths/Code/libraries_of_others/misc/SPM'));

First half of demo function


In [27]:
%%matlab
rng('default')

Nc    = 2;                                        % number of channels
Ns    = 2;                                        % number of sources

options.spatial  = 'LFP';
options.model    = 'ERP';
options.analysis = 'ERP';
M.dipfit.model   = options.model;
M.dipfit.type    = options.spatial;
M.dipfit.Nc      = Nc;
M.dipfit.Ns      = Ns;


% sspecify connectivity – reciprocal connections with condition specific
% changes in intrinsic and extrinsic connectivity
%--------------------------------------------------------------------------
A{1}    = [0 0; 1 0];
A{2}    = [0 1; 0 0];
A{3}    = [0 0; 0 0];
B{1}    = [1 0; 1 0];
C       = [1; 0];

[pE,pC] = spm_dcm_neural_priors(A,B,C,options.model);
[gE,gC] = spm_L_priors(M.dipfit);
[x,f]   = spm_dcm_x_neural(pE,options.model);

% hyperpriors (assuming a high signal to noise)
%--------------------------------------------------------------------------
hE      = 6;
hC      = 1/128;

% create model
%--------------------------------------------------------------------------
M.IS   = 'spm_gen_erp';
M.G    = 'spm_lx_erp';
M.f    = f;
M.x    = x;
M.pE   = pE;
M.pC   = pC;
M.gE   = gE;
M.gC   = gC;
M.hE   = hE;
M.hC   = hC;
M.m    = length(B);
M.n    = length(spm_vec(M.x));
M.l    = Nc;
M.ns   = 64;

% create input structure
%--------------------------------------------------------------------------
dt     = 4/1000;
pst    = (1:M.ns)*dt;
M.ons  = 64;
M.dur  = 16;
U.dt   = dt;
U.X    = [0; 1];

% specified true connectivity (P) and spatial parameters (G) – with
% condition specific effects on the intrinsic connectivity of the first
% source and its forward extrinsic connection
%--------------------------------------------------------------------------
P      = pE;
G      = gE;
P.B{1} = [-1/4 0; 1/2 0];


% generate neuronal response and data
%--------------------------------------------------------------------------
x     = spm_gen_erp(P,M,U);                 % neuronal response
L     = spm_lx_erp(G,M.dipfit);             % lead field
V     = spm_sqrtm(spm_Q(1/2,M.ns));         % square root of noise covariance
for i = 1:length(x)
    n    = exp(-hE/2)*V*randn(M.ns,Nc);     % noise
    s{i} = x{i}*L';                         % signal
    y{i} = s{i} + n;                        % data (signal plus noise)
end


% data structure specification
%--------------------------------------------------------------------------
Y.y   = y;
Y.Q   = {spm_Q(1/2,M.ns,1)};
Y.dt  = dt;
Y.pst = pst;

Push variables to python


In [31]:
%%matlab -o Y 
disp('');

In [34]:
len(Y['y'])


Out[34]:
2

Run python version

Compare matlab and python values


In [ ]:


In [ ]:


In [ ]:


In [ ]:
a

In [ ]:


In [25]:
%%matlab

% display
%--------------------------------------------------------------------------
spm_figure('Getwin','Figure 1');
subplot(2,1,1)
plot(pst,x{1},'r',pst,x{2},'b')
xlabel('time');ylabel('amplitude');
title('Hidden neuronal states','FontSize',16)



subplot(2,1,2)
plot(pst,s{1},':r',pst,s{2},':b',pst,y{1},'r',pst,y{2},'b')
xlabel('time');ylabel('amplitude');
title('Observed response','FontSize',16)


---------------------------------------------------------------------------
MatlabInterperterError                    Traceback (most recent call last)
<ipython-input-25-7f3b42da95c8> in <module>()
----> 1 get_ipython().run_cell_magic(u'matlab', u'', u"\n% display\n%--------------------------------------------------------------------------\nspm_figure('Getwin','Figure 1');\nsubplot(2,1,1)\nplot(pst,x{1},'r',pst,x{2},'b')\nxlabel('time');ylabel('amplitude');\ntitle('Hidden neuronal states','FontSize',16)\n\nsubplot(2,1,2)\nplot(pst,s{1},':r',pst,s{2},':b',pst,y{1},'r',pst,y{2},'b')\nxlabel('time');ylabel('amplitude');\ntitle('Observed response','FontSize',16)")

/home/jgriffiths/Software/anaconda2/envs/jupyter_release/lib/python2.7/site-packages/IPython/core/interactiveshell.pyc in run_cell_magic(self, magic_name, line, cell)
   2113             magic_arg_s = self.var_expand(line, stack_depth)
   2114             with self.builtin_trap:
-> 2115                 result = fn(magic_arg_s, cell)
   2116             return result
   2117 

<decorator-gen-122> in matlab(self, line, cell, local_ns)

/home/jgriffiths/Software/anaconda2/envs/jupyter_release/lib/python2.7/site-packages/IPython/core/magic.pyc in <lambda>(f, *a, **k)
    186     # but it's overkill for just that one bit of state.
    187     def magic_deco(arg):
--> 188         call = lambda f, *a, **k: f(*a, **k)
    189 
    190         if callable(arg):

/home/jgriffiths/Software/anaconda2/envs/jupyter_release/lib/python2.7/site-packages/pymatbridge/matlab_magic.pyc in matlab(self, line, cell, local_ns)
    161 
    162         try:
--> 163             result_dict = self.eval(code)
    164         except MatlabInterperterError:
    165             raise

/home/jgriffiths/Software/anaconda2/envs/jupyter_release/lib/python2.7/site-packages/pymatbridge/matlab_magic.pyc in eval(self, line)
     88 
     89         if not run_dict['success']:
---> 90             raise MatlabInterperterError(line, run_dict['content']['stdout'])
     91 
     92         # This is the matlab stdout:

MatlabInterperterError: Failed to parse and evaluate line u"\n% display\n%--------------------------------------------------------------------------\nspm_figure('Getwin','Figure 1');\nsubplot(2,1,1)\nplot(pst,x{1},'r',pst,x{2},'b')\nxlabel('time');ylabel('amplitude');\ntitle('Hidden neuronal states','FontSize',16)\n\nsubplot(2,1,2)\nplot(pst,s{1},':r',pst,s{2},':b',pst,y{1},'r',pst,y{2},'b')\nxlabel('time');ylabel('amplitude');\ntitle('Observed response','FontSize',16)".
 Matlab error message: u'Java exception occurred: \ncom.mathworks.hg.util.OutputHelperProcessingException: Problem while processing in an OutputHelper. java.lang.OutOfMemoryError: Java heap space\n\tat com.mathworks.hg.util.HGGetframeOutputHelper.generateOutput(HGGetframeOutputHelper.java:216)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.grabOnEDT(HGGetframeOutputHelper.java:363)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.access$300(HGGetframeOutputHelper.java:344)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber$1.run(HGGetframeOutputHelper.java:380)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$5$1.run(AWTUtilities.java:591)\n\tat com.mathworks.mvm.context.ThreadContext$1.call(ThreadContext.java:76)\n\tat com.mathworks.mvm.context.ThreadContext.callWithContext(ThreadContext.java:105)\n\tat com.mathworks.mvm.context.ThreadContext.runWithContext(ThreadContext.java:73)\n\tat com.mathworks.mvm.context.MvmContext.runWithContext(MvmContext.java:107)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$5.runWithOutput(AWTUtilities.java:588)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$2.watchedRun(AWTUtilities.java:475)\n\tat com.mathworks.jmi.AWTUtilities$WatchedRunnable.run(AWTUtilities.java:436)\n\tat com.mathworks.jmi.AWTUtilities$Invoker.invoke(AWTUtilities.java:490)\n\tat com.mathworks.jmi.AWTUtilities.invokeAndWait(AWTUtilities.java:304)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.grab(HGGetframeOutputHelper.java:377)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.getBufferedImage(HGGetframeOutputHelper.java:356)\n\tat com.mathworks.hg.peer.FigureClientProxyPanel.setPaintDisabled(FigureClientProxyPanel.java:66)\n\tat com.mathworks.hg.peer.PaintDisabled.setPaintDisabled(PaintDisabled.java:60)\n\tat com.mathworks.hg.peer.HeavyweightLightweightContainerFactory$FigurePanelContainerLight.disablePaint(HeavyweightLightweightContainerFactory.java:326)\n\tat com.mathworks.hg.peer.HeavyweightLightweightContainerFactory$FigurePanelContainerLight.doSetPaintDisabled(HeavyweightLightweightContainerFactory.java:363)\n\tat com.mathworks.hg.peer.HeavyweightLightweightContainerFactory$FigurePanelContainerLight.setPaintDisabled(HeavyweightLightweightContainerFactory.java:387)\n\tat sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)\n\tat sun.reflect.NativeMethodAccessorImpl.invoke(Unknown Source)\n\tat sun.reflect.DelegatingMethodAccessorImpl.invoke(Unknown Source)\n\tat java.lang.reflect.Method.invoke(Unknown Source)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$3$1.call(AWTUtilities.java:525)\n\tat com.mathworks.mvm.context.ThreadContext.callWithContext(ThreadContext.java:105)\n\tat com.mathworks.mvm.context.MvmContext.callWithContext(MvmContext.java:113)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$3.runWithOutput(AWTUtilities.java:522)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$2.watchedRun(AWTUtilities.java:475)\n\tat com.mathworks.jmi.AWTUtilities$WatchedRunnable.run(AWTUtilities.java:436)\n\tat java.awt.event.InvocationEvent.dispatch(Unknown Source)\n\tat java.awt.EventQueue.dispatchEventImpl(Unknown Source)\n\tat java.awt.EventQueue.access$200(Unknown Source)\n\tat java.awt.EventQueue$3.run(Unknown Source)\n\tat java.awt.EventQueue$3.run(Unknown Source)\n\tat java.security.AccessController.doPrivileged(Native Method)\n\tat java.security.ProtectionDomain$1.doIntersectionPrivilege(Unknown Source)\n\tat java.awt.EventQueue.dispatchEvent(Unknown Source)\n\tat java.awt.EventDispatchThread.pumpOneEventForFilters(Unknown Source)\n\tat java.awt.EventDispatchThread.pumpEventsForFilter(Unknown Source)\n\tat java.awt.EventDispatchThread.pumpEventsForHierarchy(Unknown Source)\n\tat java.awt.EventDispatchThread.pumpEvents(Unknown Source)\n\tat java.awt.EventDispatchThread.pumpEvents(Unknown Source)\n\tat java.awt.EventDispatchThread.run(Unknown Source)\nCaused by: com.mathworks.hg.util.HGGetframeOutputHelper$RasterSizeException: java.lang.OutOfMemoryError: Java heap space\n\t... 45 more\nCaused by: java.lang.OutOfMemoryError: Java heap space\n\tat java.awt.image.DataBufferInt.<init>(Unknown Source)\n\tat java.awt.image.Raster.createPackedRaster(Unknown Source)\n\tat java.awt.image.DirectColorModel.createCompatibleWritableRaster(Unknown Source)\n\tat java.awt.image.BufferedImage.<init>(Unknown Source)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$1.run(HGGetframeOutputHelper.java:176)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$5$1.run(AWTUtilities.java:591)\n\tat com.mathworks.mvm.context.ThreadContext$1.call(ThreadContext.java:76)\n\tat com.mathworks.mvm.context.ThreadContext.callWithContext(ThreadContext.java:105)\n\tat com.mathworks.mvm.context.ThreadContext.runWithContext(ThreadContext.java:73)\n\tat com.mathworks.mvm.context.MvmContext.runWithContext(MvmContext.java:107)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$5.runWithOutput(AWTUtilities.java:588)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$2.watchedRun(AWTUtilities.java:475)\n\tat com.mathworks.jmi.AWTUtilities$WatchedRunnable.run(AWTUtilities.java:436)\n\tat com.mathworks.jmi.AWTUtilities$Invoker.invoke(AWTUtilities.java:490)\n\tat com.mathworks.jmi.AWTUtilities.invokeAndWait(AWTUtilities.java:304)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper.generateOutput(HGGetframeOutputHelper.java:206)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.grabOnEDT(HGGetframeOutputHelper.java:363)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.access$300(HGGetframeOutputHelper.java:344)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber$1.run(HGGetframeOutputHelper.java:380)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$5$1.run(AWTUtilities.java:591)\n\tat com.mathworks.mvm.context.ThreadContext$1.call(ThreadContext.java:76)\n\tat com.mathworks.mvm.context.ThreadContext.callWithContext(ThreadContext.java:105)\n\tat com.mathworks.mvm.context.ThreadContext.runWithContext(ThreadContext.java:73)\n\tat com.mathworks.mvm.context.MvmContext.runWithContext(MvmContext.java:107)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$5.runWithOutput(AWTUtilities.java:588)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$2.watchedRun(AWTUtilities.java:475)\n\tat com.mathworks.jmi.AWTUtilities$WatchedRunnable.run(AWTUtilities.java:436)\n\tat com.mathworks.jmi.AWTUtilities$Invoker.invoke(AWTUtilities.java:490)\n\tat com.mathworks.jmi.AWTUtilities.invokeAndWait(AWTUtilities.java:304)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.grab(HGGetframeOutputHelper.java:377)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.getBufferedImage(HGGetframeOutputHelper.java:356)\n\tat com.mathworks.hg.peer.FigureClientProxyPanel.setPaintDisabled(FigureClientProxyPanel.java:66)'

In [19]:
% Invert model under increasing shrinkage priors on the condition specific
% change in the intrinsic (B) parameter of the first source. This range
% specified by alpha. Because the true value is non-zero, we expect  free-
% energy to decrease when the prior covariance falls to 0 and this
% parameter is effectively eliminated. The ensuing estimates are obtained
% while optimising the parameters of the spatial model and the  (neuronal)
% extrinsic parameter..
% =========================================================================


---------------------------------------------------------------------------
MatlabInterperterError                    Traceback (most recent call last)
<ipython-input-19-029d730cffbc> in <module>()
----> 1 get_ipython().run_cell_magic(u'matlab', u'', u"\n% display\n%--------------------------------------------------------------------------\nspm_figure('Getwin','Figure 1');\nsubplot(2,1,1)\nplot(pst,x{1},'r',pst,x{2},'b')\nxlabel('time');ylabel('amplitude');\ntitle('Hidden neuronal states','FontSize',16)\n\nsubplot(2,1,2)\nplot(pst,s{1},':r',pst,s{2},':b',pst,y{1},'r',pst,y{2},'b')\nxlabel('time');ylabel('amplitude');\ntitle('Observed response','FontSize',16)\n\n\n\n% Invert model under increasing shrinkage priors on the condition specific\n% change in the intrinsic (B) parameter of the first source. This range\n% specified by alpha. Because the true value is non-zero, we expect  free-\n% energy to decrease when the prior covariance falls to 0 and this\n% parameter is effectively eliminated. The ensuing estimates are obtained\n% while optimising the parameters of the spatial model and the  (neuronal)\n% extrinsic parameter..\n% =========================================================================")

/home/jgriffiths/Software/anaconda2/envs/jupyter_release/lib/python2.7/site-packages/IPython/core/interactiveshell.pyc in run_cell_magic(self, magic_name, line, cell)
   2113             magic_arg_s = self.var_expand(line, stack_depth)
   2114             with self.builtin_trap:
-> 2115                 result = fn(magic_arg_s, cell)
   2116             return result
   2117 

<decorator-gen-122> in matlab(self, line, cell, local_ns)

/home/jgriffiths/Software/anaconda2/envs/jupyter_release/lib/python2.7/site-packages/IPython/core/magic.pyc in <lambda>(f, *a, **k)
    186     # but it's overkill for just that one bit of state.
    187     def magic_deco(arg):
--> 188         call = lambda f, *a, **k: f(*a, **k)
    189 
    190         if callable(arg):

/home/jgriffiths/Software/anaconda2/envs/jupyter_release/lib/python2.7/site-packages/pymatbridge/matlab_magic.pyc in matlab(self, line, cell, local_ns)
    161 
    162         try:
--> 163             result_dict = self.eval(code)
    164         except MatlabInterperterError:
    165             raise

/home/jgriffiths/Software/anaconda2/envs/jupyter_release/lib/python2.7/site-packages/pymatbridge/matlab_magic.pyc in eval(self, line)
     88 
     89         if not run_dict['success']:
---> 90             raise MatlabInterperterError(line, run_dict['content']['stdout'])
     91 
     92         # This is the matlab stdout:

MatlabInterperterError: Failed to parse and evaluate line u"\n% display\n%--------------------------------------------------------------------------\nspm_figure('Getwin','Figure 1');\nsubplot(2,1,1)\nplot(pst,x{1},'r',pst,x{2},'b')\nxlabel('time');ylabel('amplitude');\ntitle('Hidden neuronal states','FontSize',16)\n\nsubplot(2,1,2)\nplot(pst,s{1},':r',pst,s{2},':b',pst,y{1},'r',pst,y{2},'b')\nxlabel('time');ylabel('amplitude');\ntitle('Observed response','FontSize',16)\n\n\n\n% Invert model under increasing shrinkage priors on the condition specific\n% change in the intrinsic (B) parameter of the first source. This range\n% specified by alpha. Because the true value is non-zero, we expect  free-\n% energy to decrease when the prior covariance falls to 0 and this\n% parameter is effectively eliminated. The ensuing estimates are obtained\n% while optimising the parameters of the spatial model and the  (neuronal)\n% extrinsic parameter..\n% =========================================================================".
 Matlab error message: u'Java exception occurred: \ncom.mathworks.hg.util.OutputHelperProcessingException: Problem while processing in an OutputHelper. java.lang.OutOfMemoryError: Java heap space\n\tat com.mathworks.hg.util.HGGetframeOutputHelper.generateOutput(HGGetframeOutputHelper.java:216)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.grabOnEDT(HGGetframeOutputHelper.java:363)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.access$300(HGGetframeOutputHelper.java:344)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber$1.run(HGGetframeOutputHelper.java:380)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$5$1.run(AWTUtilities.java:591)\n\tat com.mathworks.mvm.context.ThreadContext$1.call(ThreadContext.java:76)\n\tat com.mathworks.mvm.context.ThreadContext.callWithContext(ThreadContext.java:105)\n\tat com.mathworks.mvm.context.ThreadContext.runWithContext(ThreadContext.java:73)\n\tat com.mathworks.mvm.context.MvmContext.runWithContext(MvmContext.java:107)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$5.runWithOutput(AWTUtilities.java:588)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$2.watchedRun(AWTUtilities.java:475)\n\tat com.mathworks.jmi.AWTUtilities$WatchedRunnable.run(AWTUtilities.java:436)\n\tat com.mathworks.jmi.AWTUtilities$Invoker.invoke(AWTUtilities.java:490)\n\tat com.mathworks.jmi.AWTUtilities.invokeAndWait(AWTUtilities.java:304)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.grab(HGGetframeOutputHelper.java:377)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.getBufferedImage(HGGetframeOutputHelper.java:356)\n\tat com.mathworks.hg.peer.FigureClientProxyPanel.setPaintDisabled(FigureClientProxyPanel.java:66)\n\tat com.mathworks.hg.peer.PaintDisabled.setPaintDisabled(PaintDisabled.java:60)\n\tat com.mathworks.hg.peer.HeavyweightLightweightContainerFactory$FigurePanelContainerLight.disablePaint(HeavyweightLightweightContainerFactory.java:326)\n\tat com.mathworks.hg.peer.HeavyweightLightweightContainerFactory$FigurePanelContainerLight.doSetPaintDisabled(HeavyweightLightweightContainerFactory.java:363)\n\tat com.mathworks.hg.peer.HeavyweightLightweightContainerFactory$FigurePanelContainerLight.setPaintDisabled(HeavyweightLightweightContainerFactory.java:387)\n\tat sun.reflect.NativeMethodAccessorImpl.invoke0(Native Method)\n\tat sun.reflect.NativeMethodAccessorImpl.invoke(Unknown Source)\n\tat sun.reflect.DelegatingMethodAccessorImpl.invoke(Unknown Source)\n\tat java.lang.reflect.Method.invoke(Unknown Source)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$3$1.call(AWTUtilities.java:525)\n\tat com.mathworks.mvm.context.ThreadContext.callWithContext(ThreadContext.java:105)\n\tat com.mathworks.mvm.context.MvmContext.callWithContext(MvmContext.java:113)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$3.runWithOutput(AWTUtilities.java:522)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$2.watchedRun(AWTUtilities.java:475)\n\tat com.mathworks.jmi.AWTUtilities$WatchedRunnable.run(AWTUtilities.java:436)\n\tat java.awt.event.InvocationEvent.dispatch(Unknown Source)\n\tat java.awt.EventQueue.dispatchEventImpl(Unknown Source)\n\tat java.awt.EventQueue.access$200(Unknown Source)\n\tat java.awt.EventQueue$3.run(Unknown Source)\n\tat java.awt.EventQueue$3.run(Unknown Source)\n\tat java.security.AccessController.doPrivileged(Native Method)\n\tat java.security.ProtectionDomain$1.doIntersectionPrivilege(Unknown Source)\n\tat java.awt.EventQueue.dispatchEvent(Unknown Source)\n\tat java.awt.EventDispatchThread.pumpOneEventForFilters(Unknown Source)\n\tat java.awt.EventDispatchThread.pumpEventsForFilter(Unknown Source)\n\tat java.awt.EventDispatchThread.pumpEventsForHierarchy(Unknown Source)\n\tat java.awt.EventDispatchThread.pumpEvents(Unknown Source)\n\tat java.awt.EventDispatchThread.pumpEvents(Unknown Source)\n\tat java.awt.EventDispatchThread.run(Unknown Source)\nCaused by: com.mathworks.hg.util.HGGetframeOutputHelper$RasterSizeException: java.lang.OutOfMemoryError: Java heap space\n\t... 45 more\nCaused by: java.lang.OutOfMemoryError: Java heap space\n\tat java.awt.image.DataBufferInt.<init>(Unknown Source)\n\tat java.awt.image.Raster.createPackedRaster(Unknown Source)\n\tat java.awt.image.DirectColorModel.createCompatibleWritableRaster(Unknown Source)\n\tat java.awt.image.BufferedImage.<init>(Unknown Source)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$1.run(HGGetframeOutputHelper.java:176)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$5$1.run(AWTUtilities.java:591)\n\tat com.mathworks.mvm.context.ThreadContext$1.call(ThreadContext.java:76)\n\tat com.mathworks.mvm.context.ThreadContext.callWithContext(ThreadContext.java:105)\n\tat com.mathworks.mvm.context.ThreadContext.runWithContext(ThreadContext.java:73)\n\tat com.mathworks.mvm.context.MvmContext.runWithContext(MvmContext.java:107)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$5.runWithOutput(AWTUtilities.java:588)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$2.watchedRun(AWTUtilities.java:475)\n\tat com.mathworks.jmi.AWTUtilities$WatchedRunnable.run(AWTUtilities.java:436)\n\tat com.mathworks.jmi.AWTUtilities$Invoker.invoke(AWTUtilities.java:490)\n\tat com.mathworks.jmi.AWTUtilities.invokeAndWait(AWTUtilities.java:304)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper.generateOutput(HGGetframeOutputHelper.java:206)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.grabOnEDT(HGGetframeOutputHelper.java:363)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.access$300(HGGetframeOutputHelper.java:344)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber$1.run(HGGetframeOutputHelper.java:380)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$5$1.run(AWTUtilities.java:591)\n\tat com.mathworks.mvm.context.ThreadContext$1.call(ThreadContext.java:76)\n\tat com.mathworks.mvm.context.ThreadContext.callWithContext(ThreadContext.java:105)\n\tat com.mathworks.mvm.context.ThreadContext.runWithContext(ThreadContext.java:73)\n\tat com.mathworks.mvm.context.MvmContext.runWithContext(MvmContext.java:107)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$5.runWithOutput(AWTUtilities.java:588)\n\tat com.mathworks.jmi.AWTUtilities$Invoker$2.watchedRun(AWTUtilities.java:475)\n\tat com.mathworks.jmi.AWTUtilities$WatchedRunnable.run(AWTUtilities.java:436)\n\tat com.mathworks.jmi.AWTUtilities$Invoker.invoke(AWTUtilities.java:490)\n\tat com.mathworks.jmi.AWTUtilities.invokeAndWait(AWTUtilities.java:304)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.grab(HGGetframeOutputHelper.java:377)\n\tat com.mathworks.hg.util.HGGetframeOutputHelper$BufferedImageGrabber.getBufferedImage(HGGetframeOutputHelper.java:356)\n\tat com.mathworks.hg.peer.FigureClientProxyPanel.setPaintDisabled(FigureClientProxyPanel.java:66)'

In [ ]:
%%matlab

% fix all (neuronal) parameters (except those of interest)
% -------------------------------------------------------------------------
M.pC           = spm_unvec(spm_vec(pC)*0,pC);
M.pC.B{1}(1,1) = 1/8;
M.pC.B{1}(2,1) = 1/8;

% vary the prior covariance
% -------------------------------------------------------------------------
alpha = exp(-8:2);
for i = 1:length(alpha)   

    % reset shrinkage prior
    % ---------------------------------------------------------------------
    M.pC.B{1}(1,1) = alpha(i);

    % full inversion
    % ---------------------------------------------------------------------
    [Ep,Eg,Cp,Cg,S,F] = spm_nlsi_N(M,U,Y);
    Ep_all{i} = spm_vec(Ep);
    Cp_all{i} = diag(Cp);
    F_all(i)  = F;
    
end

In [ ]:
%%matlab

% Post-hoc reduction for the full (complex) model: this evaluates the free
% energy and posterior distribution over the parameters, given just the
% full prior and posterior – for any required reduced prior.
% -------------------------------------------------------------------------
for i = 1:length(alpha)
    
    rC           = M.pC;
    rC.B{1}(1,1) = alpha(i);

    [F,rEp,rCp]  = spm_log_evidence_reduce(Ep,Cp,M.pE,M.pC,M.pE,rC);
    rEp_all{i}   = spm_vec(rEp);
    rCp_all{i}   = diag(rCp);
    rF_all(i)    = F;
    
end

In [12]:
% compare full inversion and model reduction
% -------------------------------------------------------------------------
spm_figure('Getwin','Figure 2');

subplot(2,1,1)
semilogx(alpha,F_all - F_all(end),alpha,rF_all - rF_all(end))
xlabel('prior covariance');ylabel('relative log evidence');
legend('full','posthoc');
title('Full and reduced log evidence','FontSize',16)


% compare parameter estimates
% -------------------------------------------------------------------------
j     = find(spm_vec(M.pC));
pP    = spm_vec(P);
for i = 1:length(alpha)   

    fQp(i,:) = Ep_all{i}(j);
    fCp(i,:) = Cp_all{i}(j);
    rQp(i,:) = rEp_all{i}(j);
    vCp(i,:) = rCp_all{i}(j);
    
    Pp(i,:)  = pP(j);
end

subplot(2,1,2)
spm_plot_ci(rQp',vCp'),  hold on
plot(fQp,'--'),          hold on
plot(Pp,'-.'),           hold off

xlabel('prior log-covariance');ylabel('difference in log evidence');
title('MAP estimates (full -- reduced  - true -.)','FontSize',16)
set(gca,'XTickLabel',log(alpha))


return


% free energy landscape:
% Here, we evaluate the free energy, which is a functional of the data and
% conditional or posterior expectations (noting that posterior precisions
% are functions of the expectations). The free energy can be computed in a
% simple way by inverting the model using a single iteration. In what
% follows, we evaluate the free energy over a range of the two (intrinsic
% and extrinsic) coupling parameters (at the true values of the remaining
% parameters)
% =========================================================================
beta      = linspace(-1,1,32);                       % range of parameters
M.nograph = 1;
M.Nmax    = 1;
M.Gmax    = 1;
M.Hmax    = 1;

M.P   = P;
M.Q   = G;
for i = 1:length(beta)
    for j = 1:length(beta)
        
        % specify posterior expectations (and implicitly precisions)
        % -----------------------------------------------------------------
        M.P           = P;
        M.P.B{1}(1,1) = beta(i);
        M.P.B{1}(2,1) = beta(j);
        
        % evaluate free energy (without updating expectations)
        % -----------------------------------------------------------------
        [Ep,Eg,Cp,Cg,S,F] = spm_nlsi_N(M,U,Y);
        FF(i,j)           = F;
        
        % record conditional covariance at maximum
        % -----------------------------------------------------------------
        if F >= max(FF(:)), CP  = Cp; end

    end
end

% apply Occam's window
% -------------------------------------------------------------------------
FF    = FF - max(FF(:));
F     = max(FF,-64)

% Free energy landscape and associated conditional covariance: the
% conditional covariance  is the inverse precision – which is proportional
% to the curvature of the variational energy (equivalent to the curvature
% of the free energy functional of posterior or conditional expectations)
% -------------------------------------------------------------------------
spm_figure('Getwin','Figure 3');

j     = find(spm_vec(M.pC));
subplot(2,1,1)
imagesc(beta,beta,F),       hold on
contour(beta,beta,F,8,'m'), hold on
plot(pP(j(2)),pP(j(1)),'.r','MarkerSize',32), hold off
xlabel('extrinsic');ylabel('intrinsic');
title('Free-energy landscape','FontSize',16)
axis square

subplot(2,1,2)
imagesc(CP(j,j))
xlabel('parameters');ylabel('parameters');
title('Posterior covariance','FontSize',16)
axis square

In [ ]:
%%matlab