https://github.com/Data2Dynamics/d2d
Tip revision: 4c697ee160f4d91719a65cfcfe3e3b50395b6656 authored by Helge Hass on 29 July 2015, 07:45:33 UTC
new C code version
new C code version
Tip revision: 4c697ee
arSSA.m
% Stochastic Simulation Algorithm (SSA)
% see in Gillespie, The Journal of Physical Chemistry, 1977, 81(25)
%
% arSSA(nruns,scaling)
%
% nruns: number of runs [10]
% scaling: rescaling factor from species [{ones}]
% concentration to number of molecules
function arSSA(nruns, scaling)
global ar
if(~exist('nruns','var'))
nruns = 10;
end
ar.config.ssa_runs = nruns;
% update initial concentrations
arSimu(false, true);
for m=1:length(ar.model)
if(~exist('scaling','var'))
scaling{m} = ones(1,length(ar.model(m).x));
end
% flux scaling
kscaling = zeros(1, length(ar.model(m).fv));
for j=1:length(ar.model(m).fv)
qi = find(ar.model(m).N(:,j)<0);
if(~isempty(qi))
kscaling(j) = scaling{m}(qi(1));
else
qi = find(ar.model(m).N(:,j)>0);
if(~isempty(qi))
kscaling(j) = scaling{m}(qi(1));
else
error('reaction with empty N');
end
end
end
% setup arrays
for c=1:length(ar.model(m).condition)
ar.model(m).condition(c).xFineSSA = nan(length(ar.model(m).condition(c).tFine), ...
length(ar.model(m).x), nruns);
ar.model(m).condition(c).xFineSSA_lb = nan(length(ar.model(m).condition(c).tFine), ...
length(ar.model(m).x), nruns);
ar.model(m).condition(c).xFineSSA_ub = nan(length(ar.model(m).condition(c).tFine), ...
length(ar.model(m).x), nruns);
ar.model(m).condition(c).zFineSSA = nan(length(ar.model(m).condition(c).tFine), ...
length(ar.model(m).z), nruns);
ar.model(m).condition(c).zFineSSA_lb = nan(length(ar.model(m).condition(c).tFine), ...
length(ar.model(m).z), nruns);
ar.model(m).condition(c).zFineSSA_ub = nan(length(ar.model(m).condition(c).tFine), ...
length(ar.model(m).z), nruns);
if(isfield(ar.model(m).condition(c), 'tExp'))
ar.model(m).condition(c).xExpSSA = nan(length(ar.model(m).condition(c).tExp), ...
length(ar.model(m).x), nruns);
ar.model(m).condition(c).zExpSSA = nan(length(ar.model(m).condition(c).tExp), ...
length(ar.model(m).z), nruns);
end
ar.model(m).condition(c).x0_ssa = ar.model(m).condition(c).xFineSimu(1,:);
ar.model(m).condition(c).scale_x_ssa = scaling{m};
ar.model(m).condition(c).scale_v_ssa = kscaling;
end
if(isfield(ar.model(m),'data'))
for d=1:length(ar.model(m).data)
ar.model(m).data(d).yFineSSA = nan(length(ar.model(m).data(d).tFine), ...
length(ar.model(m).data(d).y), nruns);
ar.model(m).data(d).yFineSSA_lb = nan(length(ar.model(m).data(d).tFine), ...
length(ar.model(m).data(d).y), nruns);
ar.model(m).data(d).yFineSSA_ub = nan(length(ar.model(m).data(d).tFine), ...
length(ar.model(m).data(d).y), nruns);
if(isfield(ar.model(m).data(d), 'tExp'))
ar.model(m).data(d).yExpSSA = nan(length(ar.model(m).data(d).tExp), ...
length(ar.model(m).data(d).y), nruns);
end
end
end
end
% fprintf('Stochastic simulation (%i runs)...', ar.config.ssa_runs); tic;
% call mex function to simulate models
feval(ar.fkt, ar, true, false, true, true, 'condition', 'threads')
% fprintf('done (%g sec)\n', toc);