Stochastic Gillespie-type simulations using Julia
Switch branches/tags
Nothing to show
Clone or download
Latest commit f65e05e Sep 27, 2017
Permalink
Type Name Latest commit message Commit time
Failed to load latest commit information.
benchmarks Added benchmark file Jul 28, 2016
docs Fixed examples for v0.5 Jul 29, 2016
examples Update sir_tjm.jl Aug 14, 2017
notebooks Fix minor error in notebook Apr 3, 2017
paper Update paper.bib Jul 30, 2016
src Merge branch 'master' into deprecation_utf8string Sep 27, 2017
test Fixing tests Mar 28, 2017
.gitignore More Documenter tweaks Jul 22, 2016
.travis.yml Update .travis.yml Aug 11, 2017
LICENSE.md Added docstrings Feb 15, 2015
README.md Added true jump method Aug 11, 2017
REQUIRE Added true jump method Aug 11, 2017
codemeta.json Update codemeta.json Jul 30, 2016
sir.png Updated README Aug 13, 2015

README.md

Gillespie.jl

Build Status Coverage Status Gillespie Gillespie Stories in Ready status DOI

Throughput Graph

Statement of need

This is an implementation of Gillespie's direct method as well as uniformization/Jensen's method for performing stochastic simulations, which are widely used in many fields, including systems biology and epidemiology. It borrows the basic interface (although none of the code) from the R library GillespieSSA by Mario Pineda-Krch, although Gillespie.jl only implements exact methods at present, whereas GillespieSSA also includes tau-leaping, etc.. It is intended to offer performance on par with hand-coded C code; please file an issue if you find an example that is significantly slower (2 to 5 times) than C.

Installation

The stable release of Gillespie.jl can be installed from the Julia REPL using the following command.

Pkg.add("Gillespie")

The development version from this repository can be installed as follows.

Pkg.clone("https://github.com/sdwfrost/Gillespie.jl")
Pkg.build("Gillespie")

Example usage

An example of a susceptible-infected-recovered (SIR) epidemiological model is as follows.

using Gillespie
using Gadfly

function F(x,parms)
  (S,I,R) = x
  (beta,gamma) = parms
  infection = beta*S*I
  recovery = gamma*I
  [infection,recovery]
end

x0 = [999,1,0]
nu = [[-1 1 0];[0 -1 1]]
parms = [0.1/1000.0,0.01]
tf = 250.0
srand(1234)

result = ssa(x0,F,nu,parms,tf)

data = ssa_data(result)

plot_theme = Theme(
    panel_fill=colorant"white",
    default_color=colorant"black"
)
p=plot(data,
    layer(x=:time,y=:x1,Geom.step,Theme(default_color=colorant"red")),
    layer(x=:time,y=:x2,Geom.step,Theme(default_color=colorant"orange")),
    layer(x=:time,y=:x3,Geom.step,Theme(default_color=colorant"blue")),
    Guide.xlabel("Time"),
    Guide.ylabel("Number"),
    Guide.title("SSA simulation"),
    Guide.manual_color_key("Subpopulation",["S","I","R"],["red","orange","blue"]),
    plot_theme
)

SIR

Julia versions of the examples used in GillespieSSA are given in the examples directory.

Jensen's method or uniformization

The development version of Gillespie.jl includes code to simulate via uniformization (a.k.a. Jensen's method); the API is the same as for the SSA, with the addition of max_rate, the maximum rate (Float64). Optionally, another argument, thin (Bool), can be set to false to return all the jumps (including the fictitious ones), and saves a bit of time by pre-allocating the time vector. This code is under development at present, and may change. Time-varying rates can be accommodated by passing a rate function with three arguments, F(x,parms,t), where x is the discrete state, parms are the parameters, and t is the simulation time.

The true jump method

The development version of Gillespie.jl also includes code to simulate assuming time-varying rates via the true jump method; the API is the same as for the SSA, with the exception that the rate function must have three arguments, as described above.

Performance considerations

Passing functions as arguments in Julia v0.4 incurs a performance penalty. One can circumvent this by passing an immutable object, with call overloaded, as follows.

immutable G; end
call(::Type{G},x,parms) = F(x,parms)

An example of this approach is given here. This is the default behaviour in v0.5 and above.

Benchmarks

The speed of an SIR model in Gillespie.jl was compared to:

  • A version using the R package GillespieSSA
  • Handcoded versions of the SIR model in Julia, R, and Rcpp

1000 simulations were performed, and the time per simulation computed (lower is better). Benchmarks were run on a Mac Pro (Late 2013), with 3 Ghz 8-core Intel Xeon E3, 64GB 1866 Mhz RAM, running OSX v 10.11.3 (El Capitan), using Julia v0.4.5 and R v.3.3. Jupyter notebooks for Julia and R with the code and benchmarks are available as gists. A plain Julia file is also provided in the benchmarks subdirectory for ease of benchmarking locally.

Implementation Time per simulation (ms)
R (GillespieSSA) 894.25
R (handcoded) 1087.94
Rcpp (handcoded) 1.31
Julia (Gillespie.jl) 3.99
Julia (Gillespie.jl, passing object) 1.78
Julia (handcoded) 1.20

(smaller is better)

Julia performance for Gillespie.jl is much better than GillespieSSA, and close to a handcoded version in Julia (which is itself comparable to Rcpp); as compiler performance improves, the gap in performance should narrow.

Future work

Gillespie.jl is under development, and pull requests are welcome. Future enhancements include:

  • Constrained simulations (where events are forced to occur at specific times)
  • Discrete time simulation

Citation

If you use Gillespie.jl in a publication, please cite the following.

  • Frost, Simon D.W. (2016) Gillespie.jl: Stochastic Simulation Algorithm in Julia. Journal of Open Source Software 1(3) doi:0.21105/joss.00042

A Bibtex entry can be found here.