Add ModelingToolkit v10 compatibility

Project.toml

@@ -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. """
@@ -398,7 +398,7 @@ function ReactionSystem(eqs, iv, unknowns, ps;
         name = nothing,
         default_u0 = Dict(),
         default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        defaults = merge(Dict(default_u0), Dict(default_p)),
         connection_type = nothing,
         checks = true,
         networkproperties = nothing,
@@ -481,13 +481,13 @@ function ReactionSystem(eqs, iv, unknowns, ps;
     end
 
     # Creates the continuous and discrete callbacks.
-    ccallbacks = MT.SymbolicContinuousCallbacks(continuous_events)
-    dcallbacks = MT.SymbolicDiscreteCallbacks(discrete_events)
+    ccallbacks = continuous_events === nothing ? MT.SymbolicContinuousCallback[] : continuous_events
+    dcallbacks = discrete_events === nothing ? MT.SymbolicDiscreteCallback[] : 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)
+        ccallbacks, dcallbacks, metadata === nothing ? Base.ImmutableDict{Symbol,Any}() : metadata; checks = checks)
 end
 
 # Two-argument constructor (reactions/equations and time variable).

src/reactionsystem_conversions.jl

@@ -507,7 +507,7 @@ COMPLETENESS_ERROR = "A ReactionSystem must be complete before it can be convert
 
 """
 ```julia
-Base.convert(::Type{<:ODESystem},rs::ReactionSystem)
+Base.convert(::typeof(ODESystem),rs::ReactionSystem)
 ```
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.ODESystem`.
 
@@ -524,11 +524,11 @@ 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(::typeof(ODESystem), 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(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)), expand_catalyst_funs = true,
+        defaults = merge(Dict(default_u0), Dict(default_p)), expand_catalyst_funs = true,
         kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
@@ -544,7 +544,7 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs)
     ODESystem(eqs, get_iv(fullrs), us, ps;
         observed = obs,
         name,
-        defaults = _merge(defaults, defs),
+        defaults = merge(defaults, defs),
         checks,
         continuous_events = MT.get_continuous_events(fullrs),
         discrete_events = MT.get_discrete_events(fullrs),
@@ -569,7 +569,7 @@ end
 
 """
 ```julia
-Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)
+Base.convert(::typeof(NonlinearSystem),rs::ReactionSystem)
 ```
 
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.NonlinearSystem`.
@@ -592,11 +592,11 @@ 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{<:NonlinearSystem}, rs::ReactionSystem; name = nameof(rs),
+function Base.convert(::typeof(NonlinearSystem), rs::ReactionSystem; name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         remove_conserved = false, conseqs_remake_warn = true, checks = false,
         default_u0 = Dict(), default_p = Dict(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        defaults = merge(Dict(default_u0), Dict(default_p)),
         all_differentials_permitted = false, expand_catalyst_funs = true, kwargs...)
     # Error checks.
     iscomplete(rs) || error(COMPLETENESS_ERROR)
@@ -625,7 +625,7 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
     NonlinearSystem(eqs, us, ps;
         name,
         observed = obs, initialization_eqs = initeqs,
-        defaults = _merge(defaults, defs),
+        defaults = merge(defaults, defs),
         checks,
         kwargs...)
 end
@@ -661,7 +661,7 @@ end
 
 """
 ```julia
-Base.convert(::Type{<:SDESystem},rs::ReactionSystem)
+Base.convert(::typeof(SDESystem),rs::ReactionSystem)
 ```
 
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.SDESystem`.
@@ -679,11 +679,11 @@ 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(::typeof(SDESystem), 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(),
-        defaults = _merge(Dict(default_u0), Dict(default_p)),
+        defaults = merge(Dict(default_u0), Dict(default_p)),
         expand_catalyst_funs = true,
         kwargs...)
     # Error checks.
@@ -707,7 +707,7 @@ function Base.convert(::Type{<:SDESystem}, rs::ReactionSystem;
     SDESystem(eqs, noiseeqs, get_iv(flatrs), us, ps;
         observed = obs,
         name,
-        defaults = _merge(defaults, defs),
+        defaults = merge(defaults, defs),
         checks,
         continuous_events = MT.get_continuous_events(flatrs),
         discrete_events = MT.get_discrete_events(flatrs),
@@ -747,7 +747,7 @@ end
 
 """
 ```julia
-Base.convert(::Type{<:JumpSystem},rs::ReactionSystem; combinatoric_ratelaws=true)
+Base.convert(::typeof(JumpSystem),rs::ReactionSystem; combinatoric_ratelaws=true)
 ```
 
 Convert a [`ReactionSystem`](@ref) to an `ModelingToolkit.JumpSystem`.
@@ -769,10 +769,10 @@ 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(::typeof(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)), expand_catalyst_funs = true,
+        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)
@@ -814,7 +814,7 @@ function Base.convert(::Type{<:JumpSystem}, rs::ReactionSystem; name = nameof(rs
     JumpSystem(eqs, get_iv(flatrs), us, ps;
         observed = obs,
         name,
-        defaults = _merge(defaults, defs),
+        defaults = merge(defaults, defs),
         checks,
         discrete_events = MT.discrete_events(flatrs),
         continuous_events = MT.continuous_events(flatrs),
@@ -923,7 +923,7 @@ Inputs for a JumpProblem from a given `ReactionSystem`.
 # Fields
 $(FIELDS)
 """
-struct JumpInputs{S <: MT.JumpSystem, T <: SciMLBase.AbstractODEProblem}
+struct JumpInputs{S <: MT.AbstractSystem, T <: SciMLBase.AbstractODEProblem}
     """The `JumpSystem` to define the problem over"""
     sys::S
     """The problem the JumpProblem should be defined over, for example DiscreteProblem"""

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}

test/dsl/dsl_options.jl

@@ -1247,7 +1247,7 @@ let
     # The `@discrete_events` option.
     rn61 = @reaction_network rn1 begin
         @species X(t)
-        @discrete_events [X > 3.0] => [X ~ X - 1]
+        @discrete_events [X > 3.0] => [X ~ Pre(X) - 1]
     end
     rn62 = @reaction_network rn1 begin
         @species X(t)
@@ -1258,7 +1258,7 @@ let
     @test isequal(rn61, rn62)
     @test_throws Exception @eval @reaction_network begin
         @species X(t)
-        @discrete_events [X > 3.0] => [X ~ X - 1] [X < 1.0] => [X ~ X + 1]
+        @discrete_events [X > 3.0] => [X ~ Pre(X) - 1] [X < 1.0] => [X ~ X + 1]
     end
 end
 

test/miscellaneous_tests/reactionsystem_serialisation.jl

@@ -309,29 +309,29 @@ let
     ]
 
     # Prepares all events.
-    continuous_events_1 = [(A ~ t_1) => [A ~ A + 2.0, X ~ X/2]]
-    continuous_events_2 = [(A ~ t_2) => [A ~ A + 2.0, X ~ X/2]]
-    continuous_events_3 = [(A ~ t_3) => [A ~ A + 2.0, X ~ X/2]]
-    continuous_events_4 = [(A ~ t_4) => [A ~ A + 2.0, X ~ X/2]]
+    continuous_events_1 = [(A ~ t_1) => [A ~ Pre(A) + 2.0, X ~ Pre(X)/2]]
+    continuous_events_2 = [(A ~ t_2) => [A ~ Pre(A) + 2.0, X ~ Pre(X)/2]]
+    continuous_events_3 = [(A ~ t_3) => [A ~ Pre(A) + 2.0, X ~ Pre(X)/2]]
+    continuous_events_4 = [(A ~ t_4) => [A ~ Pre(A) + 2.0, X ~ Pre(X)/2]]
     discrete_events_1 = [
-        10.0 => [X2_1 ~ X2_1 + 1.0]
+        10.0 => [X2_1 ~ Pre(X2_1) + 1.0]
         [5.0, 10.0] => [b_1 ~ 2 * b_1]
-        (X > 5.0) => [X2_1 ~ X2_1 + 1.0, X ~ X - 1]
+        (X > 5.0) => [X2_1 ~ Pre(X2_1) + 1.0, X ~ Pre(X) - 1]
     ]
     discrete_events_2 = [
-        10.0 => [X2_2 ~ X2_2 + 1.0]
+        10.0 => [X2_2 ~ Pre(X2_2) + 1.0]
         [5.0, 10.0] => [b_2 ~ 2 * b_2]
-        (X > 5.0) => [X2_2 ~ X2_2 + 1.0, X ~ X - 1]
+        (X > 5.0) => [X2_2 ~ Pre(X2_2) + 1.0, X ~ Pre(X) - 1]
     ]
     discrete_events_3 = [
-        10.0 => [X2_3 ~ X2_3 + 1.0]
+        10.0 => [X2_3 ~ Pre(X2_3) + 1.0]
         [5.0, 10.0] => [b_3 ~ 2 * b_3]
-        (X > 5.0) => [X2_3 ~ X2_3 + 1.0, X ~ X - 1]
+        (X > 5.0) => [X2_3 ~ Pre(X2_3) + 1.0, X ~ Pre(X) - 1]
     ]
     discrete_events_4 = [
-        10.0 => [X2_4 ~ X2_4 + 1.0]
+        10.0 => [X2_4 ~ Pre(X2_4) + 1.0]
         [5.0, 10.0] => [b_4 ~ 2 * b_4]
-        (X > 5.0) => [X2_4 ~ X2_4 + 1.0, X ~ X - 1]
+        (X > 5.0) => [X2_4 ~ Pre(X2_4) + 1.0, X ~ Pre(X) - 1]
     ]
 
     # Creates the systems.
@@ -367,8 +367,8 @@ let
     rs = @reaction_network begin
         @equations D(V) ~ 1 - V
         @continuous_events begin
-            [X ~ 5.0] => [X ~ X + 1.0]
-            [X ~ 20.0] => [X ~ X - 1.0]
+            [X ~ 5.0] => [X ~ Pre(X) + 1.0]
+            [X ~ 20.0] => [X ~ Pre(X) - 1.0]
         end
         @discrete_events 5.0 => [d ~ d/2]
         d, X --> 0

test/reactionsystem_core/custom_crn_functions.jl

@@ -305,7 +305,7 @@ let
     rates = getfield.(jsyseqs, :rate)
     affects = getfield.(jsyseqs, :affect!)
     reqs = [ Y*X*hill(X, v, K, n), Y*X*mm(X, v, K), hillr(X, v, K, n)*Y*X, Y*X*mmr(X, v, K)]
-    affeqs = [Z ~ 1 + Z, Y ~ -1 + Y, X ~ -1 + X]
+    affeqs = [Z ~ 1 + Pre(Z), Y ~ -1 + Pre(Y), X ~ -1 + Pre(X)]
     @test all(iszero, simplify(rates .- reqs))
     @test all(aff -> isequal(aff, affeqs), affects)
 
@@ -314,7 +314,7 @@ let
     rates = getfield.(jsyseqs, :rate)
     affects = getfield.(jsyseqs, :affect!)
     reqs = [ Y*X*hill2(X, v, K, n), Y*X*mm2(X, v, K), hillr2(X, v, K, n)*Y*X, Y*X*mmr2(X, v, K)]
-    affeqs = [Z ~ 1 + Z, Y ~ -1 + Y, X ~ -1 + X]
+    affeqs = [Z ~ 1 + Pre(Z), Y ~ -1 + Pre(Y), X ~ -1 + Pre(X)]
     @test all(iszero, simplify(rates .- reqs))
     @test all(aff -> isequal(aff, affeqs), affects)
 end