
filename='/home/hespanha/matlab/SBML/SBMLmodels/curated/BIOMD0000000003.xml';

MATLAB_SBMLModel=TranslateSBML(filename);
%ViewModel(MATLAB_SBMLModel)

[net,X0]=readSBML(filename) 

Tmax=.01;

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

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

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

figure(1);clf
subplot(1,3,1)
plot(t1,mu1,'k-')
grid on
names=fun();
legend(names,'location','best','interpreter','none')


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%% SSA
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
net2stochKit(net,'ProblemDefinition',mdyn1.Mu.x0);
%type ProblemDefinition.cpp

system('make single stats');

system(sprintf('single %f single.txt',Tmax));
tx=load('single.txt');

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


system(sprintf('stats %f 1000 stats.txt',Tmax));
xfinal=load('stats.txt');

subplot(1,3,3)
hist(xfinal(:,1),20)
title(sprintf('Monte Carlo: E[%s]=%g, Std[%s]=%g',mean(xfinal(:,1)),net.species(1).id,std(xfinal(:,1)),net.species(1).id))
drawnow

