WIP: fixes for MTK changes

HISTORY.md

@@ -1,6 +1,10 @@
 # Breaking updates and feature summaries across releases
 
 ## Catalyst unreleased (master branch)
+- `ModelingToolkit.complete` now adds initial parameters corresponding to
+  initial conditions into the internal `sys.ps` field (accessed via
+  `get_ps(sys)`). For this reason `parameters(sys)` and `get_ps(sys)` are no
+  longer equivalent for a `complete`d system (as the DSL generates).
 - The Catalyst release process is changing; certain core dependencies of
   Catalyst will now be capped to ensure Catalyst releases are only installed
   with versions of dependencies for which Catalyst CI and doc build tests pass

Project.toml

@@ -56,7 +56,7 @@ JumpProcesses = "9.13.2"
 LaTeXStrings = "1.3.0"
 Latexify = "0.16.6"
 MacroTools = "0.5.5"
-ModelingToolkit = "< 9.60"
+ModelingToolkit = "9.64"
 NetworkLayout = "0.4.7"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"

src/reactionsystem.jl

@@ -71,6 +71,7 @@ end
 Base.Sort.defalg(::ReactionComplex) = Base.DEFAULT_UNSTABLE
 
 ### NetworkProperties Structure ###
+const __UNINITIALIZED_CONSERVED_CONSTS = MT.unwrap(only(@parameters __UNINITIALIZED[1]))
 
 #! format: off
 # Internal cache for various ReactionSystem calculated properties
@@ -86,6 +87,7 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     depspecs::Set{V} = Set{V}()
     conservedeqs::Vector{Equation} = Equation[]
     constantdefs::Vector{Equation} = Equation[]
+    conservedconst::BasicSymbolic{Vector{Real}} = __UNINITIALIZED_CONSERVED_CONSTS
     speciesmap::Dict{V, Int} = Dict{V, Int}()
     complextorxsmap::OrderedDict{ReactionComplex{Int}, Vector{Pair{Int, Int}}} = OrderedDict{ReactionComplex{Int},Vector{Pair{Int,Int}}}()
     complexes::Vector{ReactionComplex{Int}} = Vector{ReactionComplex{Int}}(undef, 0)
@@ -96,7 +98,6 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     linkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
     stronglinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
     terminallinkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
-
     checkedrobust::Bool = false
     robustspecies::Vector{Int} = Vector{Int}(undef, 0)
     deficiency::Int = -1
@@ -126,6 +127,7 @@ function reset!(nps::NetworkProperties{I, V}) where {I, V}
     empty!(nps.col_order)
     nps.rank = 0
     nps.nullity = 0
+    nps.conservedconst = __UNINITIALIZED_CONSERVED_CONSTS
     empty!(nps.indepspecs)
     empty!(nps.depspecs)
     empty!(nps.conservedeqs)
@@ -149,6 +151,10 @@ function reset!(nps::NetworkProperties{I, V}) where {I, V}
     nothing
 end
 
+function initialized_conserved(nps::NetworkProperties)
+    nps.conservedconst !== __UNINITIALIZED_CONSERVED_CONSTS
+end
+
 ### ReactionSystem Constructor Functions ###
 
 # Used to sort the reaction/equation vector as reactions first, equations second.
@@ -583,34 +589,61 @@ function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
     isequivalent(rn1, rn2; ignorenames = false)
 end
 
+function debug_comparer(fun, prop1, prop2, propname; debug = false)
+    if fun(prop1, prop2)
+        return true
+    else
+        debug && println("Comparison was false for property: ", propname, "\n    Found: ", prop1, " vs ", prop2)
+        return false
+    end
+end
+
 """
-    isequivalent(rn1::ReactionSystem, rn2::ReactionSystem; ignorenames = true)
+    isequivalent(rn1::ReactionSystem, rn2::ReactionSystem; ignorenames = true, 
+        debug = false)
 
-Tests whether the underlying species, parameters and reactions are the same in
-the two [`ReactionSystem`](@ref)s. Ignores the names of the systems in testing
-equality.
+Tests whether the underlying species, parameters and reactions are the same in the two
+[`ReactionSystem`](@ref)s. Ignores the names of the systems in testing equality.
 
 Notes:
-- *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
-    considered different than `(A+1)^2`.
-- Does not include `defaults` in determining equality.
+- *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be considered
+    different than `(A+1)^2`.
+- `ignorenames = false` is used when checking equality of sub and parent systems.
+- Does not check that `parent` systems are the same.
+- Pass `debug = true` to print out the field that caused the two systems to be considered
+    different.
 """
-function isequivalent(rn1::ReactionSystem, rn2::ReactionSystem; ignorenames = true)
+function isequivalent(rn1::ReactionSystem, rn2::ReactionSystem; ignorenames = true,
+        debug = false)    
+
+    # metadata type fields
     if !ignorenames
-        (nameof(rn1) == nameof(rn2)) || return false
+        debug_comparer(==, nameof(rn1), nameof(rn2), "name"; debug) || return false
     end
-
-    (get_combinatoric_ratelaws(rn1) == get_combinatoric_ratelaws(rn2)) || return false
-    isequal(get_iv(rn1), get_iv(rn2)) || return false
-    issetequal(get_sivs(rn1), get_sivs(rn2)) || return false
-    issetequal(get_unknowns(rn1), get_unknowns(rn2)) || return false
-    issetequal(get_ps(rn1), get_ps(rn2)) || return false
-    issetequal(MT.get_observed(rn1), MT.get_observed(rn2)) || return false
-    issetequal(get_eqs(rn1), get_eqs(rn2)) || return false
-
-    # subsystems
-    (length(get_systems(rn1)) == length(get_systems(rn2))) || return false
-    issetequal(get_systems(rn1), get_systems(rn2)) || return false
+    debug_comparer(==, get_combinatoric_ratelaws(rn1), get_combinatoric_ratelaws(rn2), 
+        "combinatoric_ratelaws"; debug) || return false
+    debug_comparer(==, MT.iscomplete(rn1), MT.iscomplete(rn2), "complete"; debug) || return false
+
+    # symbolic variables and parameters
+    debug_comparer(isequal, get_iv(rn1), get_iv(rn2), "ivs"; debug) || return false
+    debug_comparer(issetequal, get_sivs(rn1), get_sivs(rn2), "sivs"; debug) || return false
+    debug_comparer(issetequal, get_unknowns(rn1), get_unknowns(rn2), "unknowns"; debug) || return false
+    debug_comparer(issetequal, get_species(rn1), get_species(rn2), "species"; debug) || return false
+    debug_comparer(issetequal, get_ps(rn1), get_ps(rn2), "ps"; debug) || return false
+    debug_comparer(issetequal, MT.get_defaults(rn1), MT.get_defaults(rn2), "defaults"; debug) || return false
+
+    # equations and reactions
+    debug_comparer(issetequal, MT.get_observed(rn1), MT.get_observed(rn2), "observed"; debug) || return false
+    debug_comparer(issetequal, get_eqs(rn1), get_eqs(rn2), "eqs"; debug) || return false
+    debug_comparer(issetequal, MT.get_continuous_events(rn1), MT.get_continuous_events(rn2), "cevents"; debug) || return false
+    debug_comparer(issetequal, MT.get_discrete_events(rn1), MT.get_discrete_events(rn2), "devents"; debug) || return false
+
+    # coupled systems    
+    if (length(get_systems(rn1)) != length(get_systems(rn2))) 
+        println("Systems have different numbers of subsystems.")
+        return false
+    end
+    debug_comparer(issetequal, get_systems(rn1), get_systems(rn2), "systems"; debug) || return false
 
     true
 end

src/reactionsystem_conversions.jl

@@ -392,11 +392,9 @@ function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved
         nps = get_networkproperties(rs)
 
         # add the conservation constants as parameters and set their values
-        ps = vcat(ps, collect(eq.lhs for eq in nps.constantdefs))
+        ps = copy(ps)
+        push!(ps, nps.conservedconst)
         defs = copy(MT.defaults(rs))
-        for eq in nps.constantdefs
-            defs[eq.lhs] = eq.rhs
-        end
 
         # add the dependent species as observed
         obs = copy(MT.observed(rs))
@@ -410,11 +408,12 @@ function addconstraints!(eqs, rs::ReactionSystem, ists, ispcs; remove_conserved
     if !isempty(ceqs)
         if remove_conserved
             @info """
-                  Be careful mixing constraints and elimination of conservation laws.
-                  Catalyst does not check that the conserved equations still hold for the
-                  final coupled system of equations. Consider using `remove_conserved =
-                  false` and instead calling ModelingToolkit.structural_simplify to simplify
-                  any generated ODESystem or NonlinearSystem.
+                  Be careful mixing ODEs or algebraic equations and elimination of
+                  conservation laws. Catalyst does not check that the conserved equations
+                  still hold for the final coupled system of equations. Consider using
+                  `remove_conserved = false` and instead calling
+                  ModelingToolkit.structural_simplify to simplify any generated ODESystem or
+                  NonlinearSystem.
                   """
         end
         append!(eqs, ceqs)
@@ -503,7 +502,6 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs)
     eqs = assemble_drift(fullrs, ispcs; combinatoric_ratelaws, remove_conserved,
         include_zero_odes)
     eqs, us, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
-
     ODESystem(eqs, get_iv(fullrs), us, ps;
         observed = obs,
         name,
@@ -514,6 +512,17 @@ function Base.convert(::Type{<:ODESystem}, rs::ReactionSystem; name = nameof(rs)
         kwargs...)
 end
 
+# the following remove Initial wrapped parameters (i.e. initial conditions) from ps
+function remove_inits!(ps)
+    filter!(x -> !iscall(x) || !isa(operation(x), Initial), ps)
+    ps
+end
+
+function remove_inits(ps)
+    filter(x -> !iscall(x) || !isa(operation(x), Initial), ps)
+end
+
+
 """
 ```julia
 Base.convert(::Type{<:NonlinearSystem},rs::ReactionSystem)
@@ -555,6 +564,9 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
         as_odes = false, include_zero_odes)
     eqs, us, ps, obs, defs = addconstraints!(eqs, fullrs, ists, ispcs; remove_conserved)
 
+    # remove Initial conditions from parameters
+    remove_inits!(ps)
+
     # Throws a warning if there are differential equations in non-standard format.
     # Next, sets all differential terms to `0`.
     all_differentials_permitted || nonlinear_convert_differentials_check(rs)
@@ -719,8 +731,8 @@ function DiffEqBase.ODEProblem(rs::ReactionSystem, u0, tspan,
     # Handles potential differential algebraic equations (which requires `structural_simplify`).
     if structural_simplify
         (osys = MT.structural_simplify(osys))
-    elseif has_alg_equations(rs)
-        error("The input ReactionSystem has algebraic equations. This requires setting `structural_simplify=true` within `ODEProblem` call.")
+    # elseif has_alg_equations(rs)
+    #     error("The input ReactionSystem has algebraic equations. This requires setting `structural_simplify=true` within `ODEProblem` call.")
     else
         osys = complete(osys)
     end

test/compositional_modelling/component_based_model_creation.jl

@@ -351,7 +351,7 @@ end
 # Test throw error if there are ODE constraints and convert to NonlinearSystem.
 # Note, these can now be created.
 let
-    rn = @reaction_network rn begin
+    rn = @network_component rn begin
         @parameters k1 k2
         (k1, k2), A <--> B
     end
@@ -362,8 +362,7 @@ let
     eqs = [D(C) ~ -b * C + a * A]
     @named osys = ODESystem(eqs, t, [A, C], [a, b])
     rn2 = extend(osys, rn)
-    rn2 = complete(rn2)
-    rnodes = complete(convert(ODESystem, rn2))
+    rnodes = convert(ODESystem, complete(rn2))
 
     # Ensure right number of equations are generated.
     @variables G(t)
@@ -376,12 +375,11 @@ let
     eqs = [0 ~ -a * A + C, 0 ~ -b * C + a * A]
     @named nlsys = NonlinearSystem(eqs, [A, C], [a, b])
     rn2 = extend(nlsys, rn)
-    rn2 = complete(rn2)
-    rnodes = complete(convert(ODESystem, rn2))
-    rnnlsys = complete(convert(NonlinearSystem, rn2))
+    rn2c = complete(rn2)
+    rnodes = complete(convert(ODESystem, rn2c))
+    rnnlsys = complete(convert(NonlinearSystem, rn2c))
     @named nlsys = ODESystem(eqs, t, [A, C], [a, b])
-    rn2 = extend(nlsys, rn)
-    rn2 = complete(rn2)
+    rn2 = complete(extend(nlsys, rn))
     rnodes = convert(ODESystem, rn2)
     rnnlsys = convert(NonlinearSystem, rn2)
 end
@@ -404,7 +402,7 @@ end
 let
     @parameters b
     @species V(t) [isbcspecies = true]
-    rn = @reaction_network begin
+    rn = @network_component begin
         @parameters k
         k/$V, A + B --> C
     end

test/dsl/dsl_advanced_model_construction.jl

@@ -88,15 +88,15 @@ let
         @species A(t) B(t) C(t) D(t)
         k*$AAA, C --> D
     end
-    rn2 = ReactionSystem([Reaction(k*AAA, [C], [D])], t; name=:rn)
+    rn2 = complete(ReactionSystem([Reaction(k*AAA, [C], [D])], t; name=:rn))
     @test rn == rn2
 
     rn = @reaction_network rn begin
         @parameters k
         @species A(t) C(t) D(t)
         k, $AA + C --> D
     end
-    rn2 = ReactionSystem([Reaction(k, [AA,C], [D])], t; name=:rn)
+    rn2 = complete(ReactionSystem([Reaction(k, [AA,C], [D])], t; name=:rn))
     @test rn == rn2
 end
 let
@@ -105,9 +105,9 @@ let
         @parameters k1 k2
         (k1,k2), C + $A2 + $BB + $A2 <--> $BB + $BB
     end
-    rn2 = ReactionSystem([Reaction(k1, [C, A, B], [B], [1,2,1],[2]),
+    rn2 = complete(ReactionSystem([Reaction(k1, [C, A, B], [B], [1,2,1],[2]),
                         Reaction(k2, [B], [C, A, B], [2], [1,2,1])],
-                        t; name=:rn)
+                        t; name=:rn))
     @test rn == rn2
 end
 let
@@ -118,7 +118,7 @@ let
         @parameters α k k1 k2
         α+$kk1*$kk2*$AA, 2*$AA + B --> $AA
     end
-    rn2 = ReactionSystem([Reaction(α+kk1*kk2*AA, [A, B], [A], [2, 1], [1])], t; name=:rn)
+    rn2 = complete(ReactionSystem([Reaction(α+kk1*kk2*AA, [A, B], [A], [2, 1], [1])], t; name=:rn))
     @test rn == rn2
 end
 
@@ -324,13 +324,13 @@ end
 
 # Test species have the right metadata via the DSL.
 let
-    rn = @reaction_network begin
+    rn = @network_component begin
         k, 2*A + B --> C
     end
     @test issetequal(unknowns(rn), species(rn))
     @test all(isspecies, species(rn))
 
-    rn2 = @reaction_network begin
+    rn2 = @network_component begin
         @species A(t) = 1 B(t) = 2 [isbcspecies = true]
         k, A + 2*B --> 2*B
     end
@@ -341,6 +341,7 @@ let
     @named rn2 = extend(osys, rn2)
     rn2 = complete(rn2)
     @test issetequal(unknowns(rn2), species(rn2))
+    rn = complete(rn)
     @test all(isspecies, species(rn))
     @test Catalyst.isbc(ModelingToolkit.value(B))
     @test Catalyst.isbc(ModelingToolkit.value(A)) == false
@@ -363,7 +364,7 @@ let
     @species (X(t))[1:2] Y(t) C(t)
     rx = Reaction(k[1]*a+k[2], [X[1], X[2]], [Y, C], [1, V[1]], [V[2] * W, B])
     @named arrtest = ReactionSystem([rx], t)
-    @test arrtest == rn
+    @test complete(arrtest) == rn
 
     rn = @reaction_network twostate begin
         @parameters k[1:2]
@@ -376,5 +377,5 @@ let
     rx1 = Reaction(k[1], [X[1]], [X[2]])
     rx2 = Reaction(k[2], [X[2]], [X[1]])
     @named twostate = ReactionSystem([rx1, rx2], t)
-    @test twostate == rn
+    @test complete(twostate) == rn
 end

test/dsl/dsl_basic_model_construction.jl

@@ -23,6 +23,10 @@ function unpacksys(sys)
     get_eqs(sys), get_iv(sys), get_unknowns(sys), get_ps(sys), nameof(sys), get_systems(sys)
 end
 
+function remove_inits(ps)
+    filter(x -> !iscall(x) || !isa(operation(x), Initial), ps)
+end
+
 opname(x) = iscall(x) ? nameof(operation(x)) : nameof(x)
 alleq(xs, ys) = all(isequal(x, y) for (x, y) in zip(xs, ys))
 
@@ -49,8 +53,8 @@ function basic_test(rn, N, unknowns_syms, p_syms)
     @test length(unknowns) == length(unknowns_syms)
     @test issetequal(map(opname, unknowns), unknowns_syms)
     @test all_reactants(eqs) == Set(unknowns_syms)
-    @test length(ps) == length(p_syms)
-    @test issetequal(map(opname, ps), p_syms)
+    @test length(remove_inits(ps)) == length(p_syms)
+    @test issetequal(map(opname, remove_inits(ps)), p_syms)
 end
 
 ## Run Tests ###