change to not depend on equality after seed in stoch tests

test/reactionsystem_core/events.jl

@@ -243,10 +243,6 @@ let
     sol_dsl = solve(ODEProblem(rn_dsl, u0, tspan, ps), Tsit5())
     sol_prog = solve(ODEProblem(rn_prog, u0, tspan, ps), Tsit5())
     @test sol_dsl == sol_prog
-
-    sol_dsl = solve(SDEProblem(rn_dsl, u0, tspan, ps), ImplicitEM(); seed = 1234)
-    sol_prog = solve(SDEProblem(rn_prog, u0, tspan, ps), ImplicitEM(); seed = 1234)
-    @test sol_dsl == sol_prog
 end
 
 # Checks that misformatted events yields errors in the DSL.
@@ -313,6 +309,65 @@ end
 
 ### Additional Correctness Tests ###
 
+# Tests that events are properly triggered for SDEs.
+# Tests for continuous events, and all three types of discrete events.
+let
+    # Creates model with all types of events. The `e` parameters track whether events are triggered.
+    rn = @reaction_network begin
+        @parameters e1=0 e2=0 e3=0 e4=0
+        @continuous_events begin
+            [X ~ 1000.0] => [e1 ~ 1]
+        end
+        @discrete_events begin
+            [1.0] => [e2 ~ 1]
+            1.0 => [e3 ~ 1]
+            (Y > 1000.0) & (e4==0) => [e4 ~ 1]
+        end
+        (p,d), 0 <--> X
+        (p,d), 0 <--> Y
+    end
+
+    # Simulates the model for conditions where it *definitely* will cross `X = 1000.0`
+    u0 = [:X => 999.9, :Y => 999.9]
+    ps = [:p => 10.0, :d => 0.001]
+    sprob = SDEProblem(rn, u0, (0.0, 2.0), ps)
+    sol = solve(sprob, ImplicitEM(); seed)
+
+    # Checks that all `e` parameters have been updated properly.
+    @test sol.ps[:e1] == 1
+    @test sol.ps[:e2] == 1
+    @test sol.ps[:e3] == 1
+    @test sol.ps[:e4] == 1
+end
+
+# Tests that events are properly triggered for Jump simulations.
+# Tests for all three types of discrete events.
+let
+    # Creates model with all types of events. The `e` parameters track whether events are triggered.
+    rn = @reaction_network begin
+        @parameters e1=0 e2=0 e3=0
+        @discrete_events begin
+            [1.0] => [e1 ~ 1]
+            # 1.0 => [e2 ~ 1]
+            (X > 1000.0) & (e3==0) => [e3 ~ 1]
+        end
+        (p,d), 0 <--> X
+    end
+
+    # Simulates the model for conditions where it *definitely* will cross `X = 1000.0`
+    u0 = [:X => 999]
+    ps = [:p => 10.0, :d => 0.001]
+    dprob = DiscreteProblem(rn, u0, (0.0, 2.0), ps)
+    jprob = JumpProblem(rn, dprob, Direct(); rng)
+    sol = solve(jprob, SSAStepper(); seed)
+
+    # Checks that all `e` parameters have been updated properly.
+    # Note that periodic discrete events are currently broken for jump processes.
+    @test sol.ps[:e1] == 1
+    @test_broken sol.ps[:e2] == 1
+    @test sol.ps[:e3] == 1
+end
+
 # Compares simulations using MTK type events with those generated through callbacks.
 # Jump simulations must be handles differently (since these only accepts discrete callbacks).
 # Checks for all types of discrete callbacks, and for continuous callbacks.

test/reactionsystem_core/higher_order_reactions.jl

@@ -103,7 +103,7 @@ let
     dprob_alt2 = DiscreteProblem(u0_alt2, (0.0, 100.0), ps_alt2)
     jprob_alt2 = JumpProblem(dprob_alt2, Direct(), higher_order_network_alt2...; rng = StableRNG(1234))
 
-    # Simualtes the models.
+    # Simulates the models.
     sol_base = solve(jprob_base, SSAStepper(); seed, saveat = 1.0)
     sol_alt1 = solve(jprob_alt1, SSAStepper(); seed, saveat = 1.0)
     sol_alt2 = solve(jprob_alt2, SSAStepper(); seed, saveat = 1.0)

test/reactionsystem_core/symbolic_stoichiometry.jl

@@ -192,64 +192,68 @@ end
 # Tests symbolic stoichiometries in simulations.
 # Tests for decimal numbered symbolic stoichiometries.
 let
+    # Declares models. The references models have the `n` parameters so they can use the same 
+    # parameter vectors as the non-reference ones.
     rs_int = @reaction_network begin
-        @parameters n1::Int64 n2::Int64
-        p, 0 -->X
-        k, n1*X --> n2*Y
-        d, Y --> 0
+        @parameters n::Int64
+        (k1, k2), n*X1 <--> X2
     end
     rs_dec = @reaction_network begin
-        @parameters n1::Float64 n2::Float64
-        p, 0 -->X
-        k, n1*X --> n2*Y
-        d, Y --> 0
+        @parameters n::Float64
+        (k1, k2), n*X1 <--> X2
     end
     rs_ref_int = @reaction_network begin
-        @parameters n1::Int64 n2::Int64
-        p, 0 -->X
-        k, 3*X --> 4*Y
-        d, Y --> 0
+        @parameters n::Int64
+        (k1, k2), 3*X1 <--> X2
     end
     rs_ref_dec = @reaction_network begin
-        @parameters n1::Float64 n2::Float64
-        p, 0 -->X
-        k, 0.5*X --> 2.5*Y
-        d, Y --> 0
+        @parameters n::Float64
+        (k1, k2), 2.5*X1 <--> X2
     end
 
-    u0 = [:X => 8, :Y => 8]
-    ps_int = (:p => 2.0, :k => 0.01, :n1 => 3, :n2 => 4, :d => 0.2)
-    ps_dec = [:p => 2.0, :k => 0.01, :n1 => 0.5, :n2 => 2.5, :d => 0.2]
-    tspan = (0.0, 10.0)
+    # Set simulation settings. Initial conditions are design to start, more or less, at 
+    # steady state concentrations.
+    # Values are selected so that stochastic tests should always pass within the bounds (independent 
+    # of seed).
+    u0_int = [:X1 => 150, :X2 => 600]
+    u0_dec = [:X1 => 100, :X2 => 600]
+    tspan_det = (0.0, 1.0)
+    tspan_stoch = (0.0, 10000.0)
+    ps_int = (:k1 => 0.00001, :k2 => 0.01, :n => 3)
+    ps_dec = (:k1 => 0.00001, :k2 => 0.01, :n => 2.5)
 
     # Test ODE simulations with integer coefficients.
-    oprob_int = ODEProblem(rs_int, u0, tspan, ps_int)
-    oprob_int_ref = ODEProblem(rs_ref_int, u0, tspan, ps_int)
-    @test solve(oprob_int, Tsit5()) ≈ solve(oprob_int_ref, Tsit5())
+    oprob_int = ODEProblem(rs_int, u0_int, tspan_det, ps_int)
+    oprob_int_ref = ODEProblem(rs_ref_int, u0_int, tspan_det, ps_int)
+    @test solve(oprob_int, Tsit5())[:X1][end] ≈ solve(oprob_int_ref, Tsit5())[:X1][end]
 
     # Test ODE simulations with decimal coefficients.
-    oprob_dec = ODEProblem(rs_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false)
-    oprob_dec_ref = ODEProblem(rs_ref_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false)
-    @test solve(oprob_dec, Tsit5()) ≈ solve(oprob_dec_ref, Tsit5())
+    oprob_dec = ODEProblem(rs_dec, u0_dec, tspan_det, ps_dec; combinatoric_ratelaws = false)
+    oprob_dec_ref = ODEProblem(rs_ref_dec, u0_dec, tspan_det, ps_dec; combinatoric_ratelaws = false)
+    @test solve(oprob_dec, Tsit5())[:X1][end] ≈ solve(oprob_dec_ref, Tsit5())[:X1][end]
 
     # Test SDE simulations with integer coefficients.
-    sprob_int = SDEProblem(rs_int, u0, tspan, ps_int)
-    sprob_int_ref = SDEProblem(rs_ref_int, u0, tspan, ps_int)
-    @test solve(sprob_int, ImplicitEM(); seed) ≈ solve(sprob_int_ref, ImplicitEM(); seed)
+    sprob_int = SDEProblem(rs_int, u0_int, tspan_stoch, ps_int)
+    sprob_int_ref = SDEProblem(rs_ref_int, u0_dec, tspan_stoch, ps_int)
+    ssol_int = solve(sprob_int, ImplicitEM(); seed)
+    ssol_int_ref = solve(sprob_int_ref, ImplicitEM(); seed)
+    @test mean(ssol_int[:X1]) ≈ mean(ssol_int_ref[:X1]) atol = 2*1e0
 
     # Test SDE simulations with decimal coefficients.
-    sprob_dec = SDEProblem(rs_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false)
-    sprob_dec_ref = SDEProblem(rs_ref_dec, u0, tspan, ps_dec; combinatoric_ratelaws = false)
-    @test solve(sprob_dec, ImplicitEM(); seed) ≈ solve(sprob_dec_ref, ImplicitEM(); seed)
-
-    # Tests jump simulations with integer coefficients.
-    dprob_int = DiscreteProblem(rs_int, u0, (0.0, 100000.0), ps_int)
-    dprob_int_ref = DiscreteProblem(rs_ref_int, u0, (0.0, 100000.0), ps_int)
-    jprob_int = JumpProblem(rs_int, dprob_int, Direct(); rng)
-    jprob_int_ref = JumpProblem(rs_ref_int, dprob_int_ref, Direct(); rng)
-    sol_int = solve(jprob_int, SSAStepper(); seed)
-    sol_int_ref = solve(jprob_int_ref, SSAStepper(); seed)
-    @test mean(sol_int[:Y]) ≈ mean(sol_int_ref[:Y]) atol = 1e-2 rtol = 1e-2
+    sprob_dec = SDEProblem(rs_dec, u0_dec, tspan_stoch, ps_dec; combinatoric_ratelaws = false)
+    sprob_dec_ref = SDEProblem(rs_ref_dec, u0_dec, tspan_stoch, ps_dec; combinatoric_ratelaws = false)
+    ssol_dec = solve(sprob_dec, ImplicitEM(); seed)
+    ssol_dec_ref = solve(sprob_dec_ref, ImplicitEM(); seed)
+    @test mean(ssol_dec[:X1]) ≈ mean(ssol_dec_ref[:X1]) atol = 2*1e0
+
+    # Test Jump simulations with integer coefficients.
+    dprob_int = DiscreteProblem(rs_int, u0_int, tspan_stoch, ps_int)
+    dprob_int_ref = DiscreteProblem(rs_ref_int, u0_int, tspan_stoch, ps_int)
+    jprob_int = JumpProblem(rs_int, dprob_int, Direct(); rng, save_positions = (false, false))
+    jprob_int_ref = JumpProblem(rs_ref_int, dprob_int_ref, Direct(); rng, save_positions = (false, false))
+    jsol_int = solve(jprob_int, SSAStepper(); seed, saveat = 1.0)
+    jsol_int_ref = solve(jprob_int_ref, SSAStepper(); seed, saveat = 1.0)
+    @test mean(jsol_int[:X1]) ≈ mean(jsol_int_ref[:X1]) atol = 1e-2 rtol = 1e-2
 end
 
 # Check that jump simulations (implemented with and without symbolic stoichiometries) yield simulations

test/simulation_and_solving/simulate_SDEs.jl

@@ -128,17 +128,6 @@ let
             duW = zeros(size(rn_manual.nrp))
             rn_manual.g(duW, u0_2, ps_2, 0.0)
             @test duW ≈ g_eval(rn_catalyst, u0_1, ps_1, 0.0)
-
-            # Compares simulation with identical seed. 
-            # Cannot test the third network (requires very fiddly parameters for CLE to be defined).
-            if nameof(rn_catalyst) != :rnc9
-                sprob_1 = SDEProblem(rn_catalyst, u0_1, (0.0, 1.0), ps_1)
-                sprob_2 = SDEProblem(rn_manual.f, rn_manual.g, u0_2, (0.0, 1.0), ps_2; noise_rate_prototype = rn_manual.nrp)
-                seed =  seed = rand(rng, 1:100)
-                sol1 = solve(sprob_1, ImplicitEM(); seed)
-                sol2 = solve(sprob_2, ImplicitEM(); seed)
-                @test sol1[u0_sym] ≈ sol2.u
-            end
         end
     end
 end

test/simulation_and_solving/simulate_jumps.jl

@@ -6,6 +6,7 @@ using Catalyst, JumpProcesses, Statistics, Test
 # Sets stable rng number.
 using StableRNGs
 rng = StableRNG(12345)
+seed = rand(rng, 1:100)
 
 # Fetch test functions and networks.
 include("../test_functions.jl")
@@ -16,12 +17,17 @@ include("../test_networks.jl")
 
 # Compares jump simulations generated through hCatalyst, and manually created systems.
 let
-    # Manually declares jumps to compare Catalyst-generated jump simulations to.
+    # Declares vectors to store information about each network in. Initial conditions are approximately
+    # at the steady state.
     catalyst_networks = []
     manual_networks = []
     u0_syms = []
     ps_syms = []
+    u0s = []
+    ps = []
+    sps = []
 
+    # Manual declaration of the first network (and its simulation settings).
     rate_1_1(u, p, t) = p[1]
     rate_1_2(u, p, t) = p[2] * u[1]
     rate_1_3(u, p, t) = p[3] * u[2]
@@ -47,12 +53,16 @@ let
     jump_1_7 = ConstantRateJump(rate_1_7, affect_1_7!)
     jump_1_8 = ConstantRateJump(rate_1_8, affect_1_8!)
     jumps_1 = (jump_1_1, jump_1_2, jump_1_3, jump_1_4, jump_1_5, jump_1_6, jump_1_7,
-               jump_1_8)
+                jump_1_8)
     push!(catalyst_networks, reaction_networks_standard[5])
     push!(manual_networks, jumps_1)
     push!(u0_syms, [:X1, :X2, :X3, :X4])
     push!(ps_syms, [:p, :k1, :k2, :k3, :k4, :k5, :k6, :d])
+    push!(u0s, [:X1 => 40, :X2 => 30, :X3 => 20, :X4 => 10])
+    push!(ps, [:p => 10.0, :k1 => 1.0, :k2 => 1.0, :k3 => 1.0, :k4 => 1.0, :k5 => 1.0, :k6 => 1.0, :d => 1.0])
+    push!(sps, :X4)
 
+    # Manual declaration of the second network (and its simulation settings).
     rate_2_1(u, p, t) = p[1] / 10 + p[1] * (u[1]^p[3]) / (u[1]^p[3] + p[2]^p[3])
     rate_2_2(u, p, t) = p[4] * u[1] * u[2]
     rate_2_3(u, p, t) = p[5] * u[3]
@@ -83,7 +93,11 @@ let
     push!(manual_networks, jumps_2)
     push!(u0_syms, [:X1, :X2, :X3])
     push!(ps_syms, [:v, :K, :n, :k1, :k2, :k3, :d])
+    push!(u0s, [:X1 => 10, :X2 => 10, :X3 => 10])
+    push!(ps, [:v => 20.0, :K => 2.0, :n => 2, :k1 => 0.1, :k2 => 0.1, :k3 => 0.1, :d => 1.0])
+    push!(sps, :X3)
 
+    # Manual declaration of the third network (and its simulation settings).
     rate_3_1(u, p, t) = p[1] * binomial(u[1], 1)
     rate_3_2(u, p, t) = p[2] * binomial(u[2], 2)
     rate_3_3(u, p, t) = p[3] * binomial(u[2], 2)
@@ -107,33 +121,29 @@ let
     push!(manual_networks, jumps_3)
     push!(u0_syms, [:X1, :X2, :X3, :X4])
     push!(ps_syms, [:k1, :k2, :k3, :k4, :k5, :k6])
+    push!(u0s, [:X1 => 15, :X2 => 5, :X3 => 5, :X4 => 5])
+    push!(ps, [:k1 => 0.1, :k2 => 0.1, :k3 => 0.1, :k4 => 0.1, :k5 => 0.1, :k6 => 0.1])
+    push!(sps, :X4)
 
     # Loops through all cases, checks that identical simulations are generated with/without Catalyst.
-    for (rn_catalyst, rn_manual, u0_sym, ps_sym) in zip(catalyst_networks, manual_networks, u0_syms, ps_syms)
-        for factor in [5, 50]
-            seed = rand(rng, 1:100)
-            u0_1 = rnd_u0_Int64(rn_catalyst, rng; n = factor)
-            ps_1 = rnd_ps(rn_catalyst, rng; factor = factor/100.0)
-            dprob_1 = DiscreteProblem(rn_catalyst, u0_1, (0.0, 100.0), ps_1)
-            jprob_1 = JumpProblem(rn_catalyst, dprob_1, Direct(); rng)
-            sol1 = solve(jprob_1, SSAStepper(); seed, saveat = 1.0)
-
-            u0_2 = map_to_vec(u0_1, u0_sym)
-            ps_2 = map_to_vec(ps_1, ps_sym)
-            dprob_2 = DiscreteProblem(u0_2, (0.0, 100.0), ps_2)
-            jprob_2 = JumpProblem(dprob_2, Direct(), rn_manual...; rng)
-            sol2 = solve(jprob_2, SSAStepper(); seed, saveat = 1.0)
+    for (rn_catalyst, rn_manual, u0_sym, ps_sym, u0_1, ps_1, sp) in 
+            zip(catalyst_networks, manual_networks, u0_syms, ps_syms, u0s, ps, sps)
 
-            if nameof(rn_catalyst) == :rnh7
-                # Have spent a few hours figuring this one out. For certain seeds it actually works,
-                # but others not. This feels weird, and I didn't get any longer. I tried using non-random
-                # parameters/initial conditions, and removing the non-hill function reactions. Problem
-                # still persists.
-                @test_broken sol1[u0_sym] == sol2.u
-            else
-                @test sol1[u0_sym] == sol2.u
-            end
-        end
+        # Simulates the Catalyst-created model.
+        dprob_1 = DiscreteProblem(rn_catalyst, u0_1, (0.0, 10000.0), ps_1)
+        jprob_1 = JumpProblem(rn_catalyst, dprob_1, Direct(); rng)
+        sol1 = solve(jprob_1, SSAStepper(); seed, saveat = 1.0)
+
+        # Simulates the manually written model
+        u0_2 = map_to_vec(u0_1, u0_sym)
+        ps_2 = map_to_vec(ps_1, ps_sym)
+        dprob_2 = DiscreteProblem(u0_2, (0.0, 10000.0), ps_2)
+        jprob_2 = JumpProblem(dprob_2, Direct(), rn_manual...; rng)
+        sol2 = solve(jprob_2, SSAStepper(); seed, saveat = 1.0)
+
+        # Checks that the means are similar (the test have been check that it holds across a large 
+        # number of simulates, even without seed).
+        @test mean(sol1[sp]) ≈ mean(sol2[findfirst(u0_sym .== sp),:]) rtol = 1e-1
     end
 end