Merge pull request #940 from isaacsas/fix_parametric_stoich_tutorial

Fix parametric stoich tutorial

Author: isaacsas

docs/make.jl

@@ -41,7 +41,7 @@ makedocs(sitename = "Catalyst.jl",
     clean = true,
     pages = pages,
     pagesonly = true,
-    warnonly = [:missing_docs])
+    warnonly = true)
 
 deploydocs(repo = "github.com/SciML/Catalyst.jl.git";
     push_preview = true)

docs/src/assets/Project.toml

@@ -51,7 +51,7 @@ IncompleteLU = "0.2"
 JumpProcesses = "9.11"
 Latexify = "0.16"
 LinearSolve = "2.30"
-ModelingToolkit = "9.15"
+ModelingToolkit = "9.16.0"
 NonlinearSolve = "3.12"
 Optim = "1.9"
 Optimization = "3.25"
@@ -68,5 +68,5 @@ SpecialFunctions = "2.4"
 StaticArrays = "1.9"
 SteadyStateDiffEq = "2.2"
 StochasticDiffEq = "6.65"
-StructuralIdentifiability = "0.5.7"
-Symbolics = "5.28"
+StructuralIdentifiability = "0.5.8"
+Symbolics = "5.30.1"

docs/src/introduction_to_catalyst/introduction_to_catalyst.md

@@ -182,7 +182,7 @@ Gillespie's `Direct` method, and then solve it to generate one realization of
 the jump process:
 
 ```@example tut1
-# imports the JumpProcesses packages 
+# imports the JumpProcesses packages
 using JumpProcesses
 
 # redefine the initial condition to be integer valued
@@ -361,4 +361,4 @@ and the ODE model
 
 ---
 ## References
-[^1]: [Torkel E. Loman, Yingbo Ma, Vasily Ilin, Shashi Gowda, Niklas Korsbo, Nikhil Yewale, Chris Rackauckas, Samuel A. Isaacson, *Catalyst: Fast and flexible modeling of reaction networks*, PLOS Computational Biology (2023).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011530)
+1. [Torkel E. Loman, Yingbo Ma, Vasily Ilin, Shashi Gowda, Niklas Korsbo, Nikhil Yewale, Chris Rackauckas, Samuel A. Isaacson, *Catalyst: Fast and flexible modeling of reaction networks*, PLOS Computational Biology (2023).](https://journals.plos.org/ploscompbiol/article?id=10.1371/journal.pcbi.1011530)

docs/src/model_creation/examples/smoluchowski_coagulation_equation.md

@@ -147,6 +147,6 @@ plot!(ϕ, sol[3,:] / Nₒ, line = (:dot, 4, :purple), label = "Analytical sol--X
 
 ---
 ## References
-[^1]: [https://en.wikipedia.org/wiki/Smoluchowski\_coagulation\_equation](https://en.wikipedia.org/wiki/Smoluchowski_coagulation_equation)
-[^2]: Scott, W. T. (1968). Analytic Studies of Cloud Droplet Coalescence I, Journal of Atmospheric Sciences, 25(1), 54-65. Retrieved Feb 18, 2021, from https://journals.ametsoc.org/view/journals/atsc/25/1/1520-0469\_1968\_025\_0054\_asocdc\_2\_0\_co\_2.xml
-[^3]: Ian J. Laurenzi, John D. Bartels, Scott L. Diamond, A General Algorithm for Exact Simulation of Multicomponent Aggregation Processes, Journal of Computational Physics, Volume 177, Issue 2, 2002, Pages 418-449, ISSN 0021-9991, https://doi.org/10.1006/jcph.2002.7017.
+1. [https://en.wikipedia.org/wiki/Smoluchowski\_coagulation\_equation](https://en.wikipedia.org/wiki/Smoluchowski_coagulation_equation)
+2. Scott, W. T. (1968). Analytic Studies of Cloud Droplet Coalescence I, Journal of Atmospheric Sciences, 25(1), 54-65. Retrieved Feb 18, 2021, from https://journals.ametsoc.org/view/journals/atsc/25/1/1520-0469\_1968\_025\_0054\_asocdc\_2\_0\_co\_2.xml
+3. Ian J. Laurenzi, John D. Bartels, Scott L. Diamond, A General Algorithm for Exact Simulation of Multicomponent Aggregation Processes, Journal of Computational Physics, Volume 177, Issue 2, 2002, Pages 418-449, ISSN 0021-9991, https://doi.org/10.1006/jcph.2002.7017.

docs/src/model_creation/parametric_stoichiometry.md

@@ -9,12 +9,18 @@ is a parameter, and the number of products is the product of two parameters.
 ```@example s1
 using Catalyst, Latexify, OrdinaryDiffEq, ModelingToolkit, Plots
 revsys = @reaction_network revsys begin
+    @parameters m::Int64 n::Int64
     k₊, m*A --> (m*n)*B
     k₋, B --> A
 end
 reactions(revsys)
 ```
-Note, as always the `@reaction_network` macro defaults to setting all symbols
+Notice, as described in the [Reaction rate laws used in simulations](@ref introduction_to_catalyst_ratelaws)
+section, the default rate laws involve factorials in the stoichiometric
+coefficients. For this reason we explicitly specify `m` and `n` as integers (as
+otherwise ModelingToolkit will implicitly assume they are floating point).
+
+As always the `@reaction_network` macro defaults to setting all symbols
 neither used as a reaction substrate nor a product to be parameters. Hence, in
 this example we have two species (`A` and `B`) and four parameters (`k₊`, `k₋`,
 `m`, and `n`). In addition, the stoichiometry is applied to the rightmost symbol
@@ -36,21 +42,21 @@ We could have equivalently specified our systems directly via the Catalyst
 API. For example, for `revsys` we would could use
 ```@example s1
 t = default_t()
-@parameters k₊, k₋, m, n
+@parameters k₊ k₋ m::Int n::Int
 @species A(t), B(t)
 rxs = [Reaction(k₊, [A], [B], [m], [m*n]),
        Reaction(k₋, [B], [A])]
 revsys2 = ReactionSystem(rxs,t; name=:revsys)
 revsys2 == revsys
 ```
-which can be simplified using the `@reaction` macro to
+or
 ```@example s1
-rxs2 = [(@reaction k₊, m*A --> (m*n)*B),
+rxs2 = [(@reaction k₊, $m*A --> ($m*$n)*B),
         (@reaction k₋, B --> A)]
 revsys3 = ReactionSystem(rxs2,t; name=:revsys)
 revsys3 == revsys
 ```
-Note, the `@reaction` macro again assumes all symbols are parameters except the
+Here we interpolate in the pre-declared `m` and `n` symbolic variables using `$m` and `$n` to ensure the parameter is known to be integer-valued. The `@reaction` macro again assumes all symbols are parameters except the
 substrates or reactants (i.e. `A` and `B`). For example, in
 `@reaction k, F*A + 2(H*G+B) --> D`, the substrates are `(A,G,B)` with
 stoichiometries `(F,2*H,2)`.
@@ -62,41 +68,41 @@ osys = complete(osys)
 equations(osys)
 show(stdout, MIME"text/plain"(), equations(osys)) # hide
 ```
-Notice, as described in the [Reaction rate laws used in simulations](@ref)
-section, the default rate laws involve factorials in the stoichiometric
-coefficients. For this reason we must specify `m` and `n` as integers, and hence
-*use a tuple for the parameter mapping*
+Specifying the parameter and initial condition values,
 ```@example s1
-p  = (k₊ => 1.0, k₋ => 1.0, m => 2, n => 2)
-u₀ = [A => 1.0, B => 1.0]
+p  = (revsys.k₊ => 1.0, revsys.k₋ => 1.0, revsys.m => 2, revsys.n => 2)
+u₀ = [revsys.A => 1.0, revsys.B => 1.0]
 oprob = ODEProblem(osys, u₀, (0.0, 1.0), p)
 nothing # hide
 ```
-We can now solve and plot the system
-```@julia
+we can now solve and plot the system
+```@example s1
 sol = solve(oprob, Tsit5())
 plot(sol)
 ```
-*If we had used a vector to store parameters, `m` and `n` would be converted to
-floating point giving an error when solving the system.* **Note, currently a [bug](https://github.com/SciML/ModelingToolkit.jl/issues/2296) in ModelingToolkit has broken this example by converting to floating point when using tuple parameters, see the alternative approach below for a workaround.**
 
 An alternative approach to avoid the issues of using mixed floating point and
 integer variables is to disable the rescaling of rate laws as described in
-[Reaction rate laws used in simulations](@ref) section. This requires passing
-the `combinatoric_ratelaws=false` keyword to `convert` or to `ODEProblem` (if
-directly building the problem from a `ReactionSystem` instead of first
-converting to an `ODESystem`). For the previous example this gives the following
-(different) system of ODEs
+[Reaction rate laws used in simulations](@ref introduction_to_catalyst_ratelaws)
+section. This requires passing the `combinatoric_ratelaws=false` keyword to
+`convert` or to `ODEProblem` (if directly building the problem from a
+`ReactionSystem` instead of first converting to an `ODESystem`). For the
+previous example this gives the following (different) system of ODEs where we
+now let `m` and `n` be floating point valued parameters (the default):
 ```@example s1
+revsys = @reaction_network revsys begin
+    k₊, m*A --> (m*n)*B
+    k₋, B --> A
+end
 osys = convert(ODESystem, revsys; combinatoric_ratelaws = false)
 osys = complete(osys)
 equations(osys)
 show(stdout, MIME"text/plain"(), equations(osys)) # hide
 ```
 Since we no longer have factorial functions appearing, our example will now run
-even with floating point values for `m` and `n`:
+with `m` and `n` treated as floating point parameters:
 ```@example s1
-p  = (k₊ => 1.0, k₋ => 1.0, m => 2.0, n => 2.0)
+p  = (revsys.k₊ => 1.0, revsys.k₋ => 1.0, revsys.m => 2.0, revsys.n => 2.0)
 oprob = ODEProblem(osys, u₀, (0.0, 1.0), p)
 sol = solve(oprob, Tsit5())
 plot(sol)