start on supporting hybrid jump systems

Author: isaacsas

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"
@@ -32,14 +33,12 @@ BifurcationKit = "0f109fa4-8a5d-4b75-95aa-f515264e7665"
 CairoMakie = "13f3f980-e62b-5c42-98c6-ff1f3baf88f0"
 GraphMakie = "1ecd5474-83a3-4783-bb4f-06765db800d2"
 HomotopyContinuation = "f213a82b-91d6-5c5d-acf7-10f1c761b327"
-# StructuralIdentifiability = "220ca800-aa68-49bb-acd8-6037fa93a544"
 
 [extensions]
 CatalystBifurcationKitExtension = "BifurcationKit"
 CatalystCairoMakieExtension = "CairoMakie"
 CatalystGraphMakieExtension = "GraphMakie"
 CatalystHomotopyContinuationExtension = "HomotopyContinuation"
-# CatalystStructuralIdentifiabilityExtension = "StructuralIdentifiability"
 
 [compat]
 BifurcationKit = "0.3"
@@ -50,6 +49,7 @@ DiffEqBase = "6.83.0"
 DocStringExtensions = "0.8, 0.9"
 DynamicPolynomials = "0.5, 0.6"
 DynamicQuantities = "0.13.2, 1"
+EnumX = "1.0.4"
 GraphMakie = "0.5"
 Graphs = "1.4"
 HomotopyContinuation = "2.9"
@@ -64,7 +64,6 @@ Requires = "1.0"
 RuntimeGeneratedFunctions = "0.5.12"
 SciMLBase = "< 2.57.2"
 Setfield = "1"
-# StructuralIdentifiability = "0.5.8"
 SymbolicUtils = "2.1.2, 3.3.0"
 Symbolics = "5.30.1, 6"
 Unitful = "1.12.4"
@@ -93,7 +92,4 @@ Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
 
 [targets]
-test = ["DiffEqCallbacks", "DomainSets", "Graphviz_jll", "Logging", "NonlinearSolve", 
-        "OrdinaryDiffEq", "Pkg", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", 
-        "StableRNGs", "StaticArrays", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", 
-        "Test", "Unitful"]
+test = ["DiffEqCallbacks", "DomainSets", "Graphviz_jll", "Logging", "NonlinearSolve", "OrdinaryDiffEq", "Pkg", "Plots", "Random", "SafeTestsets", "SciMLBase", "SciMLNLSolve", "StableRNGs", "StaticArrays", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful"]

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
@@ -95,7 +95,10 @@ end
 # 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,45 @@ 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
+
+"""
+    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_conversions.jl

@@ -673,7 +673,7 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         remove_conserved = nothing, checks = false,
         default_u0 = Dict(), default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        defaults = _merge(Dict(default_u0), Dict(default_p)), physical_scales = nothing
         kwargs...)
     iscomplete(rs) || error(COMPLETENESS_ERROR)
     spatial_convert_err(rs::ReactionSystem, JumpSystem)
@@ -683,9 +683,7 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs
     flatrs = Catalyst.flatten(rs)
     error_if_constraints(JumpSystem, flatrs)
 
-    (length(MT.continuous_events(flatrs)) > 0) &&
-        (@warn "continuous_events will be dropped as they are not currently supported by JumpSystems.")
-
+    physical_scales = merge_physical_scales(reactions(rs), physical_scales)
     eqs = assemble_jumps(flatrs; combinatoric_ratelaws)
 
     # handle BC species
@@ -769,6 +767,16 @@ function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan,
         noise_rate_prototype = p_matrix, kwargs...)
 end
 
+function merge_physical_scales(rxs, physical_scales)
+    scales = get_physical_scale.(rxs)
+    if physical_scales !== nothing
+        for (idx, scale) in physical_scales
+            scales[idx] = scale
+        end
+    end
+    scales
+end
+
 """
 $(TYPEDEF)
 
@@ -817,11 +825,14 @@ plot(sol, idxs = :A)
 """
 function JumpInputs(rs::ReactionSystem, u0, tspan, p = DiffEqBase.NullParameters();
         name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        checks = false, kwargs...)
+        checks = false, physical_scales = nothing, kwargs...)
     u0map = symmap_to_varmap(rs, u0)
     pmap = symmap_to_varmap(rs, p)
-    jsys = complete(convert(JumpSystem, rs; name, combinatoric_ratelaws, checks))
-    if MT.has_variableratejumps(jsys)
+    jsys = complete(convert(JumpSystem, rs; name, combinatoric_ratelaws, checks, 
+        physical_scales))
+
+    if MT.has_variableratejumps(jsys) || MT.has_equations(jsys) || 
+            !isempty(MT.continuous_events(jsys))
         prob = ODEProblem(jsys, u0map, tspan, pmap; kwargs...)
     else
         prob = DiscreteProblem(jsys, u0map, tspan, pmap; kwargs...)