net=readNet('Dicty_Laub_Loomis.net');

%% Plot comparison of methods

Tmax=400;
X5minmax=[0,40000];
X4minmax=[0,40000];
nMC=4;

figure(1);clf
[exact,approx,MC]=plotCompare(net,2,2,Tmax,...
... %% cAMPi
    {{'mass action kin.',X5minmax},{'dm',1,'X5','k-'}},...%{'MCsample',1,'X5','c-'},
    {{'derivative matching',X5minmax},{'dm',2,'X5','b-'}},...%{'dm',3,'X5','r-'}},...
    {{'zero cumulants',X5minmax},{'zc',2,'X5','b-'}},...%,{'zc',3,'X5','r-'}},...
    {{'quasi deterministic',X5minmax},{'qd',2,'X5','b-'}}...%,{'qd',3,'X5','r-'}},...
...%    {{'distribution (Monte Carlo)',X5minmax},'MCdistribution',nMC,'X5'},...
... %% RegA
...%    {{'mass action kin.',X4minmax},{'dm',1,'X4','k-'}},...%{'MCsample',1,'X4','c-'},
...%    {{'derivative matching',X4minmax},{'dm',2,'X4','b-'}},...%,{'dm',3,'X4','r-'}},...
...%    {{'zero cumulants',X4minmax},{'zc',2,'X4','b-'}},...%,{'zc',3,'X4','r-'}},...
...%    {{'quasi deterministic',X4minmax},{'qd',2,'X4','b-'}}...%,{'qd',3,'X4','r-'}}...
...%    {{'distribution (Monte Carlo)',X4minmax},'MCdistribution',nMC,'X4'}...
);

%print -depsc cAMPi.ps

return

%% Compute closure and determine derivative matching

mdyn2dm=closureDynamics(net,2,'derMatch');
mdyn2zc=closureDynamics(net,2,'zeroCumulants');

%% Compute (exact) N derivatives of the state (for deterministic values)

mpnb=symengine;

N=2;  % # of derivatives to compute

psi=mdyn2dm.Mu.mon;

LpsiExact=sym(zeros(length(psi),N));

for i=1:N
    fprintf('Computing (exact) derivative %1... ',i);
    LpsiExact(:,i)=simplify(applyGenerator(net,psi,mpnb));

    %    LpsiExact(:,i)=simplify(subsParameters(net,LpsiExact(:,i),mpnb));

    fprintf('done\n');
    
    psi=LpsiExact(:,i);
end

%% Compute approximate N derivatives of the state (for deterministic values)

LpsiApproxDm=sym(zeros(length(psi),N));

ELpsi=mdyn2dm.approxDotMu.sym;
LpsiApproxDm(:,1)=simplify(mysubs(ELpsi,mdyn2dm.Mu.sym,mdyn2dm.Mu.mon));
for i=2:N
    fprintf('Computing approximate derivative %d... ',i);

    ELpsi=simplify(jacobian(ELpsi,mdyn2dm.Mu.sym)*mdyn2dm.dotMu.sym);
    
    % go back to monomials
    LpsiApproxDm(:,i)=simplify(mysubs(ELpsi,mdyn2dm.Mu.sym,mdyn2dm.Mu.mon));
    LpsiApproxDm(:,i)=simplify(mysubs(LpsiApproxDm(:,i),mdyn2dm.barMu.sym,mdyn2dm.barMu.mon));
    
    %    LpsiApproxDm(:,i)=simplify(subsParameters(net,LpsiApproxDm(:,i),mpnb));

    fprintf('done\n');

end

%% Derivative Error

for i=1:N
    fprintf('Derivative error %d:\n',i);
    
    disp(simplify(LpsiExact(:,i)-LpsiApproxDm(:,i)));
    
end

return

