From 5ef4f1a923724d8edc8bc5540702a00cd6a9717d Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Mon, 3 Jun 2024 17:33:55 -0400
Subject: [PATCH 1/9] init

---
 src/network_analysis.jl | 124 +++++++++++++++++++++++++++++++++-------
 1 file changed, 104 insertions(+), 20 deletions(-)

diff --git a/src/network_analysis.jl b/src/network_analysis.jl
index cc33f9b63e..93b0429d88 100644
--- a/src/network_analysis.jl
+++ b/src/network_analysis.jl
@@ -4,6 +4,8 @@
 function filter_constspecs(specs, stoich::AbstractVector{V}, smap) where {V <: Integer}
     isempty(specs) && (return Vector{Int}(), Vector{V}())
 
+    # If any species are constant, go through these manually and add their indices and 
+    # stoichiometries to `ids` and `filtered_stoich`.
     if any(isconstant, specs)
         ids = Vector{Int}()
         filtered_stoich = Vector{V}()
@@ -44,11 +46,15 @@ function reactioncomplexmap(rn::ReactionSystem)
     !isempty(nps.complextorxsmap) && return nps.complextorxsmap
     complextorxsmap = nps.complextorxsmap
 
+    # Retrieves system reactions and a map from species to their index in the species vector.
     rxs = reactions(rn)
     smap = speciesmap(rn)
     numreactions(rn) > 0 ||
         error("There must be at least one reaction to find reaction complexes.")
+
     for (i, rx) in enumerate(rxs)
+        # Create the `ReactionComplex` corresponding to the reaction's substrates. Adds it 
+        # to the reaction complex dictionary (recording it as the substrates of the i'th reaction). 
         subids, substoich = filter_constspecs(rx.substrates, rx.substoich, smap)
         subrc = sort!(ReactionComplex(subids, substoich))
         if haskey(complextorxsmap, subrc)
@@ -57,6 +63,8 @@ function reactioncomplexmap(rn::ReactionSystem)
             complextorxsmap[subrc] = [i => -1]
         end
 
+        # Create the `ReactionComplex` corresponding to the reaction's products. Adds it 
+        # to the reaction complex dictionary (recording it as the products of the i'th reaction). 
         prodids, prodstoich = filter_constspecs(rx.products, rx.prodstoich, smap)
         prodrc = sort!(ReactionComplex(prodids, prodstoich))
         if haskey(complextorxsmap, prodrc)
@@ -94,7 +102,10 @@ function reactioncomplexes(rn::ReactionSystem; sparse = false)
     isempty(get_systems(rn)) ||
         error("reactioncomplexes does not currently support subsystems.")
     nps = get_networkproperties(rn)
+
+    # If the complexes have not been cached, or the cached complexes uses a different sparsity.
     if isempty(nps.complexes) || (sparse != issparse(nps.complexes))
+        # Computes the reaction complex dictionary. Use it to create a sparse/dense matrix.
         complextorxsmap = reactioncomplexmap(rn)
         nps.complexes, nps.incidencemat = if sparse
             reactioncomplexes(SparseMatrixCSC{Int, Int}, rn, complextorxsmap)
@@ -105,8 +116,11 @@ function reactioncomplexes(rn::ReactionSystem; sparse = false)
     nps.complexes, nps.incidencemat
 end
 
+# Creates a *sparse* reaction complex matrix.
 function reactioncomplexes(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSystem,
                            complextorxsmap)
+    # Computes the I, J, and V vectors used for the sparse matrix (read about sparse matrix 
+    # representation for more information).
     complexes = collect(keys(complextorxsmap))
     Is = Int[]
     Js = Int[]
@@ -122,6 +136,7 @@ function reactioncomplexes(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSystem
     complexes, B
 end
 
+# Creates a *dense* reaction complex matrix.
 function reactioncomplexes(::Type{Matrix{Int}}, rn::ReactionSystem, complextorxsmap)
     complexes = collect(keys(complextorxsmap))
     B = zeros(Int, length(complexes), numreactions(rn))
@@ -158,6 +173,9 @@ function complexstoichmat(rn::ReactionSystem; sparse = false)
     isempty(get_systems(rn)) ||
         error("complexstoichmat does not currently support subsystems.")
     nps = get_networkproperties(rn)
+
+    # If the complexes stoichiometry matrix has not been cached, or the cached one uses a
+    # different sparsity, computes (and caches) it.
     if isempty(nps.complexstoichmat) || (sparse != issparse(nps.complexstoichmat))
         nps.complexstoichmat = if sparse
             complexstoichmat(SparseMatrixCSC{Int, Int}, rn, keys(reactioncomplexmap(rn)))
@@ -168,7 +186,10 @@ function complexstoichmat(rn::ReactionSystem; sparse = false)
     nps.complexstoichmat
 end
 
+# Creates a *sparse* reaction complex stoichiometry matrix.
 function complexstoichmat(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSystem, rcs)
+    # Computes the I, J, and V vectors used for the sparse matrix (read about sparse matrix 
+    # representation for more information).
     Is = Int[]
     Js = Int[]
     Vs = Int[]
@@ -182,6 +203,7 @@ function complexstoichmat(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSystem,
     Z = sparse(Is, Js, Vs, numspecies(rn), length(rcs))
 end
 
+# Creates a *dense* reaction complex stoichiometry matrix.
 function complexstoichmat(::Type{Matrix{Int}}, rn::ReactionSystem, rcs)
     Z = zeros(Int, numspecies(rn), length(rcs))
     for (i, rc) in enumerate(rcs)
@@ -213,6 +235,9 @@ function complexoutgoingmat(rn::ReactionSystem; sparse = false)
     isempty(get_systems(rn)) ||
         error("complexoutgoingmat does not currently support subsystems.")
     nps = get_networkproperties(rn)
+
+    # If the outgoing complexes matrix has not been cached, or the cached one uses a
+    # different sparsity, computes (and caches) it.
     if isempty(nps.complexoutgoingmat) || (sparse != issparse(nps.complexoutgoingmat))
         B = reactioncomplexes(rn, sparse = sparse)[2]
         nps.complexoutgoingmat = if sparse
@@ -224,16 +249,22 @@ function complexoutgoingmat(rn::ReactionSystem; sparse = false)
     nps.complexoutgoingmat
 end
 
+# Creates a *sparse* outgoing reaction complex stoichiometry matrix.
 function complexoutgoingmat(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSystem, B)
+    # Computes the I, J, and V vectors used for the sparse matrix (read about sparse matrix 
+    # representation for more information).
     n = size(B, 2)
     rows = rowvals(B)
     vals = nonzeros(B)
     Is = Int[]
     Js = Int[]
     Vs = Int[]
+
+    # Allocates space to the vectors (so that it is not done incrementally in the loop). 
     sizehint!(Is, div(length(vals), 2))
     sizehint!(Js, div(length(vals), 2))
     sizehint!(Vs, div(length(vals), 2))
+
     for j in 1:n
         for i in nzrange(B, j)
             if vals[i] != one(eltype(vals))
@@ -246,6 +277,7 @@ function complexoutgoingmat(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSyste
     sparse(Is, Js, Vs, size(B, 1), size(B, 2))
 end
 
+# Creates a *dense* outgoing reaction complex stoichiometry matrix.
 function complexoutgoingmat(::Type{Matrix{Int}}, rn::ReactionSystem, B)
     Δ = copy(B)
     for (I, b) in pairs(Δ)
@@ -284,10 +316,14 @@ function incidencematgraph(rn::ReactionSystem)
     nps.incidencegraph
 end
 
+# Computes the incidence graph from an *dense* incidence matrix.
 function incidencematgraph(incidencemat::Matrix{Int})
     @assert all(∈([-1, 0, 1]), incidencemat)
     n = size(incidencemat, 1)  # no. of nodes/complexes
     graph = Graphs.DiGraph(n)
+
+    # Walks through each column (corresponds to reactions). For each, find the input and output
+    # complex and add an edge representing these to the incidence graph.
     for col in eachcol(incidencemat)
         src = 0
         dst = 0
@@ -301,12 +337,16 @@ function incidencematgraph(incidencemat::Matrix{Int})
     return graph
 end
 
+# Computes the incidence graph from an *sparse* incidence matrix.
 function incidencematgraph(incidencemat::SparseMatrixCSC{Int, Int})
     @assert all(∈([-1, 0, 1]), incidencemat)
     m, n = size(incidencemat)
     graph = Graphs.DiGraph(m)
     rows = rowvals(incidencemat)
     vals = nonzeros(incidencemat)
+
+    # Loops through the (n) columns. For each column, directly find the index of the input
+    # and output complex and add an edge representing these to the incidence graph.
     for j in 1:n
         inds = nzrange(incidencemat, j)
         row = rows[inds]
@@ -387,32 +427,43 @@ rcs,incidencemat = reactioncomplexes(sir)
 ```
 """
 function deficiency(rn::ReactionSystem)
+    # Precomputes required information. `conservationlaws` caches the conservation laws in `rn`.
     nps = get_networkproperties(rn)
     conservationlaws(rn)
     r = nps.rank
     ig = incidencematgraph(rn)
     lc = linkageclasses(rn)
+
+    # Computes deficiency using its formula. Caches and returns it as output.
     nps.deficiency = Graphs.nv(ig) - length(lc) - r
     nps.deficiency
 end
 
-# Used in the subsequent function.
+# For a linkage class (set of reaction complexes that form an isolated sub-graph), and some
+# additional information of the full network, find the reactions, species, and parameters 
+# that constitute the corresponding sub-reaction network.
 function subnetworkmapping(linkageclass, allrxs, complextorxsmap, p)
+    # Finds the reactions that are part of teh sub-reaction network.
     rxinds = sort!(collect(Set(rxidx for rcidx in linkageclass
                                for rxidx in complextorxsmap[rcidx])))
-    rxs = allrxs[rxinds]
-    specset = Set(s for rx in rxs for s in rx.substrates if !isconstant(s))
-    for rx in rxs
+    newrxs = allrxs[rxinds]
+
+    # Find the species that are part of the sub-reaction network.
+    specset = Set(s for rx in newrxs for s in rx.substrates if !isconstant(s))
+    for rx in newrxs
         for product in rx.products
             !isconstant(product) && push!(specset, product)
         end
     end
-    specs = collect(specset)
+    newspecs = collect(specset)
+
+    # Find the parameters that are part of the sub-reaction network.
     newps = Vector{eltype(p)}()
-    for rx in rxs
+    for rx in newrxs
         Symbolics.get_variables!(newps, rx.rate, p)
     end
-    rxs, specs, newps   # reactions and species involved in reactions of subnetwork
+
+    newrxs, newspecs, newps
 end
 
 """
@@ -436,6 +487,8 @@ subnetworks(sir)
 """
 function subnetworks(rs::ReactionSystem)
     isempty(get_systems(rs)) || error("subnetworks does not currently support subsystems.")
+
+    # Retrieves required components. `linkageclasses` caches linkage classes in `rs`.
     lcs = linkageclasses(rs)
     rxs = reactions(rs)
     p = parameters(rs)
@@ -443,11 +496,14 @@ function subnetworks(rs::ReactionSystem)
     spatial_ivs = get_sivs(rs)
     complextorxsmap = [map(first, rcmap) for rcmap in values(reactioncomplexmap(rs))]
     subnetworks = Vector{ReactionSystem}()
+
+    # Loops through each sub-graph of connected reaction complexes. For each, create a 
+    # new `ReactionSystem` model and pushes it to the subnetworks vector.
     for i in 1:length(lcs)
-        reacs, specs, newps = subnetworkmapping(lcs[i], rxs, complextorxsmap, p)
+        newrxs, newspecs, newps = subnetworkmapping(lcs[i], rxs, complextorxsmap, p)
         newname = Symbol(nameof(rs), "_", i)
         push!(subnetworks,
-              ReactionSystem(reacs, t, specs, newps; name = newname, spatial_ivs))
+              ReactionSystem(newrxs, t, newspecs, newps; name = newname, spatial_ivs))
     end
     subnetworks
 end
@@ -475,6 +531,9 @@ function linkagedeficiencies(rs::ReactionSystem)
     lcs = linkageclasses(rs)
     subnets = subnetworks(rs)
     δ = zeros(Int, length(lcs))
+
+    # For each sub-reaction network of the reaction network, compute its deficiency. Returns
+    # the full vector of deficiencies for each sub-reaction network.
     for (i, subnet) in enumerate(subnets)
         conservationlaws(subnet)
         nps = get_networkproperties(subnet)
@@ -531,11 +590,16 @@ function isweaklyreversible(rn::ReactionSystem, subnets)
     im = get_networkproperties(rn).incidencemat
     isempty(im) &&
         error("Error, please call reactioncomplexes(rn::ReactionSystem) to ensure the incidence matrix has been cached.")
+
+    # For each sub-reaction network, caches its reaction complexes.
     sparseig = issparse(im)
     for subnet in subnets
         nps = get_networkproperties(subnet)
         isempty(nps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig)
     end
+
+    # the network is weakly reversible if each sub-network's incidenc graph is strongl 
+    # connec (i.e. each node is reachable from each node).
     all(Graphs.is_strongly_connected ∘ incidencematgraph, subnets)
 end
 
@@ -645,25 +709,43 @@ end
 
 # Used in the subsequent function.
 function cache_conservationlaw_eqs!(rn::ReactionSystem, N::AbstractMatrix, col_order)
+    # Retrieves nullity (the number of conservation laws). `r` is the rank of the netstoichmat.
     nullity = size(N, 1)
-    r = numspecies(rn) - nullity     # rank of the netstoichmat
-    sts = species(rn)
+    r = numspecies(rn) - nullity
+
+    # Creates vectors of all independent and dependent species (those that are not, and are, 
+    # eliminated by the conservation laws). Get vectors both with their indexes and the actual
+    # species symbolic variables.
+    sps = species(rn)
     indepidxs = col_order[begin:r]
-    indepspecs = sts[indepidxs]
+    indepspecs = sps[indepidxs]
     depidxs = col_order[(r + 1):end]
-    depspecs = sts[depidxs]
+    depspecs = sps[depidxs]
+
+    # Declares the conservation law parameters.
     constants = MT.unwrap.(MT.scalarize(only(
                 @parameters $(CONSERVED_CONSTANT_SYMBOL)[1:nullity] [conserved=true])))
 
+    # Computes the equations for (examples uses simple two-state system, `X1 <--> X2`):
+    # - The species eliminated through conservation laws (`conservedeqs`). E.g. `[X2 ~ X1 - Γ[1]]`.
+    # - The conserved quantity parameters (`constantdefs`). E.g. `[Γ[1] ~ X1 + X2]`.
     conservedeqs = Equation[]
     constantdefs = Equation[]
+
+    # For each conserved quantity.
     for (i, depidx) in enumerate(depidxs)
+        # Finds the coefficient (in the conservation law) of the species that is eliminated
+        # by this conservation law.
         scaleby = (N[i, depidx] != 1) ? N[i, depidx] : one(eltype(N))
         (scaleby != 0) || error("Error, found a zero in the conservation law matrix where "
-              *
-              "one was not expected.")
+            * "one was not expected.")
+
+        # Creates, for this conservation law, the sum of all independent species (weighted by
+        # the ratio between the coefficient of the species and the species which is elimianted.
         coefs = @view N[i, indepidxs]
-        terms = sum(p -> p[1] / scaleby * p[2], zip(coefs, indepspecs))
+        terms = sum((coef, sp) -> coef / scaleby * sp, zip(coefs, indepspecs))
+
+        # Computes the two equations corresponding to this conserved quantity.
         eq = depspecs[i] ~ constants[i] - terms
         push!(conservedeqs, eq)
         eq = constants[i] ~ depspecs[i] + terms
@@ -695,6 +777,8 @@ Notes:
 function conservationlaws(rs::ReactionSystem)
     nps = get_networkproperties(rs)
     !isempty(nps.conservationmat) && (return nps.conservationmat)
+
+    # If the conservation law matrix is not computed, do so and caches the result.
     nsm = netstoichmat(rs)
     nps.conservationmat = conservationlaws(nsm; col_order = nps.col_order)
     cache_conservationlaw_eqs!(rs, nps.conservationmat, nps.col_order)
@@ -708,13 +792,13 @@ Compute conserved quantities for a system with the given conservation laws.
 """
 conservedquantities(state, cons_laws) = cons_laws * state
 
-# If u0s are not given while conservation laws are present, throws an error.
+# If u0s are not given while conservation laws are present, throw an error.
 # Used in HomotopyContinuation and BifurcationKit extensions.
-# Currently only checks if any u0s are given
-# (not whether these are enough for computing conserved quantitites, this will yield a less informative error).
+# Currently only checks if any u0s are given (not whether these are enough for computing
+# conserved quantities, this will yield a less informative error).
 function conservationlaw_errorcheck(rs, pre_varmap)
     vars_with_vals = Set(p[1] for p in pre_varmap)
-    any(s -> s in vars_with_vals, species(rs)) && return
+    any(sp -> sp in vars_with_vals, species(rs)) && return
     isempty(conservedequations(Catalyst.flatten(rs))) ||
         error("The system has conservation laws but initial conditions were not provided for some species.")
 end

From 9ba924a353381250f04bcf6e13654a28a4c5bcf5 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Mon, 3 Jun 2024 17:54:35 -0400
Subject: [PATCH 2/9] up

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

diff --git a/src/network_analysis.jl b/src/network_analysis.jl
index 93b0429d88..813060bea0 100644
--- a/src/network_analysis.jl
+++ b/src/network_analysis.jl
@@ -743,7 +743,7 @@ function cache_conservationlaw_eqs!(rn::ReactionSystem, N::AbstractMatrix, col_o
         # Creates, for this conservation law, the sum of all independent species (weighted by
         # the ratio between the coefficient of the species and the species which is elimianted.
         coefs = @view N[i, indepidxs]
-        terms = sum((coef, sp) -> coef / scaleby * sp, zip(coefs, indepspecs))
+        terms = sum(coef / scaleby * sp for (coef, sp) in zip(coefs, indepspecs))
 
         # Computes the two equations corresponding to this conserved quantity.
         eq = depspecs[i] ~ constants[i] - terms

From 245b12d802bd4c84d20d532b22709e80d39dfe7a Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Tue, 4 Jun 2024 10:57:24 -0400
Subject: [PATCH 3/9] save progress

---
 src/network_analysis.jl |  2 +-
 src/reactionsystem.jl   | 64 +++++++++++++++++++++++++++++++++++++++++
 2 files changed, 65 insertions(+), 1 deletion(-)

diff --git a/src/network_analysis.jl b/src/network_analysis.jl
index 813060bea0..74245fa4a2 100644
--- a/src/network_analysis.jl
+++ b/src/network_analysis.jl
@@ -727,7 +727,7 @@ function cache_conservationlaw_eqs!(rn::ReactionSystem, N::AbstractMatrix, col_o
                 @parameters $(CONSERVED_CONSTANT_SYMBOL)[1:nullity] [conserved=true])))
 
     # Computes the equations for (examples uses simple two-state system, `X1 <--> X2`):
-    # - The species eliminated through conservation laws (`conservedeqs`). E.g. `[X2 ~ X1 - Γ[1]]`.
+    # - The species eliminated through conservation laws (`conservedeqs`). E.g. `[X2 ~ Γ[1] - X1]`.
     # - The conserved quantity parameters (`constantdefs`). E.g. `[Γ[1] ~ X1 + X2]`.
     conservedeqs = Equation[]
     constantdefs = Equation[]
diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
index edb0f9fedf..57ddadcabc 100644
--- a/src/reactionsystem.jl
+++ b/src/reactionsystem.jl
@@ -76,25 +76,89 @@ Base.Sort.defalg(::ReactionComplex) = Base.DEFAULT_UNSTABLE
 
 #! format: off
 # Internal cache for various ReactionSystem calculated properties
+# All related functionality is in the "network_analysis.jl" file. However, this must be declared
+# here as the structure is part of the `ReactionSystem` structure.
 Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Real}}
+    """Flag which is switched to `true` once any field is updated."""
     isempty::Bool = true
     netstoichmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0)
     conservationmat::Matrix{I} = Matrix{I}(undef, 0, 0)
     col_order::Vector{Int} = Int[]
+    """
+    The reaction networks *rank* (i.e. the span of the columns of its net stoichiometry matrix,
+    or its number of independent species).
+    """
     rank::Int = 0
     nullity::Int = 0
+    """
+    The set of *independent species* of the reaction system (i.e. species that will not be
+    eliminated when we eliminate the conserved quantities.
+    """
     indepspecs::Set{V} = Set{V}()
+    """
+    The set of *dependent species* of the reaction system. These species are eliminated when
+    we eliminated the conserved quantities. In the resulting `ODESystem` these become
+    observables, not unknowns.
+    """
     depspecs::Set{V} = Set{V}()
+    """
+    The equations for the (dependent) species eliminated by any conservation laws. I.e. for
+    the two simple two state system (`X1 <--> X2`) `X2` becomes a dependant species with the
+    conserved equation `X2 ~ Γ[1] - X1`.
+    """
     conservedeqs::Vector{Equation} = Equation[]
+    """
+    The equations for the conserved quantity parameters. I.e. for the two simple two state
+    system (`X1 <--> X2`) there is one conserved quantity with the equation `Γ[1] ~ X1 + X2`.
+    """
     constantdefs::Vector{Equation} = Equation[]
     speciesmap::Dict{V, Int} = Dict{V, Int}()
+    """
+    A dictionary from each reaction complex to the reactions they participate it. The value
+    mapped from each reaction complex is a pair from the reaction's index to a value which is
+    `-1` if the complex is a substrate and `+1` if the complex is a product.    
+    """
     complextorxsmap::OrderedDict{ReactionComplex{Int}, Vector{Pair{Int, Int}}} = OrderedDict{ReactionComplex{Int},Vector{Pair{Int,Int}}}()
+    """ A vector with all the reaction system's reaction complexes """
     complexes::Vector{ReactionComplex{Int}} = Vector{ReactionComplex{Int}}(undef, 0)
+    """
+    An MxN matrix where M is the number of reaction complexes and N the number of reactions.
+    Element i,j is:
+        -1 if the i'th complex is a substrate of the j'th reaction.
+        +1 if the i'th complex is a product of the j'th reaction.
+        0 if the i'th complex is not part of the j'th reaction.
+    """
     incidencemat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0)
+    """
+    An MxN matrix where M is the number of species and N the number of reaction complexes.
+    Element i,j is the coefficient of the i'th species in the j'th complex (0 entries denote
+    species that are not part of the corresponding complex). Whether the matrix is sparse
+    is designated when it is created.
+    """
     complexstoichmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0)
+    """
+    An MxN matrix where M is the number of reaction complexes and N the number of reactions.
+    Element i,j is -1 if i'th complex is a substrate of the j'th reaction (and 0 otherwise).
+    """
     complexoutgoingmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0)
+    """
+    A (directed) graph, with nodes corresponding to reaction complexes and edges to reactions.
+    There is an edge from complex i to complex j if there is a reaction converting complex
+    i to complex j.
+    """
     incidencegraph::Graphs.SimpleDiGraph{Int} = Graphs.DiGraph()
+    """
+    A vector of the connected components of the incidence graph. Each element of the
+    `linkageclasses` corresponds to a connected component, and is a vector listing all the
+    reaction complexes in that connected component.    
+    """
     linkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
+    """
+    The networks deficiency. It is computed as *n - l - r*, where *n* is the number of reaction
+    complexes, *l* is the number of linkage classes (i.e. the number of connected components
+    in the incidence graph), and *r* is the reaction networks *rank* (i.e. the span of the columns
+    of its net stoichiometry matrix, or its number of independent species).
+    """
     deficiency::Int = 0
 end
 #! format: on

From 416440001cc033408fd51a863e2f7b9bec03c053 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Tue, 4 Jun 2024 11:25:33 -0400
Subject: [PATCH 4/9] up

---
 src/network_analysis.jl |  1 +
 src/reactionsystem.jl   | 17 +++++++++++++++++
 2 files changed, 18 insertions(+)

diff --git a/src/network_analysis.jl b/src/network_analysis.jl
index 74245fa4a2..d8c533c2f8 100644
--- a/src/network_analysis.jl
+++ b/src/network_analysis.jl
@@ -692,6 +692,7 @@ conservation laws, each represented as a row in the output.
 function conservationlaws(nsm::T; col_order = nothing) where {T <: AbstractMatrix}
 
     # compute the left nullspace over the integers
+    # The `nullspace` function updates the `col_order`.
     N = MT.nullspace(nsm'; col_order)
 
     # if all coefficients for a conservation law are negative, make positive
diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
index 57ddadcabc..9e1dffa299 100644
--- a/src/reactionsystem.jl
+++ b/src/reactionsystem.jl
@@ -81,7 +81,17 @@ Base.Sort.defalg(::ReactionComplex) = Base.DEFAULT_UNSTABLE
 Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Real}}
     """Flag which is switched to `true` once any field is updated."""
     isempty::Bool = true
+    """
+    The reaction network's net stoichiometry matrix. It is a MxN matrix where M is its number of 
+    species and N its number of reaction. Element i,j is net stoichiometric change to the i'th
+    species as a result of the j'th reaction.
+    """
     netstoichmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0)
+    """
+    The reaction network's conservation law matrix. It is a MxN matrix where M is its number of 
+    conservation laws and N its number of species. Element i,j is the coefficient of species 
+    j in the i'th conservation law.
+    """
     conservationmat::Matrix{I} = Matrix{I}(undef, 0, 0)
     col_order::Vector{Int} = Int[]
     """
@@ -89,6 +99,10 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     or its number of independent species).
     """
     rank::Int = 0
+    """
+    The reaction network's *nullity* is its number of species - its rank. This is equal to its
+    number of conservation laws.
+    """
     nullity::Int = 0
     """
     The set of *independent species* of the reaction system (i.e. species that will not be
@@ -112,6 +126,9 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     system (`X1 <--> X2`) there is one conserved quantity with the equation `Γ[1] ~ X1 + X2`.
     """
     constantdefs::Vector{Equation} = Equation[]
+    """
+    A map from each species (as a symbolic variable) to its index in the species vector.
+    """
     speciesmap::Dict{V, Int} = Dict{V, Int}()
     """
     A dictionary from each reaction complex to the reactions they participate it. The value

From 8478423fca015fcf53d6f8795ee8f839ac5df628 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Tue, 4 Jun 2024 11:28:06 -0400
Subject: [PATCH 5/9] spelling fixes

---
 src/reactionsystem.jl | 20 ++++++++++----------
 1 file changed, 10 insertions(+), 10 deletions(-)

diff --git a/src/reactionsystem.jl b/src/reactionsystem.jl
index 9e1dffa299..90e833777c 100644
--- a/src/reactionsystem.jl
+++ b/src/reactionsystem.jl
@@ -82,13 +82,13 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     """Flag which is switched to `true` once any field is updated."""
     isempty::Bool = true
     """
-    The reaction network's net stoichiometry matrix. It is a MxN matrix where M is its number of 
-    species and N its number of reaction. Element i,j is net stoichiometric change to the i'th
-    species as a result of the j'th reaction.
+    The reaction network's net stoichiometry matrix. It is an MxN matrix where M is its number of 
+    species and N is its number of reactions. Element i,j is the net stoichiometric change to the
+    i'th species as a result of the j'th reaction.
     """
     netstoichmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0)
     """
-    The reaction network's conservation law matrix. It is a MxN matrix where M is its number of 
+    The reaction network's conservation law matrix. It is an MxN matrix where M is its number of 
     conservation laws and N its number of species. Element i,j is the coefficient of species 
     j in the i'th conservation law.
     """
@@ -111,18 +111,18 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     indepspecs::Set{V} = Set{V}()
     """
     The set of *dependent species* of the reaction system. These species are eliminated when
-    we eliminated the conserved quantities. In the resulting `ODESystem` these become
+    we eliminate the conserved quantities. In the resulting `ODESystem` these become
     observables, not unknowns.
     """
     depspecs::Set{V} = Set{V}()
     """
     The equations for the (dependent) species eliminated by any conservation laws. I.e. for
-    the two simple two state system (`X1 <--> X2`) `X2` becomes a dependant species with the
+    the two simple two-state system (`X1 <--> X2`) `X2` becomes a dependant species with the
     conserved equation `X2 ~ Γ[1] - X1`.
     """
     conservedeqs::Vector{Equation} = Equation[]
     """
-    The equations for the conserved quantity parameters. I.e. for the two simple two state
+    The equations for the conserved quantity parameters. I.e. for the two simple two-state
     system (`X1 <--> X2`) there is one conserved quantity with the equation `Γ[1] ~ X1 + X2`.
     """
     constantdefs::Vector{Equation} = Equation[]
@@ -139,7 +139,7 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     """ A vector with all the reaction system's reaction complexes """
     complexes::Vector{ReactionComplex{Int}} = Vector{ReactionComplex{Int}}(undef, 0)
     """
-    An MxN matrix where M is the number of reaction complexes and N the number of reactions.
+    An MxN matrix where M is the number of reaction complexes and N is the number of reactions.
     Element i,j is:
         -1 if the i'th complex is a substrate of the j'th reaction.
         +1 if the i'th complex is a product of the j'th reaction.
@@ -154,7 +154,7 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     """
     complexstoichmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0)
     """
-    An MxN matrix where M is the number of reaction complexes and N the number of reactions.
+    An MxN matrix where M is the number of reaction complexes and N is the number of reactions.
     Element i,j is -1 if i'th complex is a substrate of the j'th reaction (and 0 otherwise).
     """
     complexoutgoingmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0)
@@ -171,7 +171,7 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     """
     linkageclasses::Vector{Vector{Int}} = Vector{Vector{Int}}(undef, 0)
     """
-    The networks deficiency. It is computed as *n - l - r*, where *n* is the number of reaction
+    The network's deficiency. It is computed as *n - l - r*, where *n* is the number of reaction
     complexes, *l* is the number of linkage classes (i.e. the number of connected components
     in the incidence graph), and *r* is the reaction networks *rank* (i.e. the span of the columns
     of its net stoichiometry matrix, or its number of independent species).

From 46035327a2d7b36d2284f12f6aea01e2ba7cc605 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 13 Jul 2024 22:27:45 -0400
Subject: [PATCH 6/9] up

---
 src/network_analysis.jl | 78 +++++++++++++++++++++--------------------
 1 file changed, 40 insertions(+), 38 deletions(-)

diff --git a/src/network_analysis.jl b/src/network_analysis.jl
index 4cf0e1dbb9..a932cf9fac 100644
--- a/src/network_analysis.jl
+++ b/src/network_analysis.jl
@@ -4,8 +4,8 @@
 function filter_constspecs(specs, stoich::AbstractVector{V}, smap) where {V <: Integer}
     isempty(specs) && (return Vector{Int}(), Vector{V}())
 
-    # If any species are constant, go through these manually and add their indices and 
-    # stoichiometries to `ids` and `filtered_stoich`.
+    # if any species are constant, go through these manually and add their indices and 
+    # stoichiometries to `ids` and `filtered_stoich`
     if any(isconstant, specs)
         ids = Vector{Int}()
         filtered_stoich = Vector{V}()
@@ -46,7 +46,7 @@ function reactioncomplexmap(rn::ReactionSystem)
     !isempty(nps.complextorxsmap) && return nps.complextorxsmap
     complextorxsmap = nps.complextorxsmap
 
-    # Retrieves system reactions and a map from species to their index in the species vector.
+    # retrieves system reactions and a map from species to their index in the species vector
     rxs = reactions(rn)
     smap = speciesmap(rn)
     numreactions(rn) > 0 ||
@@ -54,7 +54,7 @@ function reactioncomplexmap(rn::ReactionSystem)
 
     for (i, rx) in enumerate(rxs)
         # Create the `ReactionComplex` corresponding to the reaction's substrates. Adds it 
-        # to the reaction complex dictionary (recording it as the substrates of the i'th reaction). 
+        # to the reaction complex dictionary (recording it as the substrates of the i'th reaction).
         subids, substoich = filter_constspecs(rx.substrates, rx.substoich, smap)
         subrc = sort!(ReactionComplex(subids, substoich))
         if haskey(complextorxsmap, subrc)
@@ -64,7 +64,7 @@ function reactioncomplexmap(rn::ReactionSystem)
         end
 
         # Create the `ReactionComplex` corresponding to the reaction's products. Adds it 
-        # to the reaction complex dictionary (recording it as the products of the i'th reaction). 
+        # to the reaction complex dictionary (recording it as the products of the i'th reaction).
         prodids, prodstoich = filter_constspecs(rx.products, rx.prodstoich, smap)
         prodrc = sort!(ReactionComplex(prodids, prodstoich))
         if haskey(complextorxsmap, prodrc)
@@ -103,7 +103,7 @@ function reactioncomplexes(rn::ReactionSystem; sparse = false)
         error("reactioncomplexes does not currently support subsystems.")
     nps = get_networkproperties(rn)
 
-    # If the complexes have not been cached, or the cached complexes uses a different sparsity.
+    # if the complexes have not been cached, or the cached complexes uses a different sparsity
     if isempty(nps.complexes) || (sparse != issparse(nps.complexes))
         # Computes the reaction complex dictionary. Use it to create a sparse/dense matrix.
         complextorxsmap = reactioncomplexmap(rn)
@@ -116,11 +116,11 @@ function reactioncomplexes(rn::ReactionSystem; sparse = false)
     nps.complexes, nps.incidencemat
 end
 
-# Creates a *sparse* reaction complex matrix.
+# creates a *sparse* reaction complex matrix
 function reactioncomplexes(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSystem,
                            complextorxsmap)
-    # Computes the I, J, and V vectors used for the sparse matrix (read about sparse matrix 
-    # representation for more information).
+    # computes the I, J, and V vectors used for the sparse matrix (read about sparse matrix 
+    # representation for more information)
     complexes = collect(keys(complextorxsmap))
     Is = Int[]
     Js = Int[]
@@ -136,7 +136,7 @@ function reactioncomplexes(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSystem
     complexes, B
 end
 
-# Creates a *dense* reaction complex matrix.
+# creates a *dense* reaction complex matrix
 function reactioncomplexes(::Type{Matrix{Int}}, rn::ReactionSystem, complextorxsmap)
     complexes = collect(keys(complextorxsmap))
     B = zeros(Int, length(complexes), numreactions(rn))
@@ -174,8 +174,8 @@ function complexstoichmat(rn::ReactionSystem; sparse = false)
         error("complexstoichmat does not currently support subsystems.")
     nps = get_networkproperties(rn)
 
-    # If the complexes stoichiometry matrix has not been cached, or the cached one uses a
-    # different sparsity, computes (and caches) it.
+    # if the complexes stoichiometry matrix has not been cached, or the cached one uses a
+    # different sparsity, computes (and caches) it
     if isempty(nps.complexstoichmat) || (sparse != issparse(nps.complexstoichmat))
         nps.complexstoichmat = if sparse
             complexstoichmat(SparseMatrixCSC{Int, Int}, rn, keys(reactioncomplexmap(rn)))
@@ -186,10 +186,10 @@ function complexstoichmat(rn::ReactionSystem; sparse = false)
     nps.complexstoichmat
 end
 
-# Creates a *sparse* reaction complex stoichiometry matrix.
+# creates a *sparse* reaction complex stoichiometry matrix
 function complexstoichmat(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSystem, rcs)
-    # Computes the I, J, and V vectors used for the sparse matrix (read about sparse matrix 
-    # representation for more information).
+    # computes the I, J, and V vectors used for the sparse matrix (read about sparse matrix 
+    # representation for more information)
     Is = Int[]
     Js = Int[]
     Vs = Int[]
@@ -203,7 +203,7 @@ function complexstoichmat(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSystem,
     Z = sparse(Is, Js, Vs, numspecies(rn), length(rcs))
 end
 
-# Creates a *dense* reaction complex stoichiometry matrix.
+# creates a *dense* reaction complex stoichiometry matrix
 function complexstoichmat(::Type{Matrix{Int}}, rn::ReactionSystem, rcs)
     Z = zeros(Int, numspecies(rn), length(rcs))
     for (i, rc) in enumerate(rcs)
@@ -236,8 +236,8 @@ function complexoutgoingmat(rn::ReactionSystem; sparse = false)
         error("complexoutgoingmat does not currently support subsystems.")
     nps = get_networkproperties(rn)
 
-    # If the outgoing complexes matrix has not been cached, or the cached one uses a
-    # different sparsity, computes (and caches) it.
+    # if the outgoing complexes matrix has not been cached, or the cached one uses a
+    # different sparsity, computes (and caches) it
     if isempty(nps.complexoutgoingmat) || (sparse != issparse(nps.complexoutgoingmat))
         B = reactioncomplexes(rn, sparse = sparse)[2]
         nps.complexoutgoingmat = if sparse
@@ -249,10 +249,10 @@ function complexoutgoingmat(rn::ReactionSystem; sparse = false)
     nps.complexoutgoingmat
 end
 
-# Creates a *sparse* outgoing reaction complex stoichiometry matrix.
+# creates a *sparse* outgoing reaction complex stoichiometry matrix
 function complexoutgoingmat(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSystem, B)
-    # Computes the I, J, and V vectors used for the sparse matrix (read about sparse matrix 
-    # representation for more information).
+    # computes the I, J, and V vectors used for the sparse matrix (read about sparse matrix 
+    # representation for more information)
     n = size(B, 2)
     rows = rowvals(B)
     vals = nonzeros(B)
@@ -260,7 +260,7 @@ function complexoutgoingmat(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSyste
     Js = Int[]
     Vs = Int[]
 
-    # Allocates space to the vectors (so that it is not done incrementally in the loop). 
+    # allocates space to the vectors (so that it is not done incrementally in the loop)
     sizehint!(Is, div(length(vals), 2))
     sizehint!(Js, div(length(vals), 2))
     sizehint!(Vs, div(length(vals), 2))
@@ -277,7 +277,7 @@ function complexoutgoingmat(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSyste
     sparse(Is, Js, Vs, size(B, 1), size(B, 2))
 end
 
-# Creates a *dense* outgoing reaction complex stoichiometry matrix.
+# creates a *dense* outgoing reaction complex stoichiometry matrix
 function complexoutgoingmat(::Type{Matrix{Int}}, rn::ReactionSystem, B)
     Δ = copy(B)
     for (I, b) in pairs(Δ)
@@ -310,7 +310,7 @@ function incidencematgraph(rn::ReactionSystem)
     nps.incidencegraph
 end
 
-# Computes the incidence graph from an *dense* incidence matrix.
+# computes the incidence graph from an *dense* incidence matrix
 function incidencematgraph(incidencemat::Matrix{Int})
     @assert all(∈([-1, 0, 1]), incidencemat)
     n = size(incidencemat, 1)  # no. of nodes/complexes
@@ -331,7 +331,7 @@ function incidencematgraph(incidencemat::Matrix{Int})
     return graph
 end
 
-# Computes the incidence graph from an *sparse* incidence matrix.
+# computes the incidence graph from an *sparse* incidence matrix
 function incidencematgraph(incidencemat::SparseMatrixCSC{Int, Int})
     @assert all(∈([-1, 0, 1]), incidencemat)
     m, n = size(incidencemat)
@@ -420,8 +420,9 @@ function terminallinkageclasses(rn::ReactionSystem)
     nps.terminallinkageclasses
 end
 
-# Helper function for terminallinkageclasses. Given a linkage class and a reaction network, say whether the linkage class is terminal, 
-# i.e. all outgoing reactions from complexes in the linkage class produce a complex also in the linkage class
+# Helper function for terminallinkageclasses. Given a linkage class and a reaction network, say 
+# whether the linkage class is terminal, i.e. all outgoing reactions from complexes in the linkage 
+# class produce a complex also in the linkage class
 function isterminal(lc::Vector, rn::ReactionSystem)
     imat = incidencemat(rn)
 
@@ -429,7 +430,8 @@ function isterminal(lc::Vector, rn::ReactionSystem)
         # Find the index of the reactant complex for a given reaction
         s = findfirst(==(-1), @view imat[:, r])
 
-        # If the reactant complex is in the linkage class, check whether the product complex is also in the linkage class. If any of them are not, return false. 
+        # If the reactant complex is in the linkage class, check whether the product complex is
+        # also in the linkage class. If any of them are not, return false. 
         if s in Set(lc)
             p = findfirst(==(1), @view imat[:, r])
             p in Set(lc) ? continue : return false
@@ -697,7 +699,7 @@ conservation laws, each represented as a row in the output.
 function conservationlaws(nsm::T; col_order = nothing) where {T <: AbstractMatrix}
 
     # compute the left nullspace over the integers
-    # The `nullspace` function updates the `col_order`.
+    # the `nullspace` function updates the `col_order`
     N = MT.nullspace(nsm'; col_order)
 
     # if all coefficients for a conservation law are negative, make positive
@@ -713,7 +715,7 @@ function conservationlaws(nsm::T; col_order = nothing) where {T <: AbstractMatri
     T(N')
 end
 
-# Used in the subsequent function.
+# used in the subsequent function
 function cache_conservationlaw_eqs!(rn::ReactionSystem, N::AbstractMatrix, col_order)
     # Retrieves nullity (the number of conservation laws). `r` is the rank of the netstoichmat.
     nullity = size(N, 1)
@@ -728,7 +730,7 @@ function cache_conservationlaw_eqs!(rn::ReactionSystem, N::AbstractMatrix, col_o
     depidxs = col_order[(r + 1):end]
     depspecs = sps[depidxs]
 
-    # Declares the conservation law parameters.
+    # declares the conservation law parameters
     constants = MT.unwrap.(MT.scalarize(only(
         @parameters $(CONSERVED_CONSTANT_SYMBOL)[1:nullity] [conserved = true])))
 
@@ -738,20 +740,20 @@ function cache_conservationlaw_eqs!(rn::ReactionSystem, N::AbstractMatrix, col_o
     conservedeqs = Equation[]
     constantdefs = Equation[]
 
-    # For each conserved quantity.
+    # for each conserved quantity
     for (i, depidx) in enumerate(depidxs)
-        # Finds the coefficient (in the conservation law) of the species that is eliminated
-        # by this conservation law.
+        # finds the coefficient (in the conservation law) of the species that is eliminated
+        # by this conservation law
         scaleby = (N[i, depidx] != 1) ? N[i, depidx] : one(eltype(N))
         (scaleby != 0) || error("Error, found a zero in the conservation law matrix where "
             * "one was not expected.")
 
-        # Creates, for this conservation law, the sum of all independent species (weighted by
-        # the ratio between the coefficient of the species and the species which is elimianted.
+        # creates, for this conservation law, the sum of all independent species (weighted by
+        # the ratio between the coefficient of the species and the species which is elimianted
         coefs = @view N[i, indepidxs]
         terms = sum(coef / scaleby * sp for (coef, sp) in zip(coefs, indepspecs))
 
-        # Computes the two equations corresponding to this conserved quantity.
+        # computes the two equations corresponding to this conserved quantity
         eq = depspecs[i] ~ constants[i] - terms
         push!(conservedeqs, eq)
         eq = constants[i] ~ depspecs[i] + terms
@@ -784,7 +786,7 @@ function conservationlaws(rs::ReactionSystem)
     nps = get_networkproperties(rs)
     !isempty(nps.conservationmat) && (return nps.conservationmat)
 
-    # If the conservation law matrix is not computed, do so and caches the result.
+    # if the conservation law matrix is not computed, do so and caches the result
     nsm = netstoichmat(rs)
     nps.conservationmat = conservationlaws(nsm; col_order = nps.col_order)
     cache_conservationlaw_eqs!(rs, nps.conservationmat, nps.col_order)

From 61fb536c87b3bb9b6a275cd89bde22cb0376b436 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 13 Jul 2024 22:32:00 -0400
Subject: [PATCH 7/9] up

---
 src/network_analysis.jl | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/network_analysis.jl b/src/network_analysis.jl
index a932cf9fac..ac8321b948 100644
--- a/src/network_analysis.jl
+++ b/src/network_analysis.jl
@@ -118,7 +118,7 @@ end
 
 # creates a *sparse* reaction complex matrix
 function reactioncomplexes(::Type{SparseMatrixCSC{Int, Int}}, rn::ReactionSystem,
-                           complextorxsmap)
+        complextorxsmap)
     # computes the I, J, and V vectors used for the sparse matrix (read about sparse matrix 
     # representation for more information)
     complexes = collect(keys(complextorxsmap))
@@ -480,7 +480,7 @@ end
 function subnetworkmapping(linkageclass, allrxs, complextorxsmap, p)
     # Finds the reactions that are part of teh sub-reaction network.
     rxinds = sort!(collect(Set(rxidx for rcidx in linkageclass
-                               for rxidx in complextorxsmap[rcidx])))
+    for rxidx in complextorxsmap[rcidx])))
     rxs = allrxs[rxinds]
     specset = Set(s for rx in rxs for s in rx.substrates if !isconstant(s))
     for rx in rxs
@@ -531,7 +531,7 @@ function subnetworks(rs::ReactionSystem)
         newrxs, newspecs, newps = subnetworkmapping(lcs[i], rxs, complextorxsmap, p)
         newname = Symbol(nameof(rs), "_", i)
         push!(subnetworks,
-              ReactionSystem(reacs, t, specs, newps; name = newname, spatial_ivs))
+            ReactionSystem(reacs, t, specs, newps; name = newname, spatial_ivs))
     end
     subnetworks
 end
@@ -745,8 +745,8 @@ function cache_conservationlaw_eqs!(rn::ReactionSystem, N::AbstractMatrix, col_o
         # finds the coefficient (in the conservation law) of the species that is eliminated
         # by this conservation law
         scaleby = (N[i, depidx] != 1) ? N[i, depidx] : one(eltype(N))
-        (scaleby != 0) || error("Error, found a zero in the conservation law matrix where "
-            * "one was not expected.")
+        (scaleby != 0) ||
+            error("Error, found a zero in the conservation law matrix where one was not expected.")
 
         # creates, for this conservation law, the sum of all independent species (weighted by
         # the ratio between the coefficient of the species and the species which is elimianted

From a4120b9797ff4cadc5f5765af71f73eee507b678 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 13 Jul 2024 23:47:51 -0400
Subject: [PATCH 8/9] merge fixes

---
 src/network_analysis.jl | 23 ++++++++++++-----------
 1 file changed, 12 insertions(+), 11 deletions(-)

diff --git a/src/network_analysis.jl b/src/network_analysis.jl
index ac8321b948..4d562159e3 100644
--- a/src/network_analysis.jl
+++ b/src/network_analysis.jl
@@ -479,11 +479,11 @@ end
 # that constitute the corresponding sub-reaction network.
 function subnetworkmapping(linkageclass, allrxs, complextorxsmap, p)
     # Finds the reactions that are part of teh sub-reaction network.
-    rxinds = sort!(collect(Set(rxidx for rcidx in linkageclass
-    for rxidx in complextorxsmap[rcidx])))
-    rxs = allrxs[rxinds]
-    specset = Set(s for rx in rxs for s in rx.substrates if !isconstant(s))
-    for rx in rxs
+    rxinds = sort!(collect(Set(
+        rxidx for rcidx in linkageclass for rxidx in complextorxsmap[rcidx])))
+    newrxs = allrxs[rxinds]
+    specset = Set(s for rx in newrxs for s in rx.substrates if !isconstant(s))
+    for rx in newrxs
         for product in rx.products
             !isconstant(product) && push!(specset, product)
         end
@@ -531,7 +531,7 @@ function subnetworks(rs::ReactionSystem)
         newrxs, newspecs, newps = subnetworkmapping(lcs[i], rxs, complextorxsmap, p)
         newname = Symbol(nameof(rs), "_", i)
         push!(subnetworks,
-            ReactionSystem(reacs, t, specs, newps; name = newname, spatial_ivs))
+            ReactionSystem(newrxs, t, newspecs, newps; name = newname, spatial_ivs))
     end
     subnetworks
 end
@@ -600,14 +600,16 @@ isweaklyreversible(rn, subnets)
 ```
 """
 function isweaklyreversible(rn::ReactionSystem, subnets)
-    im = get_networkproperties(rn).incidencemat
-    isempty(im) &&
-        error("Error, please call reactioncomplexes(rn::ReactionSystem) to ensure the incidence matrix has been cached.")
-    sparseig = issparse(im)
+    nps = get_networkproperties(rn)
+    isempty(nps.incidencemat) && reactioncomplexes(rn)
+    sparseig = issparse(nps.incidencemat)
+
     for subnet in subnets
         subnps = get_networkproperties(subnet)
         isempty(subnps.incidencemat) && reactioncomplexes(subnet; sparse = sparseig)
     end
+
+    # A network is weakly reversible if all of its subnetworks are strongly connected
     all(Graphs.is_strongly_connected ∘ incidencematgraph, subnets)
 end
 
@@ -697,7 +699,6 @@ Given the net stoichiometry matrix of a reaction system, computes a matrix of
 conservation laws, each represented as a row in the output.
 """
 function conservationlaws(nsm::T; col_order = nothing) where {T <: AbstractMatrix}
-
     # compute the left nullspace over the integers
     # the `nullspace` function updates the `col_order`
     N = MT.nullspace(nsm'; col_order)

From a0ab5da1fdf315dfcc56399bb9d3e7c337c57496 Mon Sep 17 00:00:00 2001
From: Torkel <torkel.loman@gmail.com>
Date: Sat, 13 Jul 2024 23:53:44 -0400
Subject: [PATCH 9/9] another merge thinf

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

diff --git a/src/network_analysis.jl b/src/network_analysis.jl
index 4d562159e3..74428da747 100644
--- a/src/network_analysis.jl
+++ b/src/network_analysis.jl
@@ -716,7 +716,7 @@ function conservationlaws(nsm::T; col_order = nothing) where {T <: AbstractMatri
     T(N')
 end
 
-# used in the subsequent function
+# Used in the subsequent function.
 function cache_conservationlaw_eqs!(rn::ReactionSystem, N::AbstractMatrix, col_order)
     # Retrieves nullity (the number of conservation laws). `r` is the rank of the netstoichmat.
     nullity = size(N, 1)