updates with working conservation laws

Author: TorkelE

Project.toml

@@ -44,7 +44,7 @@ HomotopyContinuation = "2.9"
 LaTeXStrings = "1.3.0"
 Latexify = "0.14, 0.15, 0.16"
 MacroTools = "0.5.5"
-ModelingToolkit = "9.6"
+ModelingToolkit = "9.7.1"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"

test/extensions/homotopy_continuation.jl

@@ -14,7 +14,8 @@ include("../test_functions.jl")
 # Tests for Symbolics initial condition input.
 # Tests for different types (Symbol/Symbolics) for parameters and initial conditions.
 # Tests that attempts to find steady states of system with conservation laws, while u0 is not provided, gives an error.
-@test_broken let # Requires https://github.com/SciML/ModelingToolkit.jl/issues/2566 to be fixed
+let 
+    # Creates the model.
     rs = @reaction_network begin
         (k1,k2), X1 <--> X2
         (k3,k4), 2X2 + X3 <--> X2_2X3
@@ -23,10 +24,12 @@ include("../test_functions.jl")
     ps = [k1 => 1.0, k2 => 2.0, k3 => 2.0, k4 => 2.0]
     u0 = [:X1 => 2.0, :X2 => 2.0, :X3 => 2.0, :X2_2X3 => 2.0]
 
-    hc_ss = hc_steady_states(rs, ps; u0=u0, show_progress=false)[1]
-    f = ODEFunction(convert(ODESystem, rs))
-    @test f(hc_ss, last.(ps), 0.0)[1] == 0.0
+    # Computes the single steady state, checks that when given to the ODE rhs, all are evaluated to 0.
+    hc_ss = hc_steady_states(rs, ps; u0=u0, show_progress=false)
+    hc_ss = Pair.(unknowns(rs), hc_ss[1])
+    @test maximum(abs.(f_eval(rs, hc_ss, ps, 0.0))) ≈ 0.0 atol=1e-12
 
+    # Checks that not giving a `u0` argument yields an error for systems with conservation laws.
     @test_throws Exception hc_steady_states(rs, ps; show_progress=false)
 end
 
@@ -55,7 +58,7 @@ end
 # Tests correctness in presence of default values.
 # Tests where some default values are overwritten with other values.
 # Tests where input ps/u0 are tuples with mixed types.
-@test_broken let # Requires https://github.com/SciML/ModelingToolkit.jl/issues/2566 to be fixed
+let
     rs_1 = @reaction_network begin
         @parameters kX1=1.0 kX2=2.0 kY1=12345.0 
         @species X1(t)=0.1 X2(t)=0.2 Y1(t)=12345.0

test/miscellaneous_tests/api.jl

@@ -533,7 +533,6 @@ end
 # Test conservation law elimination.
 # Conservation laws currently broken (you get stuck in an infinite loop in MTK or something).
 let
-    return (@test_broken false)
     rn = @reaction_network begin
         (k1, k2), A + B <--> C
         (m1, m2), D <--> E

test/miscellaneous_tests/nonlinear_solve.jl

@@ -77,8 +77,6 @@ end
 # Checks for system with conservation laws.
 # Checks using interfacing with output solution.
 let 
-    # Conservation laws currently broken (you get stuck in an infinite loop in MTK or something).
-    return (@test_broken false)
     # Creates steady state network, unpack the parameter values.
     steady_state_network_3 = @reaction_network begin
         (p,d), 0 <--> X
@@ -98,6 +96,6 @@ let
     sol2 = solve(nl_prob_2, DynamicSS(Rosenbrock23()); abstol=1e-12, reltol=1e-12)
 
     # Checks output using the ODE's drift function
-    @test f_eval([:X => sol1[X], :Y => sol1[Y], :Y2 => sol1[Y2], :XY2 => sol1[XY2]], p, 0.0) ≈ [0.0, 0.0, 0.0, 0.0] atol=1e-10
-    @test f_eval([:X => sol2[X], :Y => sol2[Y], :Y2 => sol2[Y2], :XY2 => sol2[XY2]], p, 0.0) ≈ [0.0, 0.0, 0.0, 0.0] atol=1e-10
+    @test f_eval(steady_state_network_3, [:X => sol1[X], :Y => sol1[Y], :Y2 => sol1[Y2], :XY2 => sol1[XY2]], p, 0.0) ≈ [0.0, 0.0, 0.0, 0.0] atol=1e-10
+    @test f_eval(steady_state_network_3, [:X => sol2[X], :Y => sol2[Y], :Y2 => sol2[Y2], :XY2 => sol2[XY2]], p, 0.0) ≈ [0.0, 0.0, 0.0, 0.0] atol=1e-10
 end