From c449db56b4d81b8eae8d08ec426f98f4873ce3d7 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Wed, 15 Sep 2021 16:46:50 -0400
Subject: [PATCH 01/11] updating basic API for subsystems

---
 docs/src/api/catalyst_api.md            |   5 +-
 docs/src/tutorials/advanced_examples.md |   2 +-
 docs/src/tutorials/basics.md            |   2 +-
 docs/src/tutorials/generated_systems.md |   2 +-
 docs/src/tutorials/models.md            |   2 +-
 docs/src/tutorials/using_catalyst.md    |   4 +-
 src/Catalyst.jl                         |  12 ++-
 src/networkapi.jl                       | 120 ++++++++++++++++--------
 src/reactionsystem.jl                   |   5 +-
 test/api.jl                             |   2 +-
 test/reactionsystem_components.jl       |  22 ++++-
 11 files changed, 121 insertions(+), 57 deletions(-)

diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md
index 1799efe08d..5eb03a2c9c 100644
--- a/docs/src/api/catalyst_api.md
+++ b/docs/src/api/catalyst_api.md
@@ -79,12 +79,13 @@ ReactionSystem
 ```@docs
 species
 speciesmap
-params
+ModelingToolkit.parameters
+reactionparams
 paramsmap
 reactions
 numspecies
 numparams
-numreactions
+numreactionparams
 ```
 
 ## Basic Reaction Properties
diff --git a/docs/src/tutorials/advanced_examples.md b/docs/src/tutorials/advanced_examples.md
index 8d4ae57a14..0049e7c537 100644
--- a/docs/src/tutorials/advanced_examples.md
+++ b/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())
 ```
diff --git a/docs/src/tutorials/basics.md b/docs/src/tutorials/basics.md
index 064b4c0313..87dde5a59b 100644
--- a/docs/src/tutorials/basics.md
+++ b/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
diff --git a/docs/src/tutorials/generated_systems.md b/docs/src/tutorials/generated_systems.md
index 6210ff062a..59de11d6e5 100644
--- a/docs/src/tutorials/generated_systems.md
+++ b/docs/src/tutorials/generated_systems.md
@@ -6,7 +6,7 @@ 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
+* [`parameters(rn)`](@ref) 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.
diff --git a/docs/src/tutorials/models.md b/docs/src/tutorials/models.md
index 5186ce48da..ed850fc178 100644
--- a/docs/src/tutorials/models.md
+++ b/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.
diff --git a/docs/src/tutorials/using_catalyst.md b/docs/src/tutorials/using_catalyst.md
index f563bbb98c..066f25eb62 100644
--- a/docs/src/tutorials/using_catalyst.md
+++ b/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
diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index e827acf212..d962677857 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -4,15 +4,16 @@ $(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; 
 using Symbolics
 using ModelingToolkit: Symbolic, value, istree, get_states, get_ps, get_iv, get_systems, 
                        get_eqs, get_defaults, toparam
 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
@@ -22,9 +23,12 @@ const DEFAULT_IV = (@parameters t)[1]
 
 import Base: (==), merge!, merge, hash, size, getindex, setindex, isless, Sort.defalg, length, show
 import MacroTools, LightGraphs
-const LG = LightGraphs
 import Nemo: FlintZZ, matrix, nullspace, rank
 
+# convenience shorthands for packages
+const MT = ModelingToolkit
+const LG = LightGraphs
+
 # as used in Catlab
 const USE_GV_JLL = Ref(false)
 function __init__()
@@ -54,7 +58,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
diff --git a/src/networkapi.jl b/src/networkapi.jl
index 5240e5c647..8ec65511fd 100644
--- a/src/networkapi.jl
+++ b/src/networkapi.jl
@@ -1,31 +1,61 @@
 # Functions for querying network properties.
 
+# DEPRECIATED after v9.0
+function params(network)
+    Base.depwarn("`params` is deprecated, please use `ModelingToolkit.parameters` for all system and subsystem parameters, or `reactionparams` for all parameters within system and subsystem `Reaction`s.", :params, force=true)
+    parameters(network)
+end
+
+function numparams(network)
+    Base.depwarn("`numparams` is deprecated, please use `length(ModelingToolkit.parameters)` for the total number of parameters across all systems and subsystems, or `numreactionparams` for the number of parameters within system and subsystem `Reaction`s.", :params, force=true)
+    length(parameters(network))
+end
+
 ######### Accessors: #########
 
+function filter_nonrxsys(network)
+    systems   = get_systems(network)
+    rxsystems = ReactionSystem[]
+    for sys in systems
+        (sys isa ReactionSystem) && push!(rxsystems, sys)
+    end
+    rxsystems
+end
+
+
 """
     species(network)
 
-Given a [`ReactionSystem`](@ref), return a vector of species `Variable`s.
+Given a [`ReactionSystem`](@ref), return a vector of all species defined in the system
+and any subsystems that are of type `ReactionSystem`. To get the variables in the system
+and all subsystems, including non-`ReactionSystem` subsystems, uses `states(network)`.
 
 Notes:
-- If `ModelingToolkit.get_systems(network)` is not empty, may allocate. Otherwise returns
-  `ModelingToolkit.get_states(network)`.
+- If `ModelingToolkit.get_systems(network)` is non-empty will allocate.
 """
-function species(network)
-    isempty(get_systems(network)) ? get_states(network) : states(network)
+function species(network, usenamespace=false)
+    sts     = usenamespace ? namespace_variables(get_states(network)) : get_states(network)
+    systems = filter_nonrxsys(network)
+    isempty(systems) && return sts
+    unique([sts; reduce(vcat, species.(systems,true))])
 end
 
 """
-    params(network)
+    reactionparams(network)
 
-Given a [`ReactionSystem`](@ref), return a vector of parameter `Variable`s.
+Given a [`ReactionSystem`](@ref), return a vector of all parameters defined
+within the system and any subsystems that are of type `ReactionSystem`. To get
+the parameters in the system and all subsystems, including non-`ReactionSystem`
+subsystems, use `parameters(network)`.
 
 Notes:
-- If `ModelingToolkit.get_systems(network)` is not empty, may allocate. Otherwise returns
-  `ModelingToolkit.get_ps(network)`.
+- If `ModelingToolkit.get_systems(network)` is non-empty will allocate.
 """
-function params(network)
-    isempty(get_systems(network)) ? get_ps(network) : parameters(network)
+function reactionparams(network, usenamespace=false)
+    ps = usenamespace ? namespace_parameters(get_ps(network)) : get_ps(network)
+    systems = filter_nonrxsys(network)
+    isempty(systems) && return ps
+    unique([ps; reduce(vcat, reactionparams.(systems,true))])
 end
 
 """
@@ -37,15 +67,18 @@ Notes:
 - If `ModelingToolkit.get_systems(network)` is not empty, may allocate. Otherwise returns
   `ModelingToolkit.get_eqs(network)`.
 """
-function reactions(network)
-    isempty(get_systems(network)) ? get_eqs(network) : equations(network)
+function reactions(network, usenamespace=false)
+    rxs = usenamespace ? map(eq -> namespace_equation(eq,network), get_eqs(network)) : get_eqs(network)
+    systems = filter_nonrxsys(network)
+    isempty(systems) && return rxs
+    Reaction[rxs; reduce(vcat, reactions.(systems,true), init=Reaction[])]
 end
 
 """
     speciesmap(network)
 
-Given a [`ReactionSystem`](@ref), return a Dictionary mapping from species to
-species indices. (Allocates)
+Given a [`ReactionSystem`](@ref), return a Dictionary mapping species that
+participate in `Reaction`s to their index within [`species(networw)`](@ref).
 """
 function speciesmap(network)
     Dict(S => i for (i,S) in enumerate(species(network)))
@@ -54,24 +87,26 @@ end
 """
     paramsmap(network)
 
-Given a [`ReactionSystem`](@ref), return a Dictionary mapping from parameters to
-parameter indices. (Allocates)
+Given a [`ReactionSystem`](@ref), return a Dictionary mapping from parameters that
+appear within `Reaction`s to their index within [`reactionparams(network`](@ref).
 """
 function paramsmap(network)
-    Dict(p => i for (i,p) in enumerate(params(network)))
+    Dict(p => i for (i,p) in enumerate(reactionparams(network)))
 end
 
 """
     numspecies(network)
 
-Return the number of species within the given [`ReactionSystem`](@ref).
+Return the number of species within the given [`ReactionSystem`](@ref) that
+participate in `Reaction`s. 
+
+Notes
+- If there are no subsystems this will be fast.
+- As this calls [`species`](@ref), it can be slow and will allocate if there are
+  any subsystems. 
 """
 function numspecies(network)
-    ns = length(get_states(network))
-    for sys in get_systems(network)
-        ns += numspecies(sys)
-    end
-    ns
+    length(species(network))
 end
 
 """
@@ -82,22 +117,24 @@ Return the number of reactions within the given [`ReactionSystem`](@ref).
 function numreactions(network)
     nr = length(get_eqs(network))
     for sys in get_systems(network)
-        nr += numreactions(sys)
+        (sys isa ReactionSystem) && (nr += numreactions(sys))
     end
     nr
 end
 
 """
-    numparams(network)
+    numreactionparams(network)
 
-Return the number of parameters within the given [`ReactionSystem`](@ref).
+Return the number of parameters within the given [`ReactionSystem`](@ref) that
+participate in `Reaction`s.
+
+Notes
+- If there are no subsystems this will be fast.
+- As this calls [`reactionparams`](@ref), it can be slow and will allocate
+  if there are any subsystems. 
 """
-function numparams(network)
-    np = length(get_ps(network))
-    for sys in get_systems(network)
-        np += numparams(sys)
-    end
-    np
+function numreactionparams(network)
+    length(reactionparams(network))
 end
 
 """
@@ -175,6 +212,7 @@ function substoichmat(::Type{Matrix{Int}},rn::ReactionSystem; smap=speciesmap(rn
     smat
 end
 function substoichmat(rn::ReactionSystem; sparse::Bool=false, smap=speciesmap(rn))
+    isempty(get_systems(rn)) || error("substoichmat does not currently support subsystems.")
 	sparse ? substoichmat(SparseMatrixCSC{Int,Int}, rn; smap=smap) : substoichmat(Matrix{Int}, rn; smap=smap)
 end
 
@@ -211,6 +249,7 @@ function prodstoichmat(::Type{Matrix{Int}},rn::ReactionSystem; smap=speciesmap(r
     pmat
 end
 function prodstoichmat(rn::ReactionSystem; sparse=false, smap=speciesmap(rn))
+    isempty(get_systems(rn)) || error("prodstoichmat does not currently support subsystems.")
 	sparse ? prodstoichmat(SparseMatrixCSC{Int,Int}, rn; smap=smap) : prodstoichmat(Matrix{Int}, rn; smap=smap)
 end
 
@@ -245,6 +284,7 @@ function netstoichmat(::Type{Matrix{Int}},rn::ReactionSystem; smap=speciesmap(rn
     nmat
 end
 function netstoichmat(rn::ReactionSystem; sparse=false, smap=speciesmap(rn))
+    isempty(get_systems(rn)) || error("netstoichmat does not currently support subsystems.")
 	sparse ? netstoichmat(SparseMatrixCSC{Int,Int}, rn; smap=smap) : netstoichmat(Matrix{Int}, rn; smap=smap)
 end
 
@@ -305,6 +345,7 @@ Notes:
   `reactionidx`.
 """
 function reactioncomplexmap(rn::ReactionSystem; smap=speciesmap(rn))
+    isempty(get_systems(rn)) || error("reactioncomplexmap does not currently support subsystems.")
     rxs = reactions(rn)
     numreactions(rn) > 0 || error("There must be at least one reaction to find reaction complexes.")
     complextorxsmap = OrderedDict{ReactionComplex{eltype(rxs[1].substoich)},Vector{Pair{Int,Int}}}()
@@ -378,6 +419,7 @@ B_{i j} = \begin{cases}
 """
 function reactioncomplexes(rn::ReactionSystem; sparse=false, smap=speciesmap(rn), 
                            complextorxsmap=reactioncomplexmap(rn; smap=smap))
+    isempty(get_systems(rn)) || error("reactioncomplexes does not currently support subsystems.")
 	sparse ? reactioncomplexes(SparseMatrixCSC{Int,Int}, rn, smap, complextorxsmap) :
              reactioncomplexes(Matrix{Int}, rn, smap, complextorxsmap)
 end
@@ -418,7 +460,8 @@ Notes:
 - Set sparse=true for a sparse matrix representation
 """
 function complexstoichmat(rn::ReactionSystem; sparse=false, rcs=keys(reactioncomplexmap(rn)))
-	sparse ? complexstoichmat(SparseMatrixCSC{Int,Int}, rn, rcs) : 
+    isempty(get_systems(rn)) || error("complexstoichmat does not currently support subsystems.")
+    sparse ? complexstoichmat(SparseMatrixCSC{Int,Int}, rn, rcs) : 
              complexstoichmat(Matrix{Int}, rn, rcs)
 end
 
@@ -468,7 +511,8 @@ Notes:
 - Set sparse=true for a sparse matrix representation
 """
 function complexoutgoingmat(rn::ReactionSystem; sparse=false, B=reactioncomplexes(rn,sparse=sparse)[2])
-	sparse ? complexoutgoingmat(SparseMatrixCSC{Int,Int}, rn, B) : 
+    isempty(get_systems(rn)) || error("complexoutgoingmat does not currently support subsystems.")
+    sparse ? complexoutgoingmat(SparseMatrixCSC{Int,Int}, rn, B) : 
              complexoutgoingmat(Matrix{Int}, rn, B)
 end
 
@@ -607,7 +651,7 @@ function subnetworks(rs::ReactionSystem, lcs::AbstractVector;
                   rxs = reactions(rs),
                   complextorxmap = [map(first,rcmap) for rcmap in values(reactioncomplexmap(rs))],
                   p = parameters(rs))
-
+    isempty(get_systems(rs)) || error("subnetworks does not currently support subsystems.")
     t = get_iv(rs)
     subreac = Vector{ReactionSystem}()
     for i in 1:length(lcs)
@@ -710,7 +754,7 @@ Notes:
 """
 function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
     issetequal(species(rn1), species(rn2)) || return false
-    issetequal(params(rn1), params(rn2)) || return false
+    issetequal(parameters(rn1), parameters(rn2)) || return false
     isequal(get_iv(rn1), get_iv(rn2)) || return false
     (numreactions(rn1) == numreactions(rn2)) || return false
 
@@ -896,7 +940,7 @@ the same units).
 
 """
 function validate(rx::Reaction; info::String = "")     
-    validated = ModelingToolkit._validate([rx.rate], [string(rx, ": rate")], info = info)
+    validated = MT._validate([rx.rate], [string(rx, ": rate")], info = info)
     
     subunits = isempty(rx.substrates) ? nothing : get_unit(rx.substrates[1])
     for i in 2:length(rx.substrates)
diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
index 7c1f23edaa..ec8cdd2720 100644
--- a/src/reactionsystem.jl
+++ b/src/reactionsystem.jl
@@ -203,7 +203,7 @@ struct ReactionSystem <: ModelingToolkit.AbstractTimeDependentSystem
     end
 end
 
-function ReactionSystem(eqs, iv, species, params;
+function ReactionSystem(eqs, iv, species, ps;
                         observed = [],
                         systems = [],
                         name = nothing,
@@ -214,7 +214,8 @@ function ReactionSystem(eqs, iv, species, params;
                         checks = true)
     name === nothing && throw(ArgumentError("The `name` keyword must be provided. Please consider using the `@named` macro"))
 
-    ReactionSystem(eqs, iv, species, params, observed, name, systems, defaults, connection_type, checks = checks)
+    ReactionSystem(eqs, iv, species, ps, observed, name, systems, defaults, connection_type, 
+                   checks = checks)
 end
 
 function ReactionSystem(iv; kwargs...)
diff --git a/test/api.jl b/test/api.jl
index 625963db21..a2e6bf3757 100644
--- a/test/api.jl
+++ b/test/api.jl
@@ -27,7 +27,7 @@ addspecies!(rs3, D)
 addparam!(rs3, k3)
 addparam!(rs3, k4)
 @test issetequal(species(rs3), [S, D])
-@test issetequal(params(rs3), [k3, k4])
+@test issetequal(parameters(rs3), [k3, k4])
 addreaction!(rs3, Reaction(k3, [S], [D]))
 addreaction!(rs3, Reaction(k4, [S,I], [D]))
 merge!(rs, rs3)
diff --git a/test/reactionsystem_components.jl b/test/reactionsystem_components.jl
index ea41754029..870fc6225a 100644
--- a/test/reactionsystem_components.jl
+++ b/test/reactionsystem_components.jl
@@ -101,9 +101,14 @@ sol = solve(nlprob, NewtonRaphson(), tol=1e-9)
 # adding algebraic constraints
 @parameters t, r₊, r₋, β
 @variables A(t), B(t), C(t), D(t)
-rxs = [Reaction(r₊, [A,B], [C]),
-       Reaction(r₋, [C], [A,B])]
-@named rs = ReactionSystem(rxs, t, [A,B,C], [r₊, r₋])
+rxs1 = [Reaction(r₊, [A,B], [C])]
+rxs2 = [Reaction(r₋, [C], [A,B])]
+@named rs1 = ReactionSystem(rxs1, t, [A,B,C], [r₊])
+@named rs2 = ReactionSystem(rxs2, t, [A,B,C], [r₋])
+@named rs  = extend(rs1, rs2)
+@test issetequal(states(rs), [A,B,C])
+@test issetequal(parameters(rs), [r₊,r₋])
+@test issetequal(equations(rs), union(rxs1,rxs2))
 A2 = ModelingToolkit.ParentScope(A)
 B2 = ModelingToolkit.ParentScope(B)
 nseqs = [D ~ 2*A2 + β*B2]
@@ -114,4 +119,13 @@ p  = [r₊ => 1.0, r₋ => 2.0, ns.β => 3.0]
 u₀ = [A => 1.0, B => 2.0, C => 0.0]
 oprob = ODEProblem(structural_simplify(osys), u₀, (0.0,10.0), p)
 sol = solve(oprob, Tsit5())
-@test isapprox(0, norm(sol[ns.D] .- 2*sol[A] - 3*sol[B]), atol=(100*eps())) 
\ No newline at end of file
+@test isapprox(0, norm(sol[ns.D] .- 2*sol[A] - 3*sol[B]), atol=(100*eps())) 
+
+# test API functions for composed model
+@test issetequal(species(rs), [A,B,C])
+@test issetequal(states(rs), [A,B,C,ns.D])
+@test issetequal(reactionparams(rs), [r₊,r₋])
+@test issetequal(parameters(rs), [r₊,r₋,ns.β])
+@test issetequal(reactions(rs), union(rxs1,rxs2))
+@test issetequal(filter(eq -> eq isa Reaction, equations(rs)), union(rxs1,rxs2))
+@test issetequal(filter(eq -> eq isa Equation, equations(rs)), [ModelingToolkit.namespace_equation(nseqs[1],ns)])
\ No newline at end of file

From 96cca796ce3a931788121639357b261db5172d2c Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Wed, 15 Sep 2021 16:50:25 -0400
Subject: [PATCH 02/11] more API updates

---
 src/networkapi.jl | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/src/networkapi.jl b/src/networkapi.jl
index 8ec65511fd..2d66f3fef0 100644
--- a/src/networkapi.jl
+++ b/src/networkapi.jl
@@ -736,7 +736,8 @@ function (==)(rn1::Reaction, rn2::Reaction)
     isequal(rn1.rate, rn2.rate) || return false
     issetequal(zip(rn1.substrates,rn1.substoich), zip(rn2.substrates,rn2.substoich)) || return false
     issetequal(zip(rn1.products,rn1.prodstoich), zip(rn2.products,rn2.prodstoich)) || return false
-    issetequal(rn1.netstoich, rn2.netstoich)
+    issetequal(rn1.netstoich, rn2.netstoich) || return false
+    rn1.only_use_rate == rn2.only_use_rate
 end
 
 

From b85c79296cda5438d41841c493c161f8246b39f8 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Wed, 15 Sep 2021 21:11:18 -0400
Subject: [PATCH 03/11] update API functions for subsystems

---
 HISTORY.md                              |  5 ++
 docs/src/api/catalyst_api.md            |  4 +-
 docs/src/tutorials/generated_systems.md |  3 +-
 src/Catalyst.jl                         |  2 +-
 src/networkapi.jl                       | 84 +++++++++++++------------
 test/api.jl                             |  2 +-
 6 files changed, 55 insertions(+), 45 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index 3b166dbd62..257b3dc593 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -10,6 +10,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
diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md
index 5eb03a2c9c..f9e3f65836 100644
--- a/docs/src/api/catalyst_api.md
+++ b/docs/src/api/catalyst_api.md
@@ -112,8 +112,8 @@ reactionrates
 addspecies!
 addparam!
 addreaction!
-merge!(network1::ReactionSystem, network2::ReactionSystem)
-merge(network1::ReactionSystem, network2::ReactionSystem)
+ModelingToolkit.extend
+ModelingToolkit.compose
 ```
 
 ## Network Analysis and Representations
diff --git a/docs/src/tutorials/generated_systems.md b/docs/src/tutorials/generated_systems.md
index 59de11d6e5..c889007ed1 100644
--- a/docs/src/tutorials/generated_systems.md
+++ b/docs/src/tutorials/generated_systems.md
@@ -37,4 +37,5 @@ 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!`.
+[`@add_reactions`](@ref), or composed using [`ModelingToolkit.extend`](@ref) and
+`ModelingToolkit.compose`](@ref).
diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index d962677857..e23a077f0b 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -21,7 +21,7 @@ import ModelingToolkit: check_variables, check_parameters, _iszero, _merge, chec
 const DEFAULT_IV = (@parameters t)[1]
 @reexport using ModelingToolkit
 
-import Base: (==), merge!, merge, hash, size, getindex, setindex, isless, Sort.defalg, length, show
+import Base: (==), hash, size, getindex, setindex, isless, Sort.defalg, length, show
 import MacroTools, LightGraphs
 import Nemo: FlintZZ, matrix, nullspace, rank
 
diff --git a/src/networkapi.jl b/src/networkapi.jl
index 2d66f3fef0..78f2332e9a 100644
--- a/src/networkapi.jl
+++ b/src/networkapi.jl
@@ -2,15 +2,36 @@
 
 # DEPRECIATED after v9.0
 function params(network)
-    Base.depwarn("`params` is deprecated, please use `ModelingToolkit.parameters` for all system and subsystem parameters, or `reactionparams` for all parameters within system and subsystem `Reaction`s.", :params, force=true)
+    Base.depwarn("`params` is depreciated, please use `ModelingToolkit.parameters` for all system and subsystem parameters, or `reactionparams` for all parameters within system and subsystem `Reaction`s.", :params, force=true)
     parameters(network)
 end
 
 function numparams(network)
-    Base.depwarn("`numparams` is deprecated, please use `length(ModelingToolkit.parameters)` for the total number of parameters across all systems and subsystems, or `numreactionparams` for the number of parameters within system and subsystem `Reaction`s.", :params, force=true)
+    Base.depwarn("`numparams` is depreciated, please use `length(ModelingToolkit.parameters)` for the total number of parameters across all systems and subsystems, or `numreactionparams` for the number of parameters within system and subsystem `Reaction`s.", :params, force=true)
     length(parameters(network))
 end
 
+"""
+    merge(network1::ReactionSystem, network2::ReactionSystem)
+
+Create a new network merging `network1` and `network2`.
+
+Notes:
+- Duplicate reactions between the two networks are not filtered out.
+- [`Reaction`](@ref)s are not deepcopied to minimize allocations, so the new
+  network will share underlying data arrays.
+- Subsystems are not deepcopied between the two networks and will hence be
+  shared.
+- Returns the merged network.
+"""
+function Base.merge(network1::ReactionSystem, network2::ReactionSystem)
+    Base.depwarn("`merge(sys1::ReactionSystem, sys2::ReactionNetwork)` is depreciated, please use `ModelingToolkit.extend` instead.", :merge, force=true)
+    network = make_empty_network()
+    merge!(network, network1)
+    merge!(network, network2)
+    network
+end
+
 ######### Accessors: #########
 
 function filter_nonrxsys(network)
@@ -745,28 +766,29 @@ end
     ==(rn1::ReactionSystem, rn2::ReactionSystem)
 
 Tests whether the underlying species, parameters and reactions are the same in
-the two [`ReactionSystem`](@ref)s. Ignores order network components were
-defined.
+the two [`ReactionSystem`](@ref)s. 
 
 Notes:
 - *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
     considered different than `(A+1)^2`.
-- Flattens subsystems, and hence may allocate, when checking equality.
+- Does not include `defaults` in determining equality.
 """
 function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
-    issetequal(species(rn1), species(rn2)) || return false
-    issetequal(parameters(rn1), parameters(rn2)) || return false
+    (nameof(rn1) == nameof(rn2)) || return false
     isequal(get_iv(rn1), get_iv(rn2)) || return false
-    (numreactions(rn1) == numreactions(rn2)) || return false
-
-    # the following fails for some reason, so need to use issubset
-    #issetequal(equations(rn1), equations(rn2)) || return false
-    (issubset(reactions(rn1),reactions(rn2)) && issubset(reactions(rn2),reactions(rn1))) || return false
+    issetequal(get_states(rn1), get_states(rn2)) || return false
+    issetequal(get_ps(rn1), get_ps(rn2)) || return false
+    issetequal(MT.get_observed(rn1), MT.get_observed(rn2))
+  
+    # reactions
+    # issetequal fails for some reason, so need to use issubset
+    (length(get_eqs(rn1)) == length(get_eqs(rn2))) || return false
+    issubset(get_eqs(rn1),get_eqs(rn2)) && issubset(get_eqs(rn2),get_eqs(rn1)) || return false
+
+    # subsystems
+    (length(get_systems(rn1)) == length(get_systems(rn2))) || return false
+    issetequal(get_systems(rn1), get_systems(rn2))
 
-    # BELOW SHOULD NOT BE NEEDED as species, params and equations flatten
-    #issetequal(rn1.systems, rn2.systems) || return false
-    # sys1 = rn1.systems; sys2 = rn2.systems
-    # (issubset(sys1,sys2) && issubset(sys2,sys1)) || return false
     true
 end
 
@@ -882,7 +904,6 @@ function addreaction!(network::ReactionSystem, rx::Reaction)
     length(get_eqs(network))
 end
 
-
 """
     merge!(network1::ReactionSystem, network2::ReactionSystem)
 
@@ -895,39 +916,19 @@ Notes:
 - Subsystems are not deepcopied between the two networks and will hence be
   shared.
 - Returns `network1`.
-- Does not currently handle pins.
 """
 function Base.merge!(network1::ReactionSystem, network2::ReactionSystem)
     isequal(get_iv(network1), get_iv(network2)) || 
         error("Reaction networks must have the same independent variable to be mergable.")
+    append!(get_eqs(network1), get_eqs(network2))
     union!(get_states(network1), get_states(network2))
     union!(get_ps(network1), get_ps(network2))
-    append!(get_eqs(network1), get_eqs(network2))
+    append!(get_observed(network1), get_observed(network2))    
     append!(get_systems(network1), get_systems(network2))
+    merge!(get_defaults(network1), get_defaults(network2))
     network1
 end
 
-"""
-    merge(network1::ReactionSystem, network2::ReactionSystem)
-
-Create a new network merging `network1` and `network2`.
-
-Notes:
-- Duplicate reactions between the two networks are not filtered out.
-- [`Reaction`](@ref)s are not deepcopied to minimize allocations, so the new
-  network will share underlying data arrays.
-- Subsystems are not deepcopied between the two networks and will hence be
-  shared.
-- Returns the merged network.
-- Does not currently handle pins.
-"""
-function Base.merge(network1::ReactionSystem, network2::ReactionSystem)
-    network = make_empty_network()
-    merge!(network, network1)
-    merge!(network, network2)
-    network
-end
-
 
 ###############################   units   #####################################
 
@@ -973,6 +974,9 @@ end
 Check that all species in the [`ReactionSystem`](@ref) have the same units, and
 that the rate laws of all reactions reduce to units of (species units) / (time
 units).
+
+Notes:
+- Does not check subsystems too. 
 """
 function validate(rs::ReactionSystem, info::String="")
     specs = get_states(rs)
diff --git a/test/api.jl b/test/api.jl
index a2e6bf3757..3e6b5e4f0c 100644
--- a/test/api.jl
+++ b/test/api.jl
@@ -49,7 +49,7 @@ addparam!(rs3, k3)
 addparam!(rs3, k4)
 addreaction!(rs3, Reaction(k3, [S], [D]))
 addreaction!(rs3, Reaction(k4, [S,I], [D]))
-rs4 = merge(rs, rs3)
+rs4 = extend(rs, rs3)
 @test rs2 == rs4
 
 rxs = [Reaction(k1*S, [S,I], [I], [2,3], [2]),

From 12211e567c9afbd5183800a344f59a11bd46e85e Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Wed, 15 Sep 2021 21:21:19 -0400
Subject: [PATCH 04/11] update imports

---
 src/Catalyst.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index 891af0df6f..cba1cfa36c 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -11,7 +11,7 @@ using Latexify, Requires
 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, namespace_variables, namespace_parameters
 

From 945cc668f2488d13ff030ee1799878eea93cffa2 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Wed, 15 Sep 2021 21:47:29 -0400
Subject: [PATCH 05/11] fix make model tests

---
 src/networkapi.jl          |  2 +-
 test/model_modification.jl | 18 ++++----
 test/test_networks.jl      | 88 +++++++++++++++++++-------------------
 3 files changed, 54 insertions(+), 54 deletions(-)

diff --git a/src/networkapi.jl b/src/networkapi.jl
index 28bede683b..b5c3ed3539 100644
--- a/src/networkapi.jl
+++ b/src/networkapi.jl
@@ -816,7 +816,7 @@ function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
 
     # subsystems
     (length(get_systems(rn1)) == length(get_systems(rn2))) || return false
-    issetequal(get_systems(rn1), get_systems(rn2))
+    issetequal(get_systems(rn1), get_systems(rn2)) || return false
 
     true
 end
diff --git a/test/model_modification.jl b/test/model_modification.jl
index 13cc042812..064f002221 100644
--- a/test/model_modification.jl
+++ b/test/model_modification.jl
@@ -52,7 +52,7 @@ end k0 k5 k6 k7 k8
 ### Compares test network to identical network constructed via @add_reactions ###
 identical_networks = Vector{Pair}()
 
-step_by_step_network_1 = @reaction_network begin
+step_by_step_network_1 = @reaction_network rns5 begin
     p, ∅ → X1
     (k1,k2), X1 ⟷ X2
     (k3,k4), X2 ⟷ X3
@@ -65,7 +65,7 @@ end k5 k6
 end d
 push!(identical_networks, reaction_networks_standard[5] => step_by_step_network_1)
 
-step_by_step_network_2 = @reaction_network begin
+step_by_step_network_2 = @reaction_network rns7 begin
     (p1,p2,p3), ∅ → (X1,X2,X3)
 end p1 p2 p3
 @add_reactions step_by_step_network_2 begin
@@ -75,7 +75,7 @@ end p1 p2 p3
 end k1 k2 k3 v1 K1 d1 d2 d3 d4 d5
 push!(identical_networks, reaction_networks_standard[7] => step_by_step_network_2)
 
-step_by_step_network_3 = @reaction_network begin
+step_by_step_network_3 = @reaction_network rns10 begin
     p, ∅ ⟶ X1
     (k1,k2), X1 → X2
 end p k1 k2 k3 k4
@@ -89,7 +89,7 @@ end p k1 k2 k3 k4 k5 k6
 end p k1 k2 k3 k4 k5 k6 k7 k8 d
 push!(identical_networks, reaction_networks_standard[10] => step_by_step_network_3)
 
-step_by_step_network_4 = @reaction_network begin
+step_by_step_network_4 = @reaction_network rnh7 begin
     v/10 + hill(X1,v,K,n), ∅ → X1 + X2
 end v K n k1 k2 k3 d
 @add_reactions step_by_step_network_4 begin
@@ -102,7 +102,7 @@ end k1 k2
 end k2 k3 d
 push!(identical_networks, reaction_networks_hill[7] => step_by_step_network_4)
 
-step_by_step_network_5 = @reaction_network begin
+step_by_step_network_5 = @reaction_network rnc1 begin
     (k1,k2), X1 ↔ X2
     (k3,k4), X2 ↔ X3
 end k1 k2 k3 k4 k5 k6
@@ -114,7 +114,7 @@ end k5
 end k6
 push!(identical_networks, reaction_networks_constraint[1] => step_by_step_network_5)
 
-step_by_step_network_6 = @reaction_network begin
+step_by_step_network_6 = @reaction_network rnc5 begin
     (k1,k2), X1 ↔ 2X2
 end k1 k2 k3 k4 k5 k6
 @add_reactions step_by_step_network_6 begin
@@ -125,7 +125,7 @@ end k3 k4
 end k5 k6
 push!(identical_networks, reaction_networks_constraint[5] => step_by_step_network_6)
 
-step_by_step_network_7 = @reaction_network begin
+step_by_step_network_7 = @reaction_network rnr3 begin
     k2p, Y → 0
 end k1 k2p k2pp k3p k3pp A J3 k4 m J4
 @add_reactions step_by_step_network_7 begin
@@ -140,7 +140,7 @@ end k1
 end k4 m J4
 push!(identical_networks, reaction_networks_real[3] => step_by_step_network_7)
 
-step_by_step_network_8 = @reaction_network 
+step_by_step_network_8 = @reaction_network rnw7
 @parameters k1
 addparam!(step_by_step_network_8,k1)
 @add_reactions step_by_step_network_8 begin
@@ -153,7 +153,7 @@ end k1
 end k2 k3
 push!(identical_networks, reaction_networks_weird[7] => step_by_step_network_8)
 
-step_by_step_network_9 = @reaction_network
+step_by_step_network_9 = @reaction_network rnw10
 @add_reactions step_by_step_network_9 begin
     d, 5X1 → 4X1
 end d
diff --git a/test/test_networks.jl b/test/test_networks.jl
index 7063130026..741857f8a6 100644
--- a/test/test_networks.jl
+++ b/test/test_networks.jl
@@ -10,20 +10,20 @@ reaction_networks_weird = Vector{ReactionSystem}(undef,10)
 
 
 ### Standard reaction networks. ###
-reaction_networks_standard[1] = @reaction_network begin
+reaction_networks_standard[1] = @reaction_network rns1 begin
     (p1,p2,p3), ∅ → (X1,X2,X3)
     (k1,k2), X2 ⟷ X1 + 2X3
     (k3,k4), X1 ⟷ X3
     (d1,d2,d3), (X1,X2,X3) → ∅
 end p1 p2 p3 k1 k2 k3 k4 d1 d2 d3
 
-reaction_networks_standard[2] = @reaction_network begin
+reaction_networks_standard[2] = @reaction_network rns2 begin
     mmr(X2,v1,K1), ∅ → X1
     mm(X1,v2,K2), ∅ → X2
     d, X1+X2 → ∅
 end v1 K1 v2 K2 d
 
-reaction_networks_standard[3] = @reaction_network begin
+reaction_networks_standard[3] = @reaction_network rns3 begin
     mm(X2,v1,K1), ∅ → X1
     mm(X3,v2,K2), ∅ → X2
     (k1,k2), X1 ⟷ X3
@@ -31,7 +31,7 @@ reaction_networks_standard[3] = @reaction_network begin
     d, (X1,X2,X3,X4) → ∅
 end v1 K1 v2 K2 k1 k2 k3 k4 d
 
-reaction_networks_standard[4] = @reaction_network begin
+reaction_networks_standard[4] = @reaction_network rns4 begin
     mmr(X4,v1,K1), ∅ → X1
     mmr(X1,v2,K2), ∅ → X2
     mmr(X2,v3,K3), ∅ → X3
@@ -39,7 +39,7 @@ reaction_networks_standard[4] = @reaction_network begin
     (d1,d2,d3,d4), (X1,X2,X3,X4) → ∅
 end v1 K1 v2 K2 v3 K3 v4 K4 d1 d2 d3 d4
 
-reaction_networks_standard[5] = @reaction_network begin
+reaction_networks_standard[5] = @reaction_network rns5 begin
     p, ∅ → X1
     (k1,k2), X1 ⟷ X2
     (k3,k4), X2 ⟷ X3
@@ -47,7 +47,7 @@ reaction_networks_standard[5] = @reaction_network begin
     d, X4 → ∅
 end p k1 k2 k3 k4 k5 k6 d
 
-reaction_networks_standard[6] = @reaction_network begin
+reaction_networks_standard[6] = @reaction_network rns6 begin
     (p1,p2), ∅ → (X1,X2)
     (k1,k2), 2X1 ⟷ X3
     (k3,k4), X2 + X3 ⟷ 4X4
@@ -55,21 +55,21 @@ reaction_networks_standard[6] = @reaction_network begin
     d, (X1,X2,X3,X4) → ∅
 end p1 p2 k1 k2 k3 k4 k5 k6 d
 
-reaction_networks_standard[7] = @reaction_network begin
+reaction_networks_standard[7] = @reaction_network rns7 begin
     (p1,p2,p3), ∅ → (X1,X2,X3)
     (k1,k2), X1 + 2X2 ⟷ X4
     (mm(X3,v1,K1),k3), X4 ⟷ X5
     (d1,d2,d3,d4,d5), (X1,X2,X3,X4,X5) → ∅
 end p1 p2 p3 k1 k2 k3 v1 K1 d1 d2 d3 d4 d5
 
-reaction_networks_standard[8] = @reaction_network begin
+reaction_networks_standard[8] = @reaction_network rns8 begin
     p, ∅ → 2X1
     k1, X1 → X2
     (k2, k3), X2 → X3
     d, X3 → ∅
 end p k1 k2 k3 d
 
-reaction_networks_standard[9] = @reaction_network begin
+reaction_networks_standard[9] = @reaction_network rns9 begin
     (p1,p2,p3), ∅ ⟶ (X1,X2,X3)
     (d1,d2,d3), (X1,X2,X3) ⟶ ∅
     (k1,k2), X1 + X2 ⟷ X3
@@ -77,7 +77,7 @@ reaction_networks_standard[9] = @reaction_network begin
     (k5,k6), X4 ⟷ X1 + X2
 end p1 p2 p3 d1 d2 d3 k1 k2 k3 k4 k5 k6
 
-reaction_networks_standard[10] = @reaction_network begin
+reaction_networks_standard[10] = @reaction_network rns10 begin
     p, ∅ ⟶ X1
     (k1,k2), X1 → X2
     (k3,k4), X2 → X3
@@ -90,33 +90,33 @@ end p k1 k2 k3 k4 k5 k6 k7 k8 d
 
 ### Network with hill functions ###.
 
-reaction_networks_hill[1] = @reaction_network begin
+reaction_networks_hill[1] = @reaction_network rnh1 begin
     hillr(X2,v1,K1,n1),   ∅ → X1
     hillr(X1,v2,K2,n2),   ∅ → X2
     (d1,d2),          (X1,X2) → ∅
 end v1 v2 K1 K2 n1 n2 d1 d2
 
-reaction_networks_hill[2] = @reaction_network begin
+reaction_networks_hill[2] = @reaction_network rnh2 begin
     hillr(X3,v1,K1,n1),   ∅ → X1
     hillr(X1,v2,K2,n2),   ∅ → X2
     hillr(X2,v3,K3,n3),   ∅ → X3
     (d1,d2,d3),          (X1,X2,X3) → ∅
 end v1 v2 v3 K1 K2 K3 n1 n2 n3 d1 d2 d3
 
-reaction_networks_hill[3] = @reaction_network begin
+reaction_networks_hill[3] = @reaction_network rnh3 begin
     hillr(X2,v1,K1,n1), ∅ → X1
     hill(X1,v2,K2,n2), ∅ → X2
     d, X1+X2 → ∅
 end v1 K1 n1 v2 K2 n2 d
 
-reaction_networks_hill[4] = @reaction_network begin
+reaction_networks_hill[4] = @reaction_network rnh4 begin
     hillr(X2,v1,K1,n1)*hillr(X3,v1,K1,n1), ∅ → X1
     hillr(X1,v2,K2,n2)*hillr(X3,v2,K2,n2), ∅ → X2
     hillr(X1,v3,K3,n3)*hillr(X2,v3,K3,n3), ∅ → X3
     (d1,d2,d3), (X1,X2,X3)  ⟶ ∅
 end v1 K1 n1 v2 K2 n2 v3 K3 n3 d1 d2 d3
 
-reaction_networks_hill[5] = @reaction_network begin
+reaction_networks_hill[5] = @reaction_network rnh5 begin
     hillr(X2,v1,K1,n1)*hill(X4,v1,K1,n1), ∅ → X1
     hill(X5,v2,K2,n2), ∅ → X2
     hill(X3,v3,K3,n3), ∅ → X3
@@ -128,31 +128,31 @@ reaction_networks_hill[5] = @reaction_network begin
     (d1,d2,d3,d4,d5), (X1,X2,X3,X4,X5)  ⟶ ∅
 end v1 K1 n1 v2 K2 n2 v3 K3 n3 v4 K4 n4 v5 K5 n5 k1 k2 k3 k4 k5 k6 d1 d2 d3 d4 d5
 
-reaction_networks_hill[6] = @reaction_network begin
+reaction_networks_hill[6] = @reaction_network rnh6 begin
     v/10+hill(X1,v,K,n), ∅ → X1
     d, X1 → ∅
 end v K n d
 
-reaction_networks_hill[7] = @reaction_network begin
+reaction_networks_hill[7] = @reaction_network rnh7 begin
     v/10 + hill(X1,v,K,n), ∅ → X1 + X2
     (k1,k2), X1 + X2 ↔ X3
     k3, X3 → X1
     d, (X1,X2,X3) → ∅
 end v K n k1 k2 k3 d
 
-reaction_networks_hill[8] = @reaction_network begin
+reaction_networks_hill[8] = @reaction_network rnh8 begin
     hill(X2,v1,K1,n1), ∅ → X1
     hillr(X1,v2,K2,n2)*hill(X3,v3,K3,n3), ∅ → X2
     hill(X2,v4,K4,n4), ∅ → X3
     (d1,d2,d3), (X1,X2,X3) → ∅
 end v1 K1 n1 v2 K2 n2 v3 K3 n3 v4 K4 n4 d1 d2 d3
 
-reaction_networks_hill[9] = @reaction_network begin
+reaction_networks_hill[9] = @reaction_network rnh9 begin
     hill(X1,v1,K1,n1)*hillr(X1,v2,K2,n2), ∅ → X1
     d, X1 → ∅
 end v1 K1 n1 v2 K2 n2 d
 
-reaction_networks_hill[10] = @reaction_network begin
+reaction_networks_hill[10] = @reaction_network rnh10 begin
     hill(X2,v1,K1,n1), ∅ → X1
     hillr(X4,v2,K2,n2), ∅ → X2
     (k1,k2), 2X1 + X2 ⟷ X3
@@ -164,68 +164,68 @@ end v1 K1 n1 v2 K2 n2 k1 k2 k3 k4 k5 k6 d1 d2
 
 ### Reaction networks were some linnear combination concentrations remain fixed (steady state values depends on initial conditions).
 
-reaction_networks_constraint[1] = @reaction_network begin
+reaction_networks_constraint[1] = @reaction_network rnc1 begin
     (k1,k2), X1 ↔ X2
     (k3,k4), X2 ↔ X3
     (k5,k6), X3 ↔ X1
 end k1 k2 k3 k4 k5 k6
 reaction_network_constraints[1] = [1 1 1]
 
-reaction_networks_constraint[2] = @reaction_network begin
+reaction_networks_constraint[2] = @reaction_network rnc2 begin
     (k1,k2), X1 ↔ 2X1
     (k3,k4), X1 + X2 ↔ X3
     (k5,k6), X3 ↔ X2
 end k1 k2 k3 k4 k5 k6
 reaction_network_constraints[2] = [0 1 1]
 
-reaction_networks_constraint[3] = @reaction_network begin
+reaction_networks_constraint[3] = @reaction_network rnc3 begin
     (k1,k2*X5), X1 ↔ X2
     (k3*X5,k4), X3 ↔ X4
     (p+k5*X2*X3,d), ∅ ↔ X5
 end k1 k2 k3 k4 p k5 d
 reaction_network_constraints[3] = [0 0 1 1 0; 1 1 0 0 0]
 
-reaction_networks_constraint[4] = @reaction_network begin
+reaction_networks_constraint[4] = @reaction_network rnc4 begin
     (k1,k2), X1 + X2 ↔ X3
     (mm(X3,v,K),d), ∅ ↔ X4
 end k1 k2 v K d
 reaction_network_constraints[4] = [0 1 1 0; -1 1 0 0]
 
-reaction_networks_constraint[5] = @reaction_network begin
+reaction_networks_constraint[5] = @reaction_network rnc5 begin
     (k1,k2), X1 ↔ 2X2
     (k3,k4), 2X2 ↔ 3X3
     (k5,k6), 3X3 ↔ 4X4
 end k1 k2 k3 k4 k5 k6
 reaction_network_constraints[5] = [12 6 4 3]
 
-reaction_networks_constraint[6] = @reaction_network begin
+reaction_networks_constraint[6] = @reaction_network rnc6 begin
     mmr(X1,v1,K1), X1 → X2
     mmr(X2,v2,K2), X2 → X3
     mmr(X3,v3,K3), X3 → X1
 end v1 K1 v2 K2 v3 K3
 reaction_network_constraints[6] = [1 1 1]
 
-reaction_networks_constraint[7] = @reaction_network begin
+reaction_networks_constraint[7] = @reaction_network rnc7 begin
     (k1,k2), X1 + X2 ↔ X3
     (mm(X3,v,K),d), ∅ ↔ X2
     (k3,k4), X2 ↔ X4
 end k1 k2 k3 k4 v K d
 reaction_network_constraints[7] = [1 0 1 0]
 
-reaction_networks_constraint[8] = @reaction_network begin
+reaction_networks_constraint[8] = @reaction_network rnc8 begin
     (k1,k2), X1 + X2 ↔ X3
     (mm(X3,v1,K1),mm(X4,v2,K2)), X3 ↔ X4
 end k1 k2 v1 K1 v2 K2
 reaction_network_constraints[8] = [-1 1 0 0; 0 1 1 1]
 
-reaction_networks_constraint[9] = @reaction_network begin
+reaction_networks_constraint[9] = @reaction_network rnc9 begin
     (k1,k2), X1 + X2 ↔ X3
     (k3,k4), X3 + X4 ↔ X5
     (k5,k6), X5 + X6 ↔ X7
 end k1 k2 k3 k4 k5 k6
 reaction_network_constraints[9] = [ 1 0 1 0 1 0 1; -1 1 0 0 0 0 0; 0 0 0 1 1 0 1; 0 0 0 0 0 1 1 ]
 
-reaction_networks_constraint[10] = @reaction_network begin
+reaction_networks_constraint[10] = @reaction_network rnc10 begin
     kDeg,       (w,w2,w2v,v,w2v2,vP,σB,w2σB) ⟶ ∅
     kDeg,       vPp ⟶ phos
     (kBw,kDw),  2w ⟷ w2
@@ -246,7 +246,7 @@ reaction_network_constraints[10] = [ 0 0 0 0 0 0 0 0 1 1 ]
 ### Reaction networks that are actual models that have been used ###
 
 # Brusselator.
-reaction_networks_real[1] = @reaction_network begin
+reaction_networks_real[1] = @reaction_network rnr1 begin
     A, ∅ → X
     1, 2X + Y → 3X
     B, X → Y
@@ -254,7 +254,7 @@ reaction_networks_real[1] = @reaction_network begin
 end A B;
 
 # The B.subtilis σV Lysozyme stress response.
-reaction_networks_real[2] = @reaction_network begin
+reaction_networks_real[2] = @reaction_network rnr2 begin
     v0 + hill(σ,v,K,n), ∅ → (σ+A)
     deg, (σ,A,Aσ) → ∅
     (kB,kD), A + σ ↔ Aσ
@@ -262,7 +262,7 @@ reaction_networks_real[2] = @reaction_network begin
 end v0 v K n kD kB kC deg S;
 
 # A cell cycle model
-reaction_networks_real[3] = @reaction_network begin
+reaction_networks_real[3] = @reaction_network rnr3 begin
   k1, 0 --> Y
   k2p, Y --> 0
   k2pp*P, Y --> 0
@@ -273,7 +273,7 @@ end k1 k2p k2pp k3p k3pp A J3 k4 m J4
 #reaction_networks_real[3] = cc_network
 
 # A bistable switch
-reaction_networks_real[4] = @reaction_network begin
+reaction_networks_real[4] = @reaction_network rnr4 begin
     d,    (X,Y) → ∅
     hillr(Y,v1,K1,n1), ∅ → X
     hillr(X,v2,K2,n2), ∅ → Y
@@ -281,25 +281,25 @@ end d v1 K1 n1 v2 K2 n2
 
 ### Reaction networks that contain weird functions, stuff, and other oddities ###
 
-reaction_networks_weird[1] = @reaction_network begin
+reaction_networks_weird[1] = @reaction_network rnw1 begin
     exp(p), ∅ → X
     d*exp(X), X → ∅
 end p d
 
-reaction_networks_weird[2] = @reaction_network begin
+reaction_networks_weird[2] = @reaction_network rnw2 begin
     k1, ∅ → X
     k2*log(12+X), X → Y
     k3*log(3+Y), Y → Z
     log(5,6+k4), Z → ∅
 end k1 k2 k3 k4
 
-reaction_networks_weird[3] = @reaction_network begin
+reaction_networks_weird[3] = @reaction_network rnw3 begin
     d,    (X,Y) → ∅
     hillr(hill(Y,v2,K2,n2),v1,K1,n1), ∅ → X
     hillr(hill(X,v1,K1,n1),v2,K2,n2), ∅ → Y
 end d v1 K1 n1 v2 K2 n2
 
-reaction_networks_weird[4] = @reaction_network begin
+reaction_networks_weird[4] = @reaction_network rnw4 begin
     p, ∅ → X1
     (k1,X2^X3), X1 ⟷ X2
     (X2/X4,k4), X2 ⟷ X3
@@ -307,13 +307,13 @@ reaction_networks_weird[4] = @reaction_network begin
     d*X2, X4 → ∅
 end p k1 k2 k3 k4 k5 k6 d
 
-reaction_networks_weird[5] = @reaction_network begin
+reaction_networks_weird[5] = @reaction_network rnw5 begin
     (k1*tanh(X1),k2), X1 ↔ 2X1
     (2+sin(X1),3+sin(k3)), X1 + X2 ↔ X3
     (k5,4+cos(X3+X1)+sin(X2)), X3 ↔ X2
 end k1 k2 k3 k4 k5 k6
 
-reaction_networks_weird[6] = @reaction_network begin
+reaction_networks_weird[6] = @reaction_network rnw6 begin
     (p,d), ∅ ↔ X1
     (p,d), ∅ ↔ X1
     (p,d), ∅ ↔ X1
@@ -322,14 +322,14 @@ reaction_networks_weird[6] = @reaction_network begin
     (p,d), ∅ ↔ X1
 end p d
 
-reaction_networks_weird[7] = @reaction_network begin
+reaction_networks_weird[7] = @reaction_network rnw7 begin
     k1, X1 → X2
     0, X2 → X3
     k2, X3 → X4
     k3, X4 → X5
 end k1 k2 k3
 
-reaction_networks_weird[8] = @reaction_network begin
+reaction_networks_weird[8] = @reaction_network rnw8 begin
     k1+X3, X1 → X2
     k2*X4, X2 → X3
     k3/(X5+0.01), X3 → X4
@@ -341,12 +341,12 @@ reaction_networks_weird[8] = @reaction_network begin
     X2^3+2X2^2+k9, X9 → X1
 end k1 k2 k3 k4 k5 k6 k7 k8 k9
 
-reaction_networks_weird[9] = @reaction_network begin
+reaction_networks_weird[9] = @reaction_network rnw9 begin
     p, ∅ → X1
     hill(hill(hill(hill(X1,v1,K1,n1),v2,K2,n2),v3,K3,n3),v4,K4,n4), X1 → ∅
 end p v1 K1 n1 v2 K2 n2 v3 K3 n3 v4 K4 n4
 
-reaction_networks_weird[10] = @reaction_network begin
+reaction_networks_weird[10] = @reaction_network rnw10 begin
     d, 5X1 → 4X1
 end d
 

From 16d84ed8f592a2e31954692604fb7257722d17e9 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Wed, 15 Sep 2021 21:52:16 -0400
Subject: [PATCH 06/11] fix API tests

---
 test/api.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/test/api.jl b/test/api.jl
index 1cd2117363..96aacdb63c 100644
--- a/test/api.jl
+++ b/test/api.jl
@@ -16,7 +16,7 @@ pset = Set([value(k1) => 1, value(k2) => 2])
 
 rxs2 = [Reaction(k2, [I], [R], [1], [1]),
         Reaction(k1, [S,I], [I], [1,1], [2])]
-@named rs2 = ReactionSystem(rxs2, t, [R,I,S], [k2,k1])
+rs2 = ReactionSystem(rxs2, t, [R,I,S], [k2,k1], name=:rs)
 @test rs == rs2
 
 rs3 = make_empty_network()

From e082f6cd6a75a2171a59b216b7809d4386d2aeda Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Wed, 15 Sep 2021 22:01:42 -0400
Subject: [PATCH 07/11] doc updates

---
 docs/src/api/catalyst_api.md | 20 ++++++++++--------
 src/networkapi.jl            | 40 +++++++++++++++++++++++++-----------
 2 files changed, 40 insertions(+), 20 deletions(-)

diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md
index 90187a7de9..cdeac14cd1 100644
--- a/docs/src/api/catalyst_api.md
+++ b/docs/src/api/catalyst_api.md
@@ -88,6 +88,16 @@ numparams
 numreactionparams
 ```
 
+## ModelingToolkit-Inherited Accessor Functions
+- `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.
+
 ## Basic Reaction Properties
 ```@docs
 ismassaction
@@ -99,13 +109,6 @@ 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
@@ -137,8 +140,9 @@ isweaklyreversible
 
 ## Network Comparison 
 ```@docs
-==(rn1::ReactionSystem, rn2::ReactionSystem)
 ==(rn1::Reaction, rn2::Reaction)
+isequal_without_names
+==(rn1::ReactionSystem, rn2::ReactionSystem)
 ```
 
 ## Network Visualization
diff --git a/src/networkapi.jl b/src/networkapi.jl
index b5c3ed3539..ebcca1ebfb 100644
--- a/src/networkapi.jl
+++ b/src/networkapi.jl
@@ -773,7 +773,7 @@ conservedquantities(state, cons_laws) = cons_laws * state
 ######################## reaction network operators #######################
 
 """
-    ==(rn1::Reaction, rn2::Reaction)
+    ==(rx1::Reaction, rx2::Reaction)
 
 Tests whether two [`Reaction`](@ref)s are identical.
 
@@ -782,28 +782,27 @@ Notes:
 - *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
     considered different than `(A+1)^2`.
 """
-function (==)(rn1::Reaction, rn2::Reaction)
-    isequal(rn1.rate, rn2.rate) || return false
-    issetequal(zip(rn1.substrates,rn1.substoich), zip(rn2.substrates,rn2.substoich)) || return false
-    issetequal(zip(rn1.products,rn1.prodstoich), zip(rn2.products,rn2.prodstoich)) || return false
-    issetequal(rn1.netstoich, rn2.netstoich) || return false
-    rn1.only_use_rate == rn2.only_use_rate
+function (==)(rx1::Reaction, rx2::Reaction)
+    isequal(rx1.rate, rx2.rate) || return false
+    issetequal(zip(rx1.substrates,rx1.substoich), zip(rx2.substrates,rx2.substoich)) || return false
+    issetequal(zip(rx1.products,rx1.prodstoich), zip(rx2.products,rx2.prodstoich)) || return false
+    issetequal(rx1.netstoich, rx2.netstoich) || return false
+    rx1.only_use_rate == rx2.only_use_rate
 end
 
-
 """
-    ==(rn1::ReactionSystem, rn2::ReactionSystem)
+    isequal_without_names(rn1::ReactionSystem, rn2::ReactionSystem)
 
 Tests whether the underlying species, parameters and reactions are the same in
-the two [`ReactionSystem`](@ref)s. 
+the two [`ReactionSystem`](@ref)s. Ignores the names of the systems in testing
+equality.
 
 Notes:
 - *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
     considered different than `(A+1)^2`.
 - Does not include `defaults` in determining equality.
 """
-function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
-    (nameof(rn1) == nameof(rn2)) || return false
+function isequal_without_names(rn1::ReactionSystem, rn2::ReactionSystem)
     isequal(get_iv(rn1), get_iv(rn2)) || return false
     issetequal(get_states(rn1), get_states(rn2)) || return false
     issetequal(get_ps(rn1), get_ps(rn2)) || return false
@@ -821,6 +820,23 @@ function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
     true
 end
 
+"""
+    ==(rn1::ReactionSystem, rn2::ReactionSystem)
+
+Tests whether the underlying species, parameters and reactions are the same in
+the two [`ReactionSystem`](@ref)s. Requires the systems to have the same names
+too.
+
+Notes:
+- *Does not* currently simplify rates, so a rate of `A^2+2*A+1` would be
+    considered different than `(A+1)^2`.
+- Does not include `defaults` in determining equality.
+"""
+function (==)(rn1::ReactionSystem, rn2::ReactionSystem)
+    (nameof(rn1) == nameof(rn2)) || return false
+    isequal_without_names(rn1,rn2)
+end
+
 
 ######################## functions to extend a network ####################
 

From 711b8161e86418aae43e260a43dde62d1ae7421b Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Wed, 15 Sep 2021 22:20:04 -0400
Subject: [PATCH 08/11] tweak docs

---
 docs/src/api/catalyst_api.md |  4 ++--
 src/networkapi.jl            | 22 +++++++++++-----------
 2 files changed, 13 insertions(+), 13 deletions(-)

diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md
index cdeac14cd1..aed85bb4de 100644
--- a/docs/src/api/catalyst_api.md
+++ b/docs/src/api/catalyst_api.md
@@ -79,12 +79,11 @@ ReactionSystem
 ```@docs
 species
 speciesmap
-ModelingToolkit.parameters
 reactionparams
 paramsmap
 reactions
 numspecies
-numparams
+numreactions
 numreactionparams
 ```
 
@@ -117,6 +116,7 @@ addparam!
 addreaction!
 ModelingToolkit.extend
 ModelingToolkit.compose
+merge!
 ```
 
 ## Network Analysis and Representations
diff --git a/src/networkapi.jl b/src/networkapi.jl
index ebcca1ebfb..4a4e65808c 100644
--- a/src/networkapi.jl
+++ b/src/networkapi.jl
@@ -85,8 +85,7 @@ end
 Given a [`ReactionSystem`](@ref), return a vector of all `Reactions` in the system.
 
 Notes:
-- If `ModelingToolkit.get_systems(network)` is not empty, may allocate. Otherwise returns
-  `ModelingToolkit.get_eqs(network)`.
+- If `ModelingToolkit.get_systems(network)` is not empty, will allocate.
 """
 function reactions(network, usenamespace=false)
     rxs = usenamespace ? map(eq -> namespace_equation(eq,network), get_eqs(network)) : get_eqs(network)
@@ -99,7 +98,7 @@ end
     speciesmap(network)
 
 Given a [`ReactionSystem`](@ref), return a Dictionary mapping species that
-participate in `Reaction`s to their index within [`species(networw)`](@ref).
+participate in `Reaction`s to their index within [`species(network)`](@ref).
 """
 function speciesmap(network)
     Dict(S => i for (i,S) in enumerate(species(network)))
@@ -109,7 +108,7 @@ end
     paramsmap(network)
 
 Given a [`ReactionSystem`](@ref), return a Dictionary mapping from parameters that
-appear within `Reaction`s to their index within [`reactionparams(network`](@ref).
+appear within `Reaction`s to their index within [`reactionparams(network)`](@ref).
 """
 function paramsmap(network)
     Dict(p => i for (i,p) in enumerate(reactionparams(network)))
@@ -118,8 +117,8 @@ end
 """
     numspecies(network)
 
-Return the number of species within the given [`ReactionSystem`](@ref) that
-participate in `Reaction`s. 
+Return the total number of species within the given [`ReactionSystem`](@ref) and
+subsystems that are `ReactionSystem`s. 
 
 Notes
 - If there are no subsystems this will be fast.
@@ -133,7 +132,8 @@ end
 """
     numreactions(network)
 
-Return the number of reactions within the given [`ReactionSystem`](@ref).
+Return the total number of reactions within the given [`ReactionSystem`](@ref)
+and subsystems that are `ReactionSystem`s.
 """
 function numreactions(network)
     nr = length(get_eqs(network))
@@ -146,13 +146,13 @@ end
 """
     numreactionparams(network)
 
-Return the number of parameters within the given [`ReactionSystem`](@ref) that
-participate in `Reaction`s.
+Return the total number of parameters within the given [`ReactionSystem`](@ref)
+and subsystems that are `ReactionSystem`s.
 
 Notes
 - If there are no subsystems this will be fast.
-- As this calls [`reactionparams`](@ref), it can be slow and will allocate
-  if there are any subsystems. 
+- As this calls [`reactionparams`](@ref), it can be slow and will allocate if
+  there are any subsystems. 
 """
 function numreactionparams(network)
     length(reactionparams(network))

From 3f6f4e70bfaf9aacc743ba1da6b9f9cf3c41ee52 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Thu, 16 Sep 2021 08:53:06 -0400
Subject: [PATCH 09/11] cleanup base module defs

---
 HISTORY.md                   |  4 ++--
 docs/src/api/catalyst_api.md |  2 +-
 src/Catalyst.jl              | 17 +++++++----------
 3 files changed, 10 insertions(+), 13 deletions(-)

diff --git a/HISTORY.md b/HISTORY.md
index 49e9ea90de..0313900203 100644
--- a/HISTORY.md
+++ b/HISTORY.md
@@ -1,8 +1,8 @@
 # Breaking updates and feature summaries across releases
 
 ## Catalyst unreleased (master branch) 
-- **BREAKING:** added a dependency on Nemo.jl for integer matrix linear algebra methods.
-- Added `incidencematgraph`, `linkageclasses`, `deficiency`, `subnetworks`, `linkagedeficiency`, `isreversible` and `isweaklyreversible` API functions.
+- Added `incidencematgraph`, `linkageclasses`, `deficiency`, `subnetworks`,
+  `linkagedeficiency`, `isreversible` and `isweaklyreversible` API functions.
 - Added the ability to compose `ReactionSystem`s via subsystems, and include
   either `ODESystem`s or `NonlinearSystem`s as subsystems. Note, if using
   subsystems it is not currently possible to convert to a `JumpSystem`. It is
diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md
index aed85bb4de..b4bb3a852d 100644
--- a/docs/src/api/catalyst_api.md
+++ b/docs/src/api/catalyst_api.md
@@ -116,7 +116,7 @@ addparam!
 addreaction!
 ModelingToolkit.extend
 ModelingToolkit.compose
-merge!
+merge!(network1::ReactionSystem, network2::ReactionSystem)
 ```
 
 ## Network Analysis and Representations
diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index 434385b57f..ba2083e86b 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -8,7 +8,8 @@ using SparseArrays, DiffEqBase, Reexport, DiffEqJump
 using Latexify, Requires
 
 # ModelingToolkit imports and convenience functions we use
-using ModelingToolkit; 
+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_defaults, get_observed
@@ -18,17 +19,13 @@ import ModelingToolkit: get_variables, namespace_expr, namespace_equation, get_v
 # 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
-const LG = LightGraphs
-const MT = ModelingToolkit
-
-import AbstractAlgebra; const AA = AbstractAlgebra
-
+import MacroTools, LightGraphs, AbstractAlgebra
 
+# globals for the modulate
+const LG = LightGraphs
+const AA = AbstractAlgebra
+const DEFAULT_IV = (@parameters t)[1]
 
 # as used in Catlab
 const USE_GV_JLL = Ref(false)

From 4db9ccead749704bd844675c9b22411418eeeeae Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Thu, 16 Sep 2021 09:22:45 -0400
Subject: [PATCH 10/11] update accessor docs

---
 docs/src/tutorials/generated_systems.md | 91 +++++++++++++++++++------
 1 file changed, 70 insertions(+), 21 deletions(-)

diff --git a/docs/src/tutorials/generated_systems.md b/docs/src/tutorials/generated_systems.md
index c889007ed1..5c858b00cd 100644
--- a/docs/src/tutorials/generated_systems.md
+++ b/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`.
-* [`parameters(rn)`](@ref) 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`.
+
+### [`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,9 +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 [`ModelingToolkit.extend`](@ref) and
-`ModelingToolkit.compose`](@ref).

From 179b58d50bc24d931425798b1eca43f09d8b0644 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Thu, 16 Sep 2021 09:34:23 -0400
Subject: [PATCH 11/11] update docs more

---
 docs/src/api/catalyst_api.md            | 6 ++++++
 docs/src/tutorials/generated_systems.md | 2 +-
 2 files changed, 7 insertions(+), 1 deletion(-)

diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md
index b4bb3a852d..d19dc2e8c6 100644
--- a/docs/src/api/catalyst_api.md
+++ b/docs/src/api/catalyst_api.md
@@ -76,6 +76,8 @@ ReactionSystem
 ```
 
 ## Basic System Properties
+See [The generated `ReactionSystem` and `Reaction`s](@ref) for more details
+
 ```@docs
 species
 speciesmap
@@ -88,6 +90,8 @@ 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.
@@ -96,6 +100,8 @@ numreactionparams
 - `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
diff --git a/docs/src/tutorials/generated_systems.md b/docs/src/tutorials/generated_systems.md
index 5c858b00cd..efba70c66b 100644
--- a/docs/src/tutorials/generated_systems.md
+++ b/docs/src/tutorials/generated_systems.md
@@ -1,4 +1,4 @@
-# The generated [`ReactionSystem`](@ref) and [`Reaction`](@ref)s
+# The generated `ReactionSystem` and `Reaction`s
 
 ### [`ReactionSystem`](@ref) Accessors