SciML / JumpProcesses.jl

Build and simulate jump equations like Gillespie simulations and jump diffusions with constant and state-dependent rates and mix with differential equations and scientific machine learning (SciML)
https://docs.sciml.ai/JumpProcesses/stable/
Other
135 stars 35 forks source link

EnsembleProblem of JumpProblem - All solutions are identical #184

Open TorkelE opened 3 years ago

TorkelE commented 3 years ago

Here's a minimal sample. I make a similar EnsembleProblem from both and an JumpProblem and an SDEProblem. When I solve them the SDE solutions turn out different, but the Jump ones are all identical.

using Catalyst, DifferentialEquations, Plots

rn = @reaction_network begin
    (k1,k2), X <--> Y
end k1 k2

dprob = DiscreteProblem(rn, [10,10], (0.,10.), [1.,2.])
jprob = JumpProblem(rn, dprob, Direct())    
ejprob = EnsembleProblem(jprob,prob_func=(p,i,r)->remake(p))
jsol = solve(ejprob,SSAStepper(),trajectories=4)

sprob = SDEProblem(rn, [10.,10.], (0.,10.), [1.,2.])
esprob = EnsembleProblem(sprob,prob_func=(p,i,r)->p)
ssol = solve(esprob,ImplicitEM(),trajectories=4)

plot(plot(jsol),plot(ssol),size=(900,250))

image

This is the output of Pkg.status()

  [479239e8] Catalyst v6.13.0
  [0c46a032] DifferentialEquations v6.17.1
  [91a5bcdd] Plots v1.16.6
isaacsas commented 3 years ago

If it helps in debugging this, the problem seems to go away if one doesn't use prob_func.

isaacsas commented 3 years ago

@TorkelE this seems to fix it for me, does it work for you?

ejprob = EnsembleProblem(jprob,prob_func=(p,i,r)->remake(p), safetycopy=false)

Looks like the problem is JumpProblems shouldn't make a copy of the problem unless using threads to get a continual sequence of random numbers.

isaacsas commented 3 years ago

I think this should be fixed by https://github.com/SciML/SciMLBase.jl/pull/72, but I’ll leave this open till that propagates through the ecosystem.

TorkelE commented 3 years ago

Sounds good, in the meantime

ejprob = EnsembleProblem(jprob,prob_func=(p,i,r)->remake(p), safetycopy=false)

also worked