# Increase perfomance for Gillespie Algorithm
**URL:** https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369
**Category:** Performance
**Created:** [December 1, 2021, 12:40pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369 "2021-12-01T12:40:03Z")
**Posts on this page:** 20
**Page:** 1
### Author:  [@Luis](https://discourse.julialang.org/u/Luis)
#### Post date: [December 1, 2021, 12:40pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/1 "2021-12-01T12:40:03Z")
I am relatively new to Julia and trying to implement a version of Gillespies Algorithm that can be used in population genetics.
I tried really hard to make it as fast and light weighted as possible, incorporating most of the performance tips from the Julia documentation. But despite my effort I seemed to achive the exact opposite: A very slow and memory gulping algorithm.
Unfortunately I was not even able to isolate the problem. Therefore I am very happy if anyone is interested in investing the time to look over my code and point out the memory and speed killer I implemented. Or if anyone can hand me some tools or literature that enables me to isolate the problem myselfe.
You can find all the code including some working (but slow) examples in the following git repository:
[https://github.com/roccminton/Gillespie\_Algorithm](https://github.com/roccminton/Gillespie_Algorithm)
Many thanks in advance for any help!
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [December 1, 2021, 12:41pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/2 "2021-12-01T12:41:18Z")
Unfortunately, the link leads to a 404 page
---
### Author:  [@Luis](https://discourse.julialang.org/u/Luis)
#### Post date: [December 1, 2021, 12:47pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/3 "2021-12-01T12:47:07Z")
Posting a link to a private repsoitory is not very clever. Sorry for that. I changed the visibility to public. Does the link work now.
---
### Author:  [@josuagrw](https://discourse.julialang.org/u/josuagrw)
#### Post date: [December 1, 2021, 3:04pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/4 "2021-12-01T15:04:07Z")
Yes it works now π
---
### Author:  [@Oscar\_Smith](https://discourse.julialang.org/u/Oscar_Smith)
#### Post date: [December 1, 2021, 3:34pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/5 "2021-12-01T15:34:42Z")
Can you profile this and figure out what is taking time? Debugging performance of MWEs is much easier than a full repo.
---
### Author:  [@goerch](https://discourse.julialang.org/u/goerch)
#### Post date: [December 1, 2021, 3:46pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/6 "2021-12-01T15:46:29Z")
The only example working for me on 1.8.0 seems to be `Examples\RandomMating.jl`. Is that what you are after?
---
### Author:  [@Luis](https://discourse.julialang.org/u/Luis)
#### Post date: [December 1, 2021, 4:33pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/7 "2021-12-01T16:33:56Z")
I fixed the bugs making the other examples not work, but indeed the `Examples\RandomMating.jl` is the most interesting example for me.
---
### Author:  [@Luis](https://discourse.julialang.org/u/Luis)
#### Post date: [December 1, 2021, 4:37pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/8 "2021-12-01T16:37:08Z")
Try to profile for example the function `DiploidModel.rungillespie(t,x0,model_parameter)` in the file `Examples\RandomMating.jl` (in line 28). I did as well but was unable to get helpful information out of the print output.
---
### Author:  [@goerch](https://discourse.julialang.org/u/goerch)
#### Post date: [December 1, 2021, 4:49pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/9 "2021-12-01T16:49:18Z")
After disabling `ProgressMeter` and `Plots` and putting the test into a function like so
```julia
function test()
K = 1000
b = 1.0
d = 0.9
c = (b-d)/(2K)
dni = 0.1
N = 10
model_parameter = (
birth = b,
death = d,
competition = c,
ΞΌ = dni,
Nloci = N
)
t = 0:100
x0 = DiploidModel.generatehealthypopulation(K,N)
#execute the simulation
history = DiploidModel.rungillespie(t,x0,model_parameter)
#plot simulation
#PlotFromDicts.plotmutationloadandprevalence(history)
end
```
JET reports for `@report_opt test()`
```julia
βββββ 3 possible errors found βββββ
β @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\Examples\RandomMating.jl:30 Main.DiploidModel.rungillespie(t, x0, model_parameter)
ββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\DiploidModel.jl:35 Main.DiploidModel.setupparameter(model_parameter, nβ)
βββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\DiploidModel.jl:51 Main.DiploidModel.npropagable(n0)
ββββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\DiploidModel.jl:95 Main.DiploidModel.sum(Base.Generator(#3, ps))
βββββ @ reduce.jl:549 Base.#sum#262(Base.pairs(Core.NamedTuple()), #self#, a)
ββββββ @ reduce.jl:549 Base.sum(Base.identity, a)
βββββββ @ reduce.jl:520 Base.#sum#261(Base.pairs(Core.NamedTuple()), #self#, f, a)
ββββββββ @ reduce.jl:520 Base.mapreduce(f, Base.add_sum, a)
βββββββββ @ reduce.jl:294 Base.#mapreduce#258(Base.pairs(Core.NamedTuple()), #self#, f, op, itr)
ββββββββββ @ reduce.jl:294 Base.mapfoldl(f, op, itr)
βββββββββββ @ reduce.jl:162 Base.#mapfoldl#254(Base._InitialValue(), #self#, f, op, itr)
ββββββββββββ @ reduce.jl:162 Base.mapfoldl_impl(f, op, init, itr)
βββββββββββββ @ reduce.jl:44 Base.foldl_impl(opβ², nt, itrβ²)
ββββββββββββββ @ reduce.jl:49 Base.reduce_empty_iter(op, itr)
βββββββββββββββ @ reduce.jl:370 Base.reduce_empty_iter(op, itr, Base.IteratorEltype(itr))
ββββββββββββββββ @ reduce.jl:371 Base.reduce_empty(op, Base.eltype(itr))
βββββββββββββββββ @ reduce.jl:348 Base.mapreduce_empty(#3, $(QuoteNode(Base.BottomRF{typeof(Base.add_sum)}(Base.add_sum))), _)
βββββββββββββββββ runtime dispatch detected: Base.mapreduce_empty(#3, $(QuoteNode(Base.BottomRF{typeof(Base.add_sum)}(Base.add_sum)))::Base.BottomRF{typeof(Base.add_sum)}, _::Type{Pair{Main.DiploidModel.RandomIndividual{Int64}, Int64}})
ββββββββββββββββββββββββββββββββββ
ββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\DiploidModel.jl:35 Main.DiploidModel.Gillespie.run_gillespie!(time, nβ, Main.DiploidModel.setupparameter(model_parameter, nβ), Main.DiploidModel.execute!, Main.DiploidModel.rates!, initrates, population_history)
βββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\MainFunctions.jl:17 Main.DiploidModel.Gillespie.mainiteration!(population_history, initrates, nβ, Main.DiploidModel.Gillespie.convert(Main.DiploidModel.Gillespie.Float64, Base.getindex(time, 1)), time, par, execute!, rates!)
ββββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\MainFunctions.jl:64 ct = Main.DiploidModel.Gillespie.onestep!(n0, rates, ct, step, par, ex!, r!)
βββββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\MainFunctions.jl:147 Main.DiploidModel.Gillespie.iszero(i)
ββββββ @ number.jl:42 Base.zero(x)
ββββββ runtime dispatch detected: Base.zero(x::Nothing)
ββββββββββββββββββββββ
βββββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\MainFunctions.jl:151 ex!(i, x_0, par)
ββββββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\DiploidModel.jl:194 Main.DiploidModel.death!(ps, par)
βββββββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\DiploidModel.jl:187 Base.setindex!(ps, Main.DiploidModel.-(Base.getindex(ps, Main.DiploidModel.choosefey(ps, Base.getindex(Base.getproperty(par, :popsize), 1))), Main.DiploidModel.one(Main.DiploidModel.valtype(ps))), Main.DiploidModel.choosefey(ps, Base.getindex(Base.getproperty(par, :popsize), 1)))
ββββββββ @ dict.jl:374 Base.convert(_, key0)
ββββββββ runtime dispatch detected: Base.convert(_::Type{Main.DiploidModel.RandomIndividual{Int64}}, key0::Nothing)
βββββββββββββββββββββββ
```
Does this help?
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 1, 2021, 4:55pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/10 "2021-12-01T16:55:19Z")
How does it compare to the DifferentialEquations.jl/Catalyst.jl SSAStepper?
[https://benchmarks.sciml.ai/html/Jumps/NegFeedback\_GeneExpr.html](https://benchmarks.sciml.ai/html/Jumps/NegFeedback_GeneExpr.html)
[https://catalyst.sciml.ai/dev/tutorials/using\_catalyst/](https://catalyst.sciml.ai/dev/tutorials/using_catalyst/)
It would be good to at least get a baseline.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 1, 2021, 4:56pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/11 "2021-12-01T16:56:05Z")
Or use the same baseline as these benchmarks:
[https://github.com/sdwfrost/Gillespie.jl#benchmarks](https://github.com/sdwfrost/Gillespie.jl#benchmarks)
---
### Author:  [@goerch](https://discourse.julialang.org/u/goerch)
#### Post date: [December 1, 2021, 5:06pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/12 "2021-12-01T17:06:55Z")
Iβm spotting at least one problem here:
```julia
function chooseevent(rates,total_rate)
#make it a uniform random variable in (0,total_rate)
rndm = rand(Uniform(0.0,total_rate))
#choose the rate at random
@inbounds for (index, rate) in enumerate(rates)
rndm -= rate
rndm β€ 0.0 && return index
end
end
```
`chooseevent` could return `nothing`
---
### Author:  [@isaacsas](https://discourse.julialang.org/u/isaacsas)
#### Post date: [December 1, 2021, 5:20pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/13 "2021-12-01T17:20:13Z")
@Luis can you point to any mathematical write ups of one of your models showing how they differ via the trait space from traditional Gillespie models? That might help in pointing you towards possible algorithmic improvements and/or seeing how to add your models into existing packages.
How many rates do you typically sum up each step of the algorithm to get the total rate? If it is many there are lots of possible optimized Gillespie methods one can use, typically designed around only updating individual rates that would have actually changed after the occurrence of a given jump.
---
### Author:  [@adolgert](https://discourse.julialang.org/u/adolgert)
#### Post date: [December 2, 2021, 12:10am UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/14 "2021-12-02T00:10:59Z")
Maybe I can give two pointers here, to papers that discuss what @isaacsas is describing.
- Li H, Petzold L (2006) Logarithmic direct method for discrete stochastic simulation of chemically reacting systems. Technical report. Available: [http://www.engineering.ucsb.edu/cse/Files/ldm0513.pdf](http://www.engineering.ucsb.edu/cse/Files/ldm0513.pdf). Accessed on 2012 Apr 24. - This is the standard for optimized Gillespie methods. It uses a prefix sum tree (like a Fenwick tree) to keep track of the sum of the propensities, while allowing you to modify the values of those propensities.
- [https://www.researchgate.net/publication/8379343\_Efficient\_formulation\_of\_the\_stochastic\_simulation\_algorithm\_for\_chemically\_reacting\_systems](https://www.researchgate.net/publication/8379343_Efficient_formulation_of_the_stochastic_simulation_algorithm_for_chemically_reacting_systems) - This is an overview of techniques, showing the scope of the kinds of solutions.
If you want to do a Gillespie algorithm that has an infinite number of possible reactions, with only a finite number active at any one time, then thatβs also totally possible, but it just tends to go by other names than βGillespie.β Thatβs also in the second reference, above. It entails keeping a mutable heap of the next reactions to fire and the times at which they would fire.
---
### Author:  [@Luis](https://discourse.julialang.org/u/Luis)
#### Post date: [December 2, 2021, 8:56am UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/15 "2021-12-02T08:56:20Z")
Many thanks for your effort. That is much more readable then the output I usually get!
The first error comes during the setup of the model and gets called only once, thus does not contribute too much to the problem.
Error two and three appear due to the same problem @goerch pointed out. Eventhough the function will never actually return nothing, because the handed in `total_rate` is the sum of `rates` I added another return statement at the end of the loop.
---
### Author:  [@Luis](https://discourse.julialang.org/u/Luis)
#### Post date: [December 2, 2021, 9:16am UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/16 "2021-12-02T09:16:23Z")
I am analyzing adaptive dynamics models as described in [Stochastic models for adaptive dynamics: Scaling limits and diversity](https://arxiv.org/abs/1909.02456). There usually the trait space is of infinitely large (e.g. `[-1,1]` as in the DieckmannDoeblin example) or very very big (e.g. `{0,1,2}^N` with `N >> 100` as in the DiploidModel example). And even though only a finite number of traits is alive at any time I am interested in high mutation rate regimes. There the number of different individuals alive at any time is very high.
I only came across implementations of Gillespies algorithm where there are a finite number of traits and one knows the transition rates from one trait to the other in advance (e.g. the SIR-model).
Typically I have a lot of rates to sum up to get the total rate. But unfortunately I also have to update every single rate after every single event due to the competition pressure: In the models I study individuals or die due to natural or competitive death. Every individual exceeds a small competitive pressure on every other individual. Hence if any individual exits or enters the population the death rate of every individual has to be adapted accordingly.
---
### Author:  [@Luis](https://discourse.julialang.org/u/Luis)
#### Post date: [December 3, 2021, 10:58am UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/17 "2021-12-03T10:58:11Z")
I added a SIR model to the Examples folder to compare it to the other Gillespie Algorithms around. This is the Benchmark result:
```julia
BenchmarkTools.Trial: 10000 samples with 1 evaluation.
Range (min β¦ max): 25.684 ΞΌs β¦ 2.395 ms β GC (min β¦ max): 0.00% β¦ 0.00%
Time (median): 34.669 ΞΌs β GC (median): 0.00%
Time (mean Β± Ο): 38.426 ΞΌs Β± 45.432 ΞΌs β GC (mean Β± Ο): 0.00% Β± 0.00%
βββ
ββββββ
ββββββββββ β β
ββββββββββββββββββββββββββββββββββ
βββββββββ
β
ββ
ββββ
βββββ
ββββ β
25.7 ΞΌs Histogram: log(frequency) by time 102 ΞΌs <
Memory estimate: 848 bytes, allocs estimate: 18.
```
Pleas try to do the analysis yourself, because I am not quite sure if I did it right.
---
### Author:  [@ChrisRackauckas](https://discourse.julialang.org/u/ChrisRackauckas)
#### Post date: [December 3, 2021, 11:48am UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/18 "2021-12-03T11:48:33Z")
I think this is a good clean implementation of Gillespie for small dense problems. Other implementations like the DiffEqJump ones specialize more on large problems and sparsity, while here the conditional form will only be efficient for small numbers of reactions, but would densify it.
I think one think I can see is:
[https://github.com/roccminton/Gillespie\_Algorithm/blob/master/MainFunctions.jl#L17](https://github.com/roccminton/Gillespie_Algorithm/blob/master/MainFunctions.jl#L17)
you might want to force function specialization on the higher order functions, i.e.
```julia
function run_gillespie!(time,nβ,par,execute!::F1,rates!::F2,initrates,population_history) where {F1,F2}
```
etc. going down. See if that helps out. But I think what will bite you in scaling is this:
[https://github.com/roccminton/Gillespie\_Algorithm/blob/master/Examples/SIR.jl#L22-L31](https://github.com/roccminton/Gillespie_Algorithm/blob/master/Examples/SIR.jl#L22-L31)
Itβs not necessarily bad, but as the potential `i`βs get large youβre going to have a pretty nasty conditional evaluation with bad branch prediction, so itβs a choice that will accelerate small problems at the cost of large ones.
---
### Author:  [@Luis](https://discourse.julialang.org/u/Luis)
#### Post date: [December 6, 2021, 7:39pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/19 "2021-12-06T19:39:29Z")
Thank you for the support! I will definitely check out the DiffEqJump package!
Actually for almost all my uses the `execute!` function only decides between two events (resp. `i`s), namely `birth!` and `death!`. But I belive what it slows down is the chooseevent function, when the dictionary gets big. In some applications the dictionary carries about :
[https://github.com/roccminton/Gillespie\_Algorithm/blob/8db5144e1cc07e3e8fc9a904ab65aabfe13300ed/MainFunctions.jl#L230-L241](https://github.com/roccminton/Gillespie_Algorithm/blob/8db5144e1cc07e3e8fc9a904ab65aabfe13300ed/MainFunctions.jl#L230-L241)
But I donβt know how to accelerate this. Moreover I donβt get where all the allocations come from. Especially in the `RandomMating.jl` Example. There the average dictionary has 2000-3000 entries.
---
### Author:  [@goerch](https://discourse.julialang.org/u/goerch)
#### Post date: [December 6, 2021, 7:56pm UTC](https://discourse.julialang.org/t/increase-perfomance-for-gillespie-algorithm/72369/20 "2021-12-06T19:56:29Z")
Here is the updated JET report:
```julia
βββββ 3 possible errors found βββββ
β @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\examples\randommating.jl:29 Main.DiploidModel.rungillespie(t, x0, model_parameter)
ββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\DiploidModel.jl:35 Main.DiploidModel.Gillespie.run_gillespie!(time, nβ, Main.DiploidModel.setupparameter(model_parameter, nβ), Main.DiploidModel.execute!, Main.DiploidModel.rates!, initrates, population_history)
βββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\MainFunctions.jl:17 Main.DiploidModel.Gillespie.mainiteration!(population_history, initrates, nβ, Main.DiploidModel.Gillespie.convert(Main.DiploidModel.Gillespie.Float64, Base.getindex(time, 1)), time, par, execute!, rates!)
ββββ @ C:\Users\Win10\.julia\packages\ProgressMeter\Vf8un\src\ProgressMeter.jl:934 meter#291 = ProgressMeter.Progress(ProgressMeter.length(iterable))
βββββ @ C:\Users\Win10\.julia\packages\ProgressMeter\Vf8un\src\ProgressMeter.jl:94 $(QuoteNode(ProgressMeter.var"#Progress#1#2"()))(0.1, "Progress: ", :green, ProgressMeter.stderr, ProgressMeter.nothing, ProgressMeter.BarGlyphs('|', 'β', _3, ' ', '|'), 0, 0, true, false, #self#, n)
ββββββ @ C:\Users\Win10\.julia\packages\ProgressMeter\Vf8un\src\ProgressMeter.jl:94 ProgressMeter.running_ijulia_kernel()
βββββββ @ C:\Users\Win10\.julia\packages\ProgressMeter\Vf8un\src\ProgressMeter.jl:243 Base.getproperty(ProgressMeter.Main, :IJulia)
ββββββββ @ Base.jl:31 Base.getfield(x, f)
ββββββββ variable Main.IJulia is not defined: Base.getfield(x::Module, f::Symbol)
ββββββββββββββββββββββ
ββββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\MainFunctions.jl:62 ct = Main.DiploidModel.Gillespie.onestep!(n0, rates, ct, step, par, ex!, r!)
βββββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\MainFunctions.jl:145 Main.DiploidModel.Gillespie.iszero(i)
ββββββ @ number.jl:42 Base.zero(x)
ββββββ no matching method found for call signature (Tuple{typeof(zero), Nothing}): Base.zero(x::Nothing)
ββββββββββββββββββββββ
βββββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\MainFunctions.jl:149 ex!(i, x_0, par)
ββββββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\DiploidModel.jl:194 Main.DiploidModel.death!(ps, par)
βββββββ @ C:\Users\Win10\Documents\GitHub\Gillespie_Algorithm-master\DiploidModel.jl:187 Base.setindex!(ps, Main.DiploidModel.-(Base.getindex(ps, Main.DiploidModel.choosefey(ps, Base.getindex(Base.getproperty(par, :popsize), 1))), Main.DiploidModel.one(Main.DiploidModel.valtype(ps))), Main.DiploidModel.choosefey(ps, Base.getindex(Base.getproperty(par, :popsize), 1)))
ββββββββ @ dict.jl:374 key = Base.convert(_, key0)
ββββββββ no matching method found for call signature (Tuple{typeof(convert), Type{Main.DiploidModel.RandomIndividual{Int64}}, Nothing}): key = Base.convert(_::Type{Main.DiploidModel.RandomIndividual{Int64}}, key0::Nothing)
βββββββββββββββββββββββ
```