Problems for VariableRateJumps

[isaacsas]

We need to handle problem setup for JumpSystems with VariableRate jumps. (i.e. need ODESystem and SDESystem dispatches that can take JumpSystems and construct "empty" systems). Right now only DiscreteProblems where we enforce the jumps are non-variable will work.

[ChrisRackauckas]

need ODESystem and SDESystem dispatches that can take JumpSystems and construct "empty" systems

Could you explain that part a little bit more?

[isaacsas]

For

using ModelingToolkit, DiffEqJump, OrdinaryDiffEq
@parameters β γ t
@variables S I R
rate₁   = β*S*I
affect₁ = [S ~ S - 1, I ~ I + 1]
rate₂   = γ*I*exp(-t/20)
affect₂ = [I ~ I - 1, R ~ R + 1]
j₁      = ConstantRateJump(rate₁,affect₁)
j₂      = VariableRateJump(rate₂,affect₂)
js      = JumpSystem([j₁,j₂], t, [S,I,R], [β,γ])
u₀ = [999,1,0]; p = (0.1/1000,0.01); tspan = (0.,250.)
u₀map = [S => 999, I => 1, R => 0]
parammap = [β => .1/1000, γ => .01]

How does one run a pure SSA simulation now?
One might think to do

oprob = ODEProblem((du,u,p,t) -> du .= u, u₀, tspan, p)

but then

jprob = JumpProblem(js, oprob, Direct())
sol = solve(jprob, Tsit5())

fails. Similarly though one can't do

oprob = ODEProblem(js, u₀map, tspan, parammap)

or

oprob = ODEProblem((du,u,p,t) -> du.=u, u₀map, tspan, parammap)
jprob = JumpProblem(js, oprob, Direct())
sol = solve(jprob, Tsit5())

It seems like if the input is through arrays of pairs we need a way to make dummy problems that handle the pair appropriately.

[ChrisRackauckas]

oprob = ODEProblem((du,u,p,t) -> du .= 0, u₀, tspan, p) would be needed. But yeah, we need to think about this a bit. We want JumpSystem to allow for nesting of JumpSystem for component-wise building up a system, but JumpSystem also needs to be tied to another AbstractSystem, by default a DiscreteSystem, but in general could be an ODE/SDE.

[ChrisRackauckas]

It doesn't need a whole system though, since it should just be the equations using the same variable/parameter list.

[kaandocal]

Hey, I'm currently trying to to exactly this (using the SSA with VariableRateJumps) and got stuck right here. I've tried the approaches @isaacsas pointed out, to no avail, and I haven't found anything in the documentation despite a lot of searching (the closes I got was jump-diffusion equations). Is there any way to get this to work with the current version? I'm happy to cobble together a case-specific solution as I only need it for one system, but I'm not sure how to do that... Any guidance would be much appreciated!

[isaacsas]

I think this still needs an updated version of the JumpProblem call here.

[isaacsas]

I’ll try to take a look on Thursday.

[isaacsas]

@ChrisRackauckas I think this issue is related to the ODE solve function not working when the initial values are integers.
This doesn't work

using ModelingToolkit, DiffEqJump, OrdinaryDiffEq
@parameters β γ t
@variables S I R
rate₁   = β*S*I
affect₁ = [S ~ S - 1, I ~ I + 1]
rate₂   = γ*I*exp(-t/20)
affect₂ = [I ~ I - 1, R ~ R + 1]
j₁      = ConstantRateJump(rate₁,affect₁)
j₂      = VariableRateJump(rate₂,affect₂)
js      = JumpSystem([j₁,j₂], t, [S,I,R], [β,γ])
u₀ = [999,1,0]; p = (0.1/1000,0.01); tspan = (0.,250.)
u₀map = [S => 999, I => 1, R => 0]
parammap = [β => .1/1000, γ => .01]
oprob = ODEProblem((du,u,p,t) -> du .= 0, u₀, tspan, p)
jprob = JumpProblem(js, oprob, Direct())
sol = solve(jprob, Tsit5())

with an error during initialize! at the solve call, while this seems to go through and run

using ModelingToolkit, DiffEqJump, OrdinaryDiffEq
@parameters β γ t
@variables S I R
rate₁   = β*S*I
affect₁ = [S ~ S - 1, I ~ I + 1]
rate₂   = γ*I*exp(-t/20)
affect₂ = [I ~ I - 1, R ~ R + 1]
j₁      = ConstantRateJump(rate₁,affect₁)
j₂      = VariableRateJump(rate₂,affect₂)
js      = JumpSystem([j₁,j₂], t, [S,I,R], [β,γ])
u₀ = [999.0,1.0,0.0]; p = (0.1/1000,0.01); tspan = (0.,250.)
u₀map = [S => 999.0, I => 1.0, R => 0.0]
parammap = [β => .1/1000, γ => .01]
oprob = ODEProblem((du,u,p,t) -> du .= 0, u₀, tspan, p)
jprob = JumpProblem(js, oprob, Direct())
sol = solve(jprob, Tsit5())

@kaandocal can you try something like this using floating point initial values for your problem?

[isaacsas]

Maybe we need a dedicated VRSSAStepper for pure jump problems that may have continuous components?

[ChrisRackauckas]

Is this something that can be specialized on? If the rate equations are continuously dependent and don't have an analytical solution, then you need an ODE solver.

[isaacsas]

Maybe some type of default ODEFunction we use for such problems and the ODE timesteppers can detect?