X0=5;
D0=.01;

Tmax=1000;

net=readNet('bistab2.net');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% first order model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('first-order model')
[approxDotMu,symMu,monMu,ndxMu]=closureDynamics(net,1,'derMatch','fun');

disp('Equilibrium point for X:')
numapproxDotMu=subsParameters(net,approxDotMu);
sol=solve([char(numapproxDotMu(1)),'=0,',char(numapproxDotMu(2)),'=0'],'mu_X,mu_D');
if ~isempty(sol)
  disp(sol.mu_X)
else
  disp('none')  
end  

% ODE solution
x0=double(subsmaple(monMu,{'X';'D'},{X0;D0}));

options=odeset('RelTol',1e-6,'AbsTol',1e-6,'Stats','on');
[t1,mu1]=ode23s(@(t,x)fun(x),[0,Tmax],x0,options);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% second order model
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

disp('second-order model')
[approxDotMu,symMu,monMu,ndxMu]=closureDynamics(net,2,'derMatch','fun');

% ODE solution
x0=double(subsmaple(monMu,{'X';'D'},{X0;D0}));

[t2,mu2]=ode23s(@(t,x)fun(x),[0,Tmax],x0,options);
cmu2=centerMoments(ndxMu,mu2);


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1);clf
subplot(1,3,1)
plot(t1,mu1(:,1),'k-',...
     t2,cmu2(:,1),'b-',...
     t2,[cmu2(:,1)-sqrt(cmu2(:,3)),cmu2(:,1)+sqrt(cmu2(:,3))],'b--')
grid on
names=fun();
legend([strcat(names(1),' : 1st order'),...
       strcat(names(1),' : 2nd order'),...
       strcat(names(1),'-\sigma : 2nd order'),...
       strcat(names(1),'+\sigma : 2nd order')],...
       'location','best','interpreter','none')
title(sprintf('Moment truncation: E[X]=%0.2g, Std[X]=%0.2g, E[D]=%0.2g, Std[D]=%0.2g (from 2nd trunc.)',cmu2(end,1),sqrt(cmu2(end,3)),cmu2(end,2),sqrt(cmu2(end,5))))
drawnow

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SSA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
net2stochKit(net,'ProblemDefinition',[X0,D0]);
%type ProblemDefinition.cpp

cmd='make single stats';
fprintf('executing: ''%s''...\n',cmd);
rc=system(cmd);
if rc
  fprintf('\nFailed to make stochKit executable. Will not run Monte Simulations\n');
  return
end

cmd=sprintf('rm -f single.txt;nice single %f single.txt',Tmax);
fprintf('executing: ''%s''...\n',cmd);
rc=system(cmd);
tx=load('single.txt');

subplot(1,3,2)
plot(tx(:,1),tx(:,2),'-',tx(:,1),tx(:,3),'--')
legend('X','D','location','best')
title('Monte Carlo')
drawnow

cmd=sprintf('rm -f stats.txt;nice stats %f 10000 stats.txt',Tmax);
fprintf('executing: ''%s''...\n',cmd);
rc=system(cmd);
xfinal=load('stats.txt');

subplot(1,3,3)
hist(xfinal(:,1),20)
set(gca,'yscale','log')
title(sprintf('Monte Carlo: E[X]=%0.2g, Std[X]=%0.2g, E[D]=%0.2g, Std[D]=%0.2g (from MC)',mean(xfinal(:,1)),std(xfinal(:,1)),mean(xfinal(:,2)),std(xfinal(:,2))))
drawnow
