Tmax=2000;

net=readNet('competence.net');

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

disp('first-order model')
[mdyn1,x10]=closureDynamics(net,1,'derMatch','fun');

equations=cell(subsParameters(net,mdyn.dotMu.approxDotMu));
variables=cell(mdyn.Mu.sym);
solution=solve(equations{:},variables{:});

disp('Equilibrium point for P:')
solution.mu_P

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


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

disp('second-order model')
[mdyn2,x20]=closureDynamics(net,2,'derMatch','fun');

[t2,mu2]=ode23s(@(t,x)fun(x),[0,Tmax],x20,options);
cmu2=centerMoments(ndxMu,mu2);
ndxG =[1];
ndxP =[2,8];
ndxP2=[3,11];

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

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

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SSA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
net2stochKit(net,'ProblemDefinition',[G0,P0]);
%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,5,2)
plot(tx(:,1),tx(:,ndxG(1)+1),':',tx(:,1),tx(:,ndxP(1)+1),'-',tx(:,1),tx(:,ndxP2(1)+1),'--')
legend('G','P','P2','location','best')
title('MC')
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,5,3)
hist(xfinal(:,ndxP(1)),20)
set(gca,'yscale','log')
x=1e-1:.1:max(get(gca,'xlim'));
si=sqrt(log(1+cmu2(end,ndxP(2))/cmu2(end,ndxP(1))^2));
mu=log(cmu2(end,ndxP(1)))-si/2;
hold on
plot(x,size(xfinal,1)./x/si/sqrt(2*pi).*exp(-((log(x)-mu).^2)/2/si^2))
hold off
title(sprintf('MC: E[P]=%0.2g, Std[P]=%0.2g',mean(xfinal(:,ndxP(1))),std(xfinal(:,ndxP(1)))))

subplot(1,5,4)
hist(xfinal(:,ndxP2(1)),20)
set(gca,'yscale','log')
x=1e-1:.1:max(get(gca,'xlim'));
si=sqrt(log(1+cmu2(end,ndxP2(2))/cmu2(end,ndxP2(1))^2));
mu=log(cmu2(end,ndxP2(1)))-si/2;
hold on
plot(x,size(xfinal,1)./x/si/sqrt(2*pi).*exp(-((log(x)-mu).^2)/2/si^2))
hold off
title(sprintf('MC: E[P2]=%0.2g, Std[P2]=%0.2g',mean(xfinal(:,ndxP2(1))),std(xfinal(:,ndxP2(1)))))

subplot(1,5,5)
hist(xfinal(:,ndxG(1)),20)
%set(gca,'yscale','log')
hold on
plot([0,1],size(xfinal,1)*[1-cmu2(end,ndxG(1)),cmu2(end,ndxG(1))])
hold off
title(sprintf('MC: E[G]=%0.2g, Std[G]=%0.2g',mean(xfinal(:,ndxG(1))),std(xfinal(:,ndxG(1)))))
drawnow
