Gillespie simulations: don't deepcopy(JumpProblem)

Specific Domains Modelling & Simulations
1

Example from the docs. It seems that deepcopying a JumpProblem object causes each copy to be solved with the same seed. Is this known?

using DifferentialEquations
using Catalyst
sir_model = @reaction_network begin
    β, S + I --> 2I
    ν, I --> R
end β ν

p = (0.1 / 1000, 0.01)
u₀ = [999, 1, 0]
tspan = (0.0, 250.0)
prob = DiscreteProblem(sir_model, u₀, tspan, p)
jump_prob = JumpProblem(sir_model, prob, Direct())

for _ in 1:5
    sol = solve(jump_prob, SSAStepper())
    print(sol.u[end])
end
# [0, 194, 806][0, 175, 825][1, 174, 825][0, 190, 810][0, 214, 786]
for _ in 1:5
    sol = solve(deepcopy(jump_prob), SSAStepper())
    print(sol.u[end])
end
# [1, 196, 803][1, 196, 803][1, 196, 803][1, 196, 803][1, 196, 803]
2

This is probably because pre Julia v1.7 releases use an RNG object per jump problem since the global RNG was not thread-safe.

@isaacsas is this fixed in later versions by using the new global Julia RNG?

3

Yeah, the second loop gives different results on 1.8.

@levasco what Julia version are you on?

4

Also, why do you want to deepcopy the JumpProblem? That is potentially a very expensive operation since it will deepcopy all the underlying SSA data structures and data too.

5

Thanks for the reply. I’m on 1.7.0 indeed.

As for the why, often out of habit/precaution as I’ve been surprised by Julia’s pass by reference in the past. But sometimes it’s actually needed e.g. if I update/store quantities in prob.p during integration using callbacks, I can’t reuse the same object for a new simulation.

6

Have you looked at Jump Problems · DifferentialEquations.jl? Generally I’d suggest using remake over deepcopying yourself if possible. (And if it doesn’t work right please do let us know / open a DiffEqJump issue!)

7

Also, if you are changing parameters during a simulation, make sure to look at this as you need to call a special function to ensure MassActionJumps, which Catalyst generates, are updated appropriately.

Basically, you need to always call reset_aggregated_jumps!(integrator) after changing p during a simulation to ensure the SSAs update for the changes. Otherwise they may have internal data structures that are no longer valid. Changing p between simulations you can just use remake.

8

It is strange you are seeing issues on 1.7 as that should have the new random number generator that is on 1.8 too. What DiffEqJump version are you using?

9

Just learned about remake(), seems like a much better way indeed. And thanks for the heads up regarding reset_aggregated_jumps!(), I had only used callbacks with classical ODEs so far but it was bound to bite me at some point. Version is 8.0, updating now to see if it still happens.

EDIT: Updating to 8.3 fixed it (didn’t update Julia, just the pkg)

1 Like
10

Perfect!

1 Like