Merge branch 'master' into fix_broken_doc_references

Author: TorkelE

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/introduction_to_catalyst/introduction_to_catalyst.md

@@ -181,7 +181,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
@@ -360,4 +360,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)

test/reactionsystem_core/reactionsystem.jl

@@ -159,7 +159,6 @@ end
 
 # Test with JumpSystem.
 let
-    p = rand(rng, length(k))
     @species A(t) B(t) C(t) D(t) E(t) F(t)
     rxs = [Reaction(k[1], nothing, [A]),            # 0 -> A
         Reaction(k[2], [B], nothing),            # B -> 0
@@ -193,27 +192,29 @@ let
     @test all(map(i -> typeof(equations(js)[i]) <: JumpProcesses.ConstantRateJump, cidxs))
     @test all(map(i -> typeof(equations(js)[i]) <: JumpProcesses.VariableRateJump, vidxs))
 
-    pars = rand(rng, length(k))
+    p = rand(rng, length(k))
+    pmap = parameters(js) .=> p
     u0 = rand(rng, 2:10, 6)
+    u0map = unknowns(js) .=> u0
     ttt = rand(rng)
     jumps = Vector{Union{ConstantRateJump, MassActionJump, VariableRateJump}}(undef,
                                                                               length(rxs))
 
-    jumps[1] = MassActionJump(pars[1], Vector{Pair{Int, Int}}(), [1 => 1])
-    jumps[2] = MassActionJump(pars[2], [2 => 1], [2 => -1])
-    jumps[3] = MassActionJump(pars[3], [1 => 1], [1 => -1, 3 => 1])
-    jumps[4] = MassActionJump(pars[4], [3 => 1], [1 => 1, 2 => 1, 3 => -1])
-    jumps[5] = MassActionJump(pars[5], [3 => 1], [1 => 2, 3 => -1])
-    jumps[6] = MassActionJump(pars[6], [1 => 1, 2 => 1], [1 => -1, 2 => -1, 3 => 1])
-    jumps[7] = MassActionJump(pars[7], [2 => 2], [1 => 1, 2 => -2])
-    jumps[8] = MassActionJump(pars[8], [1 => 1, 2 => 1], [2 => -1, 3 => 1])
-    jumps[9] = MassActionJump(pars[9], [1 => 1, 2 => 1], [1 => -1, 2 => -1, 3 => 1, 4 => 1])
-    jumps[10] = MassActionJump(pars[10], [1 => 2], [1 => -2, 3 => 1, 4 => 1])
-    jumps[11] = MassActionJump(pars[11], [1 => 2], [1 => -1, 2 => 1])
-    jumps[12] = MassActionJump(pars[12], [1 => 1, 2 => 3, 3 => 4],
+    jumps[1] = MassActionJump(p[1], Vector{Pair{Int, Int}}(), [1 => 1])
+    jumps[2] = MassActionJump(p[2], [2 => 1], [2 => -1])
+    jumps[3] = MassActionJump(p[3], [1 => 1], [1 => -1, 3 => 1])
+    jumps[4] = MassActionJump(p[4], [3 => 1], [1 => 1, 2 => 1, 3 => -1])
+    jumps[5] = MassActionJump(p[5], [3 => 1], [1 => 2, 3 => -1])
+    jumps[6] = MassActionJump(p[6], [1 => 1, 2 => 1], [1 => -1, 2 => -1, 3 => 1])
+    jumps[7] = MassActionJump(p[7], [2 => 2], [1 => 1, 2 => -2])
+    jumps[8] = MassActionJump(p[8], [1 => 1, 2 => 1], [2 => -1, 3 => 1])
+    jumps[9] = MassActionJump(p[9], [1 => 1, 2 => 1], [1 => -1, 2 => -1, 3 => 1, 4 => 1])
+    jumps[10] = MassActionJump(p[10], [1 => 2], [1 => -2, 3 => 1, 4 => 1])
+    jumps[11] = MassActionJump(p[11], [1 => 2], [1 => -1, 2 => 1])
+    jumps[12] = MassActionJump(p[12], [1 => 1, 2 => 3, 3 => 4],
                                [1 => -1, 2 => -3, 3 => -2, 4 => 3])
-    jumps[13] = MassActionJump(pars[13], [1 => 3, 2 => 1], [1 => -3, 2 => -1])
-    jumps[14] = MassActionJump(pars[14], Vector{Pair{Int, Int}}(), [1 => 2])
+    jumps[13] = MassActionJump(p[13], [1 => 3, 2 => 1], [1 => -3, 2 => -1])
+    jumps[14] = MassActionJump(p[14], Vector{Pair{Int, Int}}(), [1 => 2])
 
     jumps[15] = ConstantRateJump((u, p, t) -> p[15] * u[1] / (2 + u[1]),
                                  integrator -> (integrator.u[1] -= 1))
@@ -230,32 +231,34 @@ let
                                  integrator -> (integrator.u[4] -= 2; integrator.u[5] -= 1; integrator.u[6] += 2))
 
     unknownoid = Dict(unknown => i for (i, unknown) in enumerate(unknowns(js)))
-    jspmapper = ModelingToolkit.JumpSysMajParamMapper(js, pars)
+    dprob = DiscreteProblem(js, u0map, (0.0, 10.0), pmap)
+    mtkpars = dprob.p
+    jspmapper = ModelingToolkit.JumpSysMajParamMapper(js, mtkpars)
     symmaj = ModelingToolkit.assemble_maj(equations(js).x[1], unknownoid, jspmapper)
-    maj = MassActionJump(symmaj.param_mapper(pars), symmaj.reactant_stoch, symmaj.net_stoch,
+    maj = MassActionJump(symmaj.param_mapper(mtkpars), symmaj.reactant_stoch, symmaj.net_stoch,
                          symmaj.param_mapper, scale_rates = false)
     for i in midxs
-        @test_broken abs(jumps[i].scaled_rates - maj.scaled_rates[i]) < 100 * eps()
+        @test abs(jumps[i].scaled_rates - maj.scaled_rates[i]) < 100 * eps()
         @test jumps[i].reactant_stoch == maj.reactant_stoch[i]
         @test jumps[i].net_stoch == maj.net_stoch[i]
     end
     for i in cidxs
         crj = ModelingToolkit.assemble_crj(js, equations(js)[i], unknownoid)
-        @test_broken isapprox(crj.rate(u0, p, ttt), jumps[i].rate(u0, p, ttt))
-        fake_integrator1 = (u = zeros(6), p = p, t = 0.0)
-        fake_integrator2 = deepcopy(fake_integrator1)
+        @test isapprox(crj.rate(u0, mtkpars, ttt), jumps[i].rate(u0, p, ttt))
+        fake_integrator1 = (u = zeros(6), p = mtkpars, t = 0.0)
+        fake_integrator2 = (u = zeros(6), p, t = 0.0)
         crj.affect!(fake_integrator1)
         jumps[i].affect!(fake_integrator2)
-        @test fake_integrator1 == fake_integrator2
+        @test fake_integrator1.u == fake_integrator2.u
     end
     for i in vidxs
         crj = ModelingToolkit.assemble_vrj(js, equations(js)[i], unknownoid)
-        @test_broken isapprox(crj.rate(u0, p, ttt), jumps[i].rate(u0, p, ttt))
-        fake_integrator1 = (u = zeros(6), p = p, t = 0.0)
-        fake_integrator2 = deepcopy(fake_integrator1)
+        @test isapprox(crj.rate(u0, mtkpars, ttt), jumps[i].rate(u0, p, ttt))
+        fake_integrator1 = (u = zeros(6), p = mtkpars, t = 0.0)
+        fake_integrator2 = (u = zeros(6), p, t = 0.0)
         crj.affect!(fake_integrator1)
         jumps[i].affect!(fake_integrator2)
-        @test fake_integrator1 == fake_integrator2
+        @test fake_integrator1.u == fake_integrator2.u
     end
 end
 
@@ -439,7 +442,7 @@ let
         umean += sol(10.0, idxs = [B1, B2, B3, C])
     end
     umean /= Nsims
-    @test isapprox(umean[1], umean[2]; rtol = 1e-2) 
+    @test isapprox(umean[1], umean[2]; rtol = 1e-2)
     @test isapprox(umean[1], umean[3]; rtol = 1e-2)
     @test umean[4] == 10
 end