Target ModelingToolkit?

[ChrisRackauckas]

Hey Micah, Nice work. When I look at the code though:

https://github.com/mehalter/Petri.jl/blob/master/src/solvers.jl#L29-L43

I can see that the generated vector field is looping through all possible values at runtime. I think it might make sense to use some program staging here, since Petri can build some pretty big models pretty fast, and this is going to grow exponentially. If you use this process to build an MTK model though, and then have MTK output the ODE function (and some helpers like Jacobians and stuff), then the compiled versions would omit the zero'd part of these tensors, which could have a pretty wild reduction in the amount of computation that has to be done. Let me know if you want any help doing this.

[ChrisRackauckas]

Just for acceleration, modelingtoolkitize could work too, though there would be some nice things gained by targeting MTK anyways, like latex outputs.

[mehalter]

Hey Chris, thanks for the feedback! I would definitely be interested in using ModelingToolkit.jl to speed up the generation of the vector field. If you get some time, I would love some help getting started with this, its been a bit since I have looked at ModelingToolkit and am not super familiar at the moment with the current usage of it. Thanks again for your feedback and help!

[isaacsas]

This looks like a very nice project! We're actually switching DiffEqBiological over to sit on top of ModelingToolkit currently, and as part of that we've added a ReactionSystem component, here, to MT, which will be the new output from the DEBio DSL. The ReactionSystem can then be converted to ModelingToolkit-based ODE, SDE and jump models. If you read through that one file you can see how we convert a more DEBio like representation to ModelingToolkit ODEs and such -- it is a pretty minimal amount of code for ODEs/SDEs. (Jumps are more work as we try separate out into ConstantRateJumps, VariableRateJumps and MassActionJumps for performance reasons...)

It might be cool down the line to see about having conversions between the petri net representation and ReactionSystems.

[jpfairbanks]

Yeah, MT wasn't ready for us when we were using it before so we wrote this package to support the development of AlgebraicPetri.jl.

This struct seems like something we could work with.

struct ReactionSystem <: AbstractSystem
    """The reactions defining the system."""
    eqs::Vector{Reaction}
    """Independent variable (usually time)."""
    iv::Variable
    """Dependent (state) variables representing amount of each species."""
    states::Vector{Variable}
    """Parameter variables."""
    ps::Vector{Variable}
    """The name of the system"""
    name::Symbol
    """systems: The internal systems"""
    systems::Vector{ReactionSystem}
end

One of our problems with MT back in the day was support for variables indexed by integers. AlgebraicPetri contains a framework for working with systems in a way that avoids variable names, which results in systems with variables that need to be numbered or represented as x_i for i in 1:n. Is that fixed?

[jpfairbanks]

I'm also confused as to the semantics of the subsystems, how do they work?

"""
Convert a ReactionSystem to an ODESystem.
"""
function Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
    eqs = assemble_drift(rs)
    ODESystem(eqs,rs.iv,rs.states,rs.ps,name=rs.name,
              systems=convert.(ODESystem,rs.systems))
end

[isaacsas]

Do you mean something like

@variables x[1:10]

This now works in MT and you can access the variables by index, e.q. x[2].

[ChrisRackauckas]

I'm also confused as to the semantics of the subsystems, how do they work?

https://mtk.sciml.ai/dev/tutorials/ode_modeling/