Merge pull request #1257 from SciML/use_jumpinputs_in_tests

Remove` DiscreteProblem`s in tests

Author: TorkelE

test/reactionsystem_core/events.jl

@@ -113,17 +113,15 @@ let
         @test prob.ps[α] isa Int64
     end
 
-    # Handles `DiscreteProblem`s and `JumpProblem`s (these cannot contain continuous events or variables).
+    # Handles `JumpInput`s and `JumpProblem`s (these cannot contain continuous events or variables).
     discrete_events = [2.0 => [A ~ A + α]]
     @named rs_de_2 = ReactionSystem(rxs, t; discrete_events)
     rs_de_2 = complete(rs_de_2)
-    dprob = DiscreteProblem(rs_de_2, u0, (0.0, 10.0), ps)
-    jprob = JumpProblem(rs_de_2, dprob, Direct())
-    for prob in [dprob, jprob]
-        @test dprob[A] == 2
-        @test dprob.ps[α] == 1
-        @test dprob.ps[α] isa Int64
-    end
+    jin = JumpInputs(rs_de_2, u0, (0.0, 10.0), ps)
+    jprob = JumpProblem(jin)
+    @test jprob[A] == 2
+    @test jprob.ps[α] == 1
+    @test jprob.ps[α] isa Int64
 end
 
 
@@ -358,13 +356,13 @@ let
     # 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)
+    jin = JumpInputs(rn, u0, (0.0, 2.0), ps)
+    jprob = JumpProblem(jin; rng)
     sol = solve(jprob, SSAStepper(); seed)
 
     # Checks that all `e` parameters have been updated properly.
     @test sol.ps[:e1] == 1
-    @test sol.ps[:e2] == 1 
+    @test sol.ps[:e2] == 1
     @test sol.ps[:e3] == 1
 end
 
@@ -434,10 +432,10 @@ let
     # Checks for Jump simulations. (note, non-seed dependant test should be created instead)
     # Note that periodic discrete events are currently broken for jump processes (and unlikely to be fixed soon due to have events are implemented).
     callback = CallbackSet(cb_disc_1, cb_disc_2, cb_disc_3)
-    dprob = DiscreteProblem(rn, u0, tspan, ps)
-    dprob_events = DiscreteProblem(rn_dics_events, u0, tspan, ps)
-    jprob = JumpProblem(rn, dprob, Direct(); rng)
-    jprob_events = JumpProblem(rn_dics_events, dprob_events, Direct(); rng)
+    jin = JumpInputs(rn, u0, tspan, ps)
+    jin_events = JumpInputs(rn_dics_events, u0, tspan, ps)
+    jprob = JumpProblem(jin)
+    jprob_events = JumpProblem(jin_events; rng)
     sol = solve(jprob, SSAStepper(); seed, callback)
     sol_events = solve(jprob_events, SSAStepper(); seed)
     @test_broken sol == sol_events  # seems to be not identical in the sample paths

test/reactionsystem_core/higher_order_reactions.jl

@@ -89,12 +89,12 @@ let
     # Prepares JumpProblem via Catalyst.
     u0_base = rnd_u0_Int64(base_higher_order_network, rng)
     ps_base = rnd_ps(base_higher_order_network, rng)
-    dprob_base = DiscreteProblem(base_higher_order_network, u0_base, (0.0, 100.0), ps_base)
-    jprob_base = JumpProblem(base_higher_order_network, dprob_base, Direct(); rng = StableRNG(1234))
+    jin_base = JumpInputs(base_higher_order_network, u0_base, (0.0, 100.0), ps_base)
+    jprob_base = JumpProblem(jin_base; rng = StableRNG(1234))
 
     # Prepares JumpProblem partially declared manually.
-    dprob_alt1 = DiscreteProblem(higher_order_network_alt1, u0_base, (0.0, 100.0), ps_base)
-    jprob_alt1 = JumpProblem(higher_order_network_alt1, dprob_alt1, Direct(); rng = StableRNG(1234))
+    jin_alt1 = JumpInputs(higher_order_network_alt1, u0_base, (0.0, 100.0), ps_base)
+    jprob_alt1 = JumpProblem(jin_alt1; rng = StableRNG(1234))
 
     # Prepares JumpProblem via manually declared system.
     u0_alt2 = map_to_vec(u0_base, [:X1, :X2, :X3, :X4, :X5, :X6, :X7, :X8, :X9, :X10])

test/reactionsystem_core/parameter_type_designation.jl

@@ -74,13 +74,13 @@ let
     oprob = ODEProblem(rs, u0, (0.0, 1.0), p_alts[1])
     sprob = SDEProblem(rs, u0, (0.0, 1.0), p_alts[1])
     dprob = DiscreteProblem(rs, u0, (0.0, 1.0), p_alts[1])
-    jprob = JumpProblem(rs, dprob, Direct(); rng)
+    jprob = JumpProblem(JumpInputs(rs, u0, (0.0, 1.0), p_alts[1]); rng)
     nprob = NonlinearProblem(rs, u0, p_alts[1])
 
-    oinit = init(oprob, Tsit5()) 
-    sinit = init(sprob, ImplicitEM()) 
-    jinit = init(jprob, SSAStepper()) 
-    ninit = init(nprob, NewtonRaphson()) 
+    oinit = init(oprob, Tsit5())
+    sinit = init(sprob, ImplicitEM())
+    jinit = init(jprob, SSAStepper())
+    ninit = init(nprob, NewtonRaphson())
 
     osol = solve(oprob, Tsit5())
     ssol = solve(sprob, ImplicitEM(); seed)
@@ -113,7 +113,7 @@ let
         @test unwrap(mtk_struct.ps[p5]) == 3//2
         @test unwrap(mtk_struct.ps[d5]) == Float32(1.5)
     end
-    
+
     # Checks all stored variables (these should always be `Float64`).
     for mtk_struct in [oprob, sprob, dprob, jprob, nprob, oinit, sinit, jinit, ninit]
         # Checks that all variables have the correct type.

test/reactionsystem_core/reactionsystem.jl

@@ -126,10 +126,10 @@ let
     sdesys = complete(convert(SDESystem, rs))
     js = complete(convert(JumpSystem, rs))
 
-    @test ModelingToolkit.get_defaults(rs) == 
+    @test ModelingToolkit.get_defaults(rs) ==
           ModelingToolkit.get_defaults(js) == defs
 
-    # these systems add initial conditions to the defaults 
+    # these systems add initial conditions to the defaults
     @test ModelingToolkit.get_defaults(odesys) ==
           ModelingToolkit.get_defaults(sdesys)
     @test issubset(defs, ModelingToolkit.get_defaults(odesys))
@@ -551,9 +551,9 @@ let
                                    (@reaction k1, $A --> B2),
                                    (@reaction 10 * k1, ∅ --> B3)], t)
     rn = complete(rn)
-    dprob = DiscreteProblem(rn, [A => 10, C => 10, B1 => 0, B2 => 0, B3 => 0], (0.0, 10.0),
+    jin = JumpInputs(rn, [A => 10, C => 10, B1 => 0, B2 => 0, B3 => 0], (0.0, 10.0),
                             [k1 => 1.0])
-    jprob = JumpProblem(rn, dprob, Direct(); rng, save_positions = (false, false))
+    jprob = JumpProblem(jin; rng, save_positions = (false, false))
     umean = zeros(4)
     Nsims = 40000
     for i in 1:Nsims
@@ -1017,7 +1017,7 @@ let
     @test sys isa JumpSystem
     @test MT.has_equations(sys)
     @test length(massactionjumps(sys)) == 1
-    @test isempty(constantratejumps(sys)) 
+    @test isempty(constantratejumps(sys))
     @test length(variableratejumps(sys)) == 3
     @test length(odeeqs(sys)) == 4
     @test length(continuous_events(sys)) == 1
@@ -1042,7 +1042,7 @@ let
     @test sys isa JumpSystem
     @test MT.has_equations(sys)
     @test length(massactionjumps(sys)) == 1
-    @test isempty(constantratejumps(sys)) 
+    @test isempty(constantratejumps(sys))
     @test length(variableratejumps(sys)) == 2
     @test length(odeeqs(sys)) == 4
     odes = union(eqs, [D(A) ~ 0, D(B) ~ -λ*A*B, D(C) ~ 0])
@@ -1069,8 +1069,8 @@ let
     sys = jinput.sys
     @test sys isa JumpSystem
     @test MT.has_equations(sys)
-    @test isempty(massactionjumps(sys)) 
-    @test isempty(constantratejumps(sys)) 
+    @test isempty(massactionjumps(sys))
+    @test isempty(constantratejumps(sys))
     @test length(variableratejumps(sys)) == 3
     @test length(odeeqs(sys)) == 4
     odes = union(eqs, [D(A) ~ 0, D(B) ~ -λ*A*B, D(C) ~ 0])

test/reactionsystem_core/symbolic_stoichiometry.jl

@@ -19,7 +19,7 @@ include("../test_functions.jl")
 ### Base Tests ###
 
 # Checks that systems with symbolic stoichiometries, created using different approaches, are identical.
-let 
+let
     @parameters p k d::Float64 n1::Int64 n2 n3
     @species X(t) Y(t)
     rxs1 = [
@@ -71,7 +71,7 @@ begin
 end
 
 # Compares the Catalyst-generated ODE function to a manually computed ODE function.
-let 
+let
     # With combinatoric ratelaws.
     function oderhs(u, p, t)
         k,α = p
@@ -89,7 +89,7 @@ let
     end
     @test f_eval(rs, u0_1, ps_1, τ) ≈ oderhs(u0_2, ps_2, τ)
 
-    # Without combinatoric ratelaws. 
+    # Without combinatoric ratelaws.
     function oderhs_no_crl(u, p, t)
         k,α = p
         A,B,C,D = u
@@ -108,8 +108,8 @@ let
 end
 
 # Compares the Catalyst-generated SDE noise function to a manually computed SDE noise function.
-let 
-    # With combinatoric ratelaws. 
+let
+    # With combinatoric ratelaws.
     function sdenoise(u, p, t)
         k,α = p
         A,B,C,D = u
@@ -126,7 +126,7 @@ let
     end
     @test g_eval(rs, u0_1, ps_1, τ) ≈ sdenoise(u0_2, ps_2, τ)
 
-    # Without combinatoric ratelaws. 
+    # Without combinatoric ratelaws.
     function sdenoise_no_crl(u, p, t)
         k,α = p
         A,B,C,D = u
@@ -192,7 +192,7 @@ 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 
+    # 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 n::Int64
@@ -211,9 +211,9 @@ let
         (k1, k2), 2.5*X1 <--> X2
     end
 
-    # Set simulation settings. Initial conditions are design to start, more or less, at 
+    # 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 
+    # 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]
@@ -247,10 +247,10 @@ let
     @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))
+    jin_int = JumpInputs(rs_int, u0_int, tspan_stoch, ps_int)
+    jin_int_ref = JumpInputs(rs_ref_int, u0_int, tspan_stoch, ps_int)
+    jprob_int = JumpProblem(jin_int; rng, save_positions = (false, false))
+    jprob_int_ref = JumpProblem(jin_int_ref; 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
@@ -265,11 +265,11 @@ let
         @parameters n::Int64 k::Int64
         i, S + n*I --> k*I
         r, n*I --> n*R
-    end      
+    end
     sir_ref = @reaction_network begin
         i, S + I --> 2I
         r, I --> R
-    end  
+    end
 
     ps = [:i => 1e-4, :r => 1e-2, :n => 1.0, :k => 2.0]
     ps_ref = [:i => 1e-4, :r => 1e-2]
@@ -283,10 +283,10 @@ let
     @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(); rng, save_positions = (false, false)) 
-    jprob_ref = JumpProblem(sir_ref, dprob_ref, Direct(); rng, save_positions = (false, false)) 
+    jin = JumpInputs(sir, u0, tspan, ps)
+    jin_ref = JumpInputs(sir_ref, u0, tspan, ps_ref)
+    jprob = JumpProblem(jin; rng, save_positions = (false, false))
+    jprob_ref = JumpProblem(jin_ref; rng, save_positions = (false, false))
     eprob = EnsembleProblem(jprob)
     eprob_ref = EnsembleProblem(jprob_ref)
 

test/simulation_and_solving/simulate_jumps.jl

@@ -130,12 +130,12 @@ let
             zip(catalyst_networks, manual_networks, u0_syms, ps_syms, u0s, ps, sps)
 
         # 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)
+        jin_1 = JumpInputs(rn_catalyst, u0_1, (0.0, 10000.0), ps_1)
+        jprob_1 = JumpProblem(jin_1, Direct(); rng)
         sol1 = solve(jprob_1, SSAStepper(); seed, saveat = 1.0)
 
         # simulate using auto-alg
-        jprob_1b = JumpProblem(rn_catalyst, dprob_1; rng)
+        jprob_1b = JumpProblem(jin_1; rng)
         sol1b = solve(jprob_1; seed, saveat = 1.0)
         @test mean(sol1[sp]) ≈ mean(sol1b[sp]) rtol = 1e-1
 
@@ -157,8 +157,8 @@ let
     for rn in reaction_networks_all
         u0 = rnd_u0_Int64(rn, rng)
         ps = rnd_ps(rn, rng)
-        dprob = DiscreteProblem(rn, u0, (0.0, 1.0), ps)
-        jprob = JumpProblem(rn, dprob, Direct(); rng)
+        jin = JumpInputs(rn, u0, (0.0, 1.0), ps)
+        jprob = JumpProblem(jin; rng)
         @test SciMLBase.successful_retcode(solve(jprob, SSAStepper()))
     end
 end
@@ -169,8 +169,8 @@ let
         (1.2, 5), X1 ↔ X2
     end
     u0 = rnd_u0_Int64(no_param_network, rng)
-    dprob = DiscreteProblem(no_param_network, u0, (0.0, 1000.0))
-    jprob = JumpProblem(no_param_network, dprob, Direct(); rng)
+    jin = JumpInputs(no_param_network, u0, (0.0, 1000.0))
+    jprob = JumpProblem(jin; rng)
     sol = solve(jprob, SSAStepper())
     @test mean(sol[:X1]) > mean(sol[:X2])
 end

test/upstream/mtk_problem_inputs.jl

@@ -133,8 +133,8 @@ end
 # Perform jump simulations (singular and ensemble).
 let
     # Creates normal and ensemble problems.
-    base_dprob = DiscreteProblem(model, u0_alts[1], tspan, p_alts[1])
-    base_jprob = JumpProblem(model, base_dprob, Direct(); rng)
+    base_jin = JumpInputs(model, u0_alts[1], tspan, p_alts[1])
+    base_jprob = JumpProblem(base_jin; 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)
@@ -325,8 +325,8 @@ end
 # Perform jump simulations (singular and ensemble).
 let
     # Creates normal and ensemble problems.
-    base_dprob = DiscreteProblem(model_vec, u0_alts_vec[1], tspan, p_alts_vec[1])
-    base_jprob = JumpProblem(model_vec, base_dprob, Direct(); rng)
+    base_jin = JumpInputs(model_vec, u0_alts_vec[1], tspan, p_alts_vec[1])
+    base_jprob = JumpProblem(base_jin; 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)
@@ -415,7 +415,7 @@ let
     # Loops through all potential parameter sets, checking that their inputs yield errors.
     for ps in [[ps_valid]; ps_invalid], u0 in [[u0_valid]; u0s_invalid]
         # Handles all types of time-dependent systems. The `isequal` is because some case should pass.
-        for XProblem in [ODEProblem, SDEProblem, DiscreteProblem]
+        for XProblem in [ODEProblem, SDEProblem, JumpInputs]
             if isequal(ps, ps_valid) && isequal(u0, u0_valid)
                 XProblem(rn, u0, (0.0, 1.0), ps)
             else

test/upstream/mtk_structure_indexing.jl

@@ -33,7 +33,7 @@ begin
     oprob = ODEProblem(model, u0_vals, tspan, p_vals)
     sprob = SDEProblem(model,u0_vals, tspan, p_vals)
     dprob = DiscreteProblem(model, u0_vals, tspan, p_vals)
-    jprob = JumpProblem(model, deepcopy(dprob), Direct(); rng)
+    jprob = JumpProblem(JumpInputs(model, u0_vals, tspan, p_vals); rng)
     nprob = NonlinearProblem(model, u0_vals, p_vals)
     ssprob = SteadyStateProblem(model, u0_vals, p_vals)
     problems = [oprob, sprob, dprob, jprob, nprob, ssprob]
@@ -344,7 +344,7 @@ let
         ps = [k1 => 0.1, k2 => 0.2, V0 => 3.0]
         prob1 = XProblem(rs, u0, 0.001, ps; remove_conserved = true)
         Γ = prob1.f.sys.Γ
-        
+
         # Creates various `remake` version of the problem.
         prob2 = remake(prob1, u0 = [X1 => 10.0])
         prob3 = remake(prob2, u0 = [X2 => 20.0])
@@ -431,8 +431,8 @@ let
     # Creates a JumpProblem and integrator. Checks that the initial mass action rate is correct.
     u0 = [:A => 1, :B => 2, :C => 3]
     ps = [:p1 => 3.0, :p2 => 2.0]
-    dprob = DiscreteProblem(rn, u0, (0.0, 1.0), ps)
-    jprob = JumpProblem(rn, dprob, Direct())
+    jin = JumpInputs(rn, u0, (0.0, 1.0), ps)
+    jprob = JumpProblem(jin)
     jint = init(jprob, SSAStepper())
     @test jprob.massaction_jump.scaled_rates[1] == 6.0