direct.py

# A simple implementation of the Gillespie's direct method
import numpy
from functools import partial

def step(t, s, events):
    a = numpy.array([ev[0](s) for ev in events])
    atot = a.sum()
    if atot == 0.0:
        return numpy.inf, s

    r1 = numpy.random.uniform()
    r2 = numpy.random.uniform()
    tau = -numpy.log(r1) / atot

    j, a0 = 0, a[0]
    while a0 <= r2 * atot:
        j += 1
        a0 += a[j]
    return (t + tau, events[j][1](s))

def propensity(s, k, v, stoich):
    a = k * v
    for i, coef in enumerate(stoich):
        if coef < 0:
            a *= (s[i] / v) ** (-coef)
    return a

def fire(s0, stoich):
    s1 = s0.copy()
    for i, coef in enumerate(stoich):
        s1[i] += coef
    return s1


if __name__ == '__main__':
    events = []
    events.append((partial(propensity, k=1.0/30.0, v=1.0, stoich=(-1, -1, +1)),
                   partial(fire, stoich=(-1, -1, +1))))
    events.append((partial(propensity, k=1.0, v=1.0, stoich=(+1, +1, -1)),
                   partial(fire, stoich=(+1, +1, -1))))

    t, s = 0.0, numpy.array([60, 60, 0])
    for _ in range(120):
        print("{0:g} {1:s}".format(t, " ".join(["{0:g}".format(v) for v in s])))
        t, s = step(t, s, events)