% This example is discussed in the documentation of StochDynTools

clear functions
clear all

%%%%%%%%%%%%%%%%
% read .net file
%%%%%%%%%%%%%%%%
net=readNet('enzyme_with_Et_St.net');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% 1st order trunaction & latex output
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
mdyn1=closureDynamics(net,1,'derMatch');
dotMu=expand(mdyn1.approxDotMu.sym);
% dotMu=expand(subsParameters(net,dotMu));
latexmacros=mylatex()
latexMu=mylatex(mdyn1.Mu.sym,mdyn1.texrules)
latexDotMu=mylatex(dotMu,mdyn1.texrules)

%%%%%%%%%%%%%%%%%%%
% equilibrium point
%%%%%%%%%%%%%%%%%%%
muEq=momentEquilibrium(net,mdyn1);

%%%%%%%%%%%%%%%%%%%%%%
% 2nd order trunaction
%%%%%%%%%%%%%%%%%%%%%%
mdyn2=closureDynamics(net,2,'derMatch','fun');

%%%%%%%%%%%%
% simulation 
%%%%%%%%%%%%
Tmax=300;                      % maximum simulation time
[t,mu]=ode23s(@(t,x)fun(x),[0,Tmax],mdyn2.Mu.x0);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot of uncentered moments
%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(1);clf
plot(t,mu)
legend(fun(),'location','best','interpreter','none')
print -depsc uncentered.eps

%%%%%%%%%%%%%%%%%%%%%%%%%%
% plot of centered moments
%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(2);clf
h=plotCMoments(net,mdyn2,t,mu,'EF','b-','St','g-','Et','k-');
legend(h(1:3:end),'E[EF]\pm{}Std[EF]','E[St]\pm{}Std[St]','E[Et]','location','best');
print -depsc centered.eps

%%%%%%%%%%%%%%%%%%%%%%
% plot of distribution
%%%%%%%%%%%%%%%%%%%%%%
figure(3);clf
x=0:.01:4;
pdf=getDistribution(net,mdyn2,mu(end,:),'St',x,'lognormal');
plot(x,pdf)
print -depsc pdf.eps

%%%%%%%%%%%%%%%%%%%%%%%%%
% Monte Carlo Simulations
%%%%%%%%%%%%%%%%%%%%%%%%%
[Q,b,c,s,x0]=quadPropensities(net);
Q=double(subsParameters(net,Q));      % replace parameters
b=double(subsParameters(net,b));      % by their numerical values
c=double(subsParameters(net,c));
s=double(subsParameters(net,s));
nMC=1000;
Ts=(0:1:300)';
[X,Xmean,Xstd]=sampledSSA(Q,b,c,s,x0,nMC,Ts);

figure(4);clf
set(4,'PaperPositionMode','manual','papersize',[4,8.5])
orient(4,'portrait')
% plot Mean, Std.Dev. and Monte Carlos runs
StMC=reshape(X(3,:,:),nMC,[])';   % extract only St values
subplot(1,2,1)   
h=plot(Ts,Xmean(3,:)','b-',Ts,[Xmean(3,:)'-Xstd(3,:)',Xmean(3,:)'+Xstd(3,:)'],'b:',Ts,StMC(:,1:10),'k-');
axis([0,Ts(end),0,1.1*x0(3)]);
legend(h([1,2,4]),'E[St]','E[St]\pm{}Std[St]','Monte Carlo runs')

% plot distribution
subplot(1,2,2)
[n,x]=hist(StMC(end,:),1:min(20,max(StMC(end,:))));
barh(x,n)
axis([0,max(n),0,1.1*x0(3)]);
legend(sprintf('histogram at time %g',Ts(end)))
print -depsc mc.eps
