Merge pull request #1013 from isaacsas/jumpsys_auto_alg_support

add jump auto-alg support

Project.toml

@@ -52,11 +52,11 @@ DynamicQuantities = "0.13.2"
 GraphMakie = "0.5"
 Graphs = "1.4"
 HomotopyContinuation = "2.9"
-JumpProcesses = "9.3.2"
+JumpProcesses = "9.13.1"
 LaTeXStrings = "1.3.0"
 Latexify = "0.14, 0.15, 0.16"
 MacroTools = "0.5.5"
-ModelingToolkit = "9.16.0"
+ModelingToolkit = "9.30.1"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"

docs/Project.toml

@@ -52,10 +52,10 @@ GraphMakie = "0.5"
 Graphs = "1.11.1"
 HomotopyContinuation = "2.9"
 IncompleteLU = "0.2"
-JumpProcesses = "9.11"
+JumpProcesses = "9.13.1"
 Latexify = "0.16"
 LinearSolve = "2.30"
-ModelingToolkit = "9.16.0"
+ModelingToolkit = "9.30.1"
 NonlinearSolve = "3.12"
 Optim = "1.9"
 Optimization = "3.25"

src/reactionsystem_conversions.jl

@@ -42,7 +42,7 @@ function oderatelaw(rx; combinatoric_ratelaw = true)
     rl
 end
 
-# Function returning `true` for species which shouldn't change from the reactions, 
+# Function returning `true` for species which shouldn't change from the reactions,
 # including non-species variables.
 drop_dynamics(s) = isconstant(s) || isbc(s) || (!isspecies(s))
 
@@ -454,13 +454,13 @@ COMPLETENESS_ERROR = "A ReactionSystem must be complete before it can be convert
 function check_cons_warning(remove_conserved, remove_conserved_warn)
     (remove_conserved && remove_conserved_warn) || return
     @warn "You are creating a system or problem while eliminating conserved quantities. Please note,
-        due to limitations / design choices in ModelingToolkit if you use the created system to 
+        due to limitations / design choices in ModelingToolkit if you use the created system to
         create a problem (e.g. an `ODEProblem`), or are directly creating a problem, you *should not*
-        modify that problem's initial conditions for species (e.g. using `remake`). Changing initial 
-        conditions must be done by creating a new Problem from your reaction system or the 
-        ModelingToolkit system you converted it into with the new initial condition map. 
+        modify that problem's initial conditions for species (e.g. using `remake`). Changing initial
+        conditions must be done by creating a new Problem from your reaction system or the
+        ModelingToolkit system you converted it into with the new initial condition map.
         Modification of parameter values is still possible, *except* for the modification of any
-        conservation law constants ($CONSERVED_CONSTANT_SYMBOL), which is not possible. You might 
+        conservation law constants ($CONSERVED_CONSTANT_SYMBOL), which is not possible. You might
         get this warning when creating a problem directly.
 
         You can remove this warning by setting `remove_conserved_warn = false`."
@@ -568,13 +568,13 @@ function Base.convert(::Type{<:NonlinearSystem}, rs::ReactionSystem; name = name
         kwargs...)
 end
 
-# Ideally, when `ReactionSystem`s are converted to `NonlinearSystem`s, any coupled ODEs should be 
+# Ideally, when `ReactionSystem`s are converted to `NonlinearSystem`s, any coupled ODEs should be
 # on the form D(X) ~ ..., where lhs is the time derivative w.r.t. a single variable, and the rhs
 # does not contain any differentials. If this is not the case, we throw a warning to let the user
 # know that they should be careful.
 function nonlinear_convert_differentials_check(rs::ReactionSystem)
     for eq in filter(is_diff_equation, equations(rs))
-        # For each differential equation, checks in order: 
+        # For each differential equation, checks in order:
         # If there is a differential on the right hand side.
         # If the lhs is not on the form D(...).
         # If the lhs upper level function is not a differential w.r.t. time.
@@ -585,12 +585,12 @@ function nonlinear_convert_differentials_check(rs::ReactionSystem)
            !isequal(Symbolics.operation(eq.lhs), Differential(get_iv(rs))) ||
            (length(arguments(eq.lhs)) != 1) ||
            !any(isequal(arguments(eq.lhs)[1]), nonspecies(rs))
-            error("You are attempting to convert a `ReactionSystem` coupled with differential equations to a `NonlinearSystem`. However, some of these differentials are not of the form `D(x) ~ ...` where: 
+            error("You are attempting to convert a `ReactionSystem` coupled with differential equations to a `NonlinearSystem`. However, some of these differentials are not of the form `D(x) ~ ...` where:
                     (1) The left-hand side is a differential of a single variable with respect to the time independent variable, and
                     (2) The right-hand side does not contain any differentials.
                 This is generally not permitted.
 
-                If you still would like to perform this conversion, please use the `all_differentials_permitted = true` option. In this case, all differentials will be set to `0`. 
+                If you still would like to perform this conversion, please use the `all_differentials_permitted = true` option. In this case, all differentials will be set to `0`.
                 However, it is recommended to proceed with caution to ensure that the produced nonlinear equation makes sense for your intended application."
             )
         end
@@ -783,13 +783,13 @@ function DiffEqBase.DiscreteProblem(rs::ReactionSystem, u0, tspan::Tuple,
 end
 
 # JumpProblem from AbstractReactionNetwork
-function JumpProcesses.JumpProblem(rs::ReactionSystem, prob, aggregator, args...;
-        name = nameof(rs),
+function JumpProcesses.JumpProblem(rs::ReactionSystem, prob,
+        aggregator = JumpProcesses.NullAggregator(); name = nameof(rs),
         combinatoric_ratelaws = get_combinatoric_ratelaws(rs),
         checks = false, kwargs...)
     jsys = convert(JumpSystem, rs; name, combinatoric_ratelaws, checks)
     jsys = complete(jsys)
-    return JumpProblem(jsys, prob, aggregator, args...; kwargs...)
+    return JumpProblem(jsys, prob, aggregator; kwargs...)
 end
 
 # SteadyStateProblem from AbstractReactionNetwork

test/simulation_and_solving/simulate_jumps.jl

@@ -126,22 +126,27 @@ let
     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, u0_1, ps_1, sp) in 
+    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)
 
         # 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)
-
+
+        # simulate using auto-alg
+        jprob_1b = JumpProblem(rn_catalyst, dprob_1; rng)
+        sol1b = solve(jprob_1; seed, saveat = 1.0)
+        @test mean(sol1[sp]) ≈ mean(sol1b[sp]) rtol = 1e-1
+
         # 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 
+
+        # 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
@@ -162,8 +167,8 @@ end
 
 # Tests simulating a network without parameters.
 let
-    no_param_network = @reaction_network begin 
-        (1.2, 5), X1 ↔ X2 
+    no_param_network = @reaction_network begin
+        (1.2, 5), X1 ↔ X2
     end
     u0 = rnd_u0_Int64(no_param_network, rng)
     dprob = DiscreteProblem(no_param_network, u0, (0.0, 1000.0))