Handle test errors from MTKv9 update (including `ReactionSystem` completeness) #784

TorkelE · 2024-02-13T02:44:39Z

For MTK v9 all inputs to various systems have to be complete. Which means we must either add complete(::ReactionSystem) wherever these are required, or make them complete by default. Since this can be sorted out on the v14 branch even before MTK v9 I figured we could do it in a separate branch.

Currently, I make ReactionSystems complete by default. which can be changed by a DSL option/optional argument:

@reaction_network begin
    @incompelte
    (p,d), 0 <--> X
end
ReactionSystem(rxs, t; complete = false)

However, we can flip the default really easily (while also renaming @incompelte to @compelte

Deciding which to go for is still a problem. I prefer this, since it only introduces the concept of completeness when it is actually required. I think it is useful to follow how MTK does stuff, but for something which has a major impact on Catalyst (and is not just notation), I think we should do what works best, independently of MTK (whatever the case, I hope to add a "differences to MTK" section somewhere in the docs).

docs/pages.jl

TorkelE · 2024-03-23T01:15:41Z

docs/src/catalyst_applications/advanced_simulations.md

-scaling parameter*. First, we simulate a simple two-state CRN model using the
+scale the magnitude of the noise terms. Here, each reaction of the system generates a separate noise term in the CLE. If you require identical scaling for all reactions, the `@default_noise_scaling` option can be used. Else, you can supply a `noise_scaling` metadata for each individual reaction, describing how to scale the noise for that reaction.
+
+We begin with considering the first approach. First, we simulate a simple two-state CRN model using the


Additions in this file is because this PR builds upon the noise scaling one, so showing its changes as changes here. Will disappear once that gets merged.

docs/src/catalyst_functionality/compositional_modeling.md

TorkelE · 2024-03-23T01:17:36Z

docs/src/catalyst_functionality/compositional_modeling.md

 ## Compositional modeling tooling
 Catalyst supports two ModelingToolkit interfaces for composing multiple
 [`ReactionSystem`](@ref)s together into a full model. The first mechanism for
 extending a system is the `extend` command
 ```@example ex1
 using Catalyst
 basern = @reaction_network rn1 begin
+  @incomplete


Currently, the @incompelte option marks that a reaction system (created via the DSL) is not complete when created. Can change how to manage this later.

src/Catalyst.jl

src/networkapi.jl

TorkelE · 2024-03-23T01:19:19Z

src/reaction_network.jl

@@ -123,7 +123,7 @@ const forbidden_variables_error = let
 end

 # Declares the keys used for various options.
-const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables)
+const option_keys = (:species, :parameters, :variables, :ivs, :compounds, :observables, :default_noise_scaling, :incomplete)


Stuff that follows (in this file) are either the noise scaling (because this builds on that PR), or the compelteness.

test/dsl/dsl_model_construction.jl

TorkelE · 2024-03-23T01:21:22Z

test/dsl/dsl_model_construction.jl

-            @test all(abs.(f1.jac(u0, p, t) .≈ f2.jac(u0, p, t)))
-            @test all(abs.(g1(u0, p, t) .≈ g2(u0, p, t)))
+
+            @test f_eval(networks[1], u0_1, p_1, t) ≈ f_eval(networks[2], u0_2, p_2, t)


Uses the new f_eval function to evaluate the ODESystem f function.

Alsmot always, like here, f_eval is only called once for each network, do no use using a in-place version. Because most networks have different du vector sizes, reusing it across different systems is hard.

TorkelE · 2024-03-23T01:22:43Z

test/dsl/dsl_options.jl

+    rn2 = complete(rn2)
+    @test ModelingToolkit.iscomplete(rn2)
+end


Test for the new @incomplete functionality.

(need changing after final decision on completeness)

test/extensions/bifurcation_kit.jl

TorkelE · 2024-03-23T01:24:01Z

test/extensions/homotopy_continuation.jl

 ### Run Tests ###

 # Tests for network without conservation laws.
 # Tests for Symbol parameter input.
 # 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.
-let 
+@test_broken let # Requires /~https://github.com/SciML/ModelingToolkit.jl/issues/2566 to be fixed


One of the tests that are currently broken. Basically, we cannot do homotopy continuation on systems with conservation laws until that issue gets sorted out.

TorkelE · 2024-03-23T01:24:18Z

test/extensions/homotopy_continuation.jl

@@ -50,7 +53,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.
-let    
+@test_broken let # Requires /~https://github.com/SciML/ModelingToolkit.jl/issues/2566 to be fixed


Same problem as the one above.

TorkelE · 2024-03-23T01:25:48Z

test/miscellaneous_tests/nonlinear_solve.jl

-    @test isapprox(sol2[1], p[1] / p[2]; atol=1e-10)
-    @test isapprox(sol2[2]^3 / factorial(3), p[3] / p[4]; atol=1e-10)
-    @test isapprox(sol2[3], (p[5] * p[2]) / (p[6] * p[1]); atol=1e-10)
+    @test sol1[:X1] ≈ nlprob.ps[:k1] / nlprob.ps[:k2] atol=1e-10


Uses symbolic (instead of integer) indexing.

Loads of test where also written before I was familiar with nice ways to use ≈, @test, and atol, so I rewrote quite a few tests using this while I was at it.

TorkelE · 2024-03-23T01:26:46Z

test/miscellaneous_tests/nonlinear_solve.jl

 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)


Another conservation law issue, this one only appears when a XProblem is created (because these creates a MTKParameters object, which is the real culprit, not sure what the problem in these are though).

test/miscellaneous_tests/symbolic_stoichiometry.jl

TorkelE · 2024-03-23T01:28:23Z

test/miscellaneous_tests/symbolic_stoichiometry.jl

+# Gives an error. I can fix so there is not, but unsure what part of the code is what we actually are
+# testing for, and what we are not. Sam, if you have any opinions on what to preserve in this test or not,
+# tell me before I try to rewrite. 
+@test_broken let


Again, this is also broken,

test/miscellaneous_tests/units.jl

TorkelE · 2024-03-23T01:29:17Z

test/model_simulation/make_jacobian.jl

-                          -2*p[3]*u[2]*u[1] -1-p[3] * u[1] * u[1] 1;
-                          2*p[3]*u[2]*u[1] p[3]*u[1]*u[1] -1.0]
-        @test all(abs.(test_jac .- real_jac) .< 1e-9)
+        u = Dict(rnd_u0(jacobian_network_2, rng; factor))


Uses Dict to simplify evaluating real_jac in 3 lines from here.

TorkelE · 2024-03-23T01:31:48Z

test/model_simulation/simulate_ODEs.jl

+            p_1 = rnd_ps(networks[1], rng; factor)
+            (i == 3) && (p_1 =rnd_ps_Int64(networks[1], rng))
+            u0_2 = last.(u0_1)
+            p_2 = last.(p_1)


This depends on ReactionSystem not doing some weird reordering (which we do not want to promise in the future). However, here we loop through quite a few networks and just generate parameters. If we'd want to make order not matter here, we'd have to write out the vectors for each case manually, and handle them separately. I can see how to do this, but non-trivial amount of work. I'd rather do that at some later overhaul of the tests than right now (when we have tons of stuff to fix).

TorkelE · 2024-03-23T01:33:54Z

test/reactionsystem_structure/higher_order_reactions.jl

+# Currently fails because binomial only takes Int input (and X is Float64).
+# I see several solutions, but depends on whether we allow species to be specified as Int64.
+# I have marked this one as broken for now.
+@test_broken let 


This test will remain broken until we decide how to handle non-Float64 definitions of species. Which will be messy, because if we do we have to rewrite the @compound stuff to permit that as well.

Once sorted out, we can fix these tests using the best solution (which depends on whether that is allowed or not).

TorkelE · 2024-03-23T01:34:55Z

test/reactionsystem_structure/reactionsystem.jl

-    sf = SDEFunction{false}(sdesys_noise_scaling, unknowns(rs),
-                            parameters(sdesys_noise_scaling))
-    G2 = sf.g(u, p, t)
-    @test norm(G - G2) < 100 * eps()


These were moved in the noise scaling PR to the SDe simulation file (so not removed).

TorkelE · 2024-03-29T21:08:19Z

test/visualization/latexify.jl

-    end
-    @test latexify(rn; mathrm = true) != latexify(rn; mathrm = false)
-end
-


Now line 207

TorkelE · 2024-03-29T21:41:22Z

test/dsl/dsl_basics.jl

@@ -110,7 +117,8 @@ let
    @test rn == rn2
 end

-@testset "make_reaction_system can be called from another module" begin


Chanegd the @testset block to let .. end block like for all of our other tests. Not sure if there is a reason for this here, since it is all already in a testset?

test/model_simulation/u0_n_parameter_inputs.jl

isaacsas

There may be some older comments I missed deleting, so ask if any seem out of date or already settled. Sorry about that!

docs/src/catalyst_functionality/compositional_modeling.md

src/networkapi.jl

src/reaction_network.jl

src/reactionsystem.jl

isaacsas · 2024-04-02T17:36:03Z

test/reactionsystem_structure/higher_order_reactions.jl

+# Tests that Jump Systems are correct (by comparing to network with manually written higher order rates). 
+# Currently fails for reason I do not understand. Likely reason similar to the weird case in the jump tests.
+# Spent loads of time trying to figure out, the best I can get to is that it seems like the rng/seed is
+# not fully reproducible.
+let 
+    @test_broken false
+    # Declares a JumpSystem manually. Would have used Catalyst again, but `binomial` still errors when
+    # called on symbolics. For reference, here is the network as it would be created using Catalyst. 
+
+    # higher_order_network_alt2 = @reaction_network begin
+    #     p, ∅ ⟼ X1
+    #     r1 * binomial(X1, 2), 2X1 ⟾ 3X2
+    #     mm(X1, r2, K) * binomial(X2, 3), 3X2 ⟾ X3 + 2X4
+    #     r3 * binomial(X3, 1) * binomial(X4, 2), X3 + 2X4 ⟾ 3X5 + 3X6
+    #     r4 * X2 * binomial(X5, 3) * binomial(X6, 3), 3X5 + 3X6 ⟾ 3X5 + 2X7 + 4X8
+    #     r5 * binomial(X5, 3) * binomial(X7, 2) * binomial(X8, 4), 3X5 + 2X7 + 4X8 ⟾ 10X9
+    #     r6 * binomial(X9, 10), 10X9 ⟾ X10
+    #     d * binomial(X10, 2), 2X10 ⟾ ∅
+    # end


This passed before. We need to keep it working now. It is a valuable test that is pretty comprehensive...

Please open issues on MTK/Symbolics if something like binomial is now broken that used to work.

I did, I even wrote a PR to DynamicQuantities to get it solved as well, which got something, but not all the way. I will continue pushing for this though. Not sure when it will be fixed though. The @test_broken false can be removed.

If this issue is DynamicQuantites I would prefer we drop units over dropping this test and the use of binomials. We can then add units back if/when they update it.

isaacsas · 2024-04-02T17:40:24Z

test/reactionsystem_structure/higher_order_reactions.jl

+    rate1(u, p, t) = p[1]
+    rate2(u, p, t) = p[2] * binomial(u[1], 2)
+    rate3(u, p, t) = mm(u[1], p[3], p[4]) * binomial(u[2], 3)
+    rate4(u, p, t) = p[5] * binomial(u[3], 1) * binomial(u[4], 2)
+    rate5(u, p, t) = p[6] * u[2] * binomial(u[5], 3) * binomial(u[6], 3)
+    rate6(u, p, t) = p[7] * binomial(u[5], 3) * binomial(u[7], 2) * binomial(u[8], 4)
+    rate7(u, p, t) = p[8] * binomial(u[9], 10)
+    rate8(u, p, t) = p[9] * binomial(u[10], 2)
+
+    affect1!(int) = (int.u[1] += 1)
+    affect2!(int) = (int.u[1] -= 2; int.u[2] += 3;)
+    affect3!(int) = (int.u[2] -= 3; int.u[3] += 1; int.u[4] += 2;)
+    affect4!(int) = (int.u[3] -= 1; int.u[4] -= 2; int.u[5] += 3; int.u[6] += 3;)
+    affect5!(int) = (int.u[5] -= 3; int.u[6] -= 3; int.u[5] += 3; int.u[7] += 2; int.u[8] += 4;)
+    affect6!(int) = (int.u[5] -= 3; int.u[7] -= 2; int.u[8] -= 4; int.u[9] += 10;)
+    affect7!(int) = (int.u[9] -= 10; int.u[10] += 1;)
+    affect8!(int) = (int.u[10] -= 2;)
+
+    higher_order_network_alt2 = ConstantRateJump.([rate1, rate2, rate3, rate4, rate5, rate6, rate7, rate8], 
+                                [affect1!, affect2!, affect3!, affect4!, affect5!, affect6!, affect7!, affect8!])


This is a nice test too, so we should get the old Catalyst version fixed and keep your new JumpProcesses version too and compare against both going forward.

isaacsas · 2024-04-02T17:40:57Z

test/reactionsystem_structure/reactions.jl

+    @test_throws Exception Catalyst.get_noise_scaling(r1)
    @test isequal(Catalyst.get_noise_scaling(r2), η)


Do you have a has_noise_scaling function to test too (or is it tested somewhere else)?

Good catch, yes, I will add tests for this one too.

isaacsas · 2024-04-02T17:43:08Z

test/reactionsystem_structure/reactionsystem.jl

-                            parameters(sdesys_noise_scaling))
-    G2 = sf.g(u, p, t)
-    @test norm(G - G2) < 100 * eps()
+    @test_broken let # The next line throws an error.


Just trying to test the ODE rhs function gives the same values as the NonlinearProblem. We can switch to just use the function within the NonlinearProblem (I think maybe it wasn't implemented/working when we originally wrote these tests).

test/reactionsystem_structure/reactionsystem.jl

Co-authored-by: Sam Isaacson <isaacsas@users.noreply.github.com>

TorkelE changed the base branch from master to Catalyst_version_14 February 13, 2024 02:45

isaacsas mentioned this pull request Feb 28, 2024

Catalyst V14 TODO #775

Closed

47 tasks

TorkelE force-pushed the reactionsystem_completeness branch 2 times, most recently from 2614a6c to a8ed387 Compare March 19, 2024 00:16

TorkelE changed the title ~~Handle ReactionSystem completeness for MTK v9 update~~ Handle MTKv9 errors (including ReactionSystem completeness) Mar 20, 2024

TorkelE changed the title ~~Handle MTKv9 errors (including ReactionSystem completeness)~~ Handle test errors from MTKv9 update (including ReactionSystem completeness) Mar 20, 2024

TorkelE force-pushed the reactionsystem_completeness branch from dba4535 to ed4ca76 Compare March 22, 2024 21:25