WIP: MTK v10 compatibility changes

Author: ctessum

Project.toml

@@ -1,6 +1,6 @@
 name = "Catalyst"
 uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
-version = "15.0.8"
+version = "16"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -62,7 +62,7 @@ LaTeXStrings = "1.3.0"
 Latexify = "0.16.6"
 MacroTools = "0.5.5"
 Makie = "0.22.1"
-ModelingToolkit = "9.73"
+ModelingToolkit = "10"
 NetworkLayout = "0.4.7"
 Parameters = "0.12"
 Reexport = "1.0"

src/Catalyst.jl

@@ -36,7 +36,7 @@ import ModelingToolkit: get_variables, namespace_expr, namespace_equation, get_v
 
 # internal but needed ModelingToolkit functions
 import ModelingToolkit: check_variables,
-                        check_parameters, _iszero, _merge, check_units,
+                        check_parameters, _iszero, check_units,
                         get_unit, check_equations, iscomplete
 
 import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show

src/reactionsystem.jl

@@ -276,7 +276,7 @@ Notes:
   units). Unit checking can be disabled by passing the keyword argument `checks=false`.
 """
 struct ReactionSystem{V <: NetworkProperties} <:
-       MT.AbstractTimeDependentSystem
+       MT.AbstractSystem
     """The equations (reactions and algebraic/differential) defining the system."""
     eqs::Vector{CatalystEqType}
     """The Reactions defining the system. """
@@ -356,7 +356,7 @@ struct ReactionSystem{V <: NetworkProperties} <:
             check_parameters(ps, iv)
             nonrx_eqs = Equation[eq for eq in eqs if eq isa Equation]
             !isempty(nonrx_eqs) && check_equations(nonrx_eqs, iv)
-            !isempty(cevs) && check_equations(equations(cevs), iv)
+            !isnothing(cevs) && !isempty(cevs) && check_equations(equations(cevs), iv)
         end
 
         if isempty(sivs) && (checks == true || (checks & MT.CheckUnits) > 0)
@@ -480,14 +480,10 @@ function ReactionSystem(eqs, iv, unknowns, ps;
         networkproperties
     end
 
-    # Creates the continuous and discrete callbacks.
-    ccallbacks = MT.SymbolicContinuousCallbacks(continuous_events)
-    dcallbacks = MT.SymbolicDiscreteCallbacks(discrete_events)
-
     ReactionSystem(
         eqs′, rxs, iv′, sivs′, unknowns′, spcs, ps′, var_to_name, observed, name,
         systems, defaults, connection_type, nps, combinatoric_ratelaws,
-        ccallbacks, dcallbacks, metadata; checks = checks)
+        continuous_events, discrete_events, metadata; checks = checks)
 end
 
 # Two-argument constructor (reactions/equations and time variable).

src/reactionsystem_conversions.jl

@@ -386,7 +386,7 @@ function assemble_jumps(rs; combinatoric_ratelaws = true, physical_scales = noth
         error("Must have at least one reaction that will be represented as a jump when constructing a JumpSystem.")
 
     # note isvrjvec indicates which reactions are not constant rate jumps
-    # it may be that a given jump has isvrjvec[i] = true but has a physical 
+    # it may be that a given jump has isvrjvec[i] = true but has a physical
     isvrjvec = classify_vrjs(rs, physcales)
 
     rxvars = []
@@ -505,9 +505,16 @@ COMPLETENESS_ERROR = "A ReactionSystem must be complete before it can be convert
 
 ### System Conversions ###
 
+abstract type SystemType end
+
+struct ReactionRateSystem <: SystemType end
+struct ChemicalLangevinSystem <: SystemType end
+struct StochasticChemicalKineticSystem <: SystemType end
+struct NonlinearSystem <: SystemType end
+
 """
 ```julia
-Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
+Base.convert(::Type{<:System},rs::ReactionSystem)
 ```
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.ODESystem`.
 
@@ -524,7 +531,7 @@ Keyword args and default values:
   with their rational function representation when converting to another system type. Set to
   `false`` to disable.
 """
-function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs),
+function Base.convert(::Type{<:ReactionRateSystem}, rs::ReactionSystem; name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         include_zero_odes = true, remove_conserved = false, checks = false,
         default_u0 = Dict(), default_p = Dict(),
@@ -541,7 +548,7 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs)
         include_zero_odes, expand_catalyst_funs)
     eqs, us, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
 
-    ODESystem(eqs, get_iv(fullrs), us, ps;
+    System(eqs, get_iv(fullrs), us, ps;
         observed = obs,
         name,
         defaults = _merge(defaults, defs),
@@ -622,7 +629,7 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
     all_differentials_permitted || nonlinear_convert_differentials_check(rs)
     eqs = [remove_diffs(eq.lhs) ~ remove_diffs(eq.rhs) for eq in eqs]
 
-    NonlinearSystem(eqs, us, ps;
+    System(eqs, us, ps;
         name,
         observed = obs, initialization_eqs = initeqs,
         defaults = _merge(defaults, defs),
@@ -679,7 +686,7 @@ Notes:
   with their rational function representation when converting to another system type. Set to
   `false`` to disable.
 """
-function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
+function Base.convert(::Type{<:ChemicalLangevinSystem}, rs::ReactionSystem;
         name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         include_zero_odes = true, checks = false, remove_conserved = false,
         default_u0 = Dict(), default_p = Dict(),
@@ -704,7 +711,7 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
         @info "Boundary condition species detected. As constraint equations are not currently supported when converting to SDESystems, the resulting system will be undetermined. Consider using constant species instead."
     end
 
-    SDESystem(eqs, noiseeqs, get_iv(flatrs), us, ps;
+    System([eqs; noiseeqs], get_iv(flatrs), us, ps;
         observed = obs,
         name,
         defaults = _merge(defaults, defs),
@@ -728,7 +735,7 @@ Merge physical scales for a set of reactions.
 function merge_physical_scales(rxs, physical_scales, default)
     scales = get_physical_scale.(rxs)
 
-    # override metadata attached scales 
+    # override metadata attached scales
     if physical_scales !== nothing
         for (key, scale) in physical_scales
             scales[key] = scale
@@ -769,13 +776,13 @@ Notes:
   `VariableRateJump` to save the solution before and/or after the jump occurs. Defaults to
   true for both.
 """
-function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs),
+function Base.convert(::Type{<:StochasticChemicalKineticSystem}, 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)), expand_catalyst_funs = true,
         save_positions = (true, true), physical_scales = nothing, kwargs...)
     iscomplete(rs) || error(COMPLETENESS_ERROR)
-    spatial_convert_err(rs::ReactionSystem, JumpSystem)
+    spatial_convert_err(rs::ReactionSystem, StochasticChemicalKineticSystem)
     (remove_conserved !== nothing) &&
         throw(ArgumentError("Catalyst does not support removing conserved species when converting to JumpSystems."))
 
@@ -794,7 +801,7 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs
         physical_scales, save_positions)
     ists, ispcs = get_indep_sts(flatrs)
 
-    # handle coupled ODEs and BC species    
+    # handle coupled ODEs and BC species
     if (PhysicalScale.ODE in unique_scales) || has_nonreactions(flatrs)
         odeeqs = assemble_drift(flatrs, ispcs; combinatoric_ratelaws,
             remove_conserved = false, physical_scales, include_zero_odes = true)
@@ -811,7 +818,7 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs
         defs = MT.defaults(flatrs)
     end
 
-    JumpSystem(eqs, get_iv(flatrs), us, ps;
+    System(eqs, get_iv(flatrs), us, ps;
         observed = obs,
         name,
         defaults = _merge(defaults, defs),
@@ -850,8 +857,8 @@ end
 DiffEqBase.NonlinearProblem(rs::ReactionSystem, u0,
         p = DiffEqBase.NullParameters(), args...;
         name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-        remove_conserved = false, checks = false, check_length = false, 
-        structural_simplify = remove_conserved, all_differentials_permitted = false, 
+        remove_conserved = false, checks = false, check_length = false,
+        structural_simplify = remove_conserved, all_differentials_permitted = false,
         kwargs...)
 ```
 
@@ -923,7 +930,7 @@ Inputs for a JumpProblem from a given `ReactionSystem`.
 # Fields
 $(FIELDS)
 """
-struct JumpInputs{S <: MT.JumpSystem, T <: SciMLBase.AbstractODEProblem}
+struct JumpInputs{S <: MT.System, T <: SciMLBase.AbstractODEProblem}
     """The `JumpSystem` to define the problem over"""
     sys::S
     """The problem the JumpProblem should be defined over, for example DiscreteProblem"""
@@ -936,8 +943,8 @@ JumpInputs(rs::ReactionSystem, u0, tspan,
             p = DiffEqBase.NullParameters;
             name = nameof(rs),
             combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
-            checks = false, physical_scales = nothing, 
-            expand_catalyst_funs = true, 
+            checks = false, physical_scales = nothing,
+            expand_catalyst_funs = true,
             save_positions = (true, true),
             remake_warn = true, kwargs...)
 ```
@@ -988,7 +995,7 @@ function JumpInputs(rs::ReactionSystem, u0, tspan, p = DiffEqBase.NullParameters
         name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         checks = false, physical_scales = nothing, expand_catalyst_funs = true,
         save_positions = (true, true), remake_warn = true, kwargs...)
-    jsys = complete(convert(JumpSystem, rs; name, combinatoric_ratelaws, checks,
+    jsys = complete(convert(StochasticChemicalKineticSystem, rs; name, combinatoric_ratelaws, checks,
         physical_scales, expand_catalyst_funs, save_positions))
 
     if MT.has_variableratejumps(jsys) || MT.has_equations(jsys) ||
@@ -1025,7 +1032,7 @@ function DiffEqBase.DiscreteProblem(rs::ReactionSystem, u0, tspan::Tuple,
     Base.depwarn("DiscreteProblem(rn::ReactionSystem, ...) is deprecated and will be \
         removed in Catalyst 16. Use JumpInputs(rn, ...) instead.",
         :DiscreteProblem)
-    jsys = convert(JumpSystem, rs; name, combinatoric_ratelaws, checks,
+    jsys = convert(StochasticChemicalKineticSystem, rs; name, combinatoric_ratelaws, checks,
         expand_catalyst_funs)
     jsys = complete(jsys)
     return DiscreteProblem(jsys, u0, tspan, p, args...; kwargs...)
@@ -1040,7 +1047,7 @@ function JumpProcesses.JumpProblem(rs::ReactionSystem, prob::SciMLBase.AbstractD
     Base.depwarn("JumpProblem(rn::ReactionSystem, prob, ...) is \
         deprecated and will be removed in Catalyst 16. Use \
         JumpProblem(JumpInputs(rn, ...), ...) instead.", :JumpProblem)
-    jsys = convert(JumpSystem, rs; name, combinatoric_ratelaws,
+    jsys = convert(StochasticChemicalKineticSystem, rs; name, combinatoric_ratelaws,
         expand_catalyst_funs, checks)
     jsys = complete(jsys)
     return JumpProblem(jsys, prob, aggregator; kwargs...)

src/spatial_reaction_systems/lattice_reaction_systems.jl

@@ -57,7 +57,7 @@ continuous space systems with them is possible, but requires the user to determi
 (the lattice). Better support for continuous space models is a work in progress.
 - Catalyst contains extensive documentation on spatial modelling, which can be found [here](https://docs.sciml.ai/Catalyst/stable/spatial_modelling/lattice_reaction_systems/).
 """
-struct LatticeReactionSystem{Q, R, S, T} <: MT.AbstractTimeDependentSystem
+struct LatticeReactionSystem{Q, R, S, T} <: MT.AbstractSystem
     # Input values.
     """The (non-spatial) reaction system within each vertex."""
     reactionsystem::ReactionSystem{Q}

src/spatial_reaction_systems/utility.jl

@@ -40,7 +40,7 @@ function lattice_process_p(ps_in, ps_vertex_syms::Vector,
     # For uniform parameters these have size 1/(1,1). Else, they have size num_verts/(num_verts,num_verts).
     ps = lattice_process_input(ps_in, [ps_vertex_syms; ps_edge_syms])
 
-    # Split the parameter vector into one for vertex parameters and one for edge parameters. 
+    # Split the parameter vector into one for vertex parameters and one for edge parameters.
     # Next, convert their values to the correct form (vectors for vert_ps and sparse matrices for edge_ps).
     vert_ps, edge_ps = split_parameters(ps, ps_vertex_syms, ps_edge_syms)
     vert_ps = vertex_value_map(vert_ps, lrs)
@@ -164,13 +164,13 @@ function edge_value_form(values, lrs::LatticeReactionSystem, sym)
         throw(ArgumentError("Values was not provided for some edges for edge parameter $sym."))
     end
 
-    # Unlike initial conditions/vertex parameters, (unless uniform) edge parameters' values  are 
+    # Unlike initial conditions/vertex parameters, (unless uniform) edge parameters' values  are
     # always provided in the same (sparse matrix) form.
     return values
 end
 
 # Creates a map, taking each species (with transportation) to its transportation rate.
-# The species is represented by its index (in species(lrs). 
+# The species is represented by its index (in species(lrs).
 # If the rate is uniform across all edges, the transportation rate will be a size (1,1) sparse matrix.
 # Else, the rate will be a size (num_verts,num_verts) sparse matrix.
 function make_sidxs_to_transrate_map(vert_ps::Vector{Pair{R, Vector{T}}},
@@ -190,7 +190,7 @@ end
 # Computes the transport rates for all species with transportation rates. Output is a map
 # taking each species' symbolics form to its transportation rates across all edges.
 function compute_all_transport_rates(p_val_dict, lrs::LatticeReactionSystem)
-    # For all species with transportation, compute their transportation rate (across all edges). 
+    # For all species with transportation, compute their transportation rate (across all edges).
     # This is a vector, pairing each species to these rates.
     unsorted_rates = [s => compute_transport_rates(s, p_val_dict, lrs)
                       for s in spatial_species(lrs)]
@@ -199,7 +199,7 @@ function compute_all_transport_rates(p_val_dict, lrs::LatticeReactionSystem)
     return sort(unsorted_rates; by = rate -> findfirst(isequal(rate[1]), species(lrs)))
 end
 
-# For the expression describing the rate of transport (likely only a single parameter, e.g. `D`), 
+# For the expression describing the rate of transport (likely only a single parameter, e.g. `D`),
 # and the values of all our parameters, compute the transport rate(s).
 # If all parameters that the rate depends on are uniform across all edges, this becomes a length-1 vector.
 # Else it becomes a vector where each value corresponds to the rate at one specific edge.
@@ -304,10 +304,21 @@ end
 
 ### System Property Checks ###
 
-# For a Symbolic expression, and a parameter set, check if any relevant parameters have a 
+# For a Symbolic expression, and a parameter set, check if any relevant parameters have a
 # spatial component. Filters out any parameters that are edge parameters.
 function has_spatial_vertex_component(exp, ps)
     relevant_syms = Symbolics.get_variables(exp)
     value_dict = Dict(filter(p -> p[2] isa Vector, ps))
     return any(length(value_dict[sym]) > 1 for sym in relevant_syms)
 end
+
+### Utilities previously in ModelingToolkit.jl ###
+
+function todict(d)
+    eltype(d) <: Pair || throw(ArgumentError("The variable-value mapping must be a Dict."))
+    d isa Dict ? d : Dict(d)
+end
+
+_merge(d1, d2) = merge(todict(d1), todict(d2))
+
+MT.refreshed_metadata(::Nothing) = MT.MetadataT() # FIXME: Type piracy