Common model type generated by DSLs

[TorkelE]

Haven't been looking that much into ModellingToolkit for a while (my own programs haven't required it. I have also been under the impression that MT has been under rapid development, so figured I'd wait to read up on it until it started to settle down a bit). This is mostly a preemptive excuse in case what I am about to say shows ignorance of what MT is all about (but I have tried to read through what documentation is available).

Now seems like a good time to look into MT (especially if I can help with the transition of DiffEqBiological to use it). The fact that current epidemiological events have forced me to stay put in my room is also a contributing factor...

So I am under the impression that one of the functions of MT, as an intermediate representation, is to provide some kind of common interface for various modelling DSLs. So instead of each one generating all of the required functions one lets MT do it. The DSLs would try to use MT stuff to as great an extent as possible.

What I am wondering, should there be (or is there already?) some kind of common model type, which is defined by MT and generated by the DSLs? Currently, there are several functions provided by MT to help generate the fields in the types generated by the DSLs (like the reaction_network type from DiffEqBiological). There does not seem to be, however, a framework to united these different models types generated by the DSLs (as was my initial understanding).

E.g. the DiffEqBiological DSL would generate a reaction_network object, which would be a subtype of the model type implemented by MT. model would contain several fields such as ODESystem (and ODEFunction?), SDESystem etc. These are defined by MT. Next, reaction_network would contain these fields, plus some additional unique for this kind of model. Now functions like ODEProblem could take a model object as input (and DiffEqBiological would not need to define this operation). More importantly, if someone in the future makes some new cool package, working on DifferentialEquation types of models, they could make it take model as input, and separate implementations for each DSL would not be needed.

Possibly, model would need to contain fields not needed (or even, conceivable) for all DSLs. Some models might need SDEs (so these would be natural inclusion in model), but maybe in some cases, these make no sense. Thus any method acting on model would first need a check that the fields required are defined, and if not, return an error. Maybe even each DSL should define directly which fields makes no sense under that framework, and cannot be filled.

This is all rather long and rambling, and maybe not particularly useful. I am, however, still trying to understand exactly what MT is, and when and how to use it. Jacking DiffEqBiological into MT makes sense. But if MT is only an easy way to define SDESystems for a DSL, given that DiffEqBiological already has it implemented there is less to be gained by using MTs version instead (but more to be won by future DSLs which do not need to do the re-implementation). There must be something more to it?

The advantages I would hope MT could offer an implemented package like DiffEqBiological would be:

• Simpler compatibility and integration with other packages, since it could share structure with other models (as described above).
• More flexibility. E.g. one could create an initial model with DiffEqBiological, but then use MT to modify it in more detail, perhaps even combining it with non-reaction network type models.
• Possible easier (or free) integration with future functions. E.g. if someone made solver packages for spatial systems or PDEs, instead of having to manually add this to DiffEqBiological, this would be made to work on MT type stuff, and DiffEqBiological models would get access to it for free.

(This in addition to the general advantage for new DSLs to have lots of code already written. I presume we could prune down the DiffEqBiological code length quite a bit as well, which would be good)

[ChrisRackauckas]

A lot was recently added, so I plan to write a new documentation next week. But here's sneak preview in how it relates to DiffEqBio.

See the reaction system. This is the tests:

using ModelingToolkit

@parameters t a b c
@variables x(t) y(t) z(t)
rxs = [Reaction(a,[x],[y]),
       Reaction(b,[y],[x]),
       Reaction(c,[x,y],[x])]
rs = ReactionSystem(rxs,t,[x,y],[a,b])
sys = convert(ODESystem,rs)
sys = convert(SDESystem,rs)

Basically the idea with ModelingToolkit.jl is that, DSLs should be a language, but the underlying manipulations which DSLs are doing should be part of some shared framework. So it was mostly when I started doing things like specializing the (a+b)+c behavior of MTK and the sparsity stuff started happening... that's stuff that should be ModelingToolkit.jl side so that way every single DSL doesn't have to do it. For the most part, it's really just hard to maintain so many different compilers if they are all hacked to extreme performance: hacking 1 is enough maintenance. So...

But if MT is only an easy way to define SDESystems for a DSL, given that DiffEqBiological already has it implemented there is less to be gained by using MTs version instead (but more to be won by future DSLs which do not need to do the re-implementation). There must be something more to it?

In some sense it is just that, since DiffEqBiological had a bunch, and my first goal with ModelingTookit was to make sure all of those "hacks" in the DiffEqBiological lowering process got encoded into a formal system for other DSLs to use. So ModelingToolkit has all of the automatic generation of sparse Jacobians and fast code for it, etc. With MTK, DiffEqBio should be able to drop a ton of code, and really be more about "how do I translate this syntax into a ReactionSystem". I think you'll see when playing with MTK as the internals that it's quite nice too since it obeys Julia-semantics, so doing things like sharing a reaction matrix is simply reaction_matrix.^2.

You do gain a few things too already:

• Ability to define DAEs exists, so if you have a reaction system + constraints, this can auto-generate code for the solver that satisfies the constraints.
• Automatic parallelism. If you have something like @isaacsas 's 25,000 reaction system, there's an option in here multithread=true that makes the generated ODE do 25,000/threads ODEs/SDEs etc. per thread. This also applies to the sparse Jacobian calculation.
• Support for compilation that allows other types. StaticArrays work on small systems, unit checking works, etc.

So that's what we have right now, but there's a lot in the pipeline: BLT and tearing of nonlinear systems to improve numerical solver performance, automatic model order reduction so you can solve smaller approximate models, etc. So yes:

• Possible easier (or free) integration with future functions. E.g. if someone made solver packages for spatial systems or PDEs, instead of having to manually add this to DiffEqBiological, this would be made to work on MT type stuff, and DiffEqBiological models would get access to it for free.

Is probably the biggest draw.

That said, there's a few things missing. Compiling ReactionSystem to jump problems is the main one. And there isn't a handling for bifurcation problems either.

[isaacsas]

The ability for users to further modify their models with terms that aren't supported by DEBio, for example adding delay terms, combined with the broader functionality for generating sparse/optimized systems, seems to me to be the real win currently. As @ChrisRackauckas mentioned, it should also be possible to modify the add constraint macro to support DAEs now, something users have asked about before. Getting jumps out is really the only missing feature on the MT end, which seems perhaps a bit complicated due to the different types of jumps that are needed for performance reasons. (Though expression simplification would also be nice -- we did get some nice performance boosts for ODE models when manually calculating mass action derivative terms in Jacobians, eliminating extraneous multiplications and such that SymEngine generates.)

I would advocate that we don't stick with the monolithic reaction_network idea of generating everything at once for the user. Having one additional LOC to call a transform function (to ODEs, SDEs, jumps, etc), is not much more work, but I think preferable (and much more performant when studying large systems).

[TorkelE]

Thanks for the clarification, it is very helpful. Totally forgot delay terms, yes that will be a great help, constraints as well.

On the more practical side, would the goal be to let DiffEqBio generate a ReactionSystem structure, essentially using metaprogramming to transform the reaction network into

@parameters t ...
@variables ...
rxs = [Reaction(...),... 
rs = ReactionSystem(rxs,t,...,...)

and sticks rs (and possible some basic stuff) into the reaction_network output. Modelling Toolkit would then handle generating "stuff" from rs. This would tie into @isaacsas last comment, where "stuff" could be added later on. We could later turn to improving the system for manipulating reaction networks post creation.

If we are fairly certain for how the base structures should look like, it would make sense to go ahead, even if all the stuff that would act upon it isn't ready yet. Looking forward to how this will turn out!

(It seems like a later point, and I am not sure for how things will turn out, but I would likely at some point make an argument that it would be cool with a common interface for stuff generated by DiffEqBio, ParameterizedFunctions, etc.)

[isaacsas]

Yeah, this is what I was thinking. We could essentially have a DEBio reaction_network store a ReactionSystem, and our current code could be modified to just wrap the appropriate MT system generation code (for backwards compatibility, and because for now we need to keep the jump generation). Maybe this might also make the bifurcation code cleaner long term since one can directly use mathematical operations on the MT representation, avoiding much of the expression hacking we do in DEBio right now?

[TorkelE]

Possibly yes. Although it would be kinda cool to integrate into MT the ability to determine if a system is polynomial, and make such polynomial representation available. Could be useful for other stuff.

But I presume we would want the jump stuff into MT as well in the intermediate term? Should also be of general applicability. Is there anything major stopping us, or just that it is not implemented?

[isaacsas]

My understanding is that jumps are not implemented yet. I'm not sure if @ChrisRackauckas has an interface in mind for specifying jump rates and affects, that is flexible enough to handle different types of jumps, or if the plan is to also write code to analyze these expressions when generating jumps and automatically bin them into ConstantRate, VariableRate or MassAction types as appropriate.

[ChrisRackauckas]

For the general idea, check out ParameterizedFunctions.jl. It too was once SymEngine based, but now it's ModelingToolkit based. Essentially it gives a very nice interface, generates an ODESystem, and does a few MTK calls to be done. Simple and it's highly useful (it's asked about all of the time!). I'd like to get a similar thing for DiffEqBio going, since indeed the expression hacking that's going on there right now is less than optimal for future maintenance.

One of the biggest issues is that a lot is going on in the DiffEqBio backend, but it was written to get the job done instead of being something documented and given nice usability features for people to hack away at. So while it does make a sparse matrix and it works, someone really has to study the code to see why. With an MTK backend, it's just you make a matrix and call sparse, or calculate_sparse_jacobian, because these are (going to be) well-documented functions in a system that was designed for people interact with. Dealing with macros will always be messy, but I think splitting the messy parsing from the back calculations has a lot of advantages, and you can see from ParameterizedFunctions.jl that there's still a lot of parsing AST hacking that still has to happen in order to understand the syntax, but in the end the actual code generation is like 10 lines and benefits from all of the work happening here.

Yeah, this is what I was thinking. We could essentially have a DEBio reaction_network store a ReactionSystem, and our current code could be modified to just wrap the appropriate MT system generation code (for backwards compatibility, and because for now we need to keep the jump generation). Maybe this might also make the bifurcation code cleaner long term since one can directly use mathematical operations on the MT representation, avoiding much of the expression hacking we do in DEBio right now?

Yeah, I'm uncertain about what a front-end looks like, so I've just been churning through the backend. Reaction(a,[x],[y]) could be

@reaction (a) x->y

or

@reactions begin
  a, x->y
  b, y->x
  c, x+y -> x
end

could be equivalent to

rxs = [Reaction(a,[x],[y]),
       Reaction(b,[y],[x]),
       Reaction(c,[x,y],[x])]

Making it so it just generates arrays of reactions will add this composability.

Also, I don't know if I've mentioned that there's a component system. On ODEs it looks like:

using ModelingToolkit, StaticArrays, LinearAlgebra
using DiffEqBase
using Test

# Define some variables
@parameters t σ ρ β
@variables x(t) y(t) z(t)
@derivatives D'~t

eqs = [D(x) ~ σ*(y-x),
       D(y) ~ x*(ρ-z)-y,
       0 ~ x + y + β*z]

lorenz1 = ODESystem(eqs,t,[x,y,z],[σ,ρ,β],name=:lorenz1)
lorenz2 = ODESystem(eqs,t,[x,y,z],[σ,ρ,β],name=:lorenz2)

@parameters α
@variables a(t)
connnectedeqs = [D(a) ~ a*lorenz1.x + lorenz1.σ * lorenz2.y]

connected1 = ODESystem(connnectedeqs,t,[a],[α],systems=[lorenz1,lorenz2],name=:connected1)

So you can build two ODEs and then make an ODE that couples them together. The same thing is allows with ReactionSystem types (it's part of the AbstractSystem API), so you can make a Cell, and then make a model that is 20 cells with cell to cell communication at the higher level, and MTK will generate the full code there. It will soon specialize some calculations like the Jacobian on the block structure.

My understanding is that jumps are not implemented yet. I'm not sure if @ChrisRackauckas has an interface in mind for specifying jump rates and affects, that is flexible enough to handle different types of jumps, or if the plan is to also write code to analyze these expressions when generating jumps and automatically bin them into ConstantRate, VariableRate or MassAction types as appropriate.

Yes, and this is the hard part that I have avoided so far. Making it generate all of these types and figuring out the right one will just take a bit of time to parse out. For a first pass I (or whoever) might make it all ConstantRate, and then add the MassAction optimization.

[ChrisRackauckas]

BTW, the stuff for the Pantelides algorithm (https://github.com/SciML/ModelingToolkit.jl/pull/285) needs to generate connectivity graphs, so I let @YingboMa know that we want to formalize the graph interface on our systems since the jump tools use them.

(This is actually one of the main reasons I started ModelingToolkit, since we started going down this whole route of using graph algorithms at the same time I knew I'd need the same tools for DAEs, so I started seeing having 2,3,4,5,etc. implementations of the same graph tooling and...)

[isaacsas]

Graphs for jumps would be great. Shouldn't it also be possible to do within DiffEqJump using the sparsity detection features (from SparseDiffTools)? Typical dependency graphs are basically just looking at the rate functions and determining which variables they depend on, and looking at the affect functions and determining which variables they modify.

[ChrisRackauckas]

Yeah, that's more general but at the same time a bit more difficult.

BTW, one thing I was thinking about was that I think we should make a JumpSystem and types for jumps, since then ReactionSystem would be able how to create the JumpSystem, and then the JumpSystem would handle how to build Julia functions from the jumps, which are two separate things.

[isaacsas]

That sounds good, and is what I was referring to earlier. I just wasn't sure if you wanted to have separate types in MT for the jumps, or somehow have them automatically inferred based on the rate expressions when building a JumpSystem.

[isaacsas]

Why not keep DEBio for the frontend? Your @reactions is basically the @reaction_network syntax. We can then modify DEBio to generate a ReactionSystem. It already has a well-documented API for adding reactions/species/etc to an existing system, that could be easily extended to handle the new MT variables instead of Symbols/Expressions. It would also keep more of that kind of front end stuff out of MT, keeping it more compact here. I'm fine with making an eventual breaking DEBio major release that has the @reaction_network macro just return a ReactionSystem and modifies the API accordingly. (Or perhaps it should be DEBio2, a new package, though in that case there is no reason not to just put all that code here, beyond keeping MT from having a domain-specific DSL embedded in it...)

I'd also advocate that like the rest of MT, if a DSL is made here the macros should just be syntactic sugar for an underlying API that can be called programmatically (i.e. addrx!(...), addspecies!(...) or such). That was one of the things I think was problematic in DEBio, and took some work to correct.

[ChrisRackauckas]

I'd also advocate that like the rest of MT, if a DSL is made here the macros should just be syntactic sugar for an underlying API that can be called programmatically (i.e. addrx!(...), addspecies!(...) or such). That was one of the things I think was problematic in DEBio, and took some work to correct.

Indeed, I agree with this. I really just want small macros: the nice thing about MTK is directly interacting with it.

Why not keep DEBio for the frontend? Your @reactions is basically the @reaction_network syntax. We can then modify DEBio to generate a ReactionSystem. It already has a well-documented API for adding reactions/species/etc to an existing system, that could be easily extended to handle the new MT variables instead of Symbols/Expressions. It would also keep more of that kind of front end stuff out of MT, keeping it more compact here. I'm fine with making an eventual breaking DEBio major release that has the @reaction_network macro just return a ReactionSystem and modifies the API accordingly. (Or perhaps it should be DEBio2, a new package, though in that case there is no reason not to just put all that code here, beyond keeping MT from having a domain-specific DSL embedded in it...)

Yes, it just depends on how big such a frontend is. Like ParameterizedFunctions: at one point I thought maybe that DSLish kind of thing should be MTK, but no, the parsing and all of that is still a big enough job for a package. How intense is @reactions and how many features does it have? The main reason why I was mentioning that was to describe how what it lowers to could be something that allows for more composbility than what we had before. But if @reaction_network creates a ReactionSystem and since we allow making a ReactionSystem composed of ReactionSystem, that might be the right level to work at.

MTK is about the backend. I think it ends at the point where "macros start getting messy". @variable a is just an easy way to make a variable that matches its name, i.e. a = Variable(:a). I plan to do simple things like @ODESystem a = (...) which again is just ODESystem(...,name=:a). but MTK shouldn't include any more macros than what's needed for the nice syntactic sugar for using MTK: that should be left for the DSL packages. Of course, exactly where to draw that line is somewhat difficult.

DEBio definitely is beyond that line, so the question is mostly, what are the reusable backend components are resusable and what is the part where DEBio as a domain-language begins? I think MTK should have enough so that ReactionNetworkImporters directly targets MTK for example, since I think it's a little weird to have to try and make it make a Julia syntax just to put it through a macro. So MTK should having enough to be the target for an SBML reader too. And similarly, DEBio should be our Julia-based language targeting MTK, where all 3 are then using the same underlying compilation tools in harmony.

[TorkelE]

Just read through and saw the stuff about connecting systems. That is cool, and will be very useful for expanding into systems with spatial structure. Good stuff.

I think that, given that MTK contains structures like Reaction and Reaction System, it makes sense to have the basic macros like

@reaction (a) x-->y
@reactions begin
  a, x-->y
  b, y-->x
  c, x+y --> x
end

in it. These are rather simple ones anyway, just stuffing the input into the corresponding structure.

But it also makes sense to cut away larger parts, like DiffEqBio, in separate packages (and to continue having these kinds of DSLs separate). I imagine the DiffEqBio syntax will be mostly unchanged (although the mechanics behind it, and what is made and how it is stored, likely will change).

The questions are, home much of the stuff currently in DiffEqBio do we want to import to MTK? E.g.

• Should the @reaction' in MTK take arrows such as→`, there makes for pretty code?
• Should @reaction' have special routines for double-sided arrows, like↔`?
• Should @reactions allow for "compressed" reactions like

  @reactions begin
  d, (x,y,z)-->0
  end

  ?

Then there is in the other direction. MTK would allow to, e.g, adding reactions easily to an already created system. Given that DiffEqBio has a pretty notation, not all of which will be suited to export to MTK, should DiffEqBio get macros for this kind of functionality, e.g:

@add_reactions my_reaction_network begin
  d, (x,y,z)-->0
end d

? Essentially we would thus have created syntactic sugar (DiffEqBio) for syntactic sugar (MTK). Someone call a dentist.

While it is clear that there are two niches for the backend (MTK) and the frontend (DiffEqBio), I guess the problem is that there are times when one would work directly with MTK (and it thus would be nice to have some of the shortcuts implemented in DiffEqBio available there). There are also times when one would do quite extensive operations on one's model, but prefer to do it from DiffEqBio (and thus DiffEqBio could use some of the flexibility of MTK).

In principle, I'd say that anything which one would do in running code, scripts, and just implementing stuff should be mostly possible to do from DiffEqBio. Stuff which is written into functions, personal packages, and more permanent stuff, might use some MTK stuff. But the distinction might be hard to make and we will simply have to sensibly make things up as we go along.

[isaacsas]

This is a good comment; I worry about having to maintain two versions of essentially the same code, and then having to constantly update both MT and DEBio when changing or adding something.

DEBio's @reaction_network has a fair bit of symbol/expression handling that is needed to get everything out of the macro (stoichiometries, rate expressions, etc). I guess we need to see how much of that is really needed for a more minimal macro.

[TorkelE]

Making the DiffEqBio macro just generate a ReactionSystem makes this issue redundant.

[ChrisRackauckas]

Yeah, that's why I closed it. Forgot to add the message but I think we're on the same page.