Enable creation of JumpSystems from ReactionSystems

[TorkelE]

• convert(JumpSystem,rs::ReactionSystem) should now work.
• Type of Jump (Variable Rate, ConstantRate, MassAction) should now be correctly identified.

[TorkelE]

Not done, but should be in a workable state now where one can try out the code. I haven't started on tests yet.

Some irregularities:

• I had to create a separate jumpratelaw(), instead of using the oderatelaw(). The reason for this was that the later generated rates on the form A*X(t)*Y(t), these generated MethodError: objects of type Int64 are not callable when trying to solve the problem. If I replaced it instead with A*X*Y things looked like in the normal JumpSystems, and the solving worked. I am not sure what exactly is going on here, and I might have taken the wrong approach (and if not, there is likely a prettier way to solve it, but I didn't want to intrude into other function just yet).
• Relating to this, at several places I had to go from a Variable to an Operation, and so created a var2op() function (didn't find one already present). Again, I am very uncertain whenever I here have used the best approach but will leave it to someone more involved to decide.

[TorkelE]

A simple example of usage

using DiffEqJump, ModelingToolkit
@parameters t p d
@variables X(t)
rxs = [Reaction(p, nothing, [X]),
       Reaction(d, [X], nothing)
       ]
rs = ReactionSystem(rxs,t,[X],[p,d])

js = convert(JumpSystem,rs)
dprob = DiscreteProblem(js, [X => 10], (0,10.), [p =>20., d=>1.])
jprob = JumpProblem(js, dprob, Direct())
sol = solve(jprob, SSAStepper())
using Plots; plot(sol)

[ChrisRackauckas]

I had to create a separate jumpratelaw(), instead of using the oderatelaw(). The reason for this was that the later generated rates on the form AX(t)Y(t), these generated MethodError: objects of type Int64 are not callable when trying to solve the problem. If I replaced it instead with AXY things looked like in the normal JumpSystems, and the solving worked. I am not sure what exactly is going on here, and I might have taken the wrong approach (and if not, there is likely a prettier way to solve it, but I didn't want to intrude into other function just yet).

It would be nice to figure out why that's necessary. Something probably isn't using the right expression generation from variable setup: https://github.com/SciML/ModelingToolkit.jl/blob/master/src/systems/diffeqs/abstractodesystem.jl#L23-L37

[ChrisRackauckas]

I think we can make this better, but let's iteratively improve it.

[TorkelE]

There seems to be something fishy going on now. When I try the code again

using DiffEqJump, ModelingToolkit
@parameters t p d
@variables X(t)
rxs = [Reaction(p, nothing, [X]),
       Reaction(d, [X], nothing)
       ]
rs = ReactionSystem(rxs,t,[X],[p,d])

js = convert(JumpSystem,rs)
dprob = DiscreteProblem(js, [X => 10], (0,10.), [p =>20., d=>1.])
jprob = JumpProblem(js, dprob, Direct())
sol = solve(jprob, SSAStepper())
using Plots; plot(sol)

I get an

MethodError: no method matching assemble_maj(::JumpSystem, ::MassActionJump{Operation,Array{Pair{Int64,Int64},1},Array{Pair{Operation,Int64},1}}, ::Dict{Variable{Number},Int64}, ::Array{Pair{Variable{ModelingToolkit.Parameter{Number}},Float64},1})

when trying to create the JumpProblem. If I go back to the branch, before it was merged though, it is all fine. I don't know if it was something that was added between the fork and the merge which changed things up.

In general, there seems to be a problem creating creating JumpProblems from MassActionJumps now.

[ChrisRackauckas]

Yeah, something happened with the substitute stuff. @shashi could we get a SymbolicUtils.jl tag in the next hour or so?

[ChrisRackauckas]

It's UnPack.jl bounds. Things need to propagate through.

[isaacsas]

@ChrisRackauckas The JumpSystem crash seems to still be an issue.

Also, @TorkelE the SSA rate law should be different from the ODE rate law; the former should use binomials (see mass_rate_SSA in DEBio or the docs). (If I've misunderstood your code I apologize.)

[TorkelE]

Oh yeah, totally forgot about those.

[TorkelE]

Ok, I will wait with correcting until the JumpSystems run, to avoid the risk of adding errors when errors already present.

Also, @isaacsas for MassActionJumps, what should the rate be (the actual rate, or just the constant,)? That is, for A, 2X + Y --> ... the actual rate is A*binomial(X,2)*Y, but am I correct that MassActionJump deduces that by itself, and only A should be passed into it?

[isaacsas]

It should be the parameter that represents the rate constant, or any expression involving just parameters (i.e. 5*k1*k2+k3 would be fine). It just can't involve t or any other species variables. This should match the current DEBio functionality too.

[TorkelE]

Got, thanks.
(Currently when I created a MassActionJump with the input rate A for the reaction A, 2X --> the rate recorded in the structure seemed to be A/2, which got be start thinking about it, and I wanted to make sure.)

[isaacsas]

That’s what it should be!

[isaacsas]

The full rate expression will be AX(X-1)/2.

[isaacsas]

@TorkelE Just in case my answer wasn't completely clear, since I think I missed what you were actually getting at, the MassActionJump constructor, by default, adds the combinatoric scaling to the rate, see:
https://github.com/SciML/DiffEqJump.jl/blob/6f41387f472999deb54a1a6695847f19a90f11df/src/jumps.jl#L70

This is set by the constructor's scale_rates kwarg, which is true by default. I think this is the version we use through DEBio too.