update

Author: TorkelE

test/miscellaneous_tests/symbolic_stoichiometry.jl

@@ -57,7 +57,7 @@ begin
     u0_1 = (A => 3.0, B => 2.0, C => 3.0, D => 5.0)
     ps_1 = (k => 5.0, α => 2)
     u0_2 = [Int64(u[2]) for u in u0_1]
-    ps_2 = Tuple(p[2] for p in ps_1)
+    ps_2 = [p[2] for p in ps_1]
     τ = 1.5
 
     # Creates `ReactionSystem` model.
@@ -248,3 +248,46 @@ let
     jprob_int_ref = JumpProblem(rs_ref_int, dprob_int_ref, Direct(); rng)
     @test solve(jprob_int, SSAStepper(); seed) == solve(jprob_int_ref, SSAStepper(); seed)
 end
+
+# Check that jump simulations (implemented with and without symbolic stoichiometries) yield simulations
+# with identical mean number of species at the end of the simulations.
+# Also checks that ODE simulations are identical.
+let
+    # Creates the models.
+    sir = @reaction_network begin
+        @parameters n::Int64 k::Int64
+        i, S + n*I --> k*I
+        r, n*I --> n*R
+    end      
+    sir_ref = @reaction_network begin
+        i, S + I --> 2I
+        r, I --> R
+    end  
+
+    ps = [:i => 1e-4, :r => 1e-2, :n => 1.0, :k => 2.0]
+    ps_ref = [:i => 1e-4, :r => 1e-2]
+    tspan = (0.0, 250.0) # tspan[2] is selected so that it is in the middle of the outbreak peak.
+    u0 = [:S => 999.0, :I => 1.0, :R => 0.0]
+    @test issetequal(unknowns(sir), unknowns(sir_ref))
+
+    # Checks that ODE simulations are identical.
+    oprob = ODEProblem(sir, u0, tspan, ps)
+    oprob_ref = ODEProblem(sir_ref, u0, tspan, ps_ref)
+    @test solve(oprob, Tsit5()) ≈ solve(oprob_ref, Tsit5())
+
+    # Jumps. First ensemble problems for each systems is created.
+    dprob = DiscreteProblem(sir, u0, tspan, ps) 
+    dprob_ref = DiscreteProblem(sir_ref, u0, tspan, ps_ref) 
+    jprob = JumpProblem(sir, dprob, Direct(), save_positions = (false, false)) 
+    jprob_ref = JumpProblem(sir_ref, dprob_ref, Direct(), save_positions = (false, false)) 
+    eprob = EnsembleProblem(jprob)
+    eprob_ref = EnsembleProblem(jprob_ref)
+
+    # Simulates both ensemble problems. Checks that the distribution of values at the simulation ends is similar.
+    sols = solve(eprob, SSAStepper(); trajectories = 100000)
+    sols_ref = solve(eprob_ref, SSAStepper(); trajectories = 100000)
+    end_vals = [[sol[s][end] for sol in sols.u] for s in species(sir)]
+    end_vals_ref = [[sol[s][end] for sol in sols_ref.u] for s in species(sir_ref)]
+    @test mean.(end_vals_ref) ≈ mean.(end_vals) atol=1e-2 rtol = 1e-2
+    @test var.(end_vals_ref) ≈ var.(end_vals) atol=1e-2 rtol = 1e-2
+end

test/reactionsystem_structure/designating_parameter_types.jl

@@ -68,6 +68,7 @@ let
     end
 end
 
+# Test that the various structures stores the parameters using the correct type.
 # Test that the various structures stores the parameters using the correct type.
 let
     # Creates problems, integrators, and solutions.
@@ -87,8 +88,9 @@ let
     jsol = solve(jprob, SSAStepper(); seed)
     nsol = solve(nprob, NewtonRaphson())
 
-    # Checks the types of all stored parameter values.
+    # Checks all stored parameters.
     for mtk_struct in [oprob, sprob, dprob, jprob, nprob, oinit, sinit, jinit, osol, ssol, jsol, nsol]
+        # Checks that all parameters have the correct type.
         @test unwrap(mtk_struct.ps[p1]) isa Float64
         @test unwrap(mtk_struct.ps[d1]) isa Float64
         @test unwrap(mtk_struct.ps[p2]) isa Float64
@@ -99,6 +101,35 @@ let
         @test unwrap(mtk_struct.ps[d4]) isa Rational{Int64}
         @test unwrap(mtk_struct.ps[p5]) isa Rational{Int64}
         @test unwrap(mtk_struct.ps[d5]) isa Float32
+
+        # Checks that all parameters have the correct value.
+        @test unwrap(mtk_struct.ps[p1]) == 1.0
+        @test unwrap(mtk_struct.ps[d1]) == 1.0
+        @test unwrap(mtk_struct.ps[p2]) == 1.2
+        @test unwrap(mtk_struct.ps[d2]) == 1.2
+        @test unwrap(mtk_struct.ps[p3]) == 2
+        @test unwrap(mtk_struct.ps[d3]) == 2
+        @test unwrap(mtk_struct.ps[p4]) == Float32(0.5)
+        @test unwrap(mtk_struct.ps[d4]) == 1//2
+        @test unwrap(mtk_struct.ps[p5]) == 3//2
+        @test unwrap(mtk_struct.ps[d5]) == Float32(1.5)
+    end
+    
+    # Checks all stored variables
+    for mtk_struct in [oprob, sprob, dprob, jprob, nprob, oinit, sinit, jinit]
+        # Checks that all variables have the correct type.
+        @test unwrap(mtk_struct[X1]) isa Float64
+        @test unwrap(mtk_struct[X2]) isa Float64
+        @test unwrap(mtk_struct[X3]) isa Float64
+        @test unwrap(mtk_struct[X4]) isa Float64
+        @test unwrap(mtk_struct[X5]) isa Float64
+
+        # Checks that all variables have the correct value.
+        @test unwrap(mtk_struct[X1]) == 0.1
+        @test unwrap(mtk_struct[X2]) == 0.2
+        @test unwrap(mtk_struct[X3]) == 0.3
+        @test unwrap(mtk_struct[X4]) == 0.4
+        @test unwrap(mtk_struct[X5]) == 0.5
     end
 
     # Indexing currently broken for NonlinearSystem integrators (MTK intend to support this though).

test/test_functions.jl

@@ -1,7 +1,7 @@
 # Various functions that are useful for running the tests, and used across several test sets.
 
 # Fetches the (required) Random package.
-using Random
+using Random, Test
 
 ### Initial Condition/Parameter Generators ###
 
@@ -36,19 +36,31 @@ end
 ### System function evaluation ###
 
 # Evaluates the the drift function of the ODE corresponding to a reaction network.
+# Also checks that in place and out of place evaluations are identical.
 function f_eval(rs::ReactionSystem, u, p, t; combinatoric_ratelaws = true)
     prob = ODEProblem(rs, u, (0.0, 0.0), p; combinatoric_ratelaws)
-    return prob.f(prob.u0, prob.p, t)
+    du = zeros(length(u))
+    prob.f(du, prob.u0, prob.p, t)
+    @test du == prob.f(prob.u0, prob.p, t)
+    return du
 end
 
 # Evaluates the the Jacobian of the drift function of the ODE corresponding to a reaction network.
+# Also checks that in place and out of place evaluations are identical.
 function jac_eval(rs::ReactionSystem, u, p, t; combinatoric_ratelaws = true)
     prob = ODEProblem(rs, u, (0.0, 0.0), p; jac = true, combinatoric_ratelaws)
-    return prob.f.jac(prob.u0, prob.p, t)
+    J = zeros(length(u), length(u))
+    prob.f.jac(J, prob.u0, prob.p, t)
+    @test J == prob.f.jac(prob.u0, prob.p, t)
+    return J
 end
 
 # Evaluates the the diffusion function of the SDE corresponding to a reaction network.
+# Also checks that in place and out of place evaluations are identical.
 function g_eval(rs::ReactionSystem, u, p, t; combinatoric_ratelaws = true)
     prob = SDEProblem(rs, u, (0.0, 0.0), p; combinatoric_ratelaws)
-    return prob.g(prob.u0, prob.p, t)
-end
+    dW = zeros(length(u), numreactions(rs))
+    prob.g(dW, prob.u0, prob.p, t)
+    @test dW == prob.g(prob.u0, prob.p, t)
+    return dW
+end