From f31902b88bb60074db8839a2c5d1249c71037b73 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 6 Aug 2024 14:24:50 -0400 Subject: [PATCH 1/6] jump auto-alg support --- Project.toml | 2 +- docs/Project.toml | 2 +- src/reactionsystem_conversions.jl | 26 +++++++++---------- test/simulation_and_solving/simulate_jumps.jl | 17 +++++++----- 4 files changed, 26 insertions(+), 21 deletions(-) diff --git a/Project.toml b/Project.toml index 8bdaf03249..6fd0a97390 100644 --- a/Project.toml +++ b/Project.toml @@ -52,7 +52,7 @@ DynamicQuantities = "0.13.2" GraphMakie = "0.5" Graphs = "1.4" HomotopyContinuation = "2.9" -JumpProcesses = "9.3.2" +JumpProcesses = "9.13" LaTeXStrings = "1.3.0" Latexify = "0.14, 0.15, 0.16" MacroTools = "0.5.5" diff --git a/docs/Project.toml b/docs/Project.toml index c04673c790..59a4c8842d 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -52,7 +52,7 @@ GraphMakie = "0.5" Graphs = "1.11.1" HomotopyContinuation = "2.9" IncompleteLU = "0.2" -JumpProcesses = "9.11" +JumpProcesses = "9.13" Latexify = "0.16" LinearSolve = "2.30" ModelingToolkit = "9.16.0" 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..e4bcbb3cb6 100644 --- a/test/simulation_and_solving/simulate_jumps.jl +++ b/test/simulation_and_solving/simulate_jumps.jl @@ -126,24 +126,29 @@ 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) + # 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 + @test mean(sol1[sp]) ≈ mean(sol1b[sp]) rtol = 1e-1 end 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)) From a3049b51a5894f6cd118da9afa4b254742d92b27 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 6 Aug 2024 14:26:59 -0400 Subject: [PATCH 2/6] reorder test --- test/simulation_and_solving/simulate_jumps.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/test/simulation_and_solving/simulate_jumps.jl b/test/simulation_and_solving/simulate_jumps.jl index e4bcbb3cb6..2804faf99c 100644 --- a/test/simulation_and_solving/simulate_jumps.jl +++ b/test/simulation_and_solving/simulate_jumps.jl @@ -138,6 +138,8 @@ let 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) @@ -148,7 +150,6 @@ let # 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 - @test mean(sol1[sp]) ≈ mean(sol1b[sp]) rtol = 1e-1 end end From 96a1eedb5d5d59e09369358a3aa760a1d1fdafd4 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 6 Aug 2024 14:27:55 -0400 Subject: [PATCH 3/6] update --- test/simulation_and_solving/simulate_jumps.jl | 1 - 1 file changed, 1 deletion(-) diff --git a/test/simulation_and_solving/simulate_jumps.jl b/test/simulation_and_solving/simulate_jumps.jl index 2804faf99c..b03ba76cd8 100644 --- a/test/simulation_and_solving/simulate_jumps.jl +++ b/test/simulation_and_solving/simulate_jumps.jl @@ -137,7 +137,6 @@ let # 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 From 9d02d44c184a6ac55f2b6e3111a11c5d4c3f22ba Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Wed, 7 Aug 2024 15:30:49 -0400 Subject: [PATCH 4/6] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 6fd0a97390..ab490d606b 100644 --- a/Project.toml +++ b/Project.toml @@ -56,7 +56,7 @@ JumpProcesses = "9.13" 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" From 15716d4cb2bb6b054b6c60d2c16b737f65b5fd30 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Wed, 7 Aug 2024 15:31:16 -0400 Subject: [PATCH 5/6] Update Project.toml --- docs/Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/docs/Project.toml b/docs/Project.toml index 59a4c8842d..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.13" +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" From b044a769cb29826aaa14f184f77cdd811436d2ce Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Wed, 7 Aug 2024 15:31:44 -0400 Subject: [PATCH 6/6] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index ab490d606b..95b266cd27 100644 --- a/Project.toml +++ b/Project.toml @@ -52,7 +52,7 @@ DynamicQuantities = "0.13.2" GraphMakie = "0.5" Graphs = "1.4" HomotopyContinuation = "2.9" -JumpProcesses = "9.13" +JumpProcesses = "9.13.1" LaTeXStrings = "1.3.0" Latexify = "0.14, 0.15, 0.16" MacroTools = "0.5.5"