function DCM = J_dcm_model_specify(d)
% This function (1) specifies a generative model, (2) preprocess an
% individual human subject's data d, and (3) package up as a DCM.
%
% ---------- Frances Jin edited based on Karl Friston's script ------------
% October. 2021

% Set up the generative model
%==========================================================================

% dynamics and parameters of a Lorentz system (with Jacobian)
%==========================================================================
% dxdt = f(x,v)             Underlying dynamics with exogenous inputs
% s    = h(d(x)) + e        Observer model generating observed data
%--------------------------------------------------------------------------

% level 1: Thresholding function (M.h) that converts continuous states of
% psychopathology into discrete scores (here between 0 and 1) with
% thresholds parameterised by P.u. This function can be changed to apply
% multiple thresholds to generate discrete scales of the sort found in
% clinical instruments.
%--------------------------------------------------------------------------
M.h  = @(x,P)spm_phi([x(:,1) - P.u(1), x(:,2) - P.u(2), x(:,3) - P.u(3)]*8);

% level 2: the level that generates latent (psychopathology)causes y = g(x) 
% from(pathophysiological) states (x) that are subject to interventions (v).
% The first equation (M.f) models the dynamics or flow of underlying
% pathophysiology using a Lorentz system, while the second equation (M.d)
% extracts a subset of physiological states and mixes them to generate the
% psychological states (y). The parameters of these functions determine the
% underlying dynamics of the disease (P.A) and its psychopathological
% manifestation (P.B)
%--------------------------------------------------------------------------
M.f  = @(x,v,P,M)[-P.A(1) P.A(1) 0; ((1 - v*P.A(3))*P.A(2) - x(3)) -1 0; x(2) 0 -8/3]*x/P.t;
M.d  = @(x,P) x*P.B/16;

% level 3: the level that perturbed the physiological dynamics by
% intervening on the parameters governing flow. These interventions or
% exogenous inputs could be regarded as life events or pharmacological
% interventions. This exogenous input variable (M.v) is smooth and  
% fluctuating changes that are parameterised with the coefficients of a 
% discrete cosine function of time (P.v).
%--------------------------------------------------------------------------
M.v  = @(T,P) spm_dctmtx(T,numel(P.C))*P.C(:);

% Having established the structure of the generative model we can now
% specify the prior expectations over the requisite parameters; including
% the overall rate at which the dynamics are expressed (P.t). Here, P.x
% parameterises the initial states at the beginning of the timeseries.
%--------------------------------------------------------------------------
pE.A = [10 16 1];                      % parameters of pathophysiology
pE.B = [0 1 0;1 0 0;0 0 1];            % mapping to psychopathology
pE.C = [0 0 0 0 0 0 0 0];              % parameters of perturbations
pE.x = [16; 4; 32];                    % initial physiological states
pE.u = [1,1,1];                        % threshold for (three) outcomes
pE.t = 128;                            % time constant of dynamics

% The priors are completed by specifying the prior covariance, that
% determines which parameters are fixed and which are unknown and need to
% be estimated
%--------------------------------------------------------------------------
pC.A = spm_zeros(pE.A) + 1;
pC.A(1) = 0; pC.A(3) = 0;              % allow only the 2nd param to vary          
pC.B = spm_zeros(pE.B) + 1;
pC.B(3) = 0; pC.B(5) = 0; pC.B(7) = 0; % B([3, 5, 7]) are fixed based
pC.C = spm_zeros(pE.C) + 1/64;
pC.x = spm_zeros(pE.x) + 1;
pC.u = spm_zeros(pE.u) + 0;
pC.t = 1;

% model inversion with Variational Laplace
%==========================================================================

% There are two special aspects to this inversion. First, because we are
% dealing with a highly non-linear system, there can be many local maxima.
% This can be resolved using a multi-start algorithm. In other words, by
% initialising the priors in the basins of attraction of all local
% minimaone can  identify the model with the greatest free energy or
% marginal likelihood, which corresponds to the global maxima. This is 
% implemented in the model estimation step. 
% Second, because this inversion scheme uses gradient descent, it is 
% necessary to use a particular kind of feature selection (M.FS) function 
% that converts variations in the timing of the discrete changes in scores 
% into variations in amplitude. This smoothing is done automatically in 
% the same way that a link function converts a general linear model into a
% generalised linear model.

% model specification
%==========================================================================
M.Nmax = 32;                       % maximum number of iterations
M.G    = @spm_nosology_gen;        % generative function
M.FS   = @(Y)spm_conv(Y,8,0);      % feature selection  (link function)
M.pE   = pE;                       % prior expectations (parameters)
M.pC   = diag(spm_vec(pC));        % prior covariances  (parameters)
M.hE   = 2;                        % prior expectation  (log-precision)
M.hC   = 1/64;                     % prior covariances  (log-precision)

% pre-process the empirical data to eliminate NaNs if exists.
%==========================================================================

% check if NaN exists
for y = 1:size(d, 2)
    nanidx = isnan(d(:, y));
    if sum(nanidx)>0
        d(nanidx,:) = [];
    end
end
clear y nanidx
Y = d;

% Package up as a DCM
T       = size(Y,1);                % number of assessments
U       = zeros(T,1);               % Prior expectation of fluctuations
DCM      = struct();
DCM.M    = M;
DCM.xU   = U;
DCM.xY   = Y;
DCM.M.IS = DCM.M.G; % needed for spm_dcm_fit