From 8c9b4be42618ea4ac6f97e796f75be9bd6ef18fc Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Mon, 10 Jun 2024 13:35:04 -0400 Subject: [PATCH 1/7] update to master --- docs/src/model_creation/parametric_stoichiometry.md | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/docs/src/model_creation/parametric_stoichiometry.md b/docs/src/model_creation/parametric_stoichiometry.md index a0d370ef0d..1f5c2fab6c 100644 --- a/docs/src/model_creation/parametric_stoichiometry.md +++ b/docs/src/model_creation/parametric_stoichiometry.md @@ -9,6 +9,7 @@ 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 @@ -36,7 +37,7 @@ 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 @species A(t), B(t) rxs = [Reaction(k₊, [A], [B], [m], [m*n]), Reaction(k₋, [B], [A])] From 62b80042460630576b075f97f48ef536857f1474 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Mon, 10 Jun 2024 13:57:56 -0400 Subject: [PATCH 2/7] fix parametric stoich tutorial --- .../parametric_stoichiometry.md | 26 +++++++++++-------- 1 file changed, 15 insertions(+), 11 deletions(-) diff --git a/docs/src/model_creation/parametric_stoichiometry.md b/docs/src/model_creation/parametric_stoichiometry.md index 1f5c2fab6c..547fefdcb2 100644 --- a/docs/src/model_creation/parametric_stoichiometry.md +++ b/docs/src/model_creation/parametric_stoichiometry.md @@ -15,7 +15,12 @@ revsys = @reaction_network revsys begin 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) +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 @@ -37,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::Int 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)`. @@ -63,10 +68,6 @@ 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* ```@example s1 p = (k₊ => 1.0, k₋ => 1.0, m => 2, n => 2) u₀ = [A => 1.0, B => 1.0] @@ -78,8 +79,6 @@ We can now solve and plot the system 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 @@ -89,6 +88,11 @@ 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 ```@example s1 +revsys = @reaction_network revsys begin + @parameters m::Int64 n::Int64 + k₊, m*A --> (m*n)*B + k₋, B --> A +end osys = convert(ODESystem, revsys; combinatoric_ratelaws = false) osys = complete(osys) equations(osys) From 8de3c9e9108217faaef1437a51cba26f7179ca78 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Mon, 10 Jun 2024 14:17:40 -0400 Subject: [PATCH 3/7] fixes --- docs/make.jl | 2 +- docs/src/assets/Project.toml | 6 ++--- .../smoluchowski_coagulation_equation.md | 6 ++--- .../parametric_stoichiometry.md | 22 ++++++++++--------- 4 files changed, 19 insertions(+), 17 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 8595cd6ffa..32df32290d 100644 --- a/docs/make.jl +++ b/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) diff --git a/docs/src/assets/Project.toml b/docs/src/assets/Project.toml index 83ecae3cc6..b92da6ef5d 100644 --- a/docs/src/assets/Project.toml +++ b/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" diff --git a/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md b/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md index dc9427c979..d8865651b4 100644 --- a/docs/src/model_creation/examples/smoluchowski_coagulation_equation.md +++ b/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. diff --git a/docs/src/model_creation/parametric_stoichiometry.md b/docs/src/model_creation/parametric_stoichiometry.md index 547fefdcb2..2704625b25 100644 --- a/docs/src/model_creation/parametric_stoichiometry.md +++ b/docs/src/model_creation/parametric_stoichiometry.md @@ -15,7 +15,7 @@ revsys = @reaction_network revsys begin end reactions(revsys) ``` -Notice, as described in the [Reaction rate laws used in simulations](@ref) +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). @@ -68,28 +68,30 @@ osys = complete(osys) equations(osys) show(stdout, MIME"text/plain"(), equations(osys)) # hide ``` +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] 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) ``` 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 - @parameters m::Int64 n::Int64 + @parameters m n k₊, m*A --> (m*n)*B k₋, B --> A end @@ -99,7 +101,7 @@ 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) oprob = ODEProblem(osys, u₀, (0.0, 1.0), p) From 356c8b122bdee8612e23d70a42df5480f47c4a69 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Mon, 10 Jun 2024 14:20:31 -0400 Subject: [PATCH 4/7] more fixes --- docs/src/model_creation/parametric_stoichiometry.md | 1 - 1 file changed, 1 deletion(-) diff --git a/docs/src/model_creation/parametric_stoichiometry.md b/docs/src/model_creation/parametric_stoichiometry.md index 2704625b25..8fb0b53dd1 100644 --- a/docs/src/model_creation/parametric_stoichiometry.md +++ b/docs/src/model_creation/parametric_stoichiometry.md @@ -91,7 +91,6 @@ 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 - @parameters m n k₊, m*A --> (m*n)*B k₋, B --> A end From 0fd6f16c3a8b675692f3405e5e431076e97f1947 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Mon, 10 Jun 2024 14:22:39 -0400 Subject: [PATCH 5/7] restore warning level --- docs/make.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/docs/make.jl b/docs/make.jl index 32df32290d..8595cd6ffa 100644 --- a/docs/make.jl +++ b/docs/make.jl @@ -41,7 +41,7 @@ makedocs(sitename = "Catalyst.jl", clean = true, pages = pages, pagesonly = true, - warnonly = true) + warnonly = [:missing_docs]) deploydocs(repo = "github.com/SciML/Catalyst.jl.git"; push_preview = true) From 494d185c3ea4f7c7e19a19f2e499a42e91e4969f Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Mon, 10 Jun 2024 14:25:53 -0400 Subject: [PATCH 6/7] namespace mappings --- docs/src/model_creation/parametric_stoichiometry.md | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/src/model_creation/parametric_stoichiometry.md b/docs/src/model_creation/parametric_stoichiometry.md index 8fb0b53dd1..a219abfce5 100644 --- a/docs/src/model_creation/parametric_stoichiometry.md +++ b/docs/src/model_creation/parametric_stoichiometry.md @@ -70,8 +70,8 @@ show(stdout, MIME"text/plain"(), equations(osys)) # hide ``` 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 ``` @@ -102,7 +102,7 @@ show(stdout, MIME"text/plain"(), equations(osys)) # hide Since we no longer have factorial functions appearing, our example will now run 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) From c4f18d62a7a15b64a7ddce40bff8a8f9ab213d50 Mon Sep 17 00:00:00 2001 From: Sam Isaacson Date: Mon, 10 Jun 2024 14:56:13 -0400 Subject: [PATCH 7/7] final fix --- docs/make.jl | 2 +- docs/src/introduction_to_catalyst/introduction_to_catalyst.md | 4 ++-- 2 files changed, 3 insertions(+), 3 deletions(-) diff --git a/docs/make.jl b/docs/make.jl index 8595cd6ffa..32df32290d 100644 --- a/docs/make.jl +++ b/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) diff --git a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md b/docs/src/introduction_to_catalyst/introduction_to_catalyst.md index 94b36738ca..452c9cd9c6 100644 --- a/docs/src/introduction_to_catalyst/introduction_to_catalyst.md +++ b/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) \ No newline at end of file +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) \ No newline at end of file