Merge pull request #407 from isaacsas/update-api-for-subsystems

Update API and docs for subsystems

HISTORY.md

@@ -8,6 +8,11 @@
   subsystems it is not currently possible to convert to a `JumpSystem`. It is
   also not possible to include either `SDESystem`s or `JumpSystems` as
   subsystems.
+- Depreciated `merge`, use `ModelingToolkit.extend` instead.
+- Depreciated `params` and `numparams` (use `ModelingToolkit.parameters` to get
+  all parameters of a system and all subsystems, or use `reactionparams` to get
+  all parameters of a system and all `ReactionSystem` subsystems. The latter
+  correspond to those parameters used within `Reaction`s.)
 
 ## Catalyst 9.0
 *1.* **BREAKING:** `netstoichmat`, `prodstoichmat` and `substoichmat` are now

docs/src/api/catalyst_api.md

@@ -76,17 +76,33 @@ ReactionSystem
 ```
 
 ## Basic System Properties
+See [The generated `ReactionSystem` and `Reaction`s](@ref) for more details
+
 ```@docs
 species
 speciesmap
-params
+reactionparams
 paramsmap
 reactions
 numspecies
-numparams
 numreactions
+numreactionparams
 ```
 
+## ModelingToolkit-Inherited Accessor Functions
+See [The generated `ReactionSystem` and `Reaction`s](@ref) for more details
+
+- `ModelingToolkit.get_eqs(sys)`: The reactions of the system (ignores subsystems).
+- `ModelingToolkit.equations(sys)`: Collects all reactions and equations from
+  the system and all subsystems.
+- `ModelingToolkit.get_states(sys)`: The set of chemical species in the system (ignores subsystems).
+- `ModelingToolkit.states(sys)`: Collects all species and states from the system and all subsystems.
+- `ModelingToolkit.get_ps(sys)`: The parameters of the system (ignores subsystems).
+- `ModelingToolkit.parameters(sys)`: Collects all parameters from the system and all subsystems.
+- `ModelingToolkit.get_iv(sys)`: The independent variable of the system, usually time.
+- `ModelingToolkit.get_systems(sys)`: The sub-systems of `sys`.
+- `ModelingToolkit.get_defaults(sys)`: The default values for parameters and initial conditions for `sys`.
+
 ## Basic Reaction Properties
 ```@docs
 ismassaction
@@ -98,21 +114,15 @@ netstoichmat
 reactionrates
 ```
 
-## Composition and Accessor Functions for [`ReactionSystem`](@ref)s
-- `ModelingToolkit.get_eqs(sys)` or `equations(sys)`: The reactions that define the system.
-- `ModelingToolkit.get_states(sys)` or `states(sys)`: The set of chemical species in the system.
-- `ModelingToolkit.get_ps(sys)` or `parameters(sys)`: The parameters of the system.
-- `ModelingToolkit.get_iv(sys)`: The independent variable of the reaction
-  system, usually time.
-
 ## Functions to Extend a Network
 ```@docs
 @add_reactions
 addspecies!
 addparam!
 addreaction!
+ModelingToolkit.extend
+ModelingToolkit.compose
 merge!(network1::ReactionSystem, network2::ReactionSystem)
-merge(network1::ReactionSystem, network2::ReactionSystem)
 ```
 
 ## Network Analysis and Representations
@@ -136,8 +146,9 @@ isweaklyreversible
 
 ## Network Comparison 
 ```@docs
-==(rn1::ReactionSystem, rn2::ReactionSystem)
 ==(rn1::Reaction, rn2::Reaction)
+isequal_without_names
+==(rn1::ReactionSystem, rn2::ReactionSystem)
 ```
 
 ## Network Visualization

docs/src/tutorials/advanced_examples.md

@@ -21,7 +21,7 @@ tspan = (0.,4.)
 u0    = [5.]   
 osys  = convert(ODESystem, rs)
 u0map = map((x,y) -> Pair(x,y), species(rs), u0)
-pmap  = map((x,y) -> Pair(x,y), params(rs), p)
+pmap  = map((x,y) -> Pair(x,y), parameters(rs), p)
 oprob = ODEProblem(osys, u0map, tspan, pmap)
 sol   = solve(oprob, Tsit5())
 ```

docs/src/tutorials/basics.md

@@ -35,7 +35,7 @@ called `rn`). The generated `ReactionSystem` can be converted to a differential
 equation model via
 ```julia
 osys = convert(ODESystem, rn)
-oprob = ODEProblem(osys, Pair.(species(rn),u0), tspan, Pair.(params(rn),p))
+oprob = ODEProblem(osys, Pair.(species(rn),u0), tspan, Pair.(parameters(rn),p))
 ```
 or more directly via
 ```julia

docs/src/tutorials/generated_systems.md

@@ -1,18 +1,73 @@
-# The generated [`ReactionSystem`](@ref) and [`Reaction`](@ref)s
-The `@reaction_network` macro generates a [`ReactionSystem`](@ref) object, which
-has a number of fields that can be accessed directly or via the [Catalyst.jl
-API](@ref) (the recommended route). Below we list these components, with the recommended
-API method listed first:
-
-* [`species(rn)`](@ref) and `states(rn)` is a vector of all the chemical
-  species within the system, each represented as a `ModelingToolkit.Term`.
-* [`params(rn)`](@ref) and `parameters(rn)` is a vector of all the parameters
-  within the system, each represented as a `ModelingToolkit.Sym`.
-* [`reactions(rn)`](@ref) and `equations(rn)` is a vector of all the
-  `Reaction`s within the system.
-* `independent_variables(rn)` and `ModelingToolkit.get_iv(rn)` are the
-  independent variable of the system, usually `t` for time, represented as a
-  `ModelingToolkit.Sym`.
+# The generated `ReactionSystem` and `Reaction`s
+
+### [`ReactionSystem`](@ref) Accessors
+
+The `@reaction_network` macro generates a [`ReactionSystem`](@ref), an instance
+of a `ModelingToolkit.AbstractTimeDependentSystem`. It has a number of fields
+that can be accessed using the [Catalyst.jl API](@ref) or the
+[ModelingToolkit.jl Abstract System
+Interface](https://mtk.sciml.ai/dev/basics/AbstractSystem/). Below we overview
+these components.
+
+There are three basic sets of convenience accessors that will return information
+either from a top-level system, the top-level system and all sub-systems that
+are also `ReactionSystem`s (i.e. the full reaction-network), or the top-level
+system and all sub-systems (i.e. the full model). To retrieve info from just a
+base [`ReactionSystem`](@ref) `rn`, ignoring sub-systems of `rn`, one can use the ModelingToolkit
+accessors:
+
+* `get_states(rn)` is a vector that collects all the species defined within
+  `rn`. 
+* `get_ps(rn)` is a vector that collects all the parameters defined within
+  reactions in `rn`. 
+* `get_eqs(rn)` is a vector that collects all the [`Reaction`](@ref)s defined
+  within `rn`.
+* `get_iv(rn)` is the independent variable used in the system (usually `t` to
+  represent time).
+* `get_systems(rn)` is a vector of all sub-systems of `rn`.
+* `get_defaults(rn)` is a dictionary of all the default values for parameters
+  and species in `rn`.
+
+These accessors do not allocate, directly accessing internal fields of the
+`ReactionSystem`.
+
+To retrieve the corresponding information from the full reaction network
+represented by a system `rn`, which corresponds to information within both `rn`
+and all sub-systems of type `ReactionSystem`, one can call:
+
+* [`species(rn)`](@ref) is a vector collecting all the chemical species within
+  the system and any sub-systems that are also `ReactionSystems`. 
+* [`reactionparams(rn)`](@ref) is a vector of all the parameters
+  within the system and any sub-systems that are also `ReactionSystem`s.
+* [`reactions(rn)`](@ref) is a vector of all the `Reaction`s within the system
+  and any sub-systems that are also `ReactionSystem`s.
+
+These accessors will allocate unless there are no subsystems. In the latter case
+they are equivalent to the corresponding `get_*` functions.
+
+Finally, as some sub-systems may be other system types, for example specifying
+algebraic constraints with a `NonlinearSystem`, it can also be convenient to
+collect all state variables (e.g. species and algebraic variables) and such. The
+following ModelingToolkit functions provide this information
+
+* [`states(rn)`](@ref) returns all species *and variables* across the system and
+  *all sub-systems*.
+* [`parameters(rn)`](@ref) returns all parameters across the system and *all
+  sub-systems*.
+* [`equations(rn)`](@ref) returns all [`Reaction`](@ref)s and all
+  [`Equations`](@ref) defined across the system and *all sub-systems*.
+
+These accessors will allocate unless there are no subsystems. In the latter case
+they are equivalent to the corresponding `get_*` functions.
+
+Empty `ReactionSystem`s can be generated via [`make_empty_network`](@ref) or
+[`@reaction_network`](@ref) with no arguments (giving one argument to the latter
+will specify a system name). `ReactionSystem`s can be programmatically extended
+using [`addspecies!`](@ref), [`addparam!`](@ref), [`addreaction!`](@ref),
+[`@add_reactions`](@ref), or composed using [`ModelingToolkit.extend`](@ref) and
+[`ModelingToolkit.compose`](@ref).
+
+### [`Reaction`](@ref) fields
 
 Each `Reaction` within `reactions(rn)` has a number of subfields. For `rx` a
 `Reaction` we have:
@@ -24,7 +79,7 @@ Each `Reaction` within `reactions(rn)` has a number of subfields. For `rx` a
   each substrate species in `rx.substrates`.
 * `rx.prodstoich`, a vector storing the corresponding integer stoichiometry of
   each product species in `rx.products`.
-* `rx.rate`, a `Number, `ModelingToolkit.Sym` or ModelingToolkit expression
+* `rx.rate`, a `Number`, `ModelingToolkit.Sym`, or ModelingToolkit expression
   representing the reaction rate. E.g., for a reaction like `k*X, Y --> X+Y`,
   we'd have `rate = k*X`.
 * `rx.netstoich`, a vector of pairs mapping the ModelingToolkit expression for
@@ -33,8 +88,3 @@ Each `Reaction` within `reactions(rn)` has a number of subfields. For `rx` a
 * `rx.only_use_rate`, a boolean that is `true` if the reaction was made with
   non-filled arrows and should ignore mass action kinetics. `false` by default.
 
-Empty `ReactionSystem`s can be generated via [`make_empty_network`](@ref) or
-[`@reaction_network`](@ref) with no arguments (giving one argument to the latter
-will specify a system name). `ReactionSystem`s can be programmatically extended
-using [`addspecies!`](@ref), [`addparam!`](@ref), [`addreaction!`](@ref),
-[`@add_reactions`](@ref), or composed using `merge` and `merge!`.

docs/src/tutorials/models.md

@@ -25,7 +25,7 @@ sol = solve(prob, Tsit5())
 ```
 Here, the order of unknowns in `u0` and `p` matches the order that species and
 parameters first appear within the DSL. They can also be determined by examining
-the ordering within the [`species(rn)`](@ref) and [`params(rn)`](@ref) vectors,
+the ordering within the [`species(rn)`](@ref) and [`parameters(rn)`](@ref) vectors,
 or accessed more explicitly through the [`speciesmap(rn)`](@ref) and
 [`paramsmap(rn)`](@ref) dictionaries, which map the ModelingToolkit `Term`s
 and/or `Sym`s corresponding to each species or parameter to their integer id.

docs/src/tutorials/using_catalyst.md

@@ -163,7 +163,7 @@ species(repressilator)
  P₃(t)
 ```
 ```julia
-params(repressilator)
+parameters(repressilator)
 ```
 ```julia
 7-element Array{Sym{ModelingToolkit.Parameter{Real}},1}:
@@ -198,7 +198,7 @@ instead pass `odesys` directly, provided we construct mappings from each species
 to their initial value, and each parameter to their value like:
 ```julia
 u₀map  = Pair.(species(repressilator), u₀)
-pmap   = Pair.(params(repressilator), p)
+pmap   = Pair.(parameters(repressilator), p)
 oprob2 = ODEProblem(osys, u₀map, tspan, pmap)
 ```
 `oprob` and `oprob2` are functionally equivalent, each representing the same

src/Catalyst.jl

@@ -4,27 +4,28 @@ $(DocStringExtensions.README)
 module Catalyst
 
 using DocStringExtensions
-using SparseArrays, DiffEqBase, Reexport, ModelingToolkit, DiffEqJump
+using SparseArrays, DiffEqBase, Reexport, DiffEqJump
 using Latexify, Requires
 
 # ModelingToolkit imports and convenience functions we use
+using ModelingToolkit; const MT = ModelingToolkit
+@reexport using ModelingToolkit
 using Symbolics
 using ModelingToolkit: Symbolic, value, istree, get_states, get_ps, get_iv, get_systems, 
-                       get_eqs, get_defaults, toparam
+                       get_eqs, get_defaults, toparam, get_defaults, get_observed
 import ModelingToolkit: get_variables, namespace_expr, namespace_equation, get_variables!, 
-                        modified_states!, validate
+                        modified_states!, validate, namespace_variables, namespace_parameters
 
 # internal but needed ModelingToolkit functions
 import ModelingToolkit: check_variables, check_parameters, _iszero, _merge, check_units, get_unit
 
-const DEFAULT_IV = (@parameters t)[1]
-@reexport using ModelingToolkit
+import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show
+import MacroTools, LightGraphs, AbstractAlgebra
 
-import Base: (==), merge!, merge, hash, size, getindex, setindex, isless, Sort.defalg, length, show
-import MacroTools, LightGraphs
+# globals for the modulate
 const LG = LightGraphs
-
-import AbstractAlgebra; const AA = AbstractAlgebra
+const AA = AbstractAlgebra
+const DEFAULT_IV = (@parameters t)[1]
 
 # as used in Catlab
 const USE_GV_JLL = Ref(false)
@@ -55,7 +56,7 @@ include("registered_functions.jl")
 
 # functions to query network properties
 include("networkapi.jl")
-export species, params, reactions, speciesmap, paramsmap, numspecies, numreactions, numparams
+export species, reactionparams, reactions, speciesmap, paramsmap, numspecies, numreactions, numparams
 export make_empty_network, addspecies!, addparam!, addreaction!
 export dependants, dependents, substoichmat, prodstoichmat, netstoichmat
 export conservationlaws, conservedquantities