Discrete events with repeated dosage fails for JumpSystems (fine ofr e.g. ODEs)

[TorkelE]

Initially an issue from here: SciML/ModelingToolkit.jl#2611

According to Chris "It's an issue with SSAStepper not having integrator.opts.tstops for the callback"

Example:

using ModelingToolkit, JumpProcesses

@parameters p d
@variables t X(t)
rate1   = p
rate2   = X*d
affect1 = [X ~ X + 1]
affect2 = [X ~ X - 1]
j1 = ConstantRateJump(rate1, affect1)
j2 = ConstantRateJump(rate2, affect2)

# Works.
@mtkbuild js = JumpSystem([j1, j2], t, [X], [p,d])
dprob = DiscreteProblem(js, [X => 15], (0.0, 10.), [p => 2.0, d => 0.5])
jprob = JumpProblem(js, dprob, Direct())
sol = solve(jprob, SSAStepper())

# Discrete events with dosage works.
discrete_events = [2.0, 5.0] => [X ~ X + 10]
@mtkbuild js = JumpSystem([j1, j2], t, [X], [p,d]; discrete_events)
dprob = DiscreteProblem(js, [X => 15], (0.0, 10.), [p => 2.0, d => 0.5])
jprob = JumpProblem(js, dprob, Direct())
sol = solve(jprob, SSAStepper())

# Discrete events with repeated dosage fails.
discrete_events = 2.0 => [X ~ X + 10]
@mtkbuild js = JumpSystem([j1, j2], t, [X], [p,d]; discrete_events)
dprob = DiscreteProblem(js, [X => 15], (0.0, 10.), [p => 2.0, d => 0.5])
jprob = JumpProblem(js, dprob, Direct())
sol = solve(jprob, SSAStepper()) # ERROR: type NamedTuple has no field tstops

# Discrete events with works.
discrete_events = X == 10 => [X ~ X + 10]
@mtkbuild js = JumpSystem([j1, j2], t, [X], [p,d]; discrete_events)
dprob = DiscreteProblem(js, [X => 15], (0.0, 10.), [p => 2.0, d => 0.5])
jprob = JumpProblem(js, dprob, Direct())
sol = solve(jprob, SSAStepper())

For ODEs it works fine though:

# This works fine.
eqs = [Differential(t)(X) ~ p - d*X]
discrete_events = 2.0 => [X ~ X + 10]
@mtkbuild osys = ODESystem(eqs, t; discrete_events)
oprob = ODEProblem(osys, [X => 15], (0.0, 10.), [p => 2.0, d => 0.5])
solve(oprob, Tsit5())

[isaacsas]

@ChrisRackauckas sounds like maybe we should some additional API functions for querying tstops? I think SSAIntegrator going back to your original implementation has always stored them directly instead of inside the opts tuple.

[ChrisRackauckas]

Yeah it probably needs some higher level query function in DiffEqBase.

[isaacsas]

And the DiffEqCallbacks needs updates.

[ChrisRackauckas]

Yes

[isaacsas]

Fixing this is a bit more complicated as it seems DiffEqCallbacks currently relies on tstops being a specific binary heap implementation (which we don't use in JumpProcesses):

https://github.com/SciML/DiffEqCallbacks.jl/blob/1a02294bff7debdf8ecdca62b6a4cbfd52c10903/src/iterative_and_periodic.jl#L47

So it seems there also needs to be an API function to get an array representation from the tstops datastructure.

[isaacsas]

In fact, both OrdinaryDiffEq and DiffEqCallbacks seems to rely on that internal field from DataStructures.jl's heap implementation, which seems like a not so great idea in general.

[ChrisRackauckas]

I don't disagree. Someone just needs to generalize that. The issue is that the heap implementation does not expose the operation that we want here.