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

…54-967-03356142809

HISTORY.md

@@ -2,6 +2,9 @@
 
 ## Catalyst unreleased (master branch)
 
+## Catalyst 14.4
+- Symbolics 6 support.
+
 ## Catalyst 14.3
 - Support for simulating stochastic chemical kinetics models with explicitly
   time-dependent propensities (i.e. where the resulting `JumpSystem` contains

Project.toml

@@ -1,6 +1,6 @@
 name = "Catalyst"
 uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
-version = "14.3.2"
+version = "14.4"
 
 [deps]
 Combinatorics = "861a8166-3701-5b0c-9a16-15d98fcdc6aa"
@@ -48,7 +48,7 @@ Combinatorics = "1.0.2"
 DataStructures = "0.18"
 DiffEqBase = "6.83.0"
 DocStringExtensions = "0.8, 0.9"
-DynamicPolynomials = "0.5"
+DynamicPolynomials = "0.5, 0.6"
 DynamicQuantities = "0.13.2, 1"
 GraphMakie = "0.5"
 Graphs = "1.4"

docs/src/model_creation/constraint_equations.md

@@ -150,11 +150,26 @@ p = [k_on => 100.0, switch_time => 2.0, k_off => 10.0]
 ```
 Simulating our model,
 ```@example ceq3
-@named osys = ReactionSystem(rxs, t, [A, B], [k_on, k_off, switch_time]; discrete_events)
-osys = complete(osys)
+@named rs2 = ReactionSystem(rxs, t, [A, B], [k_on, k_off, switch_time]; discrete_events)
+rs2 = complete(rs2)
 
-oprob = ODEProblem(osys, u0, tspan, p)
+oprob = ODEProblem(rs2, 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 via `tstops` so that
+the ODE solver can step exactly to the specific time of our event. In the
+previous example we just manually set the numeric value of the parameter in the
+`tstops` kwarg to `solve`, however, it can often be convenient to instead get
+the value of the parameter from `oprob` and pass this numeric value. This helps
+ensure consistency between the value passed via `p` and/or symbolic defaults and
+what we pass as a `tstop` to `solve`. We can do this as
+```julia
+switch_time_val = oprob.ps[:switch_time]
+sol = solve(oprob, Tsit5(); tstops = switch_time_val)
+plot(sol)
+```
+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_simulation/examples/periodic_events_simulation.md

@@ -1,5 +1,5 @@
 # [Modelling a periodic event during ODE and jump simulations](@id periodic_event_simulation_example)
-This tutorial will describe how to simulate systems with periodic events in ODE (straightforward) and jump (slightly less, but still fairly, straightforward) simulations (SDEs use identical syntax to ODEs). We will consider a model with a [circadian rhythm](https://en.wikipedia.org/wiki/Circadian_rhythm), where a parameter represents the level of light. While outdoor light varies smoothly, in experimental settings a lamp is often simply turned on/off every 12 hours. Here we will model this toggling of the light using a periodic event that is triggered every 12 hours.
+This tutorial will describe how to simulate systems with periodic events in ODEs and jump simulations (SDEs use identical syntax). We will consider a model with a [circadian rhythm](https://en.wikipedia.org/wiki/Circadian_rhythm), where a parameter represents the level of light. While outdoor light varies smoothly, in experimental settings a lamp is often simply turned on/off every 12 hours. Here we will model this toggling of the light using a periodic event that is triggered every 12 hours.
 
 ## [Modelling a circadian periodic event in an ODE simulation](@id periodic_event_simulation_example_ode)
 We will consider a simple circadian model, consisting of a single protein ($X$), which is phosphorylated ($X \to Xᴾ$) in the presence of light ($l$). Here, the light parameter can either be $0$ (night) or $1$ (day). We can model this using a simple periodic event which switches the value of $l$ every 12 hours (here, `%` is the [remainder operator](https://docs.julialang.org/en/v1/manual/mathematical-operations/#Arithmetic-Operators)).
@@ -17,21 +17,43 @@ We can now simulate this model, observing how a 24-hour cycle is reached
 using OrdinaryDiffEq, Plots
 u0 = [:X => 150.0, :Xᴾ => 50.0]
 ps = [:kₚ => 0.1, :kᵢ => 0.1, :l => 1.0]
-oprob = ODEProblem(circadian_model, u0, (0.0, 100.0), ps)
+tspan = (0.0, 100.0)
+oprob = ODEProblem(circadian_model, u0, tspan, ps)
 sol = solve(oprob)
 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 first creating a [conditional discrete event](@ref constraint_equations_events) which is designed to trigger every 13 time units:
+We can define periodic events in an identical manner for jump simulations. Let's
+reuse our previously defined network, but now simulate it as a stochastic
+chemical kinetics jump process model
+```@example periodic_event_example
+using JumpProcesses
+u0 = [:X => 150, :Xᴾ => 50]  # define u0 as integers now
+jinput = JumpInputs(circadian_model, u0, tspan, ps)
+jprob = JumpProblem(jinput)
+jsol = solve(jprob)
+plot(jsol)
+Catalyst.PNG(plot(jsol; fmt = :png, dpi = 200)) # hide
+```
+
+Sometimes we might want to setup a more complicated condition than simply that
+our event occurs at a fixed time or with a fixed frequency. For example, suppose
+we want to skip the first event at time `t = 12` and then have the event be
+periodic after that point every 12 units of time. We can do so using a more
+general discrete callback as follows
 ```@example periodic_event_example
 circadian_model = @reaction_network begin
     @discrete_events begin
-        (t % 12 == 0) => [l ~ (l + 1)%2]
+        ((t % 12 == 0) & (t > 12)) => [l ~ (l + 1)%2]
     end
     (kₚ*l,kᵢ), X <--> Xᴾ
 end
-
+```
+Here our condition `((t % 12 == 0) & (t > 12))` determines when the event
+occurs, evaluating to `true` when `t` is a multiple of `12` that is also larger
+than `12`. We now finish specifying our model
+```@example periodic_event_example
 using JumpProcesses
 u0 = [:X => 150, :Xᴾ => 50]
 ps = [:kₚ => 0.1, :kᵢ => 0.1, :l => 1.0]
@@ -40,15 +62,20 @@ jinput = JumpInputs(circadian_model, u0, tspan, ps)
 jprob = JumpProblem(jinput)
 nothing # hide
 ```
-Next, if we simulate our model, we note that the events do not seem to be triggered
+Next, if we simulate our model, we note that the events do not seem to be
+triggered
 ```@example periodic_event_example
 sol = solve(jprob)
 plot(sol)
 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:
+The reason is that general 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:tspan[2]
+tstops = range(12.0, tspan[2]; step = 12.0)
 sol = solve(jprob; tstops)
 plot(sol)
 Catalyst.PNG(plot(sol; fmt = :png, dpi = 200)) # hide

src/dsl.jl

@@ -150,7 +150,7 @@ macro reaction_network(name::Expr, ex::Expr)
         MacroTools.striplines(ex); name = :($(esc(name.args[1])))))))
 end
 
-macro reaction_network(ex::Expr)    
+macro reaction_network(ex::Expr)
     ex = MacroTools.striplines(ex)
 
     # no name but equations: @reaction_network begin ... end ...

src/latexify_recipes.jl

@@ -28,12 +28,11 @@ const LATEX_DEFS = CatalystLatexParams()
 end
 
 function Latexify.infer_output(env, rs::ReactionSystem, args...)
-    env in [:arrows, :chem, :chemical, :arrow] && return chemical_arrows
-
-    error("The environment $env is not defined.")
-    latex_function = Latexify.get_latex_function(rs, args...)
-
-    return latex_function
+    if env in (:arrows, :chem, :chemical, :arrow)
+        return chemical_arrows
+    else
+        error("The environment $env is not defined.")
+    end
 end
 
 function processsym(s)

test/reactionsystem_core/events.jl

@@ -349,7 +349,7 @@ let
         @parameters e1=0 e2=0 e3=0
         @discrete_events begin
             [1.0] => [e1 ~ 1]
-            # 1.0 => [e2 ~ 1]
+            1.0 => [e2 ~ 1]
             (X > 1000.0) & (e3==0) => [e3 ~ 1]
         end
         (p,d), 0 <--> X
@@ -363,9 +363,8 @@ let
     sol = solve(jprob, SSAStepper(); seed)
 
     # Checks that all `e` parameters have been updated properly.
-    # Note that periodic discrete events are currently broken for jump processes (and unlikely to be fixed soon due to periodic callbacks using the internals of ODE integrator and Datastructures heap implementations).
     @test sol.ps[:e1] == 1
-    @test_broken sol.ps[:e2] == 1 # (https://github.com/SciML/JumpProcesses.jl/issues/417)
+    @test sol.ps[:e2] == 1 
     @test sol.ps[:e3] == 1
 end
 
@@ -440,8 +439,6 @@ let
     jprob = JumpProblem(rn, dprob, Direct(); rng)
     jprob_events = JumpProblem(rn_dics_events, dprob_events, Direct(); rng)
     sol = solve(jprob, SSAStepper(); seed, callback)
-    @test_broken let # (https://github.com/SciML/JumpProcesses.jl/issues/417)
-        sol_events = solve(jprob_events, SSAStepper(); seed)
-        @test sol == sol_events
-    end
-end
+    sol_events = solve(jprob_events, SSAStepper(); seed)
+    @test_broken sol == sol_events  # seems to be not identical in the sample paths
+end