function varargout = stoch(varargin)
%
% STOCH  stochastic simulation of a biochemical network
%   [Ts,Ns] = STOCH(X0,N0,K,C,D,L,TF)  performs stochastic simulation a
%   biochemical network composed of M elementary reactions, using Gillespie
%   algorithm. The input and output arguments have the following meaning:
%   
%   X0: Column vector of initial concentrations of all the species involved
%   N0: Column Vector of initial populations of all the species involved 
%   K: Row vector of rate constants of all elementary reactions
%   C: Row vector of stochastic rate constants of all elementary reactions
%   D: Stoichiometry matrix with rows correspoding to species and columns
%      corresponding to reaction channels
%   L: Stoichiometry matrix for reactants only such that L = -D.*(D < 0); 
%   TF: Final time of simulation    
%   Ts: Row vector of time points of reaction events
%   Ns: Matrix of output concentrations with a column for each time point.
% 
%   [Ts,Ns] = STOCH(X0,N0,K,C,D,L,TF,R)  performs R runs of stochastic
%   simulation. The outputs Ts and Ns are cell arrays where each cell
%   corresponds to one run. 
%  
%   [Ts,Ns,TT,NBAR] = STOCH(X0,N0,K,C,D,L,TF,R)  where R>1, also returns 
%   the ensemble of times in the row vector TT and the mean of Ns over R    
%   runs in matrix NBAR of the size of Ns.
%   
%   STOCH(X0,N0,K,C,D,L,TF,1) is the same as STOCH(X0,N0,K,C,D,L,TF).

% preliminary operations:

args = varargin;                        
if nargin<6, args{6} = 1; end           % single run by default
[n0,c,d,l,tf,runs] = args{:};           % parse inputs
oneszk = ones(size(c));
i1 = l==1;                              % reactions of type: A->B, A+B->AB
i2 = l==2;                              % reactions of type: A+A->AA
stop = tf - eps(tf);                    % simulation stop time
nOut = nargout;

% stochastic part:    

if runs < 2                                        % single run
    [yout{1:2}] = gillespie;                       % run gillispie
else                                               % multiple runs
    [tg,ng] = deal(cell(1,runs));
    for i = 1:runs, [tg{i},ng{i}] = gillespie; end % run gillispie
    yout = {tg,ng};      
    if nOut>2
        tt = unique([tg{:}]);                      % record times
        numeltt = numel(tt);                       % record populations
        nn = zeros(numel(n0), numeltt, runs);
        for i = 1:numeltt
            for j = 1:runs
                id = find(tg{j} <= tt(i),1,'last');
                nn(:,i,j) = ng{j}(:,id);
            end
        end
        yout(3:4) = {tt,mean(nn,3)};               % append population mean 
    end
end

% Output:

varargout = yout;

%   gillespie algorithm (direct method):

    function [tt,nn] = gillespie
        t = 0;                              % initial time
        n = n0;                             % initial population
        tt = [];
        nn = [];
        rand('state', sum(100*clock));      % reset random number generator
        while t <= stop
            tt = [tt t];                    % record time
            nn = [nn n];                    % record population            
            m = n(:,oneszk);                % replicate n : size(m,2) = size(k)
            b = double(~l);
            b(i1) = m(i1);                  % reactions of type: A->B, A+B->AB
            b(i2) = m(i2).*(m(i2)-1)/2;     % reactions of type: A+A->AA
            a = c.*prod(b);                 % propensity
            astr = sum(a); 
            if ~astr, break, end            % substrate utilised
            tau = -1/astr*log(rand);        % time to next reaction
            u = find(cumsum(a)>=astr*rand,1);% index of next reaction
            n = n + d(:,u);                 % update population
            t = t + tau;                    % update time
        end
        tt = [tt tf];
        nn = [nn nn(:,end)];
    end

end