Allow equations with JumpSystem / Support Jumps in System

[ChrisRackauckas]

Currently JumpSystem lives in a somewhat different world from the rest of the systems, but we should clean that up. Part of the clean up is to generalize it. Currently we have:

using ModelingToolkit, JumpProcesses
using ModelingToolkit: t_nounits as t

@parameters β γ
@variables S(t) I(t) R(t)
rate₁   = β*S*I
affect₁ = [S ~ S - 1, I ~ I + 1]
rate₂   = γ*I
affect₂ = [I ~ I - 1, R ~ R + 1]
j₁      = ConstantRateJump(rate₁,affect₁)
j₂      = ConstantRateJump(rate₂,affect₂)
j₃      = MassActionJump(2*β+γ, [R => 1], [S => 1, R => -1])
@named js      = JumpSystem([j₁,j₂,j₃], t, [S,I,R], [β,γ])

The issue is that unlike most equation arrays, the JumpSystem one only accepts jump types. That is unnecessarily restrictive, and in fact it limits us from expressing more complex systems like jump diffusions https://docs.sciml.ai/JumpProcesses/stable/tutorials/jump_diffusion/. For such equations, you need to define the system by a set of jumps, mixed with differential and algebraic equations. So what we really want is for that array to allow for any jumps + any equations.

Now when that occurs, the compilation needs to get a little bit more sophisticated. If there are algebraic equations, then the system is an ImplicitDiscreteSystem (rather than the current DiscreteSystem), and if they are all explicit then it should simplify to DiscreteSystem (though @vyudu this is one reason for removing the difference between the two). But secondly, if there are differential equations also included we need to do two things:

1. We need to validate that any differential variables are setup to have continuous states, i.e. they are not variables which are forced to be integers.
2. We should validate that any rate function of a ConstantRateJump does not contain a differential variable. That's a violation of the constant rate between jumps assumption. We can have a bypass allow_changing_constant_rates = true kwarg to allow it to build anyways, there are some use cases because the degree to which this introduces error depends on how much the rate changes between jumps, so it's not always a bad assumption. But by default we should try to point things in the right direction.
3. If there are any differential equations, the construction should require building an ODE/SDE (depending on if there are any Brownian variables).

In fact, one step we should add to this is instead just make a SymbolicJump(rate, affect), where it smartly makes it a constant rate jump or a variable rate jump depending on the types of values in the rate equation (i.e. make these jumps a MassActionJump or ConstantRateJump if all of the values in the rate function are not time-dependent functions and are not differential variables).

Now to implement this, it might make sense to do this by simply going straight to System. What we could do is instead just extend System to allow for taking in SymbolicJump/ConstantRateJump/MassActionJump/VariableRateJump in its equation list, and if it detects jumps it knows it needs to build a JumpProblem and it can do the required checks.

Note there are multiple purposes for this. One is that it would simplify the life of Catalyst.jl, as it would give it a single target. But secondly, this would also then make a very good basis on which to experiment with algorithms that move between jumps and differential equations. This has been a long time goal of me and @isaacsas, but it's always been difficult to express such algorithms. This foundation would be the perfect place to build such methods: methods which switch from discrete jumps to ODEs would just be stable transformations on such a System, where new parameters are introduced as switches and ODE equations would be introduced as well, along with callbacks.

Note that to get this out the door, it's completely fine to make compilation with differential equations in there just error. The mixture will take some work, getting the interface correct is step 1.

@TorkelE @isaacsas

[isaacsas]

Just to clarify, JumpSystems support coupled ODEs already.

[isaacsas]

(Though not in any sophisticated way that does auto analysis.)

[isaacsas]

It would be good to keep a low-level interface where the jumps can be manually specified though. What kind of overhead will first generating symbolic rates and affects for the B cell model, and then analyzing and converting them all into MassActionJumps add? In contrast, it is trivial for Catalyst to know how they should be represented and generate them as such (and it avoids any symbolic overhead since all or almost all reactions end up as MassActionJumps, which have no internal symbolic expressions). We'd also like to allow Catalyst users the option to pick the jump type, in case they for some reason want to move a jump to a different representation (which people have asked about in the past, though I can't remember why).

[ChrisRackauckas]

It would be good to keep a low-level interface where the jumps can be manually specified though.

To be clear, the plan of the System type is for it to be layers of transformations which increasingly specialize the information. So step 1 would be to do things like:

1. Identify any equations with brownian variables and separate them out to make noise equations
2. Identify any SymbolicJumps to concrete jumps
3. ...

If the user wants to bypass that detection, they could just specify the jump type. It's like today, you can bypass the detection of Brownians by directly calling SDESystem with the noise equations. However, our plan down the line is that SDESystem should just be internal, still documented but most users shouldn't have to know how to interpret their SDEs into a nondiagonal noise structure, just write the equations.

What kind of overhead will first generating symbolic rates and affects for the B cell model, and then analyzing and converting them all into MassActionJumps add?

This stuff is linear and small. There are things to worry about but this should be the bottom of the list.