clear functions

Tmax=400;

net=readNet('promoter_linear.net');

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

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

numapproxDotMu=cell(subsParameters(net,approxDotMu));
sol=solve(numapproxDotMu{:});
disp('Equilibrium point for F:')
if ~isempty(sol)
  disp(sol.mu_F)
else
  disp('none')  
end  

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

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

disp('second-order model')
[ndxMu,ndxBarMu,A,B,monMu,monBarMu,x20]=momentDynamics(net,2);
[approxDotMu,symMu,monMu2,ndxMu,monBarMu,x20]=closureDynamics(net,2,'derMatch','fun');

[monMu,B*monBarMu]
monBarMu

% ODE solution
[t2,mu2]=ode23s(@(t,x)fun(x),[0,Tmax],x20,options);

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% ODE plots
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(1);clf
subplot(1,4,1)
h1=plotCMoments(net,monMu1,t1,mu1,'F','k-');
hold on
h2=plotCMoments(net,monMu2,t2,mu2,'F','b-');
hold off
grid on
legend([h1(1),h2(1)],'1st order','2nd order','location','best')
[af,sf]=getCMoments(net,monMu,mu2(end,:),'F');
[ag,sg]=getCMoments(net,monMu,mu2(end,:),'GF');
title(sprintf('Moment truncation: E[F]\\pm{}Std[F]\nE[F]=%0.2g, Std[F]=%0.2g (2nd order)\nE[GF]=%0.2g, Std[GF]=%0.2g (2nd order)',af,sf,ag,sg));
drawnow

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SSA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
net2stochKit(net,'ProblemDefinition',x10);
%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 %0.2f single.txt',Tmax);
fprintf('executing: ''%s''...\n',cmd);
rc=system(cmd);
tx=load('single.txt');

subplot(1,4,2)
plot(tx(:,1),tx(:,2:end))
legend('GF','m','F','location','best')
title('Monte Carlo sample')
grid on
drawnow

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

subplot(1,4,3)
[n,x]=hist(xfinal(:,3),20);
bar(x,n);
%set(gca,'yscale','log')
xx=min(x):(max(x)-min(x))/1000:max(x);
pdf=getDistribution(net,monMu,mu2(end,:),'F',xx,'lognormal');
hold on
plot(xx,size(xfinal,1)*pdf*(x(2)-x(1)))
hold off
title(sprintf('Monte Carlo: F\nE[F]=%0.2g, Std[F]=%0.2g',mean(xfinal(:,3)),std(xfinal(:,3))))

subplot(1,4,4)
hist(xfinal(:,1),20)
x=0:1;
pdf=getDistribution(net,monMu,mu2(end,:),'GF',x,'binomial');
hold on
plot(x,size(xfinal,1)*pdf)
hold off
title(sprintf('Monte Carlo: GF\nE[GF]=%0.2g, Std[GF]=%0.2g',mean(xfinal(:,1)),std(xfinal(:,1))))
drawnow

