%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Stochastic simulation engine using Gillespie algorithm
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

function sim_log = ssas_gillespie(lhs, rhs, c, y0, T_MAX, dt)

    % Normal nchoosek() can't handle the case where n < k
    nchoosek2 = @(n,k)(prod((n-k+1:n))/prod(1:k)); 

    % Set up initial conditions for the simulation
    t = 0;
    Y = y0;

    sim_log = zeros(T_MAX/dt, 1 + size(Y, 1));
    sim_log(1,:) = [t, Y'];
    log_index = 2;

    % Run the simulation
    %tic
    while t < T_MAX
        % Check for extinction
        if all(Y == 0)
            disp('Reactants exhausted!');
            break;
        end

        % Compute the reaction hazards
        H = c.*prod(arrayfun(nchoosek2, repmat(Y',size(lhs, 1),1), lhs), 2);
        H_sums = cumsum(H);

        % Jump to the time of the next reaction
        t = t + exprnd(1/H_sums(end));

        % Determine which reaction occurred
        H_index = find(H_sums > H_sums(end)*rand(), 1);

        % Update the reactant quantities appropriately
        Y = Y + rhs(H_index,:)' - lhs(H_index,:)';

        % Log the new system state, if appropriate
        if (t - sim_log(log_index-1,1)) > dt
            sim_log(log_index,:) = [t, Y'];
            log_index = log_index + 1;
        end
    end
    %toc

    % Remove any leftover preallocated data points
    sim_log(log_index:end,:) = [];
end