save progress

Author: TorkelE

test/meta/mtk_problem_inputs.jl

@@ -158,6 +158,188 @@ let
     end
 end
 
+### Vector Species/Parameters Tests ###
+
+begin 
+    # Declares the model (with vector species/parameters, with/without default values, and observables).
+    @species X(t)[1:2] Y(t)[1:2] = [10.0, 20.0] XY(t)[1:2]
+    @parameters p[1:2] d[1:2] = [0.2, 0.5]
+    rxs = [
+        Reaction(p[1], [], [X[1]]),
+        Reaction(p[2], [], [X[2]]),
+        Reaction(d[1], [X[1]], []),
+        Reaction(d[2], [X[2]], []),
+        Reaction(p[1], [], [Y[1]]),
+        Reaction(p[2], [], [Y[2]]),
+        Reaction(d[1], [Y[1]], []),
+        Reaction(d[2], [Y[2]], [])
+    ]
+    observed = [XY[1] ~ X[1] + Y[1], XY[2] ~ X[2] + Y[2]]
+    @named model_vec = ReactionSystem(rxs, t; observed)
+    model_vec = complete(model_vec)
+
+    # Declares various u0 versions.
+    u0_alts_vec = [
+        # Vectors not providing default values.
+        [X => [1.0, 2.0]],
+        [X[1] => 1.0, X[2] => 2.0],
+        [model_vec.X => [1.0, 2.0]],
+        [model_vec.X[1] => 1.0, model_vec.X[2] => 2.0],
+        [:X => [1.0, 2.0]],
+        # Vectors providing default values.
+        [X => [1.0, 2.0], Y => [10.0, 20.0]],
+        [X[1] => 1.0, X[2] => 2.0, Y[1] => 10.0, Y[2] => 20.0],
+        [model_vec.X => [1.0, 2.0], model_vec.Y => [10.0, 20.0]],
+        [model_vec.X[1] => 1.0, model_vec.X[2] => 2.0, model_vec.Y[1] => 10.0, model_vec.Y[2] => 20.0],
+        [:X => [1.0, 2.0], :Y => [10.0, 20.0]],
+        # Dicts not providing default values.
+        Dict([X => [1.0, 2.0]]),
+        Dict([X[1] => 1.0, X[2] => 2.0]),
+        Dict([model_vec.X => [1.0, 2.0]]),
+        Dict([model_vec.X[1] => 1.0, model_vec.X[2] => 2.0]),
+        Dict([:X => [1.0, 2.0]]),
+        # Dicts providing default values.
+        Dict([X => [1.0, 2.0], Y => [10.0, 20.0]]),
+        Dict([X[1] => 1.0, X[2] => 2.0, Y[1] => 10.0, Y[2] => 20.0]),
+        Dict([model_vec.X => [1.0, 2.0], model_vec.Y => [10.0, 20.0]]),
+        Dict([model_vec.X[1] => 1.0, model_vec.X[2] => 2.0, model_vec.Y[1] => 10.0, model_vec.Y[2] => 20.0]),
+        Dict([:X => [1.0, 2.0], :Y => [10.0, 20.0]]),
+        # Tuples not providing default values.
+        (X => [1.0, 2.0]),
+        (X[1] => 1.0, X[2] => 2.0),
+        (model_vec.X => [1.0, 2.0]),
+        (model_vec.X[1] => 1.0, model_vec.X[2] => 2.0),
+        (:X => [1.0, 2.0]),
+        # Tuples providing default values.
+        (X => [1.0, 2.0], Y => [10.0, 20.0]),
+        (X[1] => 1.0, X[2] => 2.0, Y[1] => 10.0, Y[2] => 20.0),
+        (model_vec.X => [1.0, 2.0], model_vec.Y => [10.0, 20.0]),
+        (model_vec.X[1] => 1.0, model_vec.X[2] => 2.0, model_vec.Y[1] => 10.0, model_vec.Y[2] => 20.0),
+        (:X => [1.0, 2.0], :Y => [10.0, 20.0]),
+    ]
+
+    # Declares various ps versions.
+    p_alts_vec = [
+        # Vectors not providing default values.
+        [p => [1.0, 2.0]],
+        [p[1] => 1.0, p[2] => 2.0],
+        [model_vec.p => [1.0, 2.0]],
+        [model_vec.p[1] => 1.0, model_vec.p[2] => 2.0],
+        [:p => [1.0, 2.0]],
+        # Vectors providing default values.
+        [p => [4.0, 5.0], d => [0.2, 0.5]],
+        [p[1] => 4.0, p[2] => 5.0, d[1] => 10.0, d[2] => 20.0],
+        [model_vec.p => [4.0, 5.0], model_vec.d => [0.2, 0.5]],
+        [model_vec.p[1] => 4.0, model_vec.p[2] => 5.0, model_vec.d[1] => 10.0, model_vec.d[2] => 20.0],
+        [:p => [4.0, 5.0], :d => [0.2, 0.5]],
+        # Dicts not providing default values.
+        Dict([p => [1.0, 2.0]]),
+        Dict([p[1] => 1.0, p[2] => 2.0]),
+        Dict([model_vec.p => [1.0, 2.0]]),
+        Dict([model_vec.p[1] => 1.0, model_vec.p[2] => 2.0]),
+        Dict([:p => [1.0, 2.0]]),
+        # Dicts providing default values.
+        Dict([p => [4.0, 5.0], d => [0.2, 0.5]]),
+        Dict([p[1] => 4.0, p[2] => 5.0, d[1] => 10.0, d[2] => 20.0]),
+        Dict([model_vec.p => [4.0, 5.0], model_vec.d => [0.2, 0.5]]),
+        Dict([model_vec.p[1] => 4.0, model_vec.p[2] => 5.0, model_vec.d[1] => 10.0, model_vec.d[2] => 20.0]),
+        Dict([:p => [4.0, 5.0], :d => [0.2, 0.5]]),
+        # Tuples not providing default values.
+        (p => [1.0, 2.0]),
+        (p[1] => 1.0, p[2] => 2.0),
+        (model_vec.p => [1.0, 2.0]),
+        (model_vec.p[1] => 1.0, model_vec.p[2] => 2.0),
+        (:p => [1.0, 2.0]),
+        # Tuples providing default values.
+        (p => [4.0, 5.0], d => [0.2, 0.5]),
+        (p[1] => 4.0, p[2] => 5.0, d[1] => 10.0, d[2] => 20.0),
+        (model_vec.p => [4.0, 5.0], model_vec.d => [0.2, 0.5]),
+        (model_vec.p[1] => 4.0, model_vec.p[2] => 5.0, model_vec.d[1] => 10.0, model_vec.d[2] => 20.0),
+        (:p => [4.0, 5.0], :d => [0.2, 0.5]),
+    ]
+
+    # Declares a timespan.
+    tspan = (0.0, 10.0)
+end
+
+# Perform ODE simulations (singular and ensemble).
+let 
+    # Creates normal and ensemble problems.
+    base_oprob = ODEProblem(model_vec, u0_alts_vec[1], tspan, p_alts_vec[1])
+    base_sol = solve(base_oprob, Tsit5(); saveat = 1.0)
+    base_eprob = EnsembleProblem(base_oprob)
+    base_esol = solve(base_eprob, Tsit5(); trajectories = 2, saveat = 1.0)
+
+    # Simulates problems for all input types, checking that identical solutions are found.
+    for u0 in u0_alts_vec, p in p_alts_vec
+        oprob = remake(base_oprob; u0, p)
+        @test base_sol == solve(oprob, Tsit5(); saveat = 1.0)
+        eprob = remake(base_eprob; u0, p)
+        @test base_esol == solve(eprob, Tsit5(); trajectories = 2, saveat = 1.0)
+    end
+end
+
+# Perform SDE simulations (singular and ensemble).
+let 
+    # Creates normal and ensemble problems.
+    base_sprob = SDEProblem(model, u0_alts_vec[1], tspan, p_alts_vec[1])
+    base_sol = solve(base_sprob, ImplicitEM(); seed, saveat = 1.0)
+    base_eprob = EnsembleProblem(base_sprob)
+    base_esol = solve(base_eprob, ImplicitEM(); seed, trajectories = 2, saveat = 1.0)
+
+    # Simulates problems for all input types, checking that identical solutions are found.
+    for u0 in u0_alts_vec, p in p_alts_vec
+        sprob = remake(base_sprob; u0, p)
+        @test base_sol == solve(sprob, ImplicitEM(); seed, saveat = 1.0)
+        eprob = remake(base_eprob; u0, p)
+        @test base_esol == solve(eprob, ImplicitEM(); seed, trajectories = 2, saveat = 1.0)
+    end
+end
+
+# Perform jump simulations (singular and ensemble).
+let 
+    # Creates normal and ensemble problems.
+    base_dprob = DiscreteProblem(model, u0_alts_vec[1], tspan, p_alts_vec[1])
+    base_jprob = JumpProblem(model, base_dprob, Direct(); rng)
+    base_sol = solve(base_jprob, SSAStepper(); seed, saveat = 1.0)
+    base_eprob = EnsembleProblem(base_jprob)
+    base_esol = solve(base_eprob, SSAStepper(); seed, trajectories = 2, saveat = 1.0)
+
+    # Simulates problems for all input types, checking that identical solutions are found.
+    for u0 in u0_alts_vec, p in p_alts_vec
+        jprob = remake(base_jprob; u0, p)
+        @test base_sol == solve(base_jprob, SSAStepper(); seed, saveat = 1.0)
+        eprob = remake(base_eprob; u0, p)
+        @test base_esol == solve(eprob, SSAStepper(); seed, trajectories = 2, saveat = 1.0)
+    end
+end
+
+# Solves a nonlinear problem (EnsembleProblems are not possible for these).
+let
+    base_nlprob = NonlinearProblem(model, u0_alts_vec[1], p_alts_vec[1])
+    base_sol = solve(base_nlprob, NewtonRaphson())
+    for u0 in u0_alts_vec, p in p_alts_vec
+        nlprob = remake(base_nlprob; u0, p)
+        @test base_sol == solve(nlprob, NewtonRaphson())
+    end
+end
+
+# Perform steady state simulations (singular and ensemble).
+let 
+    # Creates normal and ensemble problems.
+    base_ssprob = SteadyStateProblem(model, u0_alts_vec[1], p_alts_vec[1])
+    base_sol = solve(base_ssprob, DynamicSS(Tsit5()))
+    base_eprob = EnsembleProblem(base_ssprob)
+    base_esol = solve(base_eprob, DynamicSS(Tsit5()); trajectories = 2)
+
+    # Simulates problems for all input types, checking that identical solutions are found.
+    for u0 in u0_alts_vec, p in p_alts_vec
+        ssprob = remake(base_ssprob; u0, p)
+        @test base_sol == solve(ssprob, DynamicSS(Tsit5()))
+        eprob = remake(base_eprob; u0, p)
+        @test base_esol == solve(eprob, DynamicSS(Tsit5()); trajectories = 2)
+    end
+end
 
 ### Checks Errors On Faulty Inputs ###
 

test/network_analysis/conservation_laws.jl

@@ -73,3 +73,25 @@ let
     @test length(cons_laws_constants) == 2
     @test count(isequal.(conserved_quantity, Num(0))) == 2
 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
+    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]] ≈ 8.0
+    @test sol[X[2]] ≈ 4.0
+end

test/reactionsystem_structure/reactionsystem.jl

@@ -279,6 +279,68 @@ 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])
+        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 for https://github.com/SciML/ModelingToolkit.jl/issues/436.