clear functions

%net=readNet('VKp244.net');
net=readNet('VKp263.net');
%net=readNet('../stochasticfocusing.net');

mdyn=momentDynamics(net,2);
disp('exact: Mu, A, B, barMu')
disp(mdyn.Mu.sym)
disp(mdyn.dotMu.A)
disp(mdyn.dotMu.B)
disp(mdyn.barMu.sym)

mdyn_dm=momentClosure(net,mdyn,'derivativeMatching');
disp('dermatch: approx barMu')
disp(mdyn_dm.approxDotMu.barMu)

mdyn_zc=momentClosure(net,mdyn,'zeroCumulants');
disp('zeroCumulants: approx barMu')
disp(mdyn_zc.approxDotMu.barMu)

mdyn_ld=momentClosure(net,mdyn,'lowDispersion');
disp('lowDispersion: approx barMu')
disp(mdyn_ld.approxDotMu.barMu)

mdyn_qd=momentClosure(net,mdyn,'quasiDeterministic');
disp('quasiDeterministic: dotPhi, approx barMu, approxdotMu.A, approxdotMu.c')
disp(mdyn_qd.dotPhi)
disp(mdyn_qd.approxDotMu.barMu)
disp(mdyn_qd.approxDotMu.A)
disp(mdyn_qd.approxDotMu.c)

mdyn_lna=momentClosure(net,mdyn,{'VanKampen','V'});
disp('linearNoiseApproximation: dotPhi, dotXi, approxDotXi.A, approxDotXi.c, approxDotMu.A, approxDotMu.c')
disp(mdyn_lna.dotPhi)
disp(mdyn_lna.Xi.sym)
disp(mdyn_lna.approxDotXi.A)
disp(mdyn_lna.approxDotXi.c)
disp(mdyn_lna.approxDotMu.A)
disp(mdyn_lna.approxDotMu.c)

disp('paused before latex outputs')
pause

fprintf('exact:\n')
fprintf('Mu      = %s\n',mylatex(mdyn.Mu.sym,mdyn.texrules))
fprintf('dotMu.A = %s\n',mylatex(mdyn.dotMu.A,mdyn.texrules))
fprintf('dotMu.B = %s\n',mylatex(mdyn.dotMu.B,mdyn.texrules))
fprintf('barMu   = %s\n',mylatex(mdyn.barMu.sym,mdyn.texrules))

fprintf('\ndermatch:\n')
fprintf('approxBarMu = %s\n',mylatex(mdyn_dm.approxDotMu.barMu,mdyn_dm.texrules))

fprintf('\nzeroCumulants:\n')
fprintf('approxBarMu = %s\n',mylatex(mdyn_zc.approxDotMu.barMu,mdyn_zc.texrules))

fprintf('\nlowDispersion:\n')
fprintf('approxBarMu = %s\n',mylatex(mdyn_ld.approxDotMu.barMu,mdyn_ld.texrules))

fprintf('\nquasiDeterministic:\n')
fprintf('Phi    = %s\n',mylatex(mdyn_qd.Phi.sym,mdyn_qd.texrules))
fprintf('dotPhi = %s\n',mylatex(mdyn_qd.dotPhi,mdyn_qd.texrules))
fprintf('approxBarMu   = %s\n',mylatex(mdyn_qd.approxDotMu.barMu,mdyn_qd.texrules))
fprintf('approxDotMu.A = %s\n',mylatex(mdyn_qd.approxDotMu.A,mdyn_qd.texrules))
fprintf('approxDotMu.c = %s\n',mylatex(mdyn_qd.approxDotMu.c,mdyn_qd.texrules))

fprintf('\nlinearNoiseApproximation:\n')
fprintf('Phi    = %s\n',mylatex(mdyn_lna.Phi.sym,mdyn_lna.texrules))
fprintf('dotPhi = %s\n',mylatex(mdyn_lna.dotPhi,mdyn_lna.texrules))
fprintf('Xi = %s\n',mylatex(mdyn_lna.Xi.sym,mdyn_lna.texrules))
fprintf('approxDotXi.A = %s\n',mylatex(mdyn_lna.approxDotXi.A,mdyn_lna.texrules))
fprintf('approxDotXi.c = %s\n',mylatex(mdyn_lna.approxDotXi.c,mdyn_lna.texrules))
fprintf('approxDotMu.A = %s\n',mylatex(mdyn_lna.approxDotMu.A,mdyn_lna.texrules))
fprintf('approxDotMu.c = %s\n',mylatex(mdyn_lna.approxDotMu.c,mdyn_lna.texrules))


Tmax=.1;
nMC=20000;

net.parameter(1).value=2;  % volume
Xminmax=[0,6];
Yminmax=[0,20];

net.parameter(1).value=20;  % volume
Xminmax=[0,26];
Yminmax=[4,35];

figure(1);clf;
papersize=2.25*[8.5 8.5];
set(1,'PaperSize',papersize,'PaperPosition',[0,0,papersize(1),papersize(2)/2])
[exact,approx,MC]=plotCompare(net,2,6,Tmax,...
    {{'derivative matching',Xminmax},{'MCsample',1,'X','c-'},{'dm',2,'X','b-'},{'dm',3,'X','r-'}},...
    {{'zero cumulants',Xminmax},{'zc',1,'X','k-'},{'zc',2,'X','b-'},{'zc',3,'X','r-'}},...
    {{'low dispersion',Xminmax},{'ld',1,'X','k-'},{'ld',2,'X','b-'},{'ld',3,'X','r-'}},...
    {{'quasi deterministic',Xminmax},{'qd',1,'X','k-'},{'qd',2,'X','b-'},{'qd',3,'X','r-'}},...
    {{'Van Kampen',Xminmax},{{'lna','V'},1,'X','k-'},{{'lna','V'},2,'X','b-'},{{'lna','V'},3,'X','r-'}},...
    {{'distribution (Monte Carlo)',Xminmax},'MCdistribution',nMC,'X'},...
    {{'derivative matching',Yminmax},{'MCsample',1,'Y','c-'},{'dm',2,'Y','b-'},{'dm',3,'Y','r-'}},...
    {{'zero cumulants',Yminmax},{'zc',1,'Y','k-'},{'zc',2,'Y','b-'},{'zc',3,'Y','r-'}},...
    {{'low dispersion',Yminmax},{'ld',1,'Y','k-'},{'ld',2,'Y','b-'},{'ld',3,'Y','r-'}},...
    {{'quasi deterministic',Yminmax},{'qd',1,'Y','k-'},{'qd',2,'Y','b-'},{'qd',3,'Y','r-'}},...
    {{'Van Kampen',Yminmax},{{'lna','V'},1,'Y','k-'},{{'lna','V'},2,'Y','b-'},{{'lna','V'},3,'Y','r-'}},...
    {{'distribution (Monte Carlo)',Yminmax},'MCdistribution',nMC,'Y'});

print('-depsc',sprintf('VanKampen_V_%d.eps',net.parameter(1).value))
