From b3d1b42906c509f061c2939a11c90d67809cfe38 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sat, 10 Aug 2024 17:02:25 -0400 Subject: [PATCH 01/19] variable rate jump support --- HISTORY.md | 36 ++++++++ README.md | 4 +- docs/src/api.md | 5 +- docs/src/index.md | 8 +- .../catalyst_for_new_julia_users.md | 10 +-- .../introduction_to_catalyst.md | 22 +++-- .../examples/basic_CRN_library.md | 40 ++++----- .../smoluchowski_coagulation_equation.md | 7 +- .../parametric_stoichiometry.md | 4 +- .../examples/periodic_events_simulation.md | 10 +-- .../simulation_introduction.md | 43 ++++++---- src/Catalyst.jl | 5 +- src/reactionsystem_conversions.jl | 83 +++++++++++++++++++ test/reactionsystem_core/reactionsystem.jl | 35 ++++++-- 14 files changed, 232 insertions(+), 80 deletions(-) diff --git a/HISTORY.md b/HISTORY.md index f9ef87e663..8cd7ddce30 100644 --- a/HISTORY.md +++ b/HISTORY.md @@ -2,6 +2,42 @@ ## Catalyst unreleased (master branch) +## Catalyst 14.3 +- Support for simulating stochastic chemical kinetics models with explicitly + time-dependent propensities (i.e. where the resulting `JumpSystem` contains + `VariableRateJump`s). As such `JumpProblem`s need to be defined over + `ODEProblem`s or `SDEProblem`s instead of `DiscreteProblem`s we have + introduced a new input struct, `JumpInputs`, that handles selecting via + analysis of the generated `JumpSystem`, i.e. one can now say + ```julia + using Catalyst, OrdinaryDiffEq, JumpProcesses, Plots + rn = @reaction_network begin + k*(1 + sin(t)), 0 --> A + end + jinput = JumpInputs(rn, [:A => 0], (0.0, 10.0), [:k => .5]) + # note that jinput.prob isa ODEProblem + jprob = JumpProblem(jinput) + sol = solve(jprob, Tsit5()) + plot(sol, idxs = :A) + + rn = @reaction_network begin + k, 0 --> A + end + jinput = JumpInputs(rn, [:A => 0], (0.0, 10.0), [:k => .5]) + # note that jinput.prob isa DiscreteProblem + jprob = JumpProblem(jinput) + sol = solve(jprob) + plot(sol, idxs = :A) + ``` + When calling solve for problems with explicit time-dependent propensities, + i.e. where `jinput.prob isa ODEProblem`, note that one must currently + explicitly select an ODE solver to handle time-stepping and integrating the + time-dependent propensities. +- Note that solutions to jump problems with explicit time-dependent + propensities, i.e. a `JumpProblem` over an `ODEProblem`, require manual + selection of the variables to plot. That is, currently `plot(sol)` will error + in this case due to limitations in the SciMLBase plot recipe. + ## Catalyst 14.2 - Support for auto-algorithm selection in `JumpProblem`s. For systems with only propensities that do not have an explicit time-dependence (i.e. that are not diff --git a/README.md b/README.md index a1cf202795..6ddd82200c 100644 --- a/README.md +++ b/README.md @@ -123,8 +123,8 @@ small network in the current example ends up being Gillespie's Direct method): # The initial conditions are now integers as we track exact populations for each species. using JumpProcesses u0_integers = [:S => 50, :E => 10, :SE => 0, :P => 0] -dprob = DiscreteProblem(model, u0_integers, tspan, ps) -jprob = JumpProblem(model, dprob) +jinput = JumpInputs(model, u0_integers, tspan, ps) +jprob = JumpProblem(jinput) jump_sol = solve(jprob) plot(jump_sol; lw = 2) ``` diff --git a/docs/src/api.md b/docs/src/api.md index fcada60612..564af6bb23 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -69,8 +69,8 @@ jumpsys = complete(jumpsys) u₀map = [S => 999, I => 1, R => 0] dprob = DiscreteProblem(jumpsys, u₀map, tspan, parammap) jprob = JumpProblem(jumpsys, dprob, Direct()) -sol = solve(jprob, SSAStepper()) -sol = solve(jprob, SSAStepper(), seed = 1234) # hide +sol = solve(jprob) +sol = solve(jprob; seed = 1234) # hide p3 = plot(sol, title = "jump") plot(p1, p2, p3; layout = (3,1)) Catalyst.PNG(plot(p1, p2, p3; layout = (3,1), fmt = :png, dpi = 200)) # hide @@ -318,6 +318,7 @@ hillar ## Transformations ```@docs Base.convert +JumpInputs ModelingToolkit.structural_simplify set_default_noise_scaling ``` diff --git a/docs/src/index.md b/docs/src/index.md index 06d93dca9f..5e03723904 100644 --- a/docs/src/index.md +++ b/docs/src/index.md @@ -134,10 +134,10 @@ The same model can be used as input to other types of simulations. E.g. here we # The initial conditions are now integers as we track exact populations for each species. using JumpProcesses u0_integers = [:S => 50, :E => 10, :SE => 0, :P => 0] -dprob = DiscreteProblem(model, u0_integers, tspan, ps) -jprob = JumpProblem(model, dprob, Direct()) -jump_sol = solve(jprob, SSAStepper()) -jump_sol = solve(jprob, SSAStepper(); seed = 1234) # hide +jinput = JumpInputs(model, u0_integers, tspan, ps) +jprob = JumpProblem(jinput) +jump_sol = solve(jprob) +jump_sol = solve(jprob; seed = 1234) # hide plot(jump_sol; lw = 2) ``` diff --git a/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md b/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md index 061c33f93e..872bcc27c7 100644 --- a/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md +++ b/docs/src/introduction_to_catalyst/catalyst_for_new_julia_users.md @@ -166,16 +166,14 @@ params = [:b => 0.2, :k => 1.0] nothing # hide ``` -Previously we have bundled this information into an `ODEProblem` (denoting a deterministic *ordinary differential equation*). Now we wish to simulate our model as a jump process (where each reaction event corresponds to a discrete change in the state of the system). We do this by first creating a `DiscreteProblem`, and then using this as an input to a `JumpProblem`. +Previously we have bundled this information into an `ODEProblem` (denoting a deterministic *ordinary differential equation*). Now we wish to simulate our model as a jump process (where each reaction event corresponds to a discrete change in the state of the system). We do this by first processing the inputs to work in a jump model -- an extra step needed for jump models that can be avoided for ODE/SDE models -- and then creating a `JumpProblem` from the inputs: ```@example ex2 using JumpProcesses # hide -dprob = DiscreteProblem(sir_model, u0, tspan, params) -jprob = JumpProblem(sir_model, dprob) +jinput = JumpInputs(sir_model, u0, tspan, params) +jprob = JumpProblem(jinput) nothing # hide ``` -Again, the order in which the inputs are given to the `DiscreteProblem` and the `JumpProblem` is important. - -Finally, we can simulate our model using the `solve` function, and plot the solution using the `plot` function. +Finally, we can now simulate our model using the `solve` function, and plot the solution using the `plot` function. ```@example ex2 sol = solve(jprob) sol = solve(jprob; seed=1234) # hide diff --git a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md index c2e21c3d63..11addc6057 100644 --- a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md +++ b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md @@ -202,11 +202,12 @@ using JumpProcesses # redefine the initial condition to be integer valued u₀map = [:m₁ => 0, :m₂ => 0, :m₃ => 0, :P₁ => 20, :P₂ => 0, :P₃ => 0] -# next we create a discrete problem to encode that our species are integer-valued: -dprob = DiscreteProblem(rn, u₀map, tspan, pmap) +# next we process the inputs for the jump problem +jinputs = JumpInputs(rn, u₀map, tspan, pmap) -# now, we create a JumpProblem, and specify Gillespie's Direct Method as the solver: -jprob = JumpProblem(rn, dprob, Direct()) +# now, we create a JumpProblem, and let a solver be chosen for us automatically +# in this case Gillespie's Direct Method will be selected +jprob = JumpProblem(jinputs) # now, let's solve and plot the jump process: sol = solve(jprob) @@ -215,10 +216,17 @@ Catalyst.PNG(plot(sol; fmt = :png, dpi = 200)) # hide ``` We see that oscillations remain, but become much noisier. Note, in constructing -the `JumpProblem` we could have used any of the SSAs that are part of JumpProcesses -instead of the `Direct` method, see the list of SSAs (i.e., constant rate jump -aggregators) in the +the `JumpProblem` we could have specified any of the SSAs that are part of +JumpProcesses instead of letting JumpProcesses auto-select a solver, see the +list of SSAs (i.e., constant rate jump aggregators) in the [documentation](https://docs.sciml.ai/JumpProcesses/stable/jump_types/#Jump-Aggregators-for-Exact-Simulation). +For example, to choose the `SortingDirect` method we would instead say +```@example tut1 +jprob = JumpProblem(jinputs, SortingDirect()) +sol = solve(jprob) +plot(sol) +Catalyst.PNG(plot(sol; fmt = :png, dpi = 200)) # hide +``` Common questions that arise in using the JumpProcesses SSAs (i.e. Gillespie methods) are collated in the [JumpProcesses FAQ](https://docs.sciml.ai/JumpProcesses/stable/faq/). diff --git a/docs/src/model_creation/examples/basic_CRN_library.md b/docs/src/model_creation/examples/basic_CRN_library.md index 4198fbf07f..618e195b1c 100644 --- a/docs/src/model_creation/examples/basic_CRN_library.md +++ b/docs/src/model_creation/examples/basic_CRN_library.md @@ -34,10 +34,10 @@ nothing # hide Next, a stochastic chemical kinetics-based jump simulation: ```@example crn_library_birth_death using JumpProcesses -dprob = DiscreteProblem(bd_process, u0, tspan, ps) -jprob = JumpProblem(bd_process, dprob, Direct()) -jsol = solve(jprob, SSAStepper()) -jsol = solve(jprob, SSAStepper(); seed = 12) # hide +jinput = JumpInputs(bd_process, u0, tspan, ps) +jprob = JumpProblem(jinput) +jsol = solve(jprob) +jsol = solve(jprob; seed = 12) # hide nothing # hide ``` Finally, we plot the results: @@ -106,10 +106,10 @@ ssol = solve(sprob, STrapezoid()) ssol = solve(sprob, STrapezoid(); seed = 12) # hide using JumpProcesses -dprob = DiscreteProblem(mm_system, u0, tspan, ps) -jprob = JumpProblem(mm_system, dprob, Direct()) -jsol = solve(jprob, SSAStepper()) -jsol = solve(jprob, SSAStepper(), seed = 12) # hide +jinput = JumpInputs(mm_system, u0, tspan, ps) +jprob = JumpProblem(jinput) +jsol = solve(jprob) +jsol = solve(jprob; seed = 12) # hide using Plots oplt = plot(osol; title = "Reaction rate equation (ODE)") @@ -145,15 +145,15 @@ plot(osol; title = "Reaction rate equation (ODE)", size=(800,350)) Next, we perform 3 different Jump simulations. Note that for the stochastic model, the occurrence of an outbreak is not certain. Rather, there is a possibility that it fizzles out without a noteworthy peak. ```@example crn_library_sir using JumpProcesses -dprob = DiscreteProblem(sir_model, u0, tspan, ps) -jprob = JumpProblem(sir_model, dprob, Direct()) +jinput = JumpInputs(sir_model, u0, tspan, ps) +jprob = JumpProblem(jinput) -jsol1 = solve(jprob, SSAStepper()) -jsol2 = solve(jprob, SSAStepper()) -jsol3 = solve(jprob, SSAStepper()) -jsol1 = solve(jprob, SSAStepper(), seed = 1) # hide -jsol2 = solve(jprob, SSAStepper(), seed = 2) # hide -jsol3 = solve(jprob, SSAStepper(), seed = 3) # hide +jsol1 = solve(jprob) +jsol2 = solve(jprob) +jsol3 = solve(jprob) +jsol1 = solve(jprob, seed = 1) # hide +jsol2 = solve(jprob, seed = 2) # hide +jsol3 = solve(jprob, seed = 3) # hide jplt1 = plot(jsol1; title = "Outbreak") jplt2 = plot(jsol2; title = "Outbreak") @@ -242,10 +242,10 @@ ps = [:v₀ => 0.1, :v => 2.0, :K => 10.0, :n => 2, :d => 0.1] oprob = ODEProblem(sa_loop, u0, tspan, ps) osol = solve(oprob) -dprob = DiscreteProblem(sa_loop, u0, tspan, ps) -jprob = JumpProblem(sa_loop, dprob, Direct()) -jsol = solve(jprob, SSAStepper()) -jsol = solve(jprob, SSAStepper(), seed = 12) # hide +jinput = JumpInputs(sa_loop, u0, tspan, ps) +jprob = JumpProblem(jinput) +jsol = solve(jprob) +jsol = solve(jprob, seed = 12) # hide fplt = plot(osol; lw = 3, label = "Reaction rate equation (ODE)") plot!(fplt, jsol; lw = 3, label = "Stochastic chemical kinetics (Jump)", yguide = "X", size = (800,350)) diff --git a/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md b/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md index d8865651b4..6b0ab3de76 100644 --- a/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md +++ b/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md @@ -104,10 +104,9 @@ rs = complete(rs) We now convert the [`ReactionSystem`](@ref) into a `ModelingToolkit.JumpSystem`, and solve it using Gillespie's direct method. For details on other possible solvers (SSAs), see the [DifferentialEquations.jl](https://docs.sciml.ai/DiffEqDocs/stable/types/jump_types/) documentation ```@example smcoag1 # solving the system -jumpsys = complete(convert(JumpSystem, rs)) -dprob = DiscreteProblem(rs, u₀map, tspan) -jprob = JumpProblem(jumpsys, dprob, Direct(), save_positions = (false, false)) -jsol = solve(jprob, SSAStepper(), saveat = tspan[2] / 30) +jinputs = JumpInputs(rs, u₀map, tspan) +jprob = JumpProblem(jinputs, Direct(); save_positions = (false, false)) +jsol = solve(jprob; saveat = tspan[2] / 30) nothing #hide ``` Lets check the results for the first three polymers/cluster sizes. We compare to the analytical solution for this system: diff --git a/docs/src/model_creation/parametric_stoichiometry.md b/docs/src/model_creation/parametric_stoichiometry.md index 47ca2e2d7f..0eac502aae 100644 --- a/docs/src/model_creation/parametric_stoichiometry.md +++ b/docs/src/model_creation/parametric_stoichiometry.md @@ -179,7 +179,7 @@ u₀ = symmap_to_varmap(jsys, [:G₊ => 1, :G₋ => 0, :P => 1]) tspan = (0., 6.0) # time interval to solve over dprob = DiscreteProblem(jsys, u₀, tspan, p) jprob = JumpProblem(jsys, dprob, Direct()) -sol = solve(jprob, SSAStepper()) +sol = solve(jprob) plot(sol.t, sol[jsys.P], legend = false, xlabel = "time", ylabel = "P(t)") ``` To double check our results are consistent with MomentClosure.jl, let's @@ -192,7 +192,7 @@ function getmean(jprob, Nsims, tv) Pmean = zeros(length(tv)) @variables P(t) for n in 1:Nsims - sol = solve(jprob, SSAStepper()) + sol = solve(jprob) Pmean .+= sol(tv, idxs=P) end Pmean ./= Nsims diff --git a/docs/src/model_simulation/examples/periodic_events_simulation.md b/docs/src/model_simulation/examples/periodic_events_simulation.md index b9bbfd9e76..ab091fac76 100644 --- a/docs/src/model_simulation/examples/periodic_events_simulation.md +++ b/docs/src/model_simulation/examples/periodic_events_simulation.md @@ -23,7 +23,7 @@ plot(sol) ``` ## [Modelling a circadian periodic event in a jump simulation](@id periodic_event_simulation_example_ode) -While defining periodic events is easy for ODE and SDE simulations, due to how events are internally implemented, these cannot currently be used for jump simulations. Instead, there is a workaround which includes firs creating a [conditional discrete event](@ref constraint_equations_events) which is designed to trigger every 13 time units: +While defining periodic events is easy for ODE and SDE simulations, due to how events are internally implemented, these cannot currently be used for jump simulations. Instead, there is a workaround which includes first creating a [conditional discrete event](@ref constraint_equations_events) which is designed to trigger every 13 time units: ```@example periodic_event_example circadian_model = @reaction_network begin @discrete_events begin @@ -35,20 +35,20 @@ end using JumpProcesses u0 = [:X => 150, :Xᴾ => 50] ps = [:kₚ => 0.1, :kᵢ => 0.1, :l => 1.0] -dprob = DiscreteProblem(circadian_model, u0, (0.0, 100.0), ps) -jprob = JumpProblem(circadian_model, dprob, Direct()) +jinput = JumpInputs(circadian_model, u0, (0.0, 100.0), ps) +jprob = JumpProblem(jinput) nothing # hide ``` Next, if we simulate our model, we note that the events do not seem to be triggered ```@example periodic_event_example -sol = solve(jprob, SSAStepper()) +sol = solve(jprob) plot(sol) plot(sol, density = 10000, fmt = :png) # hide ``` The reason is that discrete callbacks' conditions are only checked at the end of each simulation time step, and the probability that these will coincide with the callback's trigger times ($t = 12, 24, 36, ...$) is infinitely small. Hence, we must directly instruct our simulation to stop at these time points. We can do this using the `tstops` argument: ```@example periodic_event_example tstops = 12.0:12.0:dprob.tspan[2] -sol = solve(jprob, SSAStepper(); tstops) +sol = solve(jprob; tstops) plot(sol) plot(sol, density = 10000, fmt = :png) # hide ``` diff --git a/docs/src/model_simulation/simulation_introduction.md b/docs/src/model_simulation/simulation_introduction.md index f2aafd8ec7..f417d83794 100644 --- a/docs/src/model_simulation/simulation_introduction.md +++ b/docs/src/model_simulation/simulation_introduction.md @@ -1,5 +1,5 @@ # [Model Simulation Introduction](@id simulation_intro) -Catalyst's core functionality is the creation of *chemical reaction network* (CRN) models that can be simulated using ODE, SDE, and jump simulations. How such simulations are carried out has already been described in [Catalyst's introduction](@ref introduction_to_catalyst). This page provides a deeper introduction, giving some additional background and introducing various simulation-related options. +Catalyst's core functionality is the creation of *chemical reaction network* (CRN) models that can be simulated using ODE, SDE, and jump simulations. How such simulations are carried out has already been described in [Catalyst's introduction](@ref introduction_to_catalyst). This page provides a deeper introduction, giving some additional background and introducing various simulation-related options. Here we will focus on the basics, with other sections of the simulation documentation describing various specialised features, or giving advice on performance. Anyone who plans on using Catalyst's simulation functionality extensively is recommended to also read the documentation on [solution plotting](@ref simulation_plotting), and on how to [interact with simulation problems, integrators, and solutions](@ref simulation_structure_interfacing). Anyone with an application for which performance is critical should consider reading the corresponding page on performance advice for [ODEs](@ref ode_simulation_performance) or [SDEs](@ref sde_simulation_performance). @@ -133,7 +133,7 @@ Here follows a list of solver options which might be of interest to the user. - `adaptive`: Toggles adaptive time stepping for valid methods. Default to `true`. - `dt`: For non-adaptive simulations, sets the step size (also sets the initial step size for adaptive methods). - `saveat`: Determines the time points at which the simulation is saved. E.g. for `saveat = 2.0` the simulation is saved every second time unit. If not given, the solution is saved after each time step. -- `save_idxs`: Provides a vector of species whose values should be saved during the simulation. E.g. for `save_idxs = [:X1]`, only the value of species $X1$ is saved. +- `save_idxs`: Provides a vector of species whose values should be saved during the simulation. E.g. for `save_idxs = [:X1]`, only the value of species $X1$ is saved. - `maxiters`: The maximum number of time steps of the simulation. If this number is reached, the simulation is terminated. - `seed`: Sets a seed for stochastic simulations. Stochastic simulations with the same seed generate identical results. @@ -156,7 +156,7 @@ ps = Dict([:k1 => 2.0, :k2 => 5.0]) oprob = ODEProblem(two_state_model, u0, tspan, ps) nothing # hide ``` -The forms used for `u0` and `ps` does not need to be the same (but can e.g. be a vector and a tuple). +The forms used for `u0` and `ps` does not need to be the same (but can e.g. be a vector and a tuple). !!! note It is possible to [designate specific types for parameters](@ref dsl_advanced_options_parameter_types). When this is done, the tuple form for providing parameter values should be preferred. @@ -186,7 +186,7 @@ sol = solve(sprob, STrapezoid()) sol = solve(sprob, STrapezoid(); seed = 123) # hide plot(sol) ``` -we can see that while this simulation (unlike the ODE ones) exhibits some fluctuations. +we can see that while this simulation (unlike the ODE ones) exhibits some fluctuations. !!! note Unlike for ODE and jump simulations, there are no good heuristics for automatically selecting suitable SDE solvers. Hence, for SDE simulations a solver must be provided. `STrapezoid` will work for a large number of cases. When this is not the case, however, please check the list of [available SDE solvers](https://docs.sciml.ai/DiffEqDocs/stable/solvers/sde_solve/) for a suitable alternative (making sure to select one compatible with non-diagonal noise and the [Ito interpretation]https://en.wikipedia.org/wiki/It%C3%B4_calculus). @@ -292,7 +292,7 @@ While the `@default_noise_scaling` option is unavailable for [programmatically c Catalyst uses the [JumpProcesses.jl](https://github.com/SciML/JumpProcesses.jl) package to perform jump simulations. This section provides a brief introduction, with [JumpProcesses's documentation](https://docs.sciml.ai/JumpProcesses/stable/) providing a more extensive description. -Jump simulations are performed using so-called `JumpProblem`s. Unlike ODEs and SDEs (for which the corresponding problem types can be created directly), jump simulations require first creating an intermediary `DiscreteProblem`. In this example, we first declare our two-state model and its initial conditions, time span, and parameter values. +Jump simulations are performed using so-called `JumpProblem`s. Unlike ODEs and SDEs (for which the corresponding problem types can be created directly), jump simulations require first processing inputs into a correct format creating an intermediary `JumpInputs`. In this example, we first declare our two-state model and its initial conditions, time span, and parameter values. ```@example simulation_intro_jumps using Catalyst, JumpProcesses, Plots two_state_model = @reaction_network begin @@ -306,19 +306,19 @@ nothing # hide !!! note Since jump simulations typically simulate the integer copy-numbers of each species present in the system, we designate our initial conditions for jump simulations as integers. Decimal-numbered initial conditions (and thus jump simulations) are, however, also possible. While ODE and SDE simulations accept integer initial conditions, these will be converted to decimal numbers. -Next, we bundle these into a `DiscreteProblem` (similarly to how `ODEProblem`s and `SDEProblem`s are created): +Next, we process these into a `JumpInputs`: ```@example simulation_intro_jumps -dprob = DiscreteProblem(two_state_model, u0, tspan, ps) +jinput = JumpInputs(two_state_model, u0, tspan, ps) nothing # hide ``` This is then used as input to a `JumpProblem`. The `JumpProblem` also requires the CRN model and [an aggregator](@ref simulation_intro_jumps_solver_designation) as input. ```@example simulation_intro_jumps -jprob = JumpProblem(two_state_model, dprob, Direct()) +jprob = JumpProblem(jinput) nothing # hide ``` The `JumpProblem` can now be simulated using `solve` (just like any other problem type). ```@example simulation_intro_jumps -sol = solve(jprob, SSAStepper()) +sol = solve(jprob) nothing # hide ``` If we plot the solution we can see how the system's state does not change continuously, but instead in discrete jumps (due to the occurrence of the individual reactions of the system). @@ -328,27 +328,34 @@ plot(sol) ``` ### [Designating aggregators and simulation methods for jump simulations](@id simulation_intro_jumps_solver_designation) -Jump simulations (just like ODEs and SDEs) are performed using solver methods. Unlike ODEs and SDEs, jump simulations are carried out by two different types of methods acting in tandem. First, an *aggregator* method is used to (after each reaction) determine the time to, and type of, the next reaction. Next, a simulation method is used to actually carry out the simulation. +Jump simulations (just like ODEs and SDEs) are performed using stochastic simulation algorithms (SSAs) to generate exact samples of the underlying jump process. In JumpProcesses.jl and Catalyst, we call SSAs *aggregators*. These methods determine the time until, and type of, the next reaction in a system. A separate time-stepping method is then used to actually step from one reaction instance to the next. -Several different aggregators are available (a full list is provided [here](https://docs.sciml.ai/JumpProcesses/stable/jump_types/#Jump-Aggregators-for-Exact-Simulation)). To designate a specific one, provide it as the third argument to the `JumpProblem`. E.g. to designate that the sorting direct method (`SortingDirect`) should be used, use: +Several different aggregators are available (a full list is provided [here](https://docs.sciml.ai/JumpProcesses/stable/jump_types/#Jump-Aggregators-for-Exact-Simulation)). To designate a specific one, provide it as the second argument to the `JumpProblem`. E.g. to designate that the sorting direct method (`SortingDirect`) should be used, use: ```@example simulation_intro_jumps -jprob = JumpProblem(two_state_model, dprob, SortingDirect()) +jprob = JumpProblem(jinput, SortingDirect()) nothing # hide ``` -Especially for large systems, the choice of aggregator is relevant to simulation performance. - -Next, a simulation method can be provided (like for ODEs and SDEs) as the second argument to `solve`. Currently, the only relevant solver is `SSAStepper()` (which is the one used throughout Catalyst's documentation). Other choices are primarily relevant to combined ODE/SDE + jump simulations, or inexact simulations. These situations are described in more detail [here](https://docs.sciml.ai/JumpProcesses/stable/jump_solve/). +Especially for large systems, the choice of aggregator can dramatically impact +simulation performance. ### [Jump simulations where some rate depends on time](@id simulation_intro_jumps_variableratejumps) -For some models, the rate of some reactions depend on time. E.g. consider the following [circadian model](https://en.wikipedia.org/wiki/Circadian_rhythm), where the production rate of some protein ($P$) depends on a sinusoid function: +For some models, the rate terms of reactions may explicitly depend on time. E.g. consider the following [circadian clock (inspired) model](https://en.wikipedia.org/wiki/Circadian_rhythm), where the production rate of some protein ($P$) depends on a sinusoid function: ```@example simulation_intro_jumps circadian_model = @reaction_network begin A*(sin(2π*f*t - ϕ)+1)/2, 0 --> P d, P --> 0 end ``` -This type of model will generate so called *variable rate jumps*. Simulation of such model is non-trivial (and Catalyst currently lacks a good interface for this). A detailed description of how to carry out jump simulations for models with time-dependant rates can be found [here](https://docs.sciml.ai/JumpProcesses/stable/tutorials/simple_poisson_process/#VariableRateJumps-for-processes-that-are-not-constant-between-jumps). - +This type of model will generate so called *variable rate jumps* (`VariableRateJump`s in JumpProcesses.jl). Such models can be simulated in Catalyst too, but note that now a method for time-stepping the solver must be provided to `solve`. Here ODE solvers should be given as they are used to handle integrating the explicitly time-dependent propensities for problems with variable rates, i.e. the proceeding example can be solved like +```@example simulation_intro_jumps +u0map = [:P => 0] +pmap = [:f => 1.0, :A => 2.0, :ϕ => 0.0, :d => 1.0] +tspan = (0.0, 24.0) +jinputs = JumpInputs(circadian_model, u0map, tspan, pmap) +jprob = JumpProblem(jinputs) +sol = solve(jprob, Tsit5()) # use the Tsit5 ODE solver to time-step +plot(sol; idxs = :P, lw = 2) +``` --- ## [Citation](@id simulation_intro_citation) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index c82373e234..341132e0f0 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -112,9 +112,8 @@ export params, numparams # Conversions of the `ReactionSystem` structure. include("reactionsystem_conversions.jl") -export ODEProblem, - SDEProblem, JumpProblem, NonlinearProblem, DiscreteProblem, - SteadyStateProblem +export ODEProblem, SDEProblem, JumpProblem, NonlinearProblem, DiscreteProblem, + SteadyStateProblem, JumpInputs export ismassaction, oderatelaw, jumpratelaw export symmap_to_varmap diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl index 6458b6c572..90ef509245 100644 --- a/src/reactionsystem_conversions.jl +++ b/src/reactionsystem_conversions.jl @@ -769,6 +769,81 @@ function DiffEqBase.SDEProblem(rs::ReactionSystem, u0, tspan, noise_rate_prototype = p_matrix, kwargs...) end +""" +$(TYPEDEF) + +Inputs for a JumpProblem from a given `ReactionSystem`. + +# Fields +$(FIELDS) +""" +struct JumpInputs{S <: MT.JumpSystem, T <: SciMLBase.AbstractODEProblem} + """The `JumpSystem` to define the problem over""" + sys::S + """The problem the JumpProblem should be defined over, for example DiscreteProblem""" + prob::T +end + + +""" + jumpinput = JumpInputs(rs::ReactionSystem, u0map, tspan, + pmap = DiffEqBase.NullParameters; + name = nameof(rs), + combinatoric_ratelaws = get_combinatoric_ratelaws(rs), + checks = false, kwargs...) + + Constructs the input to build a JumpProblem for the given reaction system. + +Example: +```julia +using Catalyst, OrdinaryDiffEq, JumpProcesses, Plots +rn = @reaction_network begin + k*(1 + sin(t)), 0 --> A +end +jinput = JumpInputs(rn, [:A => 0], (0.0, 10.0), [:k => .5]) +@assert jinput.prob isa ODEProblem +jprob = JumpProblem(jinput) +sol = solve(jprob, Tsit5()) +plot(sol, idxs = :A) + +rn = @reaction_network begin + k, 0 --> A +end +jinput = JumpInputs(rn, [:A => 0], (0.0, 10.0), [:k => .5]) +@assert jinput.prob isa DiscreteProblem +jprob = JumpProblem(jinput) +sol = solve(jprob) +plot(sol, idxs = :A) +``` +""" +function JumpInputs(rs::ReactionSystem, u0, tspan, p = DiffEqBase.NullParameters(); + name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), + checks = false, kwargs...) + + u0map = symmap_to_varmap(rs, u0) + pmap = symmap_to_varmap(rs, p) + jsys = complete(convert(JumpSystem, rs; name, combinatoric_ratelaws, checks)) + if MT.has_variableratejumps(jsys) + prob = ODEProblem(jsys, u0map, tspan, pmap; kwargs...) + else + prob = DiscreteProblem(jsys, u0map, tspan, pmap; kwargs...) + end + JumpInputs(jsys, prob) +end + +function Base.summary(io::IO, jinputs::JumpInputs) + type_color, no_color = SciMLBase.get_colorizers(io) + print(io, + type_color, nameof(typeof(jinputs)), + no_color, " storing" , "\n", + no_color, " JumpSystem: ", type_color, nameof(jinputs.sys), "\n", + no_color, " Problem type: ", type_color, nameof(typeof(jinputs.prob))) +end + +function Base.show(io::IO, mime::MIME"text/plain", jinputs::JumpInputs) + summary(io, jinputs) +end + # DiscreteProblem from AbstractReactionNetwork function DiffEqBase.DiscreteProblem(rs::ReactionSystem, u0, tspan::Tuple, p = DiffEqBase.NullParameters(), args...; @@ -792,6 +867,14 @@ function JumpProcesses.JumpProblem(rs::ReactionSystem, prob, return JumpProblem(jsys, prob, aggregator; kwargs...) end +# JumpProblem for JumpInputs +function JumpProcesses.JumpProblem(jinputs::JumpInputs, + agg::JumpProcesses.AbstractAggregatorAlgorithm = JumpProcesses.NullAggregator(); + kwargs...) + JumpProblem(jinputs.sys, jinputs.prob, agg; kwargs...) +end + + # SteadyStateProblem from AbstractReactionNetwork function DiffEqBase.SteadyStateProblem(rs::ReactionSystem, u0, p = DiffEqBase.NullParameters(), args...; diff --git a/test/reactionsystem_core/reactionsystem.jl b/test/reactionsystem_core/reactionsystem.jl index 1ba9142bf1..685b373ffe 100644 --- a/test/reactionsystem_core/reactionsystem.jl +++ b/test/reactionsystem_core/reactionsystem.jl @@ -268,7 +268,7 @@ end # Checks that it works for programmatic/dsl-based modelling. # Checks that all forms of model input (parameter/initial condition and vector/non-vector) are # handled properly. -let +let # Declares programmatic model. @parameters p[1:2] k d1 d2 @species (X(t))[1:2] Y1(t) Y2(t) @@ -284,7 +284,7 @@ let # Declares DSL-based model. rs_dsl = @reaction_network rs begin - @parameters p[1:2] k d1 d2 + @parameters p[1:2] k d1 d2 @species (X(t))[1:2] Y1(t) Y2(t) (p[1],p[2]), 0 --> (X[1],X[2]) k, (X[1],X[2]) --> (Y1,Y2) @@ -632,15 +632,15 @@ end let # Stores a parameter, a species, and a variable (with identical names) in different variables. x_p = let - only(@parameters x) + only(@parameters x) end x_sp = let - only(@species x(t)) + only(@species x(t)) end x_v = let - only(@variables x(t)) + only(@variables x(t)) end - + # Checks that creating systems with different in combination produces errors. # Currently broken on MTK, potentially fix in Catalyst once sorted out there (https://github.com/SciML/ModelingToolkit.jl/issues/2883). @parameters d @@ -991,7 +991,7 @@ let Catalyst.reset!(nps) δ = deficiency(rn) @test size(nps.incidencemat) == (3,3) - + Catalyst.reset!(nps) δ_l = linkagedeficiencies(rn) @test size(nps.incidencemat) == (3,3) @@ -1004,3 +1004,24 @@ let weakrev = isweaklyreversible(rn, sns) @test size(nps.incidencemat) == (3,3) end + +# test JumpInputs function auto problem selection +let + rn = @reaction_network begin + k*(1 + sin(t)), 0 --> A + end + jinput = JumpInputs(rn, [:A => 0], (0.0, 10.0), [:k => .5]) + @test jinput.prob isa ODEProblem + jprob = JumpProblem(jinput; rng) + sol = solve(jprob, Tsit5()) + @test sol(10.0; idxs = :A) > 0 + + rn = @reaction_network begin + k, 0 --> A + end + jinput = JumpInputs(rn, [:A => 0], (0.0, 10.0), [:k => .5]) + @test jinput.prob isa DiscreteProblem + jprob = JumpProblem(jinput; rng) + sol = solve(jprob) + @test sol(10.0; idxs = :A) > 0 +end \ No newline at end of file From 1f3e37b1eda42638feb22ba8d1b588531d732582 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sat, 10 Aug 2024 17:26:02 -0400 Subject: [PATCH 02/19] fix doc errs --- .../model_creation/constraint_equations.md | 13 +++++----- .../examples/periodic_events_simulation.md | 9 +++---- .../ode_simulation_performance.md | 24 +++++++++---------- .../simulation_introduction.md | 1 + 4 files changed, 25 insertions(+), 22 deletions(-) diff --git a/docs/src/model_creation/constraint_equations.md b/docs/src/model_creation/constraint_equations.md index b8468dc62a..35097cabd7 100644 --- a/docs/src/model_creation/constraint_equations.md +++ b/docs/src/model_creation/constraint_equations.md @@ -23,7 +23,7 @@ There are several ways we can create our Catalyst model with the two reactions and ODE for $V(t)$. One approach is to use compositional modeling, create separate `ReactionSystem`s and `ODESystem`s with their respective components, and then extend the `ReactionSystem` with the `ODESystem`. Let's begin by -creating these two systems. +creating these two systems. Here, to create differentials with respect to time (for our differential equations), we must import the time differential operator from Catalyst. We do this through `D = default_time_deriv()`. Here, `D(V)` denotes the differential of the variable `V` with respect to time. @@ -129,9 +129,10 @@ rs = complete(rs) oprob = ODEProblem(rs, [], (0.0, 10.0)) sol = solve(oprob, Tsit5()) -plot(sol; plotdensity = 1000) +plot(sol) +Catalyst.PNG(plot(sol; fmt = :png, dpi = 200)) # hide ``` -We can also model discrete events. Similar to our example with continuous events, we start by creating reaction equations, parameters, variables, and unknowns. +We can also model discrete events. Similar to our example with continuous events, we start by creating reaction equations, parameters, variables, and unknowns. ```@example ceq3 t = default_t() @parameters k_on switch_time k_off @@ -139,7 +140,7 @@ t = default_t() rxs = [(@reaction k_on, A --> B), (@reaction k_off, B --> A)] ``` -Now we add an event such that at time `t` (`switch_time`), `k_on` is set to zero. +Now we add an event such that at time `t` (`switch_time`), `k_on` is set to zero. ```@example ceq3 discrete_events = (t == switch_time) => [k_on ~ 0.0] @@ -147,7 +148,7 @@ u0 = [:A => 10.0, :B => 0.0] tspan = (0.0, 4.0) p = [k_on => 100.0, switch_time => 2.0, k_off => 10.0] ``` -Simulating our model, +Simulating our model, ```@example ceq3 @named osys = ReactionSystem(rxs, t, [A, B], [k_on, k_off, switch_time]; discrete_events) osys = complete(osys) @@ -156,4 +157,4 @@ oprob = ODEProblem(osys, u0, tspan, p) sol = solve(oprob, Tsit5(); tstops = 2.0) plot(sol) ``` -Note that for discrete events we need to set a stop time, `tstops`, so that the ODE solver can step exactly to the specific time of our event. For a detailed discussion on how to directly use the lower-level but more flexible DifferentialEquations.jl event/callback interface, see the [tutorial](https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Event-handling-using-callbacks) on event handling using callbacks. \ No newline at end of file +Note that for discrete events we need to set a stop time, `tstops`, so that the ODE solver can step exactly to the specific time of our event. For a detailed discussion on how to directly use the lower-level but more flexible DifferentialEquations.jl event/callback interface, see the [tutorial](https://docs.sciml.ai/Catalyst/stable/catalyst_applications/advanced_simulations/#Event-handling-using-callbacks) on event handling using callbacks. \ No newline at end of file diff --git a/docs/src/model_simulation/examples/periodic_events_simulation.md b/docs/src/model_simulation/examples/periodic_events_simulation.md index ab091fac76..0b6a9126a5 100644 --- a/docs/src/model_simulation/examples/periodic_events_simulation.md +++ b/docs/src/model_simulation/examples/periodic_events_simulation.md @@ -35,7 +35,8 @@ end using JumpProcesses u0 = [:X => 150, :Xᴾ => 50] ps = [:kₚ => 0.1, :kᵢ => 0.1, :l => 1.0] -jinput = JumpInputs(circadian_model, u0, (0.0, 100.0), ps) +tspan = (0.0, 100.0) +jinput = JumpInputs(circadian_model, u0, tspan, ps) jprob = JumpProblem(jinput) nothing # hide ``` @@ -43,14 +44,14 @@ Next, if we simulate our model, we note that the events do not seem to be trigge ```@example periodic_event_example sol = solve(jprob) plot(sol) -plot(sol, density = 10000, fmt = :png) # hide +Catalyst.PNG(plot(sol; fmt = :png, dpi = 200)) # hide ``` The reason is that discrete callbacks' conditions are only checked at the end of each simulation time step, and the probability that these will coincide with the callback's trigger times ($t = 12, 24, 36, ...$) is infinitely small. Hence, we must directly instruct our simulation to stop at these time points. We can do this using the `tstops` argument: ```@example periodic_event_example -tstops = 12.0:12.0:dprob.tspan[2] +tstops = 12.0:12.0:tspan[2] sol = solve(jprob; tstops) plot(sol) -plot(sol, density = 10000, fmt = :png) # hide +Catalyst.PNG(plot(sol; fmt = :png, dpi = 200)) # hide ``` ## [Plotting the light level](@id periodic_event_simulation_plotting_light) diff --git a/docs/src/model_simulation/ode_simulation_performance.md b/docs/src/model_simulation/ode_simulation_performance.md index 11934526ff..0dbe9877a0 100644 --- a/docs/src/model_simulation/ode_simulation_performance.md +++ b/docs/src/model_simulation/ode_simulation_performance.md @@ -31,7 +31,7 @@ oprob = ODEProblem(brusselator, u0, tspan, ps) sol1 = solve(oprob, Tsit5()) plot(sol1) -plot(sol1, plotdensity = 1000, fmt = :png) # hide +Catalyst.PNG(plot(sol1; fmt = :png, dpi = 200)) # hide ``` We get a warning, indicating that the simulation was terminated. Furthermore, the resulting plot ends at $t ≈ 12$, meaning that the simulation was not completed (as the simulation's endpoint is $t = 20$). Indeed, we can confirm this by checking the *return code* of the solution object: @@ -42,11 +42,11 @@ Next, we instead try the `Rodas5P` solver (which is designed for stiff problems) ```@example ode_simulation_performance_1 sol2 = solve(oprob, Rodas5P()) plot(sol2) -``` +``` This time the simulation was successfully completed, which can be confirmed by checking the return code: ```@example ode_simulation_performance_1 sol2.retcode -``` +``` Generally, ODE solvers can be divided into [*explicit* and *implicit* solvers](https://en.wikipedia.org/wiki/Explicit_and_implicit_methods). Roughly, explicit solvers are better for non-stiff problems, with implicit solvers being required for stiff problems. While we could use implicit solvers for all problems (to guarantee successful simulations irrespective of stiffness), these are typically less performant (as compared to the explicit solvers) on equations that are non-stiff. An important property of implicit solvers is that they require the *computation of a [Jacobian](https://en.wikipedia.org/wiki/Jacobian_matrix_and_determinant)* as part of their routine. This means that the various options for efficient Jacobian computation [described later in this tutorial](@ref ode_simulation_performance_jacobian) are only relevant to implicit solvers. @@ -68,12 +68,12 @@ ps = [:p => 1.0, :d => 0.2] oprob = ODEProblem(bd_model, u0, tspan, ps) solve(oprob, Tsit5()) nothing # hide -``` +``` If no solver argument is provided to `solve`, one is automatically selected: ```@example ode_simulation_performance_2 solve(oprob) nothing # hide -``` +``` While the default choice is typically enough for most single simulations, if performance is important, it can be worthwhile exploring the available solvers to find one that is especially suited for the given problem. A complete list of possible ODE solvers, with advice on optimal selection, can be found [here](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/). This section will give some general advice. The most important part of solver selection is to select one appropriate for [the problem's stiffness](@ref ode_simulation_performance_stiffness). Generally, the `Tsit5` solver is good for non-stiff problems, and `Rodas5P` for stiff problems. For large stiff problems (with many species), `FBDF` can be a good choice. We can illustrate the impact of these choices by simulating our birth-death process using the `Tsit5`, `Vern7` (an explicit solver yielding [low error in the solution](@ref ode_simulation_performance_error)), `Rodas5P`, and `FBDF` solvers (benchmarking their respective performance using [BenchmarkTools.jl](https://github.com/JuliaCI/BenchmarkTools.jl)): @@ -93,7 +93,7 @@ Numerical ODE simulations [approximate an ODEs' continuous solutions as discrete @btime solve(oprob, Tsit5(); abstol=1e-12, reltol=1e-12) ``` It should be noted, however, that the result of the second simulation is a lot more accurate. Thus, ODE solver performance cannot be determined solely from run time, but is a composite of run -time and error. Benchmarks comparing solver performance (by plotting the run time as a function of the error) for various CRN models can be found in the [SciMLBenchmarks repository](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/Bio/BCR/). +time and error. Benchmarks comparing solver performance (by plotting the run time as a function of the error) for various CRN models can be found in the [SciMLBenchmarks repository](https://docs.sciml.ai/SciMLBenchmarksOutput/stable/Bio/BCR/). Generally, whether solution error is a consideration depends on the application. If you want to compute the trajectory of an expensive space probe as it is sent from Earth, to slingshot Jupiter, and then reach Pluto a few years later, ensuring a minimal error will be essential. However, if you want to simulate a simple CRN to determine whether it oscillates for a given parameter set, a small error will not constitute a problem. An important aspect with regard to error is that it affects the selection of the optimal solver. E.g. if tolerance is low (generating larger errors) the `Rosenbrock23` method performs well for small, stiff, problems (again, more details can be found in [OrdinaryDiffEq's documentation](https://docs.sciml.ai/DiffEqDocs/stable/solvers/ode_solve/)). @@ -136,14 +136,14 @@ using LinearSolve solve(oprob, Rodas5P(linsolve = KLUFactorization())) nothing # hide ``` -A full list of potential linear solvers can be found [here](https://docs.sciml.ai/LinearSolve/dev/solvers/solvers/#Full-List-of-Methods). Typically, the default choice performs well. +A full list of potential linear solvers can be found [here](https://docs.sciml.ai/LinearSolve/dev/solvers/solvers/#Full-List-of-Methods). Typically, the default choice performs well. A unique approach to the linear solvers is to use a *matrix-free Newton-Krylov method*. These do not actually compute the Jacobian, but rather *the effect of multiplying it with a vector*. They are typically advantageous for large systems (with large Jacobians), and can be designated using the `KrylovJL_GMRES` linear solver: ```@example ode_simulation_performance_3 solve(oprob, Rodas5P(linsolve = KrylovJL_GMRES())) nothing # hide ``` -Since these methods do not depend on a Jacobian, certain Jacobian options (such as [computing it symbolically](@ref ode_simulation_performance_symbolic_jacobian)) are irrelevant to them. +Since these methods do not depend on a Jacobian, certain Jacobian options (such as [computing it symbolically](@ref ode_simulation_performance_symbolic_jacobian)) are irrelevant to them. ### [Designating preconditioners for Jacobian-free linear solvers](@id ode_simulation_performance_preconditioners) When an implicit method solves a linear equation through an (iterative) matrix-free Newton-Krylov method, the rate of convergence depends on the numerical properties of the matrix defining the linear system. To speed up convergence, a [*preconditioner*](https://en.wikipedia.org/wiki/Preconditioner) can be applied to both sides of the linear equation, attempting to create an equation that converges faster. Preconditioners are only relevant when using matrix-free Newton-Krylov methods. @@ -188,11 +188,11 @@ oprob = ODEProblem(rs, u0, (0.0, 10.0), ps; remove_conserved = true) sol = solve(oprob) nothing # hide ``` -Conservation law elimination is not expected to ever impact performance negatively; it simply results in a (possibly) lower-dimensional system of ODEs to solve. However, eliminating conserved species may have minimal performance benefits; it is model-dependent whether elimination results in faster ODE solving times and/or increased solution accuracy. +Conservation law elimination is not expected to ever impact performance negatively; it simply results in a (possibly) lower-dimensional system of ODEs to solve. However, eliminating conserved species may have minimal performance benefits; it is model-dependent whether elimination results in faster ODE solving times and/or increased solution accuracy. ## [Parallelisation on CPUs and GPUs](@id ode_simulation_performance_parallelisation) -Whenever an ODE is simulated a large number of times (e.g. when investigating its behaviour for different parameter values), the best way to improve performance is to [parallelise the simulation over multiple processing units](https://en.wikipedia.org/wiki/Parallel_computing). Indeed, an advantage of the Julia programming language is that it was designed after the advent of parallel computing, making it well-suited for this task. Roughly, parallelisation can be divided into parallelisation on [CPUs](https://en.wikipedia.org/wiki/Central_processing_unit) and on [GPUs](https://en.wikipedia.org/wiki/General-purpose_computing_on_graphics_processing_units). CPU parallelisation is most straightforward, while GPU parallelisation requires specialised ODE solvers (which Catalyst have access to). +Whenever an ODE is simulated a large number of times (e.g. when investigating its behaviour for different parameter values), the best way to improve performance is to [parallelise the simulation over multiple processing units](https://en.wikipedia.org/wiki/Parallel_computing). Indeed, an advantage of the Julia programming language is that it was designed after the advent of parallel computing, making it well-suited for this task. Roughly, parallelisation can be divided into parallelisation on [CPUs](https://en.wikipedia.org/wiki/Central_processing_unit) and on [GPUs](https://en.wikipedia.org/wiki/General-purpose_computing_on_graphics_processing_units). CPU parallelisation is most straightforward, while GPU parallelisation requires specialised ODE solvers (which Catalyst have access to). Both CPU and GPU parallelisation require first building an `EnsembleProblem` (which defines the simulations you wish to perform) and then supplying this with the correct parallelisation options. `EnsembleProblem`s have [previously been introduced in Catalyst's documentation](@ref ensemble_simulations) (but in the context of convenient bundling of similar simulations, rather than to improve performance), with a more throughout description being found in [OrdinaryDiffEq's documentation](https://docs.sciml.ai/DiffEqDocs/stable/features/ensemble/#ensemble). Finally, a general documentation of parallel computing in Julia is available [here](https://docs.julialang.org/en/v1/manual/parallel-computing/). @@ -232,7 +232,7 @@ and output the `ODEProblem` simulated in the i'th simulation. Let us assume that we wish to simulate our model 100 times, for $kP = 0.01, 0.02, ..., 0.99, 1.0$. We define our `prob_func` using [`remake`](@ref simulation_structure_interfacing_problems_remake): ```@example ode_simulation_performance_4 -function prob_func(prob, i, repeat) +function prob_func(prob, i, repeat) return remake(prob; p = [:kP => 0.01*i]) end nothing # hide @@ -315,7 +315,7 @@ nothing # hide ``` When we declare our `prob_func` and `EnsembleProblem` we need to ensure that the updated `ODEProblem` uses `Float32`: ```@example ode_simulation_performance_5 -function prob_func(prob, i, repeat) +function prob_func(prob, i, repeat) return remake(prob; p = [:kP => 0.0001f0*i]) end eprob = EnsembleProblem(oprob; prob_func = prob_func) diff --git a/docs/src/model_simulation/simulation_introduction.md b/docs/src/model_simulation/simulation_introduction.md index f417d83794..1b899be539 100644 --- a/docs/src/model_simulation/simulation_introduction.md +++ b/docs/src/model_simulation/simulation_introduction.md @@ -348,6 +348,7 @@ end ``` This type of model will generate so called *variable rate jumps* (`VariableRateJump`s in JumpProcesses.jl). Such models can be simulated in Catalyst too, but note that now a method for time-stepping the solver must be provided to `solve`. Here ODE solvers should be given as they are used to handle integrating the explicitly time-dependent propensities for problems with variable rates, i.e. the proceeding example can be solved like ```@example simulation_intro_jumps +using OrdinaryDiffEq u0map = [:P => 0] pmap = [:f => 1.0, :A => 2.0, :ϕ => 0.0, :d => 1.0] tspan = (0.0, 24.0) From 865ca2bf76e5007c529cee76420d4ca0c001dc55 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sat, 10 Aug 2024 17:33:27 -0400 Subject: [PATCH 03/19] format --- docs/src/model_simulation/simulation_introduction.md | 2 +- src/reactionsystem_conversions.jl | 5 +---- 2 files changed, 2 insertions(+), 5 deletions(-) diff --git a/docs/src/model_simulation/simulation_introduction.md b/docs/src/model_simulation/simulation_introduction.md index 1b899be539..e2c453708b 100644 --- a/docs/src/model_simulation/simulation_introduction.md +++ b/docs/src/model_simulation/simulation_introduction.md @@ -311,7 +311,7 @@ Next, we process these into a `JumpInputs`: jinput = JumpInputs(two_state_model, u0, tspan, ps) nothing # hide ``` -This is then used as input to a `JumpProblem`. The `JumpProblem` also requires the CRN model and [an aggregator](@ref simulation_intro_jumps_solver_designation) as input. +This is then used as input to a `JumpProblem`: ```@example simulation_intro_jumps jprob = JumpProblem(jinput) nothing # hide diff --git a/src/reactionsystem_conversions.jl b/src/reactionsystem_conversions.jl index 90ef509245..cc143c3c54 100644 --- a/src/reactionsystem_conversions.jl +++ b/src/reactionsystem_conversions.jl @@ -784,7 +784,6 @@ struct JumpInputs{S <: MT.JumpSystem, T <: SciMLBase.AbstractODEProblem} prob::T end - """ jumpinput = JumpInputs(rs::ReactionSystem, u0map, tspan, pmap = DiffEqBase.NullParameters; @@ -819,7 +818,6 @@ plot(sol, idxs = :A) function JumpInputs(rs::ReactionSystem, u0, tspan, p = DiffEqBase.NullParameters(); name = nameof(rs), combinatoric_ratelaws = get_combinatoric_ratelaws(rs), checks = false, kwargs...) - u0map = symmap_to_varmap(rs, u0) pmap = symmap_to_varmap(rs, p) jsys = complete(convert(JumpSystem, rs; name, combinatoric_ratelaws, checks)) @@ -835,7 +833,7 @@ function Base.summary(io::IO, jinputs::JumpInputs) type_color, no_color = SciMLBase.get_colorizers(io) print(io, type_color, nameof(typeof(jinputs)), - no_color, " storing" , "\n", + no_color, " storing", "\n", no_color, " JumpSystem: ", type_color, nameof(jinputs.sys), "\n", no_color, " Problem type: ", type_color, nameof(typeof(jinputs.prob))) end @@ -874,7 +872,6 @@ function JumpProcesses.JumpProblem(jinputs::JumpInputs, JumpProblem(jinputs.sys, jinputs.prob, agg; kwargs...) end - # SteadyStateProblem from AbstractReactionNetwork function DiffEqBase.SteadyStateProblem(rs::ReactionSystem, u0, p = DiffEqBase.NullParameters(), args...; From 8858fd077fb319208347eecd5025177a18afa3f3 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sun, 11 Aug 2024 17:38:25 -0400 Subject: [PATCH 04/19] add accuracy test --- Project.toml | 4 +- docs/Project.toml | 4 +- test/simulation_and_solving/simulate_jumps.jl | 48 ++++++++++++++++++- 3 files changed, 50 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index f57db0db45..958b3164ae 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.13.1" +JumpProcesses = "9.13.2" LaTeXStrings = "1.3.0" Latexify = "0.16.5" MacroTools = "0.5.5" -ModelingToolkit = "9.30.1" +ModelingToolkit = "9.32" Parameters = "0.12" Reexport = "0.2, 1.0" Requires = "1.0" diff --git a/docs/Project.toml b/docs/Project.toml index 7acc721352..3adf7702cd 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.1" +JumpProcesses = "9.13.2" Latexify = "0.16.5" LinearSolve = "2.30" -ModelingToolkit = "9.30.1" +ModelingToolkit = "9.32" NonlinearSolve = "3.12" Optim = "1.9" Optimization = "3.25" diff --git a/test/simulation_and_solving/simulate_jumps.jl b/test/simulation_and_solving/simulate_jumps.jl index b03ba76cd8..e457867514 100644 --- a/test/simulation_and_solving/simulate_jumps.jl +++ b/test/simulation_and_solving/simulate_jumps.jl @@ -163,8 +163,6 @@ let end end -### Other Tests ### - # Tests simulating a network without parameters. let no_param_network = @reaction_network begin @@ -176,3 +174,49 @@ let sol = solve(jprob, SSAStepper()) @test mean(sol[:X1]) > mean(sol[:X2]) end + +# test accuracy of a mixed system that includes variable rate jumps +# this should also help with indirectly testing dep graphs are setup ok +let + rn = @reaction_network gene_model begin + α*(1 + sin(t)), D --> D + P + μ*(1 + cos(t)), P --> ∅ + k₊, D + P --> D⁻ + k₋, D⁻ --> D + P + α, D2 --> D2 + P2 + μ, P2 --> P3 + end + u0map = [rn.D => 1.0, rn.P => 0.0, rn.D⁻ => 0.0, rn.D2 => 1.0, rn.P2 => 0.0, + rn.P3 => 0.0] + pmap = [rn.α => 10.0, rn.μ => 1.0, rn.k₊ => 1.0, rn.k₋ => 2.0] + tspan = (0.0, 25.0) + jinput = JumpInputs(rn, u0map, tspan, pmap) + + # the direct method needs no dep graphs so is good as a baseline for comparison + jprobdm = JumpProblem(jinput, Direct(); save_positions = (false, false), rng) + jprobsd = JumpProblem(jinput, SortingDirect(); save_positions = (false, false), rng) + @test issetequal(jprobsd.discrete_jump_aggregation.dep_gr, [[1,2],[2]]) + jprobrssa = JumpProblem(jinput, RSSA(); save_positions = (false, false), rng) + @test issetequal(jprobrssa.discrete_jump_aggregation.vartojumps_map, [[],[],[],[1],[2],[]]) + @test issetequal(jprobrssa.discrete_jump_aggregation.jumptovars_map, [[5],[5,6]]) + N = 1000 # number of simulations to run + function getmean(N, prob) + m1 = 0.0 + m2 = 0.0 + for _ in 1:N + jsol = solve(prob, Tsit5(); saveat = tspan[2]) + m1 += jsol(prob.prob.tspan[2]; idxs = :P) + m2 += jsol(prob.prob.tspan[2]; idxs = :P3) + end + m1 /= N + m2 /= N + return m1, m2 + end + means1 = zeros(2) + means2 = zeros(2) + for (i,prob) in enumerate((jprobdm, jprobsd)) # skip rssa due JumpProcesses #439 bug + means1[i],means2[i] = getmean(N, prob) + end + @test (means1[1] - means1[2]) < .1 * means1[1] + @test (means2[1] - means2[2]) < .1 * means2[1] +end \ No newline at end of file From c6a84ac5cc17a97202142de3ebd2fb9c2a9f77a8 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sun, 11 Aug 2024 19:57:11 -0400 Subject: [PATCH 05/19] Update simulate_jumps.jl --- test/simulation_and_solving/simulate_jumps.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/test/simulation_and_solving/simulate_jumps.jl b/test/simulation_and_solving/simulate_jumps.jl index e457867514..a5d91d758b 100644 --- a/test/simulation_and_solving/simulate_jumps.jl +++ b/test/simulation_and_solving/simulate_jumps.jl @@ -1,7 +1,7 @@ ### Prepares Tests ### # Fetch packages. -using Catalyst, JumpProcesses, Statistics, Test +using Catalyst, JumpProcesses, OrdinaryDiffEq, Statistics, Test # Sets stable rng number. using StableRNGs @@ -219,4 +219,4 @@ let end @test (means1[1] - means1[2]) < .1 * means1[1] @test (means2[1] - means2[2]) < .1 * means2[1] -end \ No newline at end of file +end From 36c10b8661f44054571a7a70361ed7fd1315189f Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Sun, 11 Aug 2024 21:09:00 -0400 Subject: [PATCH 06/19] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 958b3164ae..46da48167e 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Catalyst" uuid = "479239e8-5488-4da2-87a7-35f2df7eef83" -version = "14.2" +version = "14.3" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" From 739051b14f09d0a34b7334e52ddf2c0d07a26c4b Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 13 Aug 2024 14:30:24 -0400 Subject: [PATCH 07/19] cleanup metadata test and fix error mesg bug --- src/reaction.jl | 21 ++++++++++----------- 1 file changed, 10 insertions(+), 11 deletions(-) diff --git a/src/reaction.jl b/src/reaction.jl index 165bbeff37..8f85d77343 100644 --- a/src/reaction.jl +++ b/src/reaction.jl @@ -351,7 +351,7 @@ MT.is_alg_equation(rx::Reaction) = false """ get_symbolics(set, rx::Reaction) -Returns all symbolic variables that are part of a reaction. This includes all variables +Returns all symbolic variables that are part of a reaction. This includes all variables encountered in: - Rates. - Among substrates and products. @@ -365,7 +365,7 @@ end """ get_variables!(set, rx::Reaction) -Adds all symbolic variables that are part of a reaction to set. This includes all variables +Adds all symbolic variables that are part of a reaction to set. This includes all variables encountered in: - Rates. - Among substrates and products. @@ -412,7 +412,7 @@ function MT.modified_unknowns!(munknowns, rx::Reaction, sts::AbstractVector) munknowns end -### `Reaction`-specific Functions ### +### `Reaction`-specific Functions ### """ isbcbalanced(rx::Reaction) @@ -495,12 +495,11 @@ getmetadata(reaction, :description) ``` """ function getmetadata(reaction::Reaction, md_key::Symbol) - if !hasmetadata(reaction, md_key) - error("The reaction does not have a metadata field $md_key. It does have the following metadata fields: $(keys(getmetadata_dict(reaction))).") - end metadata = getmetadata_dict(reaction) - return metadata[findfirst(isequal(md_key, entry[1]) - for entry in getmetadata_dict(reaction))][2] + idx = findfirst(isequal(md_key, entry[1]) for entry in metadata) + (idx === nothing) && + error("The reaction does not have a metadata field $md_key. It does have the following metadata fields: $(first.(values(metadata))).") + return metadata[idx][2] end ### Implemented Reaction Metadata ### @@ -622,10 +621,10 @@ getmisc(reaction) Notes: - The `misc` field can contain any valid Julia structure. This mean that Catalyst cannot check it -for symbolic variables that are added here. This means that symbolic variables (e.g. parameters of -species) that are stored here are not accessible to Catalyst. This can cause troubles when e.g. +for symbolic variables that are added here. This means that symbolic variables (e.g. parameters of +species) that are stored here are not accessible to Catalyst. This can cause troubles when e.g. creating a `ReactionSystem` programmatically (in which case any symbolic variables stored in the -`misc` metadata field should also be explicitly provided to the `ReactionSystem` constructor). +`misc` metadata field should also be explicitly provided to the `ReactionSystem` constructor). """ function getmisc(reaction::Reaction) From ad7c3859106b53dba9eb74e09b5b7d20506b0edc Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 13 Aug 2024 14:36:59 -0400 Subject: [PATCH 08/19] fix math descript typos --- docs/src/introduction_to_catalyst/math_models_intro.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/introduction_to_catalyst/math_models_intro.md b/docs/src/introduction_to_catalyst/math_models_intro.md index 1357178f56..0a2647a799 100644 --- a/docs/src/introduction_to_catalyst/math_models_intro.md +++ b/docs/src/introduction_to_catalyst/math_models_intro.md @@ -28,11 +28,11 @@ with $\alpha^k = (\alpha_1^k,\dots,\alpha_M^k)$ its substrate stoichiometry vect As explained in [the Catalyst introduction](@ref introduction_to_catalyst), for a mass action reaction where the preceding reaction has a fixed rate constant, $k$, this function would be the rate law ```math -a_k(\mathbf{X}(t)) = k \prod_{m=1}^M \frac{(X_m(t))^{\sigma_m^k}}{\sigma_m^k!}, +a_k(\mathbf{X}(t)) = k \prod_{m=1}^M \frac{(X_m(t))^{\alpha_m^k}}{\alpha_m^k!}, ``` for RRE ODE and CLE SDE models, and the propensity function ```math -a_k(\mathbf{X}(t)) = k \prod_{m=1}^M \frac{X_m(t) (X_m(t)-1) \dots (X_m(t)-\sigma_m^k+1)}{\sigma_m^k!}, +a_k(\mathbf{X}(t)) = k \prod_{m=1}^M \frac{X_m(t) (X_m(t)-1) \dots (X_m(t)-\alpha_m^k+1)}{\alpha_m^k!}, ``` for stochastic chemical kinetics jump process models. @@ -41,7 +41,7 @@ For the reaction $2A + B \overset{k}{\to} 3 C$ we would have ```math \mathbf{X}(t) = (A(t), B(t), C(t)) ``` -with $\sigma_1 = 2$, $\sigma_2 = 1$, $\sigma_3 = 0$, $\beta_1 = 0$, $\beta_2 = +with $\alpha_1 = 2$, $\alpha_2 = 1$, $\alpha_3 = 0$, $\beta_1 = 0$, $\beta_2 = 0$, $\beta_3 = 3$, $\nu_1 = -2$, $\nu_2 = -1$, and $\nu_3 = 3$. For an ODE/SDE model we would have the rate law ```math From b0cf75ed998819a59de01c3ae941d2e374b8021f Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Tue, 13 Aug 2024 15:43:26 -0400 Subject: [PATCH 09/19] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 46da48167e..5d9e1ae58d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Catalyst" uuid = "479239e8-5488-4da2-87a7-35f2df7eef83" -version = "14.3" +version = "14.3.1" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa" From e693e7573465b5a3b30a5998d3384d034fcd40ae Mon Sep 17 00:00:00 2001 From: Torkel Date: Wed, 14 Aug 2024 10:39:07 -0400 Subject: [PATCH 10/19] init --- .../homotopy_continuation_extension.jl | 8 ++++---- src/Catalyst.jl | 2 +- src/chemistry_functionality.jl | 2 +- src/dsl.jl | 2 +- src/latexify_recipes.jl | 2 +- src/registered_functions.jl | 2 +- test/dsl/dsl_options.jl | 4 ++-- test/miscellaneous_tests/compound_macro.jl | 10 +++++----- 8 files changed, 16 insertions(+), 16 deletions(-) diff --git a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl index 916a124423..b0093e5a03 100644 --- a/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl +++ b/ext/CatalystHomotopyContinuationExtension/homotopy_continuation_extension.jl @@ -70,10 +70,10 @@ end function ___make_int_exps(expr) !iscall(expr) && return expr if (operation(expr) == ^) - if isinteger(arguments(expr)[2]) - return arguments(expr)[1]^Int64(arguments(expr)[2]) + if isinteger(sorted_arguments(expr)[2]) + return sorted_arguments(expr)[1]^Int64(sorted_arguments(expr)[2]) else - error("An non integer ($(arguments(expr)[2])) was found as a variable exponent. Non-integer exponents are not supported for homotopy continuation based steady state finding.") + error("An non integer ($(sorted_arguments(expr)[2])) was found as a variable exponent. Non-integer exponents are not supported for homotopy continuation based steady state finding.") end end end @@ -83,7 +83,7 @@ function remove_denominators(expr) s_expr = simplify_fractions(expr) !iscall(expr) && return expr if operation(s_expr) == / - return remove_denominators(arguments(s_expr)[1]) + return remove_denominators(sorted_arguments(s_expr)[1]) end if operation(s_expr) == + return sum(remove_denominators(arg) for arg in arguments(s_expr)) diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 341132e0f0..242323a0e2 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -23,7 +23,7 @@ using RuntimeGeneratedFunctions RuntimeGeneratedFunctions.init(@__MODULE__) import Symbolics: BasicSymbolic -using Symbolics: iscall +using Symbolics: iscall, sorted_arguments using ModelingToolkit: Symbolic, value, get_unknowns, get_ps, get_iv, get_systems, get_eqs, get_defaults, toparam, get_var_to_name, get_observed, getvar diff --git a/src/chemistry_functionality.jl b/src/chemistry_functionality.jl index ff3b663fe6..1dd993ff16 100644 --- a/src/chemistry_functionality.jl +++ b/src/chemistry_functionality.jl @@ -106,7 +106,7 @@ function make_compound(expr) # Expression which when evaluated gives a vector with all the ivs of the components. ivs_get_expr = :(unique(reduce( - vcat, [arguments(ModelingToolkit.unwrap(comp)) + vcat, [sorted_arguments(ModelingToolkit.unwrap(comp)) for comp in $components]))) # Creates the found expressions that will create the compound species. diff --git a/src/dsl.jl b/src/dsl.jl index 7f5e2d6125..c7e020230e 100644 --- a/src/dsl.jl +++ b/src/dsl.jl @@ -776,7 +776,7 @@ function read_observed_options(options, species_n_vars_declared, ivs_sorted) dep_var_expr = :(filter(!MT.isparameter, Symbolics.get_variables($(obs_eq.args[3])))) ivs_get_expr = :(unique(reduce( - vcat, [arguments(MT.unwrap(dep)) + vcat, [sorted_arguments(MT.unwrap(dep)) for dep in $dep_var_expr]))) ivs_get_expr_sorted = :(sort($(ivs_get_expr); by = iv -> findfirst(MT.getname(iv) == ivs for ivs in $ivs_sorted))) diff --git a/src/latexify_recipes.jl b/src/latexify_recipes.jl index d1921056d7..a39ed97f22 100644 --- a/src/latexify_recipes.jl +++ b/src/latexify_recipes.jl @@ -37,7 +37,7 @@ function Latexify.infer_output(env, rs::ReactionSystem, args...) end function processsym(s) - args = Symbolics.sorted_arguments(s) + args = sorted_arguments(s) name = MT.getname(s) if length(args) <= 1 var = value(Symbolics.variable(name)) diff --git a/src/registered_functions.jl b/src/registered_functions.jl index 2c6bbebad4..1109e32100 100644 --- a/src/registered_functions.jl +++ b/src/registered_functions.jl @@ -140,7 +140,7 @@ end # it in its expanded form. If not, returns the input expression. function expand_catalyst_function(expr) is_catalyst_function(expr) || (return expr) - args = arguments(expr) + args = sorted_arguments(expr) if operation(expr) == Catalyst.mm return args[2] * args[1] / (args[1] + args[3]) elseif operation(expr) == Catalyst.mmr diff --git a/test/dsl/dsl_options.jl b/test/dsl/dsl_options.jl index 791bdca75f..2aecc8fffb 100644 --- a/test/dsl/dsl_options.jl +++ b/test/dsl/dsl_options.jl @@ -575,8 +575,8 @@ let end end V,W = getfield.(observed(rn), :lhs) - @test isequal(arguments(ModelingToolkit.unwrap(V)), Any[Catalyst.get_iv(rn), Catalyst.get_sivs(rn)[1], Catalyst.get_sivs(rn)[2]]) - @test isequal(arguments(ModelingToolkit.unwrap(W)), Any[Catalyst.get_iv(rn), Catalyst.get_sivs(rn)[2]]) + @test isequal(Symbolics.sorted_arguments(ModelingToolkit.unwrap(V)), Any[Catalyst.get_iv(rn), Catalyst.get_sivs(rn)[1], Catalyst.get_sivs(rn)[2]]) + @test isequal(Symbolics.sorted_arguments(ModelingToolkit.unwrap(W)), Any[Catalyst.get_iv(rn), Catalyst.get_sivs(rn)[2]]) end # Checks that metadata is written properly. diff --git a/test/miscellaneous_tests/compound_macro.jl b/test/miscellaneous_tests/compound_macro.jl index a65d809973..485b6795e8 100644 --- a/test/miscellaneous_tests/compound_macro.jl +++ b/test/miscellaneous_tests/compound_macro.jl @@ -121,11 +121,11 @@ let @compound SO2(t,x,y) ~ S + 2O # Checks they have the correct independent variables. - @test issetequal(arguments(ModelingToolkit.unwrap(CO2)), [t]) - @test issetequal(arguments(ModelingToolkit.unwrap(NH4)), [x]) - @test issetequal(arguments(ModelingToolkit.unwrap(H2O)), [t, x]) - @test issetequal(arguments(ModelingToolkit.unwrap(PH4)), [t, x]) - @test issetequal(arguments(ModelingToolkit.unwrap(SO2)), [t, x, y]) + @test issetequal(Symbolics.sorted_arguments(ModelingToolkit.unwrap(CO2)), [t]) + @test issetequal(Symbolics.sorted_arguments(ModelingToolkit.unwrap(NH4)), [x]) + @test issetequal(Symbolics.sorted_arguments(ModelingToolkit.unwrap(H2O)), [t, x]) + @test issetequal(Symbolics.sorted_arguments(ModelingToolkit.unwrap(PH4)), [t, x]) + @test issetequal(Symbolics.sorted_arguments(ModelingToolkit.unwrap(SO2)), [t, x, y]) end ### Interpolation Tests ### From f2657a162d54839e023999858d44a673d440ae4e Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 16 Aug 2024 11:28:17 -0400 Subject: [PATCH 11/19] metadata as dispatch from SymbolicUtils --- Project.toml | 2 + src/Catalyst.jl | 1 + src/reaction.jl | 36 +++++++++++++--- test/dsl/dsl_advanced_model_construction.jl | 28 ++++++------ test/reactionsystem_core/reaction.jl | 48 +++++++++++---------- 5 files changed, 72 insertions(+), 43 deletions(-) diff --git a/Project.toml b/Project.toml index 5d9e1ae58d..66cddb75ea 100644 --- a/Project.toml +++ b/Project.toml @@ -23,6 +23,7 @@ RuntimeGeneratedFunctions = "7e49a35a-f44a-4d26-94aa-eba1b4ca6b47" SciMLBase = "0bca4576-84f4-4d90-8ffe-ffa030f20462" Setfield = "efcf1570-3423-57d1-acb7-fd33fddbac46" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +SymbolicUtils = "d1185830-fcd6-423d-90d6-eec64667417b" Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7" Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d" @@ -65,6 +66,7 @@ SciMLBase = "2.46" Setfield = "1" StructuralIdentifiability = "0.5.8" Symbolics = "5.30.1" +SymbolicUtils = "2.1.2" Unitful = "1.12.4" julia = "1.10" diff --git a/src/Catalyst.jl b/src/Catalyst.jl index 242323a0e2..2c5c6003a2 100644 --- a/src/Catalyst.jl +++ b/src/Catalyst.jl @@ -45,6 +45,7 @@ import DataStructures: OrderedDict, OrderedSet import Parameters: @with_kw_noshow import Symbolics: occursin, wrap import Symbolics.RewriteHelpers: hasnode, replacenode +import SymbolicUtils: getmetadata, hasmetadata, setmetadata # globals for the modulate function default_time_deriv() diff --git a/src/reaction.jl b/src/reaction.jl index 8f85d77343..f4208f06dd 100644 --- a/src/reaction.jl +++ b/src/reaction.jl @@ -461,9 +461,10 @@ function getmetadata_dict(reaction::Reaction) end """ -hasmetadata(reaction::Reaction, md_key::Symbol) + SymbolicUtils.hasmetadata(reaction::Reaction, md_key::Symbol) -Checks if a `Reaction` have a certain metadata field. If it does, returns `true` (else returns `false`). +Checks if a `Reaction` have a certain metadata field. If it does, returns `true` (else +returns `false`). Arguments: - `reaction`: The reaction for which we wish to check if a specific metadata field exist. @@ -475,14 +476,15 @@ reaction = @reaction k, 0 --> X, [description="Production reaction"] hasmetadata(reaction, :description) ``` """ -function hasmetadata(reaction::Reaction, md_key::Symbol) +function SymbolicUtils.hasmetadata(reaction::Reaction, md_key::Symbol) return any(isequal(md_key, entry[1]) for entry in getmetadata_dict(reaction)) end """ -getmetadata(reaction::Reaction, md_key::Symbol) + SymbolicUtils.getmetadata(reaction::Reaction, md_key::Symbol) -Retrieves a certain metadata value from a `Reaction`. If the metadata does not exist, throws an error. +Retrieves a certain metadata value from a `Reaction`. If the metadata does not exist, throws +an error. Arguments: - `reaction`: The reaction for which we wish to retrieve a specific metadata value. @@ -494,7 +496,7 @@ reaction = @reaction k, 0 --> X, [description="Production reaction"] getmetadata(reaction, :description) ``` """ -function getmetadata(reaction::Reaction, md_key::Symbol) +function SymbolicUtils.getmetadata(reaction::Reaction, md_key::Symbol) metadata = getmetadata_dict(reaction) idx = findfirst(isequal(md_key, entry[1]) for entry in metadata) (idx === nothing) && @@ -502,7 +504,27 @@ function getmetadata(reaction::Reaction, md_key::Symbol) return metadata[idx][2] end -### Implemented Reaction Metadata ### +""" + SymbolicUtils.setmetadata(rx::Reaction, key::Symbol, val) + +Sets the metadata with key `key` to the value `val`, overwriting if already present or +adding if not present. + +Arguments: +- `rx`: The reaction to add the metadata too. +- `key`: `Symbol` representing the metadata's key (i.e. name). +- `val`: value for the metadata. +""" +function SymbolicUtils.setmetadata(rx::Reaction, key::Symbol, val) + mdvec = getmetadata_dict(rx) + idx = findfirst(isequal(key, first(md)) for md in mdvec) + if idx === nothing + push!(mdvec, key => val) + else + mdvec[idx] = key => val + end + nothing +end # Noise scaling. """ diff --git a/test/dsl/dsl_advanced_model_construction.jl b/test/dsl/dsl_advanced_model_construction.jl index 5638c9173c..b988adcb95 100644 --- a/test/dsl/dsl_advanced_model_construction.jl +++ b/test/dsl/dsl_advanced_model_construction.jl @@ -184,20 +184,20 @@ let # Checks DSL reactions are correct. rxs = reactions(rs) @test isequal([r1, r2, r3], rxs) - @test isequal(Catalyst.getmetadata_dict(r1), Catalyst.getmetadata_dict(rxs[1])) - @test isequal(Catalyst.getmetadata_dict(r2), Catalyst.getmetadata_dict(rxs[2])) - @test isequal(Catalyst.getmetadata_dict(r3), Catalyst.getmetadata_dict(rxs[3])) + @test isequal(getmetadata_dict(r1), getmetadata_dict(rxs[1])) + @test isequal(getmetadata_dict(r2), getmetadata_dict(rxs[2])) + @test isequal(getmetadata_dict(r3), getmetadata_dict(rxs[3])) # Checks that accessor functions works on the DSL. - @test Catalyst.hasmetadata(rxs[1], :noise_scaling) - @test !Catalyst.hasmetadata(rxs[1], :md_1) - @test !Catalyst.hasmetadata(rxs[2], :noise_scaling) - @test Catalyst.hasmetadata(rxs[2], :md_1) - @test !Catalyst.hasmetadata(rxs[3], :noise_scaling) - @test !Catalyst.hasmetadata(rxs[3], :md_1) + @test getmetadata(rxs[1], :noise_scaling) + @test !getmetadata(rxs[1], :md_1) + @test !getmetadata(rxs[2], :noise_scaling) + @test getmetadata(rxs[2], :md_1) + @test !getmetadata(rxs[3], :noise_scaling) + @test !getmetadata(rxs[3], :md_1) - @test isequal(Catalyst.getmetadata(rxs[1], :noise_scaling), η) - @test isequal(Catalyst.getmetadata(rxs[2], :md_1), 1.0) + @test isequal(getmetadata(rxs[1], :noise_scaling), η) + @test isequal(getmetadata(rxs[2], :md_1), 1.0) # Test that metadata works for @reaction macro. rx1 = @reaction k, 2X --> X2, [noise_scaling=$η] @@ -205,9 +205,9 @@ let rx3 = @reaction k, 2X --> X2 @test isequal([rx1, rx2, rx3], rxs) - @test isequal(Catalyst.getmetadata_dict(rx1), Catalyst.getmetadata_dict(rxs[1])) - @test isequal(Catalyst.getmetadata_dict(rx2), Catalyst.getmetadata_dict(rxs[2])) - @test isequal(Catalyst.getmetadata_dict(rx3), Catalyst.getmetadata_dict(rxs[3])) + @test isequal(getmetadata_dict(rx1), getmetadata_dict(rxs[1])) + @test isequal(getmetadata_dict(rx2), getmetadata_dict(rxs[2])) + @test isequal(getmetadata_dict(rx3), getmetadata_dict(rxs[3])) end # Checks that repeated metadata throws errors. diff --git a/test/reactionsystem_core/reaction.jl b/test/reactionsystem_core/reaction.jl index 36d5d946d0..ac5bfc58c6 100644 --- a/test/reactionsystem_core/reaction.jl +++ b/test/reactionsystem_core/reaction.jl @@ -132,11 +132,15 @@ let metadata = [:noise_scaling => 0.0] r = Reaction(k, [X], [X2], [2], [1]; metadata=metadata) - @test Catalyst.getmetadata_dict(r) == [:noise_scaling => 0.0] - @test Catalyst.hasmetadata(r, :noise_scaling) - @test !Catalyst.hasmetadata(r, :nonexisting_metadata) - @test Catalyst.getmetadata(r, :noise_scaling) == 0.0 - @test_throws Exception Catalyst.getmetadata(r, :misc) + @test getmetadata_dict(r) == [:noise_scaling => 0.0] + @test getmetadata(r, :noise_scaling) + @test !getmetadata(r, :nonexisting_metadata) + @test getmetadata(r, :noise_scaling) == 0.0 + @test_throws Exception getmetadata(r, :misc) + setmetadata(r, :test_metadata, 1234) + @test getmetadata(r, :test_metadata) == 1234 + setmetadata(r, :test_metadata, 1111) + @test getmetadata(r, :test_metadata) == 1111 metadata_repeated = [:noise_scaling => 0.0, :noise_scaling => 1.0, :metadata_entry => "unused"] @test_throws Exception Reaction(k, [X], [X2], [2], [1]; metadata=metadata_repeated) @@ -152,8 +156,8 @@ let r2 = Reaction(k, [X], [X2], [2], [1]; metadata=metadata) @test isequal(r1, r2) - @test Catalyst.getmetadata_dict(r1) == Pair{Symbol,Any}[] - @test !Catalyst.hasmetadata(r1, :md) + @test getmetadata_dict(r1) == Pair{Symbol,Any}[] + @test !getmetadata(r1, :md) end # Tests creation. @@ -172,21 +176,21 @@ let push!(metadata, :md_6 => (0.1, 2.0)) r = Reaction(k, [X], [X2], [2], [1]; metadata=metadata) - @test Catalyst.getmetadata_dict(r) isa Vector{Pair{Symbol,Any}} - @test Catalyst.hasmetadata(r, :md_1) - @test Catalyst.hasmetadata(r, :md_2) - @test Catalyst.hasmetadata(r, :md_3) - @test Catalyst.hasmetadata(r, :md_4) - @test Catalyst.hasmetadata(r, :md_5) - @test Catalyst.hasmetadata(r, :md_6) - @test !Catalyst.hasmetadata(r, :md_8) - - @test isequal(Catalyst.getmetadata(r, :md_1), 1.0) - @test isequal(Catalyst.getmetadata(r, :md_2), false) - @test isequal(Catalyst.getmetadata(r, :md_3), "Hello world") - @test isequal(Catalyst.getmetadata(r, :md_4), :sym) - @test isequal(Catalyst.getmetadata(r, :md_5), X + X2^k -1) - @test isequal(Catalyst.getmetadata(r, :md_6), (0.1, 2.0)) + @test getmetadata_dict(r) isa Vector{Pair{Symbol,Any}} + @test getmetadata(r, :md_1) + @test getmetadata(r, :md_2) + @test getmetadata(r, :md_3) + @test getmetadata(r, :md_4) + @test getmetadata(r, :md_5) + @test getmetadata(r, :md_6) + @test !getmetadata(r, :md_8) + + @test isequal(getmetadata(r, :md_1), 1.0) + @test isequal(getmetadata(r, :md_2), false) + @test isequal(getmetadata(r, :md_3), "Hello world") + @test isequal(getmetadata(r, :md_4), :sym) + @test isequal(getmetadata(r, :md_5), X + X2^k -1) + @test isequal(getmetadata(r, :md_6), (0.1, 2.0)) end # Tests the noise scaling metadata. From 23eab99779dbc371b4f7ee1e77f1a454c2b4f892 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 16 Aug 2024 11:31:58 -0400 Subject: [PATCH 12/19] update comment --- src/reaction.jl | 2 ++ 1 file changed, 2 insertions(+) diff --git a/src/reaction.jl b/src/reaction.jl index f4208f06dd..e6202a73ac 100644 --- a/src/reaction.jl +++ b/src/reaction.jl @@ -526,6 +526,8 @@ function SymbolicUtils.setmetadata(rx::Reaction, key::Symbol, val) nothing end +### Catalyst Defined Reaction Metadata ### + # Noise scaling. """ hasnoisescaling(reaction::Reaction) From c501104f77f08a88c986689a4827c1018c448a30 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 16 Aug 2024 11:34:52 -0400 Subject: [PATCH 13/19] fix --- test/dsl/dsl_advanced_model_construction.jl | 12 ++++++------ test/reactionsystem_core/reaction.jl | 6 +++--- 2 files changed, 9 insertions(+), 9 deletions(-) diff --git a/test/dsl/dsl_advanced_model_construction.jl b/test/dsl/dsl_advanced_model_construction.jl index b988adcb95..e8e73b51fe 100644 --- a/test/dsl/dsl_advanced_model_construction.jl +++ b/test/dsl/dsl_advanced_model_construction.jl @@ -184,9 +184,9 @@ let # Checks DSL reactions are correct. rxs = reactions(rs) @test isequal([r1, r2, r3], rxs) - @test isequal(getmetadata_dict(r1), getmetadata_dict(rxs[1])) - @test isequal(getmetadata_dict(r2), getmetadata_dict(rxs[2])) - @test isequal(getmetadata_dict(r3), getmetadata_dict(rxs[3])) + @test isequal(Catalyst.getmetadata_dict(r1), Catalyst.getmetadata_dict(rxs[1])) + @test isequal(Catalyst.getmetadata_dict(r2), Catalyst.getmetadata_dict(rxs[2])) + @test isequal(Catalyst.getmetadata_dict(r3), Catalyst.getmetadata_dict(rxs[3])) # Checks that accessor functions works on the DSL. @test getmetadata(rxs[1], :noise_scaling) @@ -205,9 +205,9 @@ let rx3 = @reaction k, 2X --> X2 @test isequal([rx1, rx2, rx3], rxs) - @test isequal(getmetadata_dict(rx1), getmetadata_dict(rxs[1])) - @test isequal(getmetadata_dict(rx2), getmetadata_dict(rxs[2])) - @test isequal(getmetadata_dict(rx3), getmetadata_dict(rxs[3])) + @test isequal(Catalyst.getmetadata_dict(rx1), Catalyst.getmetadata_dict(rxs[1])) + @test isequal(Catalyst.getmetadata_dict(rx2), Catalyst.getmetadata_dict(rxs[2])) + @test isequal(Catalyst.getmetadata_dict(rx3), Catalyst.getmetadata_dict(rxs[3])) end # Checks that repeated metadata throws errors. diff --git a/test/reactionsystem_core/reaction.jl b/test/reactionsystem_core/reaction.jl index ac5bfc58c6..171be69a05 100644 --- a/test/reactionsystem_core/reaction.jl +++ b/test/reactionsystem_core/reaction.jl @@ -132,7 +132,7 @@ let metadata = [:noise_scaling => 0.0] r = Reaction(k, [X], [X2], [2], [1]; metadata=metadata) - @test getmetadata_dict(r) == [:noise_scaling => 0.0] + @test Catalyst.getmetadata_dict(r) == [:noise_scaling => 0.0] @test getmetadata(r, :noise_scaling) @test !getmetadata(r, :nonexisting_metadata) @test getmetadata(r, :noise_scaling) == 0.0 @@ -156,7 +156,7 @@ let r2 = Reaction(k, [X], [X2], [2], [1]; metadata=metadata) @test isequal(r1, r2) - @test getmetadata_dict(r1) == Pair{Symbol,Any}[] + @test Catalyst.getmetadata_dict(r1) == Pair{Symbol,Any}[] @test !getmetadata(r1, :md) end @@ -176,7 +176,7 @@ let push!(metadata, :md_6 => (0.1, 2.0)) r = Reaction(k, [X], [X2], [2], [1]; metadata=metadata) - @test getmetadata_dict(r) isa Vector{Pair{Symbol,Any}} + @test Catalyst.getmetadata_dict(r) isa Vector{Pair{Symbol,Any}} @test getmetadata(r, :md_1) @test getmetadata(r, :md_2) @test getmetadata(r, :md_3) From ccd296e035c36964228ea8792e62094e0737d4d6 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 16 Aug 2024 11:38:19 -0400 Subject: [PATCH 14/19] fix --- test/dsl/dsl_advanced_model_construction.jl | 12 ++++++------ test/reactionsystem_core/reaction.jl | 20 ++++++++++---------- 2 files changed, 16 insertions(+), 16 deletions(-) diff --git a/test/dsl/dsl_advanced_model_construction.jl b/test/dsl/dsl_advanced_model_construction.jl index e8e73b51fe..c4e46bb521 100644 --- a/test/dsl/dsl_advanced_model_construction.jl +++ b/test/dsl/dsl_advanced_model_construction.jl @@ -189,12 +189,12 @@ let @test isequal(Catalyst.getmetadata_dict(r3), Catalyst.getmetadata_dict(rxs[3])) # Checks that accessor functions works on the DSL. - @test getmetadata(rxs[1], :noise_scaling) - @test !getmetadata(rxs[1], :md_1) - @test !getmetadata(rxs[2], :noise_scaling) - @test getmetadata(rxs[2], :md_1) - @test !getmetadata(rxs[3], :noise_scaling) - @test !getmetadata(rxs[3], :md_1) + @test hasmetadata(rxs[1], :noise_scaling) + @test !hasmetadata(rxs[1], :md_1) + @test !hasmetadata(rxs[2], :noise_scaling) + @test hasmetadata(rxs[2], :md_1) + @test !hasmetadata(rxs[3], :noise_scaling) + @test !hasmetadata(rxs[3], :md_1) @test isequal(getmetadata(rxs[1], :noise_scaling), η) @test isequal(getmetadata(rxs[2], :md_1), 1.0) diff --git a/test/reactionsystem_core/reaction.jl b/test/reactionsystem_core/reaction.jl index 171be69a05..3d3b5396f7 100644 --- a/test/reactionsystem_core/reaction.jl +++ b/test/reactionsystem_core/reaction.jl @@ -133,8 +133,8 @@ let r = Reaction(k, [X], [X2], [2], [1]; metadata=metadata) @test Catalyst.getmetadata_dict(r) == [:noise_scaling => 0.0] - @test getmetadata(r, :noise_scaling) - @test !getmetadata(r, :nonexisting_metadata) + @test hasmetadata(r, :noise_scaling) + @test !hasmetadata(r, :nonexisting_metadata) @test getmetadata(r, :noise_scaling) == 0.0 @test_throws Exception getmetadata(r, :misc) setmetadata(r, :test_metadata, 1234) @@ -157,7 +157,7 @@ let @test isequal(r1, r2) @test Catalyst.getmetadata_dict(r1) == Pair{Symbol,Any}[] - @test !getmetadata(r1, :md) + @test !hasmetadata(r1, :md) end # Tests creation. @@ -177,13 +177,13 @@ let r = Reaction(k, [X], [X2], [2], [1]; metadata=metadata) @test Catalyst.getmetadata_dict(r) isa Vector{Pair{Symbol,Any}} - @test getmetadata(r, :md_1) - @test getmetadata(r, :md_2) - @test getmetadata(r, :md_3) - @test getmetadata(r, :md_4) - @test getmetadata(r, :md_5) - @test getmetadata(r, :md_6) - @test !getmetadata(r, :md_8) + @test hasmetadata(r, :md_1) + @test hasmetadata(r, :md_2) + @test hasmetadata(r, :md_3) + @test hasmetadata(r, :md_4) + @test hasmetadata(r, :md_5) + @test hasmetadata(r, :md_6) + @test !hasmetadata(r, :md_8) @test isequal(getmetadata(r, :md_1), 1.0) @test isequal(getmetadata(r, :md_2), false) From 204ce6a08a0f29d4f84f6f7f0916986d760390ce Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 16 Aug 2024 12:21:07 -0400 Subject: [PATCH 15/19] try latest Symbolics and SymbolicUtils --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 66cddb75ea..995ce2f3ea 100644 --- a/Project.toml +++ b/Project.toml @@ -65,8 +65,8 @@ RuntimeGeneratedFunctions = "0.5.12" SciMLBase = "2.46" Setfield = "1" StructuralIdentifiability = "0.5.8" -Symbolics = "5.30.1" -SymbolicUtils = "2.1.2" +Symbolics = "5.30.1, 6" +SymbolicUtils = "2.1.2, 3.1.2" Unitful = "1.12.4" julia = "1.10" From 6a881596e815f33075e35a7664b595ea5ead98ef Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 16 Aug 2024 12:21:53 -0400 Subject: [PATCH 16/19] doc project too --- docs/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index 3adf7702cd..b67138775e 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -73,4 +73,4 @@ StaticArrays = "1.9" SteadyStateDiffEq = "2.2" StochasticDiffEq = "6.65" StructuralIdentifiability = "0.5.8" -Symbolics = "5.30.1" +Symbolics = "6" From 168a6abe66459fc6be592830f0ff7f2f4f7afdb3 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 16 Aug 2024 12:35:08 -0400 Subject: [PATCH 17/19] don't update Symbolics and SymbolicUtils --- docs/Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/Project.toml b/docs/Project.toml index b67138775e..3adf7702cd 100644 --- a/docs/Project.toml +++ b/docs/Project.toml @@ -73,4 +73,4 @@ StaticArrays = "1.9" SteadyStateDiffEq = "2.2" StochasticDiffEq = "6.65" StructuralIdentifiability = "0.5.8" -Symbolics = "6" +Symbolics = "5.30.1" From c2304834931b172722f356845c904bfa3b9fa9ed Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 16 Aug 2024 12:35:31 -0400 Subject: [PATCH 18/19] don't update Symbolics and SymbolicUtils --- Project.toml | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/Project.toml b/Project.toml index 995ce2f3ea..66cddb75ea 100644 --- a/Project.toml +++ b/Project.toml @@ -65,8 +65,8 @@ RuntimeGeneratedFunctions = "0.5.12" SciMLBase = "2.46" Setfield = "1" StructuralIdentifiability = "0.5.8" -Symbolics = "5.30.1, 6" -SymbolicUtils = "2.1.2, 3.1.2" +Symbolics = "5.30.1" +SymbolicUtils = "2.1.2" Unitful = "1.12.4" julia = "1.10" From feb7e5a5a375c42cc738cb59a71477910a4324e2 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Fri, 16 Aug 2024 13:52:40 -0400 Subject: [PATCH 19/19] Update Project.toml --- Project.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/Project.toml b/Project.toml index 66cddb75ea..3f65c5fb4f 100644 --- a/Project.toml +++ b/Project.toml @@ -1,6 +1,6 @@ name = "Catalyst" uuid = "479239e8-5488-4da2-87a7-35f2df7eef83" -version = "14.3.1" +version = "14.3.2" [deps] Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"