Merge pull request #469 from isaacsas/use_symmaps_in_docs

Use symmaps in Problems and explain in the docs

Author: isaacsas

HISTORY.md

@@ -1,10 +1,47 @@
 # Breaking updates and feature summaries across releases
 
 ## Catalyst unreleased (master branch) 
-- Added `symmap_to_varmap` and `setdefaults!` to allow setting initial conditions and parameter values using symbols.
+
+## Catalyst 10.4
+- Added `symmap_to_varmap`, `setdefaults!`, and updated all `*Problem(rn,...)`
+  calls to allow setting initial conditions and parameter values using symbol
+  maps. See the [Catalyst API](https://catalyst.sciml.ai/dev/) for details.
+  These allow using regular Julia `Symbols` to specify parameter values and
+  initial conditions. i.e. to set defaults we can do
+  ```julia
+  rn = @reaction_network begin
+      α, S + I --> 2I
+      β, I --> R
+  end α β
+  setdefaults!(rn, [:S => 999.0, :I => 1.0, :R => 0.0, :α => 1e-4, :β => .01])
+  op    = ODEProblem(rn, [], (0.0,250.0), [])
+  sol   = solve(op, Tsit5()) 
+  ```
+  To explicitly pass initial conditions and parameters using symbols we can do
+  ```julia
+  rn = @reaction_network begin
+      α, S + I --> 2I
+      β, I --> R
+  end α β
+  u0 = [:S => 999.0, :I => 1.0, :R => 0.0]
+  p  = (:α => 1e-4, :β => .01)
+  op    = ODEProblem(rn, u0, (0.0,250.0), p)
+  sol   = solve(op, Tsit5())  
+  ```
+  In each case ModelingToolkit symbolic variables can be used instead of
+  `Symbol`s, e.g.
+  ```julia
+  @parameters α β
+  @variables t S(t) I(t) R(t)
+  setdefaults!(rn, [S => 999.0, I => 1.0, R => 0.0, α => 1e-4, β => .01])
+  ```
 
 ## Catalyst 10.3
-- **BREAKING:** The order of the parameters in the `ReactionSystem`'s `.ps` field has been changed (only when created through the `@reaction_network` macro). Previously they were ordered according to the order with which they appeared in the macro. Now they are ordered according the to order with which they appeard after the `end` part. E.g. in
+- **BREAKING:** The order of the parameters in the `ReactionSystem`'s `.ps`
+  field has been changed (only when created through the `@reaction_network`
+  macro). Previously they were ordered according to the order with which they
+  appeared in the macro. Now they are ordered according the to order with which
+  they appeard after the `end` part. E.g. in
   ```julia
   rn = @reaction_network begin
     (p,d), 0 <--> X

README.md

@@ -85,17 +85,18 @@ the current master branch.
 #### Gillespie Simulations of Michaelis-Menten Enzyme Kinetics
 
 ```julia
+using Catalyst, Plots, DiffEqJump
 rs = @reaction_network begin
   c1, S + E --> SE
   c2, SE --> S + E
   c3, SE --> P + E
 end c1 c2 c3
-p = (0.00166,0.0001,0.1)   # [c1,c2,c3]
+p     = (:c1 => 0.00166, :c2 => 0.0001, :c3 => 0.1)   
 tspan = (0., 100.)
-u0 = [301., 100., 0., 0.]  # [S,E,SE,P]
+u0    = [:S => 301, :E => 100, :SE => 0, :P => 0]  
 
 # solve JumpProblem
-dprob = DiscreteProblem(rs, species(rs) .=> u0, tspan, parameters(rs) .=> p)
+dprob = DiscreteProblem(rs, u0, tspan, p)
 jprob = JumpProblem(rs, dprob, Direct())
 jsol = solve(jprob, SSAStepper())
 plot(jsol,lw=2,title="Gillespie: Michaelis-Menten Enzyme Kinetics")
@@ -106,16 +107,16 @@ plot(jsol,lw=2,title="Gillespie: Michaelis-Menten Enzyme Kinetics")
 #### Adaptive SDEs for A Birth-Death Process
 
 ```julia
-using Catalyst, Plots, StochasticDiffEq, DiffEqJump
+using Catalyst, Plots, StochasticDiffEq
 rs = @reaction_network begin
   c1, X --> 2X
   c2, X --> 0
   c3, 0 --> X
 end c1 c2 c3
-p = (1.0,2.0,50.) # [c1,c2,c3]
+p     = (:c1 => 1.0, :c2 => 2.0, :c3 => 50.) 
 tspan = (0.,10.)
-u0 = [5.]         # [X]
-sprob = SDEProblem(rs, species(rs) .=> u0, tspan, parameters(rs) .=> p)
+u0    = [:X => 5.]        
+sprob = SDEProblem(rs, u0, tspan, p)
 ssol  = solve(sprob, LambaEM(), reltol=1e-3)
 plot(ssol,lw=2,title="Adaptive SDE: Birth-Death Process")
 ```

docs/src/api/catalyst_api.md

@@ -185,7 +185,6 @@ reorder_states!
 addparam!
 addreaction!
 setdefaults!
-symmap_to_varmap
 ModelingToolkit.extend
 ModelingToolkit.compose
 Catalyst.flatten
@@ -260,3 +259,8 @@ ModelingToolkit.structural_simplify
 validate(rx::Reaction; info::String = "")
 validate(rs::ReactionSystem, info::String="")
 ```
+
+## Utility Functions
+```@docs
+symmap_to_varmap
+```

docs/src/faqs.md

@@ -13,6 +13,67 @@ osys = convert(ODESystem, rn; combinatoric_ratelaws=false)
 Disabling these rescalings should work for all conversions of `ReactionSystem`s
 to other `ModelingToolkit.AbstractSystem`s.
 
+## How to set default values for initial conditions and parameters?
+When directly constructing a `ReactionSystem` these can be passed to the
+constructor, and allow solving the system without needing initial condition or
+parameter vectors in the generated problem. For example
+```julia
+using Catalyst, Plots, OrdinaryDiffEq
+@parameters β ν
+@variables t S(t) I(t) R(t)
+rx1 = Reaction(β, [S,I], [I], [1,1], [2])
+rx2 = Reaction(ν, [I], [R])
+defs = [β => 1e-4, ν => .01, S => 999.0, I => 1.0, R => 0.0]
+@named sir = ReactionSystem([rx1,rx2],t; defaults=defs)
+oprob = ODEProblem(sir, [], (0.0,250.0))
+sol = solve(oprob, Tsit5())
+plot(sol)
+```
+alternatively we could also have said
+```julia
+@parameters β=1e-4 ν=.01
+@variables t S(t)=999.0 I(t)=1.0 R(t)=0.0
+rx1 = Reaction(β, [S,I], [I], [1,1], [2])
+rx2 = Reaction(ν, [I], [R])
+@named sir = ReactionSystem([rx1,rx2],t)
+```
+
+The `@reaction_network` macro does not currently provide a way to specify
+default values, however, they can be added after creating the system via the
+`setdefaults!` command, like
+```julia
+sir = @reaction_network sir begin
+    β, S + I --> 2I
+    ν, I --> R
+end β ν
+setdefaults!(sir, [:β => 1e-4, :ν => .01, :S => 999.0, :I => 1.0, :R => 0.0])
+```
+
+## How to specify initial conditions and parameters values for `ODEProblem` and other problem types?
+To explicitly pass initial conditions and parameters we can use mappings from
+Julia `Symbol`s corresponding to each variable/parameter to values, or from
+ModelingToolkit symbolic variables to each variable/parameter. Using `Symbol`s
+we have
+```julia
+rn = @reaction_network begin
+    α, S + I --> 2I
+    β, I --> R
+end α β
+u0 = [:S => 999.0, :I => 1.0, :R => 0.0]
+p  = (:α => 1e-4, :β => .01)
+op  = ODEProblem(rn, u0, (0.0,250.0), p)
+sol = solve(op, Tsit5())  
+```
+while using ModelingToolkit symbolic variables we have
+```julia
+@parameters α β
+@variables t S(t) I(t) R(t)
+u0 = [S => 999.0, I => 1.0, R => 0.0]
+p  = (α => 1e-4, β => .01)
+op  = ODEProblem(rn, u0, (0.0,250.0), p)
+sol = solve(op, Tsit5())  
+```
+
 ## How to modify generated ODEs?
 Conversion to other `ModelingToolkit.AbstractSystem`s allows the possibility to
 modify the system with further terms that are difficult to encode as a chemical

docs/src/index.md

@@ -80,11 +80,11 @@ which in Jupyter notebooks will give the figure
 
 To generate and solve a mass action ODE version of the model we use
 ```julia
-using DiffEqBase, OrdinaryDiffEq
-p     = [.1/1000, .01]           # [α,β]
+using OrdinaryDiffEq
+p     = [:α => .1/1000, :β => .01]
 tspan = (0.0,250.0)
-u0    = [999.0,1.0,0.0]          # [S,I,R] at t=0
-op    = ODEProblem(rn, species(rn) .=> u0, tspan, parameters(rn) .=> p)
+u0    = [:S => 999.0, :I => 1.0, :R => 0.0]
+op    = ODEProblem(rn, u0, tspan, p)
 sol   = solve(op, Tsit5())       # use Tsit5 ODE solver
 ```
 which we can plot as

docs/src/tutorials/advanced_examples.md

@@ -8,20 +8,21 @@ Jacobians with or without threading and/or parallelization), creating LaTeX
 representations for systems, etc.
 
 Note, when generating problems from other system types, `u0` and `p` must
-provide vectors of `Pair`s that map each species or parameter to their numerical
-value. E.g., for the Michaelis-Menten example above we'd use
+provide vectors, tuples or dictionaries of `Pair`s that map each the symbolic
+variables for each species or parameter to their numerical value. E.g., for the
+Michaelis-Menten example above we'd use
 ```julia
 rs = @reaction_network begin
   c1, X --> 2X
   c2, X --> 0
   c3, 0 --> X
 end c1 c2 c3
-p     = (1.0,2.0,50.)
+p     = (:c1 => 1.0, :c2 => 2.0, :c3 => 50.)
+pmap  = symmap_to_varmap(rs,p)   # convert Symbol map to symbolic variable map
 tspan = (0.,4.)
-u0    = [5.]   
+u0    = [:X => 5.]   
+u0map = symmap_to_varmap(rs,u0)  # convert Symbol map to symbolic variable map
 osys  = convert(ODESystem, rs)
-u0map = map((x,y) -> Pair(x,y), species(rs), u0)
-pmap  = map((x,y) -> Pair(x,y), parameters(rs), p)
 oprob = ODEProblem(osys, u0map, tspan, pmap)
 sol   = solve(oprob, Tsit5())
 ```

docs/src/tutorials/basic_examples.md

@@ -8,9 +8,9 @@ rs = @reaction_network begin
   c2, X --> 0
   c3, 0 --> X
 end c1 c2 c3
-p = (1.0,2.0,50.) # [c1,c2,c3]
+p = (:c1 => 1.0, :c2 => 2.0, :c3 => 50.) 
 tspan = (0.,4.)
-u0 = [5.]         # [X]
+u0 = [:X => 5.]
 
 # solve ODEs
 oprob = ODEProblem(rs, u0, tspan, p)
@@ -25,7 +25,7 @@ sprob = SDEProblem(rs, u0, tspan, p)
 ssol  = solve(sprob, EM(), dt=.01)
 
 # solve JumpProblem
-u0 = [5]
+u0 = [:X => 5]
 dprob = DiscreteProblem(rs, u0, tspan, p)
 jprob = JumpProblem(rs, dprob, Direct())
 jsol = solve(jprob, SSAStepper())
@@ -39,16 +39,16 @@ rs = @reaction_network begin
   c2, SE --> S + E
   c3, SE --> P + E
 end c1 c2 c3
-p = (0.00166,0.0001,0.1)   # [c1,c2,c3]
+p = (:c1 => 0.00166, :c2 => 0.0001, :c3 => 0.1) 
 tspan = (0., 100.)
-u0 = [301., 100., 0., 0.]  # [S,E,SE,P]
+u0 = [:S => 301., :E => 100., :SE => 0., :P => 0.]  
 
 # solve ODEs
 oprob = ODEProblem(rs, u0, tspan, p)
 osol  = solve(oprob, Tsit5())
 
 # solve JumpProblem
-u0 = [301, 100, 0, 0] 
+u0 = [:S => 301, :E => 100, :SE => 0, :P => 0] 
 dprob = DiscreteProblem(rs, u0, tspan, p)
 jprob = JumpProblem(rs, dprob, Direct())
 jsol = solve(jprob, SSAStepper())

docs/src/tutorials/dsl.md

@@ -35,16 +35,42 @@ these are the same as for normal variables in Julia.
 
 The chemical reaction model is generated by the `@reaction_network` macro and
 stored in the `rn` variable (a normal Julia variable, which does not need to be
-called `rn`). The generated `ReactionSystem` can be converted to a differential
-equation model via
+called `rn`). It corresponds to a [`ReactionSystem`](@ref), a symbolic
+representation of the chemical network.  The generated `ReactionSystem` can be
+converted to a symbolic differential equation model via
 ```julia
 osys  = convert(ODESystem, rn)
-oprob = ODEProblem(osys, species(rn) .=> u0, tspan, parameters(rn) .=> p)
 ```
-or more directly via
+
+We can then convert the symbolic ODE model into a compiled, optimized
+representation for use in the SciML ODE solvers by constructing an `ODEProblem`.
+Creating an `ODEProblem` also requires our specifying the initial conditions for
+the model. We do this by creating a mapping from each symbolic variable
+representing a chemical species to its initial value
+```julia
+# define the symbolic variables
+@variables t, X(t), Y(t), Z(t), XY(t), Z1(t), Z2(t)
+
+# create the mapping
+u0 = [X => 1.0, Y => 1.0, XY => 1.0, Z1 => 1.0, Z2 => 1.0]
+```
+Alternatively, we can create a mapping use Julia `Symbol`s for each variable, and then convert them to a mapping involving symbolic variables like
+```julia
+u0 = symmap_to_varmap(rn, [:X => 1.0, :Y => 1.0, :XY => 1.0, :Z1 => 1.0, :Z2 => 1.0]) 
+```
+Given the mapping, we can then create an `ODEProblem` from our symbolic `ODESystem`
 ```julia
-oprob = ODEProblem(rn, species(rn) .=> u0, tspan, parameters(rn) .=> p)
+oprob = ODEProblem(osys, u0, tspan, [])
 ```
+
+Catalyst provides a shortcut to avoid having to explicitly `convert` to an
+`ODESystem` and/or use `symmap_to_varmap`, allowing direct construction 
+of the `ODEProblem` like
+```julia
+u0 = [:X => 1.0, :Y => 1.0, :XY => 1.0, :Z1 => 1.0, :Z2 => 1.0]
+oprob = ODEProblem(rn, u0, tspan, [])
+```
+
 For more detailed examples, see the [Basic Chemical Reaction Network
 Examples](@ref). The generated differential equations use the law of mass
 action. For the above example, the ODEs are then