Merge pull request #1121 from isaacsas/pdmps

Support Jump-ODE models

Author: isaacsas

.github/workflows/Test.yml

@@ -30,6 +30,7 @@ jobs:
           - "pre"
         group:
           - Core
+          - Hybrid
           - IO
           - Spatial
           - Extensions

Project.toml

@@ -9,6 +9,7 @@ DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 DynamicPolynomials = "7c1d4256-1411-5781-91ec-d7bc3513ac07"
 DynamicQuantities = "06fc5a27-2a28-4c7c-a15d-362465fb6821"
+EnumX = "4e289a0a-7415-4d19-859d-a7e5c4648b56"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 JumpProcesses = "ccbc3e58-028d-4f4c-8cd5-9ae44345cda5"
 LaTeXStrings = "b964fa9f-0449-5b57-a5c2-d3ea65f4040f"
@@ -51,8 +52,8 @@ Combinatorics = "1.0.2"
 DataStructures = "0.18"
 DiffEqBase = "6.165.0"
 DocStringExtensions = "0.8, 0.9"
-DynamicPolynomials = "0.5, 0.6"
-DynamicQuantities = "0.13.2, 1"
+DynamicPolynomials = "0.6"
+DynamicQuantities = "1"
 GraphMakie = "0.5"
 Graphs = "1.4"
 HomotopyContinuation = "2.9"
@@ -64,7 +65,7 @@ Makie = "0.22.1"
 ModelingToolkit = "9.69"
 NetworkLayout = "0.4.7"
 Parameters = "0.12"
-Reexport = "0.2, 1.0"
+Reexport = "1.0"
 Requires = "1.0"
 RuntimeGeneratedFunctions = "0.5.12"
 SciMLBase = "2.77"

src/Catalyst.jl

@@ -4,7 +4,7 @@ $(DocStringExtensions.README)
 module Catalyst
 
 using DocStringExtensions
-using SparseArrays, DiffEqBase, Reexport, Setfield
+using SparseArrays, DiffEqBase, Reexport, Setfield, EnumX
 using LaTeXStrings, Latexify, Requires
 using LinearAlgebra, Combinatorics
 using JumpProcesses: JumpProcesses, JumpProblem,
@@ -14,7 +14,7 @@ using JumpProcesses: JumpProcesses, JumpProblem,
 # ModelingToolkit imports and convenience functions we use
 using ModelingToolkit
 const MT = ModelingToolkit
-using DynamicQuantities#, Unitful # Having Unitful here as well currently gives an error.
+using DynamicQuantities #, Unitful # Having Unitful here as well currently gives an error.
 
 @reexport using ModelingToolkit
 using Symbolics
@@ -78,7 +78,10 @@ const forbidden_symbols_error = union(Set([:im, :nothing, CONSERVED_CONSTANT_SYM
 # The `Reaction` structure and its functions.
 include("reaction.jl")
 export isspecies
-export Reaction
+export Reaction, PhysicalScale
+
+# Union type for `Reaction`s and `Equation`s.
+const CatalystEqType = Union{Reaction, Equation}
 
 # The `ReactionSystem` structure and its functions.
 include("reactionsystem.jl")

src/reaction.jl

@@ -244,9 +244,6 @@ function Reaction(rate, subs, prods; kwargs...)
     Reaction(rate, subs, prods, sstoich, pstoich; kwargs...)
 end
 
-# Union type for `Reaction`s and `Equation`s.
-const CatalystEqType = Union{Reaction, Equation}
-
 ### Base Function Dispatches ###
 
 # Used by `Base.show`.
@@ -679,6 +676,48 @@ function getmisc(reaction::Reaction)
     end
 end
 
+############## Metadata for the mathematical type of a reaction ##############
+
+"""
+    @enumx PhysicalScale
+
+EnumX instance representing the physical scale of a reaction. 
+
+Notes: The following values are possible:
+- `Auto`: (DEFAULT) Lets Catalyst decide at the time of system conversion and/or
+  problem generation at what physical scale to represent the reaction.
+- `ODE`: The reaction is to be treated via an ordinary differential equation term.
+- `SDE`: The reaction is to be treated via a stochastic differential equation (CLE) term.
+- `Jump`: The reaction is to be treated via a jump process (stochastic chemical kinetics)
+  term, letting Catalyst decide the specific jump type.
+- `VariableRateJump`: The reaction is to be treated as a jump process (stochastic chemical
+  kinetics) term, specifically assigning it to a VariableRateJump.
+"""
+@enumx PhysicalScale begin
+    Auto             # the default that lets Catalyst decide 
+    ODE
+    SDE
+    Jump             # lets Catalyst decide the jump type
+    VariableRateJump # forces a VariableRateJump
+end
+
+const JUMP_SCALES = (PhysicalScale.Jump, PhysicalScale.VariableRateJump)
+const NON_CONSTANT_JUMP_SCALES = (PhysicalScale.ODE, PhysicalScale.SDE, PhysicalScale.VariableRateJump)
+
+"""
+    has_physical_scale(rx::Reaction)
+
+Returns `true` if the input reaction has the `physical_scale` metadata field assigned, 
+else `false`.
+"""
+function has_physical_scale(rx::Reaction)
+    return hasmetadata(rx, :physical_scale)
+end
+
+function get_physical_scale(rx::Reaction)
+    return has_physical_scale(rx) ? getmetadata(rx, :physical_scale) : PhysicalScale.Auto
+end
+
 ### Units Handling ###
 
 """

src/reactionsystem.jl

@@ -871,6 +871,16 @@ function numreactions(network)
     nr
 end
 
+"""
+    has_nonreactions(network)
+
+Check if the given `network` has any non-reaction equations such as ODEs or algebraic
+equations.
+"""
+function has_nonreactions(network)
+    numreactions(network) != length(equations(network))
+end
+
 """
     nonreactions(network)