diff --git a/Project.toml b/Project.toml index 8bdaf03249..95b266cd27 100644 --- a/Project.toml +++ b/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" diff --git a/docs/Project.toml b/docs/Project.toml index c04673c790..c466d8a534 100644 --- a/docs/Project.toml +++ b/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" diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl index f2c83507fc..6458b6c572 100644 --- a/src/reactionsystem_conversions.jl +++ b/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 diff --git a/test/simulation_and_solving/simulate_jumps.jl b/test/simulation_and_solving/simulate_jumps.jl index 370b51036e..b03ba76cd8 100644 --- a/test/simulation_and_solving/simulate_jumps.jl +++ b/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))