function [pE,pC] = spm_hdm_priors(m,h)
% returns priors for a hemodynamic dynmaic causal model
% FORMAT [pE,pC] = spm_hdm_priors(m,[h])
% m.. - number of inputs
% h   - number of hemodynamic modes (default = 2)
%
% pE  - prior expectations
% pC  - prior covariances
%
% (5) biophysical parameters
%    P(1) - signal decay      - d(ds/dt)/ds)
%    P(2) - autoregulation    - d(ds/dt)/df)
%    P(3) - transit time                (t0)
%    P(4) - exponent for Fout(v)     (alpha)
%    P(5) - resting oxygen extraction   (E0)
%
% plus (m) efficacy priors
%    P(6) - ....
%
%___________________________________________________________________________
% @(#)spm_hdm_priors.m	2.6 Karl Friston 03/02/18


% default 2 - hemodynamic [eigen]modes
%---------------------------------------------------------------------------
if nargin < 2
	h = 2;
end

% biophysical parameters with prior expectation and
%---------------------------------------------------------------------------
pE    = [   0.65      0.41      0.98      0.32      0.34  ];

% covariance restricted to h modes (v) scaled by eigenvales (e) {see below)
%---------------------------------------------------------------------------
v     = [  -0.0535    0.0095   -0.1117   -0.0040   -0.0026
	   -0.0604   -0.0319    0.0430   -0.0077    0.0026
	    0.1116   -0.0347   -0.2539   -0.0169   -0.0115
	    0.1985    0.1698    0.4984   -0.4493    0.4434
	    0.0029    0.2081    1.9582   -0.5209   -1.1634]';

e     = [   2.1225    1.2006    0.3519    0.0039    0.0012];

% set variance of minor modes to zero
%---------------------------------------------------------------------------
i     = (h + 1):5;
e(i)  = 0;
pC    = v*diag(e)*v';

% append m efficacy priors
%---------------------------------------------------------------------------
pE    = [pE(:); zeros(m,1)];
pC    = blkdiag(pC,eye(m));

return


% NOTES: sample covariances from Friston et al (2000)
%---------------------------------------------------------------------------
qC    = [   0.0150    0.0052    0.0283    0.0002   -0.0027
	    0.0052    0.0020    0.0104    0.0004   -0.0013
	    0.0283    0.0104    0.0568    0.0010   -0.0069
	    0.0002    0.0004    0.0010    0.0013   -0.0010
	   -0.0027   -0.0013   -0.0069   -0.0010    0.0024];


% NOTES: Reduce rank of prior covariances for computational expediancy
%---------------------------------------------------------------------------

% assume independent priors in parameter space
%---------------------------------------------------------------------------
qC    = diag(diag(qC));


% model specification (single node DCM)
%---------------------------------------------------------------------------
M.f   = 'spm_fx_HRF';
M.g  = 'spm_lambda_HRF';
M.x   = [0 1 1 1]';
M.pE  = [pE(1:5) 1];
M.m   = 1;
M.n   = 4;
M.l   = 1;
M.N   = 32;
M.dt  = 1/2;

% compute partial derivatives w.r.t. hemodynamic parameters [J] dy(t)/dp
%---------------------------------------------------------------------------
P     = M.pE;
p     = length(P);
dp    = 1e-6;
[k J] = spm_nlsi(M);
for i = 1:5
	M.pE    = P;
	M.pE(i) = M.pE(i) + dp;
	[k q]   = spm_nlsi(M);
	Jq(:,i) = (q - J)/dp;
end

% implied covariance of impulse response
%---------------------------------------------------------------------------
Cq    = Jq*qC*Jq';

% reduce to h hemodynamic modes in measurement space
%---------------------------------------------------------------------------
[v e] = spm_svd(Cq);
e     = diag(e);
v     = pinv(Jq)*v;
qC    = v*diag(e)*v';


% NOTES: graphics - eigenvalues of qC
%---------------------------------------------------------------------------
subplot(2,2,1)
bar(e)
xlabel('eigen mode')
title('eigenvalue')
set(gca,'XLim',[0 6])
axis square
grid on

% graphics - response differentials
%---------------------------------------------------------------------------
subplot(2,2,2)
plot([1:M.N]*M.dt,Jq*v(:,1),[1:M.N]*M.dt,Jq*v(:,2),'-.')
xlabel('PST {secs}')
title('hemodynamic modes')
axis square
grid on
