clear functions

net=readNet('stochasticfocusing.net');

Tmax=4;
X1minmax=[0,4];
X2minmax=[0,22];
X3minmax=[0,3];
X4minmax=[880,1080];
nMC=400;

figure(1);clf
[exact,approx,MC]=plotCompare(net,4,6,Tmax,...
    {{'derivative matching',X1minmax},{'MCsample',1,'X1','c-'},{'dm',2,'X1','b-'}},...
    {{'zero cumulants',X1minmax},{'zc',1,'X1','k-'},{'zc',2,'X1','b-'}},...
    {{'low dispersion',X1minmax},{'ld',1,'X1','k-'},{'ld',2,'X1','b-'}},...
    {{'quasi deterministic',X1minmax},{'qd',1,'X1','k-'},{'qd',2,'X1','b-'}},...
    {{'linear noise approx.',X1minmax},{{'lna','V'},1,'X1','k-'},{{'lna','V'},2,'X1','b-'}},... 
    {{'distribution (Monte Carlo)',X1minmax},'MCdistribution',nMC,'X1'},...
    {{'derivative matching',X2minmax},{'MCsample',1,'X2','c-'},{'dm',2,'X2','b-'}},...
    {{'zero cumulants',X2minmax},{'zc',1,'X2','k-'},{'zc',2,'X2','b-'}},...
    {{'low dispersion',X2minmax},{'ld',1,'X2','k-'},{'ld',2,'X2','b-'}},...
    {{'quasi deterministic',X2minmax},{'qd',1,'X2','k-'},{'qd',2,'X2','b-'}},...
    {{'linear noise approx.',X2minmax},{{'lna','V'},1,'X2','k-'},{{'lna','V'},2,'X2','b-'}},... 
    {{'distribution (Monte Carlo)',X2minmax},'MCdistribution',nMC,'X2'},...
    {{'derivative matching',X3minmax},{'MCsample',1,'X3','c-'},{'dm',2,'X3','b-'}},...
    {{'zero cumulants',X3minmax},{'zc',1,'X3','k-'},{'zc',2,'X3','b-'}},...
    {{'low dispersion',X3minmax},{'ld',1,'X3','k-'},{'ld',2,'X3','b-'}},...
    {{'quasi deterministic',X3minmax},{'qd',1,'X3','k-'},{'qd',2,'X3','b-'}},...
    {{'linear noise approx.',X3minmax},{{'lna','V'},1,'X3','k-'},{{'lna','V'},2,'X3','b-'}},... 
    {{'distribution (Monte Carlo)',X3minmax},'MCdistribution',nMC,'X3'},...
    {{'derivative matching',X4minmax},{'MCsample',1,'X4','c-'},{'dm',2,'X4','b-'}},...
    {{'zero cumulants',X4minmax},{'zc',1,'X4','k-'},{'zc',2,'X4','b-'}},...
    {{'low dispersion',X4minmax},{'ld',1,'X4','k-'},{'ld',2,'X4','b-'}},...
    {{'quasi deterministic',X4minmax},{'qd',1,'X4','k-'},{'qd',2,'X4','b-'}},...
    {{'linear noise approx.',X4minmax},{{'lna','V'},1,'X4','k-'},{{'lna','V'},2,'X4','b-'}},... 
    {{'distribution (Monte Carlo)',X4minmax},'MCdistribution',nMC,'X4'});

return

mdyn=closureDynamics(net,2,'ld');
[muEq,phiEq,xiEq]=momentEquilibrium(net,mdyn,1,0,0);
for k=1:length(muEq)
  fprintf('Solution %d/%d\n',k,length(muEq));  
  for i=1:length(muEq{k})
    fprintf('    %s\t=%s\n',char(mdyn.Mu.sym(i)),char(muEq{k}(i)))
  end  
end  

