Merge branch 'master' into compathelper/new_version/2024-08-12-00-21-…

…56-903-00340686961

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

Project.toml

@@ -1,6 +1,6 @@
 name = "Catalyst"
 uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
-version = "14.2"
+version = "14.3.2"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -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"
 
@@ -52,11 +53,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"
@@ -65,6 +66,7 @@ SciMLBase = "2.46"
 Setfield = "1"
 StructuralIdentifiability = "0.5.8"
 Symbolics = "5.30.1, 6"
+SymbolicUtils = "2.1.2"
 Unitful = "1.12.4"
 julia = "1.10"
 

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)
 ```

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"

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
 ```

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)
 ```
 

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

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/).

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

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,25 +129,26 @@ 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
 @species A(t) B(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]
 
 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. 
+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.

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))

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:

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