Merge pull request #956 from SciML/new_default_prob_update_tests

[v14 - Awaiting Catalyst update] New default prob update tests

test/network_analysis/conservation_laws.jl

@@ -230,6 +230,55 @@ let
     @test all(sol2[osys.X1 + osys.X2] .== 4.0)
 end
 
+# Tests system problem updating when conservation laws are eliminated.
+# Checks that the correct values are used after the conservation law species are updated.
+# Here is an issue related to the broken tests: https://github.com/SciML/Catalyst.jl/issues/952
+let
+    # Create model and fetch the conservation parameter (Γ).
+    t = default_t()
+    @parameters k1 k2
+    @species X1(t) X2(t)
+    rxs = [
+        Reaction(k1, [X1], [X2]),
+        Reaction(k2, [X2], [X1])
+    ]
+    @named rs = ReactionSystem(rxs, t)
+    osys = convert(ODESystem, complete(rs); remove_conserved = true, remove_conserved_warn = false)
+    osys = complete(osys)
+    @unpack Γ = osys
+
+    # Creates an `ODEProblem`.
+    u0 = [X1 => 1.0, X2 => 2.0]
+    ps = [k1 => 0.1, k2 => 0.2]
+    oprob = ODEProblem(osys, u0, (0.0, 1.0), ps)
+
+    # Check `ODEProblem` content.
+    @test oprob[X1] == 1.0
+    @test oprob[X2] == 2.0
+    @test oprob.ps[k1] == 0.1
+    @test oprob.ps[k2] == 0.2
+    @test oprob.ps[Γ[1]] == 3.0
+
+    # Currently, any kind of updating of species or the conservation parameter(s) is not possible.
+
+    # Update problem parameters using `remake`.
+    oprob_new = remake(oprob; p = [k1 => 0.3, k2 => 0.4])
+    @test oprob_new.ps[k1] == 0.3
+    @test oprob_new.ps[k2] == 0.4
+    integrator = init(oprob_new, Tsit5())
+    @test integrator.ps[k1] == 0.3
+    @test integrator.ps[k2] == 0.4
+
+    # Update problem parameters using direct indexing.
+    oprob.ps[k1] = 0.5
+    oprob.ps[k2] = 0.6
+    @test oprob.ps[k1] == 0.5
+    @test oprob.ps[k2] == 0.6
+    integrator = init(oprob, Tsit5())
+    @test integrator.ps[k1] == 0.5
+    @test integrator.ps[k2] == 0.6
+end
+
 ### Other Tests ###
 
 # Checks that `JumpSystem`s with conservation laws cannot be generated.
@@ -288,4 +337,27 @@ let
         @test_nowarn XProblem(rn, u0, ps; remove_conserved_warn = false)
         @test_nowarn XProblem(rn, u0, ps; remove_conserved = true, remove_conserved_warn = false)
     end
-end
+end
+
+# Conservation law simulations for vectorised species.
+let 
+    # Prepares the model.
+    t = default_t()
+    @species X(t)[1:2]
+    @parameters k[1:2]
+    rxs = [
+        Reaction(k[1], [X[1]], [X[2]]),
+        Reaction(k[2], [X[2]], [X[1]])
+    ]
+    @named rs = ReactionSystem(rxs, t)
+    rs = complete(rs)
+
+    # Checks that simulation reaches known equilibrium.
+    @test_broken false # Currently broken on MTK .
+    # u0 = [:X => [3.0, 9.0]]
+    # ps = [:k => [1.0, 2.0]]
+    # oprob = ODEProblem(rs, u0, (0.0, 1000.0), ps; remove_conserved = true)
+    # sol = solve(oprob, Vern7())
+    # @test sol[X[1]][end] ≈ 8.0
+    # @test sol[X[2]][end] ≈ 4.0
+end

test/reactionsystem_core/reactionsystem.jl

@@ -262,6 +262,70 @@ let
     end
 end
 
+### Nich Model Declarations ###
+
+# Checks model with vector species and parameters.
+# Checks that it works for programmatic/dsl-based modelling.
+# Checks that all forms of model input (parameter/initial condition and vector/non-vector) are
+# handled properly.
+let 
+    # Declares programmatic model.
+    @parameters p[1:2] k d1 d2
+    @species (X(t))[1:2] Y1(t) Y2(t)
+    rxs = [
+        Reaction(p[1], [], [X[1]]),
+        Reaction(p[2], [], [X[2]]),
+        Reaction(k, [X[1]], [Y1]),
+        Reaction(k, [X[2]], [Y2]),
+        Reaction(d1, [Y1], []),
+        Reaction(d2, [Y2], []),
+    ]
+    rs_prog = complete(ReactionSystem(rxs, t; name = :rs))
+
+    # Declares DSL-based model.
+    rs_dsl = @reaction_network rs begin
+        @parameters p[1:2] k d1 d2 
+        @species (X(t))[1:2] Y1(t) Y2(t)
+        (p[1],p[2]), 0 --> (X[1],X[2])
+        k, (X[1],X[2]) --> (Y1,Y2)
+        (d1,d2), (Y1,Y2) --> 0
+    end
+
+    # Checks equivalence.
+    rs_dsl == rs_prog
+
+    # Creates all possible initial conditions and parameter values.
+    u0_alts = [
+        [X => [2.0, 5.0], Y1 => 0.2, Y2 => 0.5],
+        [X[1] => 2.0, X[2] => 5.0, Y1 => 0.2, Y2 => 0.5],
+        [rs_dsl.X => [2.0, 5.0], rs_dsl.Y1 => 0.2, rs_dsl.Y2 => 0.5],
+        [rs_dsl.X[1] => 2.0, X[2] => 5.0, rs_dsl.Y1 => 0.2, rs_dsl.Y2 => 0.5],
+        [:X => [2.0, 5.0], :Y1 => 0.2, :Y2 => 0.5]
+    ]
+    ps_alts = [
+        [p => [1.0, 10.0], d1 => 5.0, d2 => 4.0, k => 2.0],
+        [p[1] => 1.0, p[2] => 10.0, d1 => 5.0, d2 => 4.0, k => 2.0],
+        [rs_dsl.p => [1.0, 10.0], rs_dsl.d1 => 5.0, rs_dsl.d2 => 4.0, rs_dsl.k => 2.0],
+        [rs_dsl.p[1] => 1.0, p[2] => 10.0, rs_dsl.d1 => 5.0, rs_dsl.d2 => 4.0, rs_dsl.k => 2.0],
+        [:p => [1.0, 10.0], :d1 => 5.0, :d2 => 4.0, :k => 2.0]
+    ]
+
+    # Loops through all inputs and check that the correct steady state is reached
+    # Target steady state: (X1, X2, Y1, Y2) = (p1/k, p2/k, p1/d1, p2/d2).
+    # Technically only one model needs to be check. However, "equivalent" models in MTK can still
+    # have slight differences, so checking for both here to be certain.
+    for rs in [rs_prog, rs_dsl]
+        oprob = ODEProblem(rs, u0_alts[1], (0.0, 10000.), ps_alts[1])
+        @test_broken false # Cannot currently `remake` this problem/
+        # for rs in [rs_prog, rs_dsl], u0 in u0_alts, p in ps_alts
+        #     oprob_remade = remake(oprob; u0, p)
+        #     sol = solve(oprob_remade, Vern7(); abstol = 1e-8, reltol = 1e-8)
+        #     @test sol[end] ≈ [0.5, 5.0, 0.2, 2.5]
+        # end
+    end
+end
+
+### Other Tests ###
 
 ### Test Show ###
 
@@ -751,4 +815,4 @@ end
 # the code might need adaptation to take the updated reaction system into account.
 let
     @test_nowarn Catalyst.reactionsystem_uptodate_check()
-end
+end

test/upstream/mtk_problem_inputs.jl

@@ -158,7 +158,6 @@ let
     end
 end
 
-
 ### Checks Errors On Faulty Inputs ###
 
 # Checks various erroneous problem inputs, ensuring that these throw errors.