diff --git a/test/reactionsystem_core/events.jl b/test/reactionsystem_core/events.jl index 5a5c4b509f..ab15dd68ae 100644 --- a/test/reactionsystem_core/events.jl +++ b/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. diff --git a/test/reactionsystem_core/higher_order_reactions.jl b/test/reactionsystem_core/higher_order_reactions.jl index f575cd5f55..7922905792 100644 --- a/test/reactionsystem_core/higher_order_reactions.jl +++ b/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) diff --git a/test/reactionsystem_core/symbolic_stoichiometry.jl b/test/reactionsystem_core/symbolic_stoichiometry.jl index 745a9c50af..797a78f239 100644 --- a/test/reactionsystem_core/symbolic_stoichiometry.jl +++ b/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 diff --git a/test/simulation_and_solving/simulate_SDEs.jl b/test/simulation_and_solving/simulate_SDEs.jl index 4b1c935435..5684db4132 100644 --- a/test/simulation_and_solving/simulate_SDEs.jl +++ b/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 diff --git a/test/simulation_and_solving/simulate_jumps.jl b/test/simulation_and_solving/simulate_jumps.jl index 6425758b89..9c8c02db8a 100644 --- a/test/simulation_and_solving/simulate_jumps.jl +++ b/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