
% Load All parameters

reactiondata; 
pars = {NAV,k,c,M,p,l,n0,S0,nu,tf,runs}; % parameters cell array

% Run stochastic simulation and solve GMA 

cells = stoch_gma(pars);

% extract cell arrays

t = cells{1}; 
n = cells{2}; 
t_gma = cells{3};  
n_gma = cells{4};  
mean_n = cells{5};  
t_gilsp = cells{7};    

% Plotting results:

% First plotting the single run stochastic simulation.

figure

hstair = stairs(t{1},n{1}(:,[1 4])); % plot first realisation
set(hstair(1),'LineWidth',1.5,'LineStyle','-','Color','blue'); % Making one line slightly thicker for b/w printing.
set(hstair(2),'LineWidth',1,'LineStyle','-','Color','red');
hold on
hplot = plot(t_gma,n_gma(:,1),'b-.',t_gma,n_gma(:,4),'r:');      % plot GMA's prediction
set(hplot(2),'LineWidth',2)
hold off
l_hdl = legend([hstair;hplot],'{\langle}#{{\it\fontname{Times}S_1}\rangle} (stochastic)','{\langle}#{{\it\fontname{Times}S_4}\rangle} (stochastic)','#{\fontname{Times}\itS_1} (GMA)','#{\fontname{Times}\itS_4} (GMA)',0);

xlabel(['time  (sec) '],'FontSize',8)
ylabel('count of molecules','FontSize',8)
%title('\fontsize{8} single stochastic realization vs GMA solution')

axis([0 tf 0 max(n0)+1])

% Shifting the axis within the figure

posi = get(gca,'Position');
set(gca,'Position',posi.*[1 1.15 1 1],'FontSize',8);

% positioning the legend:

l_posi= get(l_hdl,'Position');
set(l_hdl,'Position',[(posi(3)-l_posi(3))*0.96 posi(2)+0.5*posi(4)-0.5*l_posi(4) l_posi(3) l_posi(4)])

% printing adjustments

set(gcf, 'PaperType', 'A4', 'PaperOrientation', 'portrait', 'PaperPositionMode', 'manual', 'PaperUnits', 'centimeters', 'PaperPosition', [2 2 8 6.3]);

textobj = findobj('type', 'text');
set(textobj, 'fontunits', 'points', 'fontsize', 8);

% Saving the figure into a .eps file:

print  -depsc 'single_run2.eps'

% Now plotting the average of multiple runs:

figure

hstair = stairs(t_gilsp,mean_n(:,[1 4])); % plot mean over realisations

set(hstair(1),'LineWidth',1.5,'LineStyle','-','Color','blue'); % Making one line slightly thicker for b/w printing.
set(hstair(2),'LineWidth',1,'LineStyle','-','Color','red');
hold on
hplot = plot(t_gma,n_gma(:,1),'b-.',t_gma,n_gma(:,4),'r:');       % plot GMA's prediction
set(hplot(2),'LineWidth',1.5)
hold off
xlabel('time  (sec)','FontSize',8)
ylabel('count of molecules','FontSize',8)
%title(['\fontsize{8} mean over ',int2str(runs),' realizations vs GMA solution'])
l_hdl = legend([hstair;hplot],'{\langle}#{{\it\fontname{Times}S_1}\rangle} (stochastic)','{\langle}#{{\it\fontname{Times}S_4}\rangle} (stochastic)','#{\fontname{Times}\itS_1} (GMA)','#{\fontname{Times}\itS_4} (GMA)',0);

axis([0 tf 0 max(n0)+1])

% Shifting the axis within the figure

posi = get(gca,'Position');
set(gca,'Position',posi.*[1 1.15 1 1],'FontSize',8);

% positioning the legend:

l_posi= get(l_hdl,'Position');
set(l_hdl,'Position',[(posi(3)-l_posi(3))*0.96 posi(2)+0.5*posi(4)-0.5*l_posi(4) l_posi(3) l_posi(4)])

% printing adjustments

set(gcf, 'PaperType', 'A4', 'PaperOrientation', 'portrait', 'PaperPositionMode', 'manual', 'PaperUnits', 'centimeters', 'PaperPosition', [2 2 8 6.3]);

textobj = findobj('type', 'text');
set(textobj, 'fontunits', 'points', 'fontsize', 8);

% Printing the figure into an .eps file:

print -depsc 'multiple_run2.eps'