From cad63af4d92a5ea536a13930afee350234c2bc56 Mon Sep 17 00:00:00 2001
From: Nikhil Yewale <yewale.nikhil65@gmail.com>
Date: Thu, 26 Aug 2021 22:23:28 +0530
Subject: [PATCH 01/10] adding SparseArrays dependency

adding SparseArrays dependency for sparse representations of `stoich` matrices
---
 src/Catalyst.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index 376c17aa7f..051a70cc64 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -5,7 +5,7 @@ module Catalyst
 
 using DocStringExtensions
 using DiffEqBase, Reexport, ModelingToolkit, DiffEqJump
-
+using SparseArrays
 # ModelingToolkit imports and convenience functions we use
 using ModelingToolkit: Symbolic, value, istree, get_states, get_ps, get_iv, get_systems, 
                        get_eqs, get_defaults, toparam

From 815655158a94a29b4638421a9de138447e45d69b Mon Sep 17 00:00:00 2001
From: Nikhil Yewale <yewale.nikhil65@gmail.com>
Date: Thu, 26 Aug 2021 22:30:41 +0530
Subject: [PATCH 02/10] updating `stoichmat` functions and complex network
 functions

updating `stoichmat` functions and complex network functions for sparse representation
---
 src/networkapi.jl | 160 +++++++++++++++++++++++++++++++++++++++++-----
 1 file changed, 145 insertions(+), 15 deletions(-)

diff --git a/src/networkapi.jl b/src/networkapi.jl
index 33f71e22ee..2eda2933ae 100644
--- a/src/networkapi.jl
+++ b/src/networkapi.jl
@@ -135,12 +135,27 @@ end
 
 
 """
-    substoichmat(rn; smap=speciesmap(rn))
+    substoichmat(rn,sparsity=false; smap=speciesmap(rn))
 
 Returns the substrate stoichiometry matrix
+
+Note:
+- Set sparsity=true for sparse representation
 """
-function substoichmat(rn; smap=speciesmap(rn))
-    smat = zeros(Int,(numreactions(rn),numspecies(rn)))
+function substoichmat(::Type{SparseMatrixCSC}, rn::ReactionSystem; smap=speciesmap(rn))
+    Is=Int64[];  Js=Int64[];  Vs=Int64[];
+    for (k,rx) in enumerate(reactions(rn))
+        stoich = rx.substoich
+        for (i,sub) in enumerate(rx.substrates)
+            push!(Is,k)
+            push!(Js, smap[sub])
+            push!(Vs, stoich[i])
+        end
+    end
+    smat = sparse(Is,Js,Vs,numreactions(rn),numspecies(rn))
+end
+function substoichmat(::Type{Matrix},rn::ReactionSystem; smap=speciesmap(rn))
+    smat = zeros(Int,numreactions(rn),numspecies(rn))
     for (k,rx) in enumerate(reactions(rn))
         stoich = rx.substoich
         for (i,sub) in enumerate(rx.substrates)
@@ -149,13 +164,32 @@ function substoichmat(rn; smap=speciesmap(rn))
     end
     smat
 end
+function substoichmat(rn::ReactionSystem, sparsity::Bool=false; smap=speciesmap(rn))
+	sparsity ? substoichmat(SparseMatrixCSC, rn; smap=smap) : substoichmat(Matrix, rn; smap=smap)
+end
+
 
 """
-    prodstoichmat(rn; smap=speciesmap(rn))
+    prodstoichmat(rn,sparsity=false; smap=speciesmap(rn))
 
 Returns the product stoichiometry matrix
+
+Note:
+- Set sparsity=true for sparse representation
 """
-function prodstoichmat(rn; smap=speciesmap(rn))
+function prodstoichmat(::Type{SparseMatrixCSC},rn::ReactionSystem; smap=speciesmap(rn))
+    Is=Int64[];  Js=Int64[];  Vs=Int64[];
+    for (k,rx) in enumerate(reactions(rn))
+        stoich = rx.prodstoich
+        for (i,prod) in enumerate(rx.products)
+			push!(Is,k)
+			push!(Js, smap[prod])
+			push!(Vs, stoich[i])
+        end
+    end
+    smat = sparse(Is,Js,Vs,numreactions(rn),numspecies(rn))
+end
+function prodstoichmat(::Type{Matrix},rn::ReactionSystem; smap=speciesmap(rn))
     pmat = zeros(Int,(numreactions(rn),numspecies(rn)))
     for (k,rx) in enumerate(reactions(rn))
         stoich = rx.prodstoich
@@ -165,14 +199,31 @@ function prodstoichmat(rn; smap=speciesmap(rn))
     end
     pmat
 end
+function prodstoichmat(rn::ReactionSystem, sparsity::Bool=false; smap=speciesmap(rn))
+	sparsity ? prodstoichmat(SparseMatrixCSC, rn; smap=smap) : prodstoichmat(Matrix, rn; smap=smap)
+end
 
 
 """
-    netstoichmat(rn; smap=speciesmap(rn))
+    netstoichmat(rn,sparsity=false; smap=speciesmap(rn))
 
 Returns the net stoichiometry matrix
+
+Note:
+- Set sparsity=true for sparse representation
 """
-function netstoichmat(rn; smap=speciesmap(rn))
+function netstoichmat(::Type{SparseMatrixCSC},rn::ReactionSystem; smap=speciesmap(rn))
+    Is=Int64[];  Js=Int64[];  Vs=Int64[];
+    for (k,rx) in pairs(reactions(rn))
+        for (spec,coef) in rx.netstoich
+			push!(Is,k)
+			push!(Js, smap[spec])
+			push!(Vs, coef)
+        end
+    end
+    nmat = sparse(Is,Js,Vs,numreactions(rn),numspecies(rn))
+end
+function netstoichmat(::Type{Matrix},rn::ReactionSystem; smap=speciesmap(rn))
     nmat = zeros(Int,(numreactions(rn),numspecies(rn)))
     for (k,rx) in pairs(reactions(rn))
         for (spec,coef) in rx.netstoich
@@ -181,6 +232,9 @@ function netstoichmat(rn; smap=speciesmap(rn))
     end
     nmat
 end
+function netstoichmat(rn::ReactionSystem, sparsity::Bool=false; smap=speciesmap(rn))
+	sparsity ? netstoichmat(SparseMatrixCSC, rn; smap=smap) : netstoichmat(Matrix, rn; smap=smap)
+end
 
 
 ######################## reaction complexes and reaction rates ###############################
@@ -226,8 +280,9 @@ Base.setindex!(rc::ReactionComplex, t::ReactionComplexElement, i...) =
 Base.isless(a::ReactionComplexElement, b::ReactionComplexElement) = isless(a.speciesid, b.speciesid)
 Base.Sort.defalg(::ReactionComplex{T}) where {T <: Integer} = Base.DEFAULT_UNSTABLE
 
+
 """
-    reactioncomplexes(network, smap=speciesmap(rn))
+    reactioncomplexes(network,sparsity=false; smap=speciesmap(rn))
 
 Calculate the reaction complexes and complex incidence matrix for the given [`ReactionSystem`](@ref). 
 
@@ -238,8 +293,9 @@ Notes:
     Bᵢⱼ = -1, if the i'th complex is the substrate of the j'th reaction,
            1, if the i'th complex is the product of the j'th reaction,
            0, otherwise
+- Set sparsity=true for sparse representation of incidence matrix
 """
-function reactioncomplexes(rn; smap=speciesmap(rn))
+function reactioncomplexes(::Type{SparseMatrixCSC},rn::ReactionSystem; smap=speciesmap(rn))
     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}}}()
@@ -260,7 +316,42 @@ function reactioncomplexes(rn; smap=speciesmap(rn))
             complextorxsmap[prodrc] = [i => 1]
         end
     end
-    
+
+    complexes = collect(keys(complextorxsmap))
+
+    Is=Int64[];  Js=Int64[];  Vs=Int64[];
+	for (i,c) in enumerate(complexes)
+        for (j,σ) in complextorxsmap[c]
+			push!(Is, i)
+			push!(Js, j)
+			push!(Vs, σ)
+        end
+    end
+	B = sparse(Is,Js,Vs,length(complexes), numreactions(rn))
+    complexes,B
+end
+function reactioncomplexes(::Type{Matrix},rn::ReactionSystem; smap=speciesmap(rn))
+    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}}}()
+    for (i,rx) in enumerate(rxs)
+        reactantids = isempty(rx.substrates) ? Vector{Int}() : [smap[sub] for sub in rx.substrates]
+        subrc = sort!(ReactionComplex(reactantids, copy(rx.substoich)))
+        if haskey(complextorxsmap, subrc)
+            push!(complextorxsmap[subrc], i => -1)
+        else
+            complextorxsmap[subrc] = [i => -1]
+        end
+
+        productids = isempty(rx.products) ? Vector{Int}() : [smap[prod] for prod in rx.products]
+        prodrc = sort!(ReactionComplex(productids, copy(rx.prodstoich)))
+        if haskey(complextorxsmap, prodrc)
+            push!(complextorxsmap[prodrc], i => 1)
+        else
+            complextorxsmap[prodrc] = [i => 1]
+        end
+    end
+
     complexes = collect(keys(complextorxsmap))
     B = zeros(Int64, length(complexes), numreactions(rn));
     for (i,c) in enumerate(complexes)
@@ -270,6 +361,10 @@ function reactioncomplexes(rn; smap=speciesmap(rn))
     end
     complexes,B
 end
+function reactioncomplexes(rn::ReactionSystem,sparsity::Bool=false; smap=speciesmap(rn))
+	sparsity ? reactioncomplexes(SparseMatrixCSC,rn;smap=smap) : reactioncomplexes(Matrix,rn;smap=smap)
+end
+
 
 """
     reaction_rates(network)
@@ -282,15 +377,29 @@ end
 
 
 """
-    complexstoichmat(network; rcs=reactioncomplexes(rn)[1]))
+    complexstoichmat(network,sparsity=false; rcs=reactioncomplexes(rn)[1]))
 
 Given a [`ReactionSystem`](@ref) and vector of reaction complexes, return a
 matrix with positive entries of size num_of_species x num_of_complexes, where
 the non-zero positive entries in the kth column denote stoichiometric
 coefficients of the species participating in the kth reaction complex.
+
+Note:
+- Set sparsity=true for sparse representation
 """
-function complexstoichmat(rn; rcs=reactioncomplexes(rn)[1])
-    Z = zeros(Int64, numspecies(rn), length(rcs));
+function complexstoichmat(::Type{SparseMatrixCSC},rn::ReactionSystem; rcs=reactioncomplexes(rn)[1])
+    Is=Int64[];  Js=Int64[];  Vs=Int64[];
+    for (i,rc) in enumerate(rcs)
+        for rcel in rc
+			push!(Is,rcel.speciesid)
+			push!(Js, i)
+			push!(Vs, rcel.speciesstoich)
+        end
+    end
+    Z = sparse(Is,Js,Vs, numspecies(rn), length(rcs))
+end
+function complexstoichmat(::Type{Matrix},rn::ReactionSystem; rcs=reactioncomplexes(rn)[1])
+    Z=zeros(Int64, numspecies(rn), length(rcs))
     for (i,rc) in enumerate(rcs)
         for rcel in rc
             Z[rcel.speciesid,i] = rcel.speciesstoich
@@ -298,9 +407,12 @@ function complexstoichmat(rn; rcs=reactioncomplexes(rn)[1])
     end
     Z
 end
+function complexstoichmat(rn::ReactionSystem, sparsity::Bool=false; rcs=reactioncomplexes(rn,sparsity)[1])
+	sparsity ? complexstoichmat(SparseMatrixCSC,rn;rcs = rcs) : complexstoichmat(Matrix,rn;rcs=rcs)
+end
 
 """
-    complexoutgoingmat(network; B=reactioncomplexes(rn)[2])
+    complexoutgoingmat(network,sparsity=false; B=reactioncomplexes(rn)[2])
 
 Given a [`ReactionSystem`](@ref) and complex incidence matrix, B, return a matrix
 of size num_of_complexes x num_of_reactions.
@@ -309,14 +421,32 @@ Notes
 - The complex outgoing matrix, Δ, is defined by 
     Δᵢⱼ = 0,    if Bᵢⱼ = 1 
     Δᵢⱼ = Bᵢⱼ,  otherwise
+- Set sparsity=true for sparse representation
 """
-function complexoutgoingmat(rn; B=reactioncomplexes(rn)[2])
+function complexoutgoingmat(::Type{SparseMatrixCSC}; B=reactioncomplexes(rn)[2])
+    Δ = copy(B)
+	n = size(Δ,2)
+	rows = rowvals(Δ)
+	vals = nonzeros(Δ)
+	for j = 1:n
+	   for i in nzrange(Δ, j)
+	      if vals[i] == 1
+			  Δ[rows[i],j] = 0
+		  end
+	   end
+	end
+    dropzeros(Δ)
+end
+function complexoutgoingmat(::Type{Matrix}; B=reactioncomplexes(rn)[2])
     Δ = copy(B)
     for (I,b) in pairs(Δ)
         (b == 1) && (Δ[I] = 0)
     end
     Δ
 end
+function complexoutgoingmat(rn::ReactionSystem, sparsity::Bool=false; B=reactioncomplexes(rn,sparsity)[2])
+	sparsity ? complexoutgoingmat(SparseMatrixCSC;B = B) : complexoutgoingmat(Matrix;B=B)
+end
 
 
 ################################################################################################

From d38423610b919129c4a74db4545300569e320fb3 Mon Sep 17 00:00:00 2001
From: Nikhil Yewale <yewale.nikhil65@gmail.com>
Date: Thu, 26 Aug 2021 22:36:01 +0530
Subject: [PATCH 03/10] updating tests for sparse representation

updating tests for sparse representation of `stoichmat` and complexesmat functions
---
 test/api.jl | 44 +++++++++++++++++++++++---------------------
 1 file changed, 23 insertions(+), 21 deletions(-)

diff --git a/test/api.jl b/test/api.jl
index 4b09be7157..3bc59d1399 100644
--- a/test/api.jl
+++ b/test/api.jl
@@ -77,8 +77,8 @@ smat = [1 2 0;
         0 3 0]
 pmat = [0 2 0;
         1 0 2]
-@test all(smat .== substoichmat(rnmat))
-@test all(pmat .== prodstoichmat(rnmat))
+@test all(smat .== substoichmat(rnmat) .== Array(substoichmat(rnmat,true)))
+@test all(pmat .== prodstoichmat(rnmat) .== Array(prodstoichmat(rnmat,true)))
               
 ############## testing newly added intermediate complexes reaction networks##############
 rns  = Vector{ReactionSystem}(undef,6)
@@ -103,9 +103,9 @@ B = [-1 0 0 0;
       0 0 0 1]
 Δ = [-1 0 0 0; 0 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 0; 0 0 0 -1; 0 0 0 0]
 rcs,B2 = reactioncomplexes(rns[1])
-@test B == B2
-@test Z == complexstoichmat(rns[1])
-@test Δ == complexoutgoingmat(rns[1])
+@test B == B2 == Array(reactioncomplexes(rns[1],true)[2])
+@test Z == complexstoichmat(rns[1]) == Array(complexstoichmat(rns[1],true))
+@test Δ == complexoutgoingmat(rns[1]) == Array(complexoutgoingmat(rns[1],true))
 
 # mass-action rober
 rns[2] = @reaction_network begin
@@ -123,9 +123,9 @@ B = [-1 0 0;
       0 0 1]
 Δ = [-1 0 0; 0 0 0; 0 -1 0; 0 0 -1; 0 0 0]
 rcs,B2 = reactioncomplexes(rns[2])
-@test B == B2
-@test Z == complexstoichmat(rns[2])
-@test Δ == complexoutgoingmat(rns[2])
+@test B == B2 == Array(reactioncomplexes(rns[2],true)[2])
+@test Z == complexstoichmat(rns[2]) == Array(complexstoichmat(rns[2],true))
+@test Δ == complexoutgoingmat(rns[2]) == Array(complexoutgoingmat(rns[2],true))
 
 
 #  some rational functions as rates
@@ -145,9 +145,9 @@ B = [-1 0 0 1;
       0 0 0 -1]
 Δ = [-1 0 0 0; 0 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 -1]
 rcs,B2 = reactioncomplexes(rns[3])
-@test B == B2
-@test Z == complexstoichmat(rns[3])
-@test Δ == complexoutgoingmat(rns[3])
+@test B == B2 == Array(reactioncomplexes(rns[3],true)[2])
+@test Z == complexstoichmat(rns[3]) == Array(complexstoichmat(rns[3],true))
+@test Δ == complexoutgoingmat(rns[3]) ==  Array(complexoutgoingmat(rns[3],true))
 
 # repressilator
 rns[4]  = @reaction_network begin
@@ -184,9 +184,9 @@ B = [-1 -1 -1 1 -1 1 -1 1 -1 0 0 0 1 1 1;
        0 0 0 0 0 0 0 -1 0 0 0 -1 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1]
 rcs,B2 = reactioncomplexes(rns[4])
-@test B == B2
-@test Z == complexstoichmat(rns[4])
-@test Δ == complexoutgoingmat(rns[4])
+@test B == B2 == Array(reactioncomplexes(rns[4],true)[2])
+@test Z == complexstoichmat(rns[4]) == Array(complexstoichmat(rns[4],true))
+@test Δ == complexoutgoingmat(rns[4]) == Array(complexoutgoingmat(rns[4],true))
 
 #brusselator
 rns[5] = @reaction_network begin
@@ -204,9 +204,9 @@ B = [-1 0 0 1;
       0 0 1 0]
 Δ = [-1 0 0 0; 0 0 -1 -1; 0 -1 0 0; 0 0 0 0; 0 0 0 0]
 rcs,B2 = reactioncomplexes(rns[5])
-@test B == B2
-@test Z == complexstoichmat(rns[5])
-@test Δ == complexoutgoingmat(rns[5])
+@test B == B2 == Array(reactioncomplexes(rns[5],true)[2])
+@test Z == complexstoichmat(rns[5]) == Array(complexstoichmat(rns[5],true))
+@test Δ == complexoutgoingmat(rns[5]) == Array(complexoutgoingmat(rns[5],true))
 
 # some rational functions as rates
 rns[6] = @reaction_network begin
@@ -232,9 +232,9 @@ B = [-1 1 0 0 0 0;
       0 0 0 0 1 -1]
 Δ = [-1 0 0 0 0 0; 0 -1 0 0 0 0; 0 0 -1 0 0 0; 0 0 0 -1 0 0; 0 0 0 0 -1 0; 0 0 0 0 0 -1]
 rcs,B2 = reactioncomplexes(rns[6])
-@test all(B .== B2)
-@test Z == complexstoichmat(rns[6]; rcs=rcs)
-@test Δ == complexoutgoingmat(rns[6], B=B2)
+@test B == B2 == Array(reactioncomplexes(rns[6],true)[2])
+@test Z == complexstoichmat(rns[6]) == Array(complexstoichmat(rns[6],true))
+@test Δ == complexoutgoingmat(rns[6]) == Array(complexoutgoingmat(rns[6],true))
                      
 
 reaction_networks_standard = Vector{ReactionSystem}(undef,10)
@@ -420,8 +420,10 @@ end d v1 K1 n1 v2 K2 n2
 myrn = [reaction_networks_standard;reaction_networks_hill;reaction_networks_real]
 for i in 1:length(myrn)
     local rcs,B = reactioncomplexes(myrn[i])
+    @test B == Array(reactioncomplexes(myrn[i],true)[2])
     local Z = complexstoichmat(myrn[i])
-    @test Z*B == (netstoichmat(myrn[i])')
+    @test Z == Array(complexstoichmat(myrn[i],true))
+    @test Z*B == (netstoichmat(myrn[i])') == Array(netstoichmat(myrn[i],true)')
 end
 
 

From 85d04a2b387f701e77e09651823aa612e293ae76 Mon Sep 17 00:00:00 2001
From: Nikhil Yewale <yewale.nikhil65@gmail.com>
Date: Thu, 26 Aug 2021 22:38:45 +0530
Subject: [PATCH 04/10] update project.toml for sparsearrays

update project.toml for sparsearrays
---
 Project.toml | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/Project.toml b/Project.toml
index 375973fc68..92d6345c81 100644
--- a/Project.toml
+++ b/Project.toml
@@ -14,6 +14,7 @@ ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
@@ -43,6 +44,7 @@ SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
+SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [targets]
-test = ["Graphviz_jll", "LinearAlgebra", "OrdinaryDiffEq", "Random", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful"]
+test = ["Graphviz_jll", "LinearAlgebra", "OrdinaryDiffEq", "Random", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful","SparseArrays"]

From 0a1566db0fa2e9fb8a07e94407d7349dfc76ea9d Mon Sep 17 00:00:00 2001
From: Nikhil Yewale <yewale.nikhil65@gmail.com>
Date: Thu, 26 Aug 2021 22:39:55 +0530
Subject: [PATCH 05/10] add using SparseArrays

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

diff --git a/test/api.jl b/test/api.jl
index 3bc59d1399..302b9cf562 100644
--- a/test/api.jl
+++ b/test/api.jl
@@ -1,5 +1,5 @@
 using Catalyst, DiffEqBase, ModelingToolkit, Test
-
+using SparseArrays
 using ModelingToolkit: value
 
 @parameters t k1 k2

From 754813205e8e6ebce867ad00690fb8cf7ac82e4d Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 31 Aug 2021 10:30:54 -0400
Subject: [PATCH 06/10] updateand transpose basic stoichmats

---
 Project.toml      |  3 +-
 src/Catalyst.jl   |  4 +--
 src/networkapi.jl | 89 ++++++++++++++++++++++++-----------------------
 3 files changed, 49 insertions(+), 47 deletions(-)

diff --git a/Project.toml b/Project.toml
index 92d6345c81..45ce0fd386 100644
--- a/Project.toml
+++ b/Project.toml
@@ -44,7 +44,6 @@ SteadyStateDiffEq = "9672c7b4-1e72-59bd-8a11-6ac3964bc41f"
 StochasticDiffEq = "789caeaf-c7a9-5a7d-9973-96adeb23e2a0"
 Test = "8dfed614-e22c-5e08-85e1-65c5234f0b40"
 Unitful = "1986cc42-f94f-5a68-af5c-568840ba703d"
-SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [targets]
-test = ["Graphviz_jll", "LinearAlgebra", "OrdinaryDiffEq", "Random", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful","SparseArrays"]
+test = ["Graphviz_jll", "LinearAlgebra", "OrdinaryDiffEq", "Random", "SafeTestsets", "StableRNGs", "Statistics", "SteadyStateDiffEq", "StochasticDiffEq", "Test", "Unitful"]
diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index 051a70cc64..41ebe97b57 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -4,8 +4,8 @@ $(DocStringExtensions.README)
 module Catalyst
 
 using DocStringExtensions
-using DiffEqBase, Reexport, ModelingToolkit, DiffEqJump
-using SparseArrays
+using SparseArrays, DiffEqBase, Reexport, ModelingToolkit, DiffEqJump
+
 # ModelingToolkit imports and convenience functions we use
 using ModelingToolkit: Symbolic, value, istree, get_states, get_ps, get_iv, get_systems, 
                        get_eqs, get_defaults, toparam
diff --git a/src/networkapi.jl b/src/networkapi.jl
index 2eda2933ae..0298d7973f 100644
--- a/src/networkapi.jl
+++ b/src/networkapi.jl
@@ -135,105 +135,108 @@ end
 
 
 """
-    substoichmat(rn,sparsity=false; smap=speciesmap(rn))
+    substoichmat(rn; sparse=false, smap=speciesmap(rn))
 
-Returns the substrate stoichiometry matrix
+Returns the substrate stoichiometry matrix, S, with Sᵢⱼ the stoichiometric
+coefficient of the ith substrate within the jth reaction.
 
 Note:
-- Set sparsity=true for sparse representation
+- Set sparse=true for a sparse matrix representation
 """
-function substoichmat(::Type{SparseMatrixCSC}, rn::ReactionSystem; smap=speciesmap(rn))
-    Is=Int64[];  Js=Int64[];  Vs=Int64[];
+function substoichmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem; smap=speciesmap(rn))
+    Is=Int[];  Js=Int[];  Vs=Int[];
     for (k,rx) in enumerate(reactions(rn))
         stoich = rx.substoich
         for (i,sub) in enumerate(rx.substrates)
-            push!(Is,k)
-            push!(Js, smap[sub])
+            push!(Js, k)
+            push!(Is, smap[sub])
             push!(Vs, stoich[i])
         end
     end
-    smat = sparse(Is,Js,Vs,numreactions(rn),numspecies(rn))
+    sparse(Is,Js,Vs,numspecies(rn),numreactions(rn))
 end
-function substoichmat(::Type{Matrix},rn::ReactionSystem; smap=speciesmap(rn))
-    smat = zeros(Int,numreactions(rn),numspecies(rn))
+function substoichmat(::Type{Matrix{Int}},rn::ReactionSystem; smap=speciesmap(rn))
+    smat = zeros(Int, numspecies(rn), numreactions(rn))
     for (k,rx) in enumerate(reactions(rn))
         stoich = rx.substoich
         for (i,sub) in enumerate(rx.substrates)
-            smat[k,smap[sub]] = stoich[i]
+            smat[smap[sub],k] = stoich[i]
         end
     end
     smat
 end
-function substoichmat(rn::ReactionSystem, sparsity::Bool=false; smap=speciesmap(rn))
-	sparsity ? substoichmat(SparseMatrixCSC, rn; smap=smap) : substoichmat(Matrix, rn; smap=smap)
+function substoichmat(rn::ReactionSystem; sparse::Bool=false; smap=speciesmap(rn))
+	sparse ? substoichmat(SparseMatrixCSC{Int,Int}, rn; smap=smap) : substoichmat(Matrix{Int}, rn; smap=smap)
 end
 
 
 """
-    prodstoichmat(rn,sparsity=false; smap=speciesmap(rn))
+    prodstoichmat(rn; sparse=false, smap=speciesmap(rn))
 
-Returns the product stoichiometry matrix
+Returns the product stoichiometry matrix, P, with Pᵢⱼ the stoichiometric
+coefficient of the ith product within the jth reaction.
 
 Note:
-- Set sparsity=true for sparse representation
+- Set sparse=true for a sparse matrix representation
 """
-function prodstoichmat(::Type{SparseMatrixCSC},rn::ReactionSystem; smap=speciesmap(rn))
-    Is=Int64[];  Js=Int64[];  Vs=Int64[];
+function prodstoichmat(::Type{SparseMatrixCSC{Int}}, rn::ReactionSystem; smap=speciesmap(rn))
+    Is=Int[];  Js=Int[];  Vs=Int[];
     for (k,rx) in enumerate(reactions(rn))
         stoich = rx.prodstoich
         for (i,prod) in enumerate(rx.products)
-			push!(Is,k)
-			push!(Js, smap[prod])
+			push!(Js, k)
+			push!(Is, smap[prod])
 			push!(Vs, stoich[i])
         end
     end
-    smat = sparse(Is,Js,Vs,numreactions(rn),numspecies(rn))
+    sparse(Is,Js,Vs,numspecies(rn),numreactions(rn))
 end
-function prodstoichmat(::Type{Matrix},rn::ReactionSystem; smap=speciesmap(rn))
-    pmat = zeros(Int,(numreactions(rn),numspecies(rn)))
+function prodstoichmat(::Type{Matrix{Int}},rn::ReactionSystem; smap=speciesmap(rn))
+    pmat = zeros(Int, numspecies(rn), numreactions(rn))
     for (k,rx) in enumerate(reactions(rn))
         stoich = rx.prodstoich
         for (i,prod) in enumerate(rx.products)
-            pmat[k,smap[prod]] = stoich[i]
+            pmat[smap[prod],k] = stoich[i]
         end
     end
     pmat
 end
-function prodstoichmat(rn::ReactionSystem, sparsity::Bool=false; smap=speciesmap(rn))
-	sparsity ? prodstoichmat(SparseMatrixCSC, rn; smap=smap) : prodstoichmat(Matrix, rn; smap=smap)
+function prodstoichmat(rn::ReactionSystem; sparse=false, smap=speciesmap(rn))
+	sparse ? prodstoichmat(SparseMatrixCSC{Int,Int}, rn; smap=smap) : prodstoichmat(Matrix{Int}, rn; smap=smap)
 end
 
 
 """
     netstoichmat(rn,sparsity=false; smap=speciesmap(rn))
 
-Returns the net stoichiometry matrix
+Returns the net stoichiometry matrix, N, with Nᵢⱼ the net stoichiometric
+coefficient of the ith species within the jth reaction.
 
 Note:
-- Set sparsity=true for sparse representation
+- Set sparse=true for a sparse matrix representation
 """
-function netstoichmat(::Type{SparseMatrixCSC},rn::ReactionSystem; smap=speciesmap(rn))
-    Is=Int64[];  Js=Int64[];  Vs=Int64[];
+function netstoichmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem; smap=speciesmap(rn))
+    Is=Int[];  Js=Int[];  Vs=Int[];
     for (k,rx) in pairs(reactions(rn))
         for (spec,coef) in rx.netstoich
-			push!(Is,k)
-			push!(Js, smap[spec])
+			push!(Js, k)
+			push!(Is, smap[spec])
 			push!(Vs, coef)
         end
     end
-    nmat = sparse(Is,Js,Vs,numreactions(rn),numspecies(rn))
+    sparse(Is,Js,Vs,numspecies(rn),numreactions(rn))
 end
-function netstoichmat(::Type{Matrix},rn::ReactionSystem; smap=speciesmap(rn))
-    nmat = zeros(Int,(numreactions(rn),numspecies(rn)))
+function netstoichmat(::Type{Matrix{Int}},rn::ReactionSystem; smap=speciesmap(rn))
+    nmat = zeros(Int,numspecies(rn),numreactions(rn))
     for (k,rx) in pairs(reactions(rn))
         for (spec,coef) in rx.netstoich
-            nmat[k,smap[spec]] = coef
+            nmat[smap[spec],k] = coef
         end
     end
     nmat
 end
-function netstoichmat(rn::ReactionSystem, sparsity::Bool=false; smap=speciesmap(rn))
-	sparsity ? netstoichmat(SparseMatrixCSC, rn; smap=smap) : netstoichmat(Matrix, rn; smap=smap)
+function netstoichmat(rn::ReactionSystem; sparse=false, smap=speciesmap(rn))
+	sparse ? netstoichmat(SparseMatrixCSC{Int,Int}, rn; smap=smap) : netstoichmat(Matrix{Int}, rn; smap=smap)
 end
 
 
@@ -295,7 +298,7 @@ Notes:
            0, otherwise
 - Set sparsity=true for sparse representation of incidence matrix
 """
-function reactioncomplexes(::Type{SparseMatrixCSC},rn::ReactionSystem; smap=speciesmap(rn))
+function reactioncomplexes(::Type{SparseMatrixCSC{Int,Int}},rn::ReactionSystem; smap=speciesmap(rn))
     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}}}()
@@ -327,10 +330,10 @@ function reactioncomplexes(::Type{SparseMatrixCSC},rn::ReactionSystem; smap=spec
 			push!(Vs, σ)
         end
     end
-	B = sparse(Is,Js,Vs,length(complexes), numreactions(rn))
+	B = sparse(Is,Js,Vs,length(complexes),numreactions(rn))
     complexes,B
 end
-function reactioncomplexes(::Type{Matrix},rn::ReactionSystem; smap=speciesmap(rn))
+function reactioncomplexes(::Type{Matrix{Int}},rn::ReactionSystem; smap=speciesmap(rn))
     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}}}()
@@ -361,8 +364,8 @@ function reactioncomplexes(::Type{Matrix},rn::ReactionSystem; smap=speciesmap(rn
     end
     complexes,B
 end
-function reactioncomplexes(rn::ReactionSystem,sparsity::Bool=false; smap=speciesmap(rn))
-	sparsity ? reactioncomplexes(SparseMatrixCSC,rn;smap=smap) : reactioncomplexes(Matrix,rn;smap=smap)
+function reactioncomplexes(rn::ReactionSystem; sparse=false; smap=speciesmap(rn))
+	sparse ? reactioncomplexes(SparseMatrixCSC{Int,Int},rn;smap=smap) : reactioncomplexes(Matrix{Int},rn;smap=smap)
 end
 
 

From 636e88f6f9f21779cc2d50f0e9764c349195c4f5 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 31 Aug 2021 11:46:12 -0400
Subject: [PATCH 07/10] update complex mat funs

---
 docs/src/api/catalyst_api.md |   1 +
 src/networkapi.jl            | 127 ++++++++++++++++++-----------------
 2 files changed, 66 insertions(+), 62 deletions(-)

diff --git a/docs/src/api/catalyst_api.md b/docs/src/api/catalyst_api.md
index 293b1412e4..c582f8222b 100644
--- a/docs/src/api/catalyst_api.md
+++ b/docs/src/api/catalyst_api.md
@@ -121,6 +121,7 @@ conservationlaws
 conservedquantities
 ReactionComplexElement
 ReactionComplex
+reactioncomplexmap
 reactioncomplexes
 complexstoichmat
 complexoutgoingmat
diff --git a/src/networkapi.jl b/src/networkapi.jl
index 0298d7973f..4214cb0337 100644
--- a/src/networkapi.jl
+++ b/src/networkapi.jl
@@ -283,22 +283,20 @@ Base.setindex!(rc::ReactionComplex, t::ReactionComplexElement, i...) =
 Base.isless(a::ReactionComplexElement, b::ReactionComplexElement) = isless(a.speciesid, b.speciesid)
 Base.Sort.defalg(::ReactionComplex{T}) where {T <: Integer} = Base.DEFAULT_UNSTABLE
 
-
 """
-    reactioncomplexes(network,sparsity=false; smap=speciesmap(rn))
+    reactioncomplexmap(rn::ReactionSystem; smap=speciesmap(rn))
 
-Calculate the reaction complexes and complex incidence matrix for the given [`ReactionSystem`](@ref). 
+Find each [`ReactionComplex`](@ref) within the specified system, constructing a
+mapping from the complex to vectors that indicate which reactions it appears in
+as substrates and products.
 
 Notes:
-- returns a pair of a vector of [`ReactionComplex`](@ref)s and the complex incidence matrix.
-- An empty [`ReactionComplex`](@ref) denotes the null (∅) state (from reactions like ∅ -> A or A -> ∅).
-- The complex incidence matrix, B, is number of complexes by number of reactions with
-    Bᵢⱼ = -1, if the i'th complex is the substrate of the j'th reaction,
-           1, if the i'th complex is the product of the j'th reaction,
-           0, otherwise
-- Set sparsity=true for sparse representation of incidence matrix
+- Each [`ReactionComplex`](@ref) is mapped to a vector of pairs, with each pair
+  having the form `reactionidx => ± 1`, where `1` indicates the complex appears
+  as a substrate and `+1` as a product in the reaction with integer label
+  `reactionidx`.
 """
-function reactioncomplexes(::Type{SparseMatrixCSC{Int,Int}},rn::ReactionSystem; smap=speciesmap(rn))
+function reactioncomplexmap(rn::ReactionSystem; smap=speciesmap(rn))
     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}}}()
@@ -319,10 +317,28 @@ function reactioncomplexes(::Type{SparseMatrixCSC{Int,Int}},rn::ReactionSystem;
             complextorxsmap[prodrc] = [i => 1]
         end
     end
+    complextorxsmap
+end
 
-    complexes = collect(keys(complextorxsmap))
 
-    Is=Int64[];  Js=Int64[];  Vs=Int64[];
+"""
+    reactioncomplexes(network; sparse=false, smap=speciesmap(rn))
+
+Calculate the reaction complexes and complex incidence matrix for the given [`ReactionSystem`](@ref). 
+
+Notes:
+- returns a pair of a vector of [`ReactionComplex`](@ref)s and the complex incidence matrix.
+- An empty [`ReactionComplex`](@ref) denotes the null (∅) state (from reactions like ∅ -> A or A -> ∅).
+- The complex incidence matrix, B, is number of complexes by number of reactions with
+    Bᵢⱼ = -1, if the i'th complex is the substrate of the j'th reaction,
+           1, if the i'th complex is the product of the j'th reaction,
+           0, otherwise
+- Set sparse=true for a sparse matrix representation of the incidence matrix
+"""
+function reactioncomplexes(::Type{SparseMatrixCSC{Int,Int}},rn::ReactionSystem; 
+                            smap=speciesmap(rn), complextorxsmap=reactioncomplexmap(rn; smap=smap))
+    complexes = collect(keys(complextorxsmap))
+    Is=Int[];  Js=Int[];  Vs=Int[];
 	for (i,c) in enumerate(complexes)
         for (j,σ) in complextorxsmap[c]
 			push!(Is, i)
@@ -333,30 +349,10 @@ function reactioncomplexes(::Type{SparseMatrixCSC{Int,Int}},rn::ReactionSystem;
 	B = sparse(Is,Js,Vs,length(complexes),numreactions(rn))
     complexes,B
 end
-function reactioncomplexes(::Type{Matrix{Int}},rn::ReactionSystem; smap=speciesmap(rn))
-    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}}}()
-    for (i,rx) in enumerate(rxs)
-        reactantids = isempty(rx.substrates) ? Vector{Int}() : [smap[sub] for sub in rx.substrates]
-        subrc = sort!(ReactionComplex(reactantids, copy(rx.substoich)))
-        if haskey(complextorxsmap, subrc)
-            push!(complextorxsmap[subrc], i => -1)
-        else
-            complextorxsmap[subrc] = [i => -1]
-        end
-
-        productids = isempty(rx.products) ? Vector{Int}() : [smap[prod] for prod in rx.products]
-        prodrc = sort!(ReactionComplex(productids, copy(rx.prodstoich)))
-        if haskey(complextorxsmap, prodrc)
-            push!(complextorxsmap[prodrc], i => 1)
-        else
-            complextorxsmap[prodrc] = [i => 1]
-        end
-    end
-
+function reactioncomplexes(::Type{Matrix{Int}},rn::ReactionSystem; 
+                            smap=speciesmap(rn), complextorxsmap=reactioncomplexmap(rn; smap=smap))
     complexes = collect(keys(complextorxsmap))
-    B = zeros(Int64, length(complexes), numreactions(rn));
+    B = zeros(Int, length(complexes), numreactions(rn));
     for (i,c) in enumerate(complexes)
         for (j,σ) in complextorxsmap[c]
             B[i,j] = σ
@@ -365,7 +361,7 @@ function reactioncomplexes(::Type{Matrix{Int}},rn::ReactionSystem; smap=speciesm
     complexes,B
 end
 function reactioncomplexes(rn::ReactionSystem; sparse=false; smap=speciesmap(rn))
-	sparse ? reactioncomplexes(SparseMatrixCSC{Int,Int},rn;smap=smap) : reactioncomplexes(Matrix{Int},rn;smap=smap)
+	sparse ? reactioncomplexes(SparseMatrixCSC{Int,Int}, rn; smap=smap) : reactioncomplexes(Matrix{Int}, rn; smap=smap)
 end
 
 
@@ -380,7 +376,7 @@ end
 
 
 """
-    complexstoichmat(network,sparsity=false; rcs=reactioncomplexes(rn)[1]))
+    complexstoichmat(network; sparse=false, rcs=keys(reactioncomplexmap(rn)))
 
 Given a [`ReactionSystem`](@ref) and vector of reaction complexes, return a
 matrix with positive entries of size num_of_species x num_of_complexes, where
@@ -388,21 +384,23 @@ the non-zero positive entries in the kth column denote stoichiometric
 coefficients of the species participating in the kth reaction complex.
 
 Note:
-- Set sparsity=true for sparse representation
+- `rcs` correspond to an iterable of the `ReactionComplexes`, i.e.
+  `rcs=keys(reactioncomplexmap(rn))` or `reactioncomplexes(rn)[1]`.
+- Set sparse=true for a sparse matrix representation
 """
-function complexstoichmat(::Type{SparseMatrixCSC},rn::ReactionSystem; rcs=reactioncomplexes(rn)[1])
-    Is=Int64[];  Js=Int64[];  Vs=Int64[];
+function complexstoichmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem; rcs=keys(reactioncomplexmap(rn)))
+    Is=Int[];  Js=Int[];  Vs=Int[];
     for (i,rc) in enumerate(rcs)
         for rcel in rc
-			push!(Is,rcel.speciesid)
+			push!(Is, rcel.speciesid)
 			push!(Js, i)
 			push!(Vs, rcel.speciesstoich)
         end
     end
     Z = sparse(Is,Js,Vs, numspecies(rn), length(rcs))
 end
-function complexstoichmat(::Type{Matrix},rn::ReactionSystem; rcs=reactioncomplexes(rn)[1])
-    Z=zeros(Int64, numspecies(rn), length(rcs))
+function complexstoichmat(::Type{Matrix{Int}}, rn::ReactionSystem; rcs=keys(reactioncomplexmap(rn)))
+    Z=zeros(Int, numspecies(rn), length(rcs))
     for (i,rc) in enumerate(rcs)
         for rcel in rc
             Z[rcel.speciesid,i] = rcel.speciesstoich
@@ -410,45 +408,50 @@ function complexstoichmat(::Type{Matrix},rn::ReactionSystem; rcs=reactioncomplex
     end
     Z
 end
-function complexstoichmat(rn::ReactionSystem, sparsity::Bool=false; rcs=reactioncomplexes(rn,sparsity)[1])
-	sparsity ? complexstoichmat(SparseMatrixCSC,rn;rcs = rcs) : complexstoichmat(Matrix,rn;rcs=rcs)
+function complexstoichmat(rn::ReactionSystem; sparse=false, rcs=keys(reactioncomplexmap(rn)))
+	sparse ? complexstoichmat(SparseMatrixCSC{Int,Int}, rn; rcs=rcs) : complexstoichmat(Matrix{Int}, rn; rcs=rcs)
 end
 
 """
-    complexoutgoingmat(network,sparsity=false; B=reactioncomplexes(rn)[2])
+    complexoutgoingmat(network; sparse=false, B=reactioncomplexes(rn)[2])
 
 Given a [`ReactionSystem`](@ref) and complex incidence matrix, B, return a matrix
-of size num_of_complexes x num_of_reactions.
+of size num_of_complexes x num_of_reactions that identifies substrate complexes.
 
 Notes
 - The complex outgoing matrix, Δ, is defined by 
     Δᵢⱼ = 0,    if Bᵢⱼ = 1 
     Δᵢⱼ = Bᵢⱼ,  otherwise
-- Set sparsity=true for sparse representation
+- Set sparse=true for a sparse matrix representation
 """
-function complexoutgoingmat(::Type{SparseMatrixCSC}; B=reactioncomplexes(rn)[2])
-    Δ = copy(B)
-	n = size(Δ,2)
-	rows = rowvals(Δ)
-	vals = nonzeros(Δ)
+function complexoutgoingmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem; B=reactioncomplexes(rn,sparse=true)[2])    
+    n = size(B,2)
+	rows = rowvals(B)
+	vals = nonzeros(B)
+    Is = Int[]; Js = Int[]; Vs = Int[]
+    sizehint!(Is, div(length(vals),2))
+    sizehint!(Js, div(length(vals),2))
+    sizehint!(Vs, div(length(vals),2))
 	for j = 1:n
-	   for i in nzrange(Δ, j)
-	      if vals[i] == 1
-			  Δ[rows[i],j] = 0
-		  end
+	   for i in nzrange(B, j)
+	      if vals[i] != one(eltype(vals)) 
+            push!(Is, rows[i])
+            push!(Js, j) 
+            push!(Vs, vals[i])
+          end
 	   end
 	end
-    dropzeros(Δ)
+    sparse(Is,Js,Vs,size(B,1),size(B,2))
 end
-function complexoutgoingmat(::Type{Matrix}; B=reactioncomplexes(rn)[2])
+function complexoutgoingmat(::Type{Matrix{Int}}, rn::ReactionSystem; B=reactioncomplexes(rn,sparse=false)[2])
     Δ = copy(B)
     for (I,b) in pairs(Δ)
         (b == 1) && (Δ[I] = 0)
     end
     Δ
 end
-function complexoutgoingmat(rn::ReactionSystem, sparsity::Bool=false; B=reactioncomplexes(rn,sparsity)[2])
-	sparsity ? complexoutgoingmat(SparseMatrixCSC;B = B) : complexoutgoingmat(Matrix;B=B)
+function complexoutgoingmat(rn::ReactionSystem; sparse=false, B=reactioncomplexes(rn,sparse=sparse)[2])
+	sparse ? complexoutgoingmat(SparseMatrixCSC{Int,Int}, rn; B=B) : complexoutgoingmat(Matrix{Int}, rn; B=B)
 end
 
 

From 218bf11dca2c150fd16f4002d5da5f56b528d018 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 31 Aug 2021 12:02:04 -0400
Subject: [PATCH 08/10] update tests

---
 src/networkapi.jl |  5 +++--
 test/api.jl       | 48 +++++++++++++++++++++++------------------------
 2 files changed, 27 insertions(+), 26 deletions(-)

diff --git a/src/networkapi.jl b/src/networkapi.jl
index 4214cb0337..885d2829c1 100644
--- a/src/networkapi.jl
+++ b/src/networkapi.jl
@@ -465,11 +465,12 @@ 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::AbstractMatrix)::Matrix
-    n_reac, n_spec = size(nsm)
+    n_spec,n_reac = size(nsm)
     
     # We basically have to compute the left null space of the matrix
     # over the integers; this is best done using its Smith Normal Form.
-    nsm_conv = AbstractAlgebra.matrix(AbstractAlgebra.ZZ, nsm)
+    # note, we transpose as this was written when netstoichmat was reac by spec
+    nsm_conv = AbstractAlgebra.matrix(AbstractAlgebra.ZZ, nsm')
     S, T, U = AbstractAlgebra.snf_with_transform(nsm_conv)
     
     # Zero columns of S (which occur after nonzero columns in SNF)
diff --git a/test/api.jl b/test/api.jl
index 302b9cf562..b1c9647df8 100644
--- a/test/api.jl
+++ b/test/api.jl
@@ -77,8 +77,8 @@ smat = [1 2 0;
         0 3 0]
 pmat = [0 2 0;
         1 0 2]
-@test all(smat .== substoichmat(rnmat) .== Array(substoichmat(rnmat,true)))
-@test all(pmat .== prodstoichmat(rnmat) .== Array(prodstoichmat(rnmat,true)))
+@test smat == substoichmat(rnmat) == Matrix(substoichmat(rnmat, sparse=true))
+@test pmat == prodstoichmat(rnmat) == Matrix(prodstoichmat(rnmat, sparse=true))
               
 ############## testing newly added intermediate complexes reaction networks##############
 rns  = Vector{ReactionSystem}(undef,6)
@@ -103,9 +103,9 @@ B = [-1 0 0 0;
       0 0 0 1]
 Δ = [-1 0 0 0; 0 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 0; 0 0 0 -1; 0 0 0 0]
 rcs,B2 = reactioncomplexes(rns[1])
-@test B == B2 == Array(reactioncomplexes(rns[1],true)[2])
-@test Z == complexstoichmat(rns[1]) == Array(complexstoichmat(rns[1],true))
-@test Δ == complexoutgoingmat(rns[1]) == Array(complexoutgoingmat(rns[1],true))
+@test B == B2 == Matrix(reactioncomplexes(rns[1], sparse=true)[2])
+@test Z == complexstoichmat(rns[1]) == Matrix(complexstoichmat(rns[1], sparse=true))
+@test Δ == complexoutgoingmat(rns[1]) == Matrix(complexoutgoingmat(rns[1], sparse=true))
 
 # mass-action rober
 rns[2] = @reaction_network begin
@@ -123,9 +123,9 @@ B = [-1 0 0;
       0 0 1]
 Δ = [-1 0 0; 0 0 0; 0 -1 0; 0 0 -1; 0 0 0]
 rcs,B2 = reactioncomplexes(rns[2])
-@test B == B2 == Array(reactioncomplexes(rns[2],true)[2])
-@test Z == complexstoichmat(rns[2]) == Array(complexstoichmat(rns[2],true))
-@test Δ == complexoutgoingmat(rns[2]) == Array(complexoutgoingmat(rns[2],true))
+@test B == B2 == Matrix(reactioncomplexes(rns[2], sparse=true)[2])
+@test Z == complexstoichmat(rns[2]) == Matrix(complexstoichmat(rns[2], sparse=true))
+@test Δ == complexoutgoingmat(rns[2]) == Matrix(complexoutgoingmat(rns[2], sparse=true))
 
 
 #  some rational functions as rates
@@ -145,9 +145,9 @@ B = [-1 0 0 1;
       0 0 0 -1]
 Δ = [-1 0 0 0; 0 0 0 0; 0 -1 0 0; 0 0 -1 0; 0 0 0 -1]
 rcs,B2 = reactioncomplexes(rns[3])
-@test B == B2 == Array(reactioncomplexes(rns[3],true)[2])
-@test Z == complexstoichmat(rns[3]) == Array(complexstoichmat(rns[3],true))
-@test Δ == complexoutgoingmat(rns[3]) ==  Array(complexoutgoingmat(rns[3],true))
+@test B == B2 == Matrix(reactioncomplexes(rns[3], sparse=true)[2])
+@test Z == complexstoichmat(rns[3]) == Matrix(complexstoichmat(rns[3], sparse=true))
+@test Δ == complexoutgoingmat(rns[3]) ==  Matrix(complexoutgoingmat(rns[3], sparse=true))
 
 # repressilator
 rns[4]  = @reaction_network begin
@@ -184,9 +184,9 @@ B = [-1 -1 -1 1 -1 1 -1 1 -1 0 0 0 1 1 1;
        0 0 0 0 0 0 0 -1 0 0 0 -1 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0;
        0 0 0 0 0 0 0 0 0 0 0 0 -1 0 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 -1 0; 0 0 0 0 0 0 0 0 0 0 0 0 0 0 -1]
 rcs,B2 = reactioncomplexes(rns[4])
-@test B == B2 == Array(reactioncomplexes(rns[4],true)[2])
-@test Z == complexstoichmat(rns[4]) == Array(complexstoichmat(rns[4],true))
-@test Δ == complexoutgoingmat(rns[4]) == Array(complexoutgoingmat(rns[4],true))
+@test B == B2 == Matrix(reactioncomplexes(rns[4], sparse=true)[2])
+@test Z == complexstoichmat(rns[4]) == Matrix(complexstoichmat(rns[4], sparse=true))
+@test Δ == complexoutgoingmat(rns[4]) == Matrix(complexoutgoingmat(rns[4], sparse=true))
 
 #brusselator
 rns[5] = @reaction_network begin
@@ -204,9 +204,9 @@ B = [-1 0 0 1;
       0 0 1 0]
 Δ = [-1 0 0 0; 0 0 -1 -1; 0 -1 0 0; 0 0 0 0; 0 0 0 0]
 rcs,B2 = reactioncomplexes(rns[5])
-@test B == B2 == Array(reactioncomplexes(rns[5],true)[2])
-@test Z == complexstoichmat(rns[5]) == Array(complexstoichmat(rns[5],true))
-@test Δ == complexoutgoingmat(rns[5]) == Array(complexoutgoingmat(rns[5],true))
+@test B == B2 == Matrix(reactioncomplexes(rns[5], sparse=true)[2])
+@test Z == complexstoichmat(rns[5]) == Matrix(complexstoichmat(rns[5], sparse=true))
+@test Δ == complexoutgoingmat(rns[5]) == Matrix(complexoutgoingmat(rns[5], sparse=true))
 
 # some rational functions as rates
 rns[6] = @reaction_network begin
@@ -232,15 +232,15 @@ B = [-1 1 0 0 0 0;
       0 0 0 0 1 -1]
 Δ = [-1 0 0 0 0 0; 0 -1 0 0 0 0; 0 0 -1 0 0 0; 0 0 0 -1 0 0; 0 0 0 0 -1 0; 0 0 0 0 0 -1]
 rcs,B2 = reactioncomplexes(rns[6])
-@test B == B2 == Array(reactioncomplexes(rns[6],true)[2])
-@test Z == complexstoichmat(rns[6]) == Array(complexstoichmat(rns[6],true))
-@test Δ == complexoutgoingmat(rns[6]) == Array(complexoutgoingmat(rns[6],true))
+@test B == B2 == Matrix(reactioncomplexes(rns[6], sparse=true)[2])
+@test Z == complexstoichmat(rns[6]) == Matrix(complexstoichmat(rns[6], sparse=true))
+@test Δ == complexoutgoingmat(rns[6]) == Matrix(complexoutgoingmat(rns[6], sparse=true))
                      
 
 reaction_networks_standard = Vector{ReactionSystem}(undef,10)
 reaction_networks_hill = Vector{ReactionSystem}(undef,10)
 reaction_networks_real = Vector{ReactionSystem}(undef,3)
-### some more networks to test whether Z*B == netstoichmat(rn)' or not. ###
+### some more networks to test whether Z*B == netstoichmat(rn) or not. ###
 reaction_networks_standard[1] = @reaction_network begin
     (p1,p2,p3), ∅ → (X1,X2,X3)
     (k1,k2), X2 ⟷ X1 + 2X3
@@ -420,10 +420,10 @@ end d v1 K1 n1 v2 K2 n2
 myrn = [reaction_networks_standard;reaction_networks_hill;reaction_networks_real]
 for i in 1:length(myrn)
     local rcs,B = reactioncomplexes(myrn[i])
-    @test B == Array(reactioncomplexes(myrn[i],true)[2])
+    @test B == Matrix(reactioncomplexes(myrn[i], sparse=true)[2])
     local Z = complexstoichmat(myrn[i])
-    @test Z == Array(complexstoichmat(myrn[i],true))
-    @test Z*B == (netstoichmat(myrn[i])') == Array(netstoichmat(myrn[i],true)')
+    @test Z == Matrix(complexstoichmat(myrn[i], sparse=true))
+    @test Z*B == netstoichmat(myrn[i]) == Matrix(netstoichmat(myrn[i], sparse=true))
 end
 
 

From f3f7eb888a66971917fd9402dbe3962c63dfa880 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 31 Aug 2021 12:46:39 -0400
Subject: [PATCH 09/10] fix tests

---
 src/networkapi.jl |  6 +++---
 test/api.jl       | 10 ++++++----
 2 files changed, 9 insertions(+), 7 deletions(-)

diff --git a/src/networkapi.jl b/src/networkapi.jl
index 885d2829c1..4436e1453b 100644
--- a/src/networkapi.jl
+++ b/src/networkapi.jl
@@ -165,7 +165,7 @@ function substoichmat(::Type{Matrix{Int}},rn::ReactionSystem; smap=speciesmap(rn
     end
     smat
 end
-function substoichmat(rn::ReactionSystem; sparse::Bool=false; smap=speciesmap(rn))
+function substoichmat(rn::ReactionSystem; sparse::Bool=false, smap=speciesmap(rn))
 	sparse ? substoichmat(SparseMatrixCSC{Int,Int}, rn; smap=smap) : substoichmat(Matrix{Int}, rn; smap=smap)
 end
 
@@ -179,7 +179,7 @@ coefficient of the ith product within the jth reaction.
 Note:
 - Set sparse=true for a sparse matrix representation
 """
-function prodstoichmat(::Type{SparseMatrixCSC{Int}}, rn::ReactionSystem; smap=speciesmap(rn))
+function prodstoichmat(::Type{SparseMatrixCSC{Int,Int}}, rn::ReactionSystem; smap=speciesmap(rn))
     Is=Int[];  Js=Int[];  Vs=Int[];
     for (k,rx) in enumerate(reactions(rn))
         stoich = rx.prodstoich
@@ -360,7 +360,7 @@ function reactioncomplexes(::Type{Matrix{Int}},rn::ReactionSystem;
     end
     complexes,B
 end
-function reactioncomplexes(rn::ReactionSystem; sparse=false; smap=speciesmap(rn))
+function reactioncomplexes(rn::ReactionSystem; sparse=false, smap=speciesmap(rn))
 	sparse ? reactioncomplexes(SparseMatrixCSC{Int,Int}, rn; smap=smap) : reactioncomplexes(Matrix{Int}, rn; smap=smap)
 end
 
diff --git a/test/api.jl b/test/api.jl
index b1c9647df8..42ac64718c 100644
--- a/test/api.jl
+++ b/test/api.jl
@@ -73,10 +73,12 @@ rnmat = @reaction_network begin
        β, 3I --> 2R + S
    end α β
 
-smat = [1 2 0;
-        0 3 0]
-pmat = [0 2 0;
-        1 0 2]
+smat = [1 0;
+        2 3;
+        0 0]
+pmat = [0 1;
+        2 0;
+        0 2]
 @test smat == substoichmat(rnmat) == Matrix(substoichmat(rnmat, sparse=true))
 @test pmat == prodstoichmat(rnmat) == Matrix(prodstoichmat(rnmat, sparse=true))
               

From 6a2b79eafac8376586c5483491d8980d5da7dd17 Mon Sep 17 00:00:00 2001
From: Sam Isaacson <isaacsas@users.noreply.github.com>
Date: Tue, 31 Aug 2021 12:48:27 -0400
Subject: [PATCH 10/10] export reactioncomplexmap

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

diff --git a/src/Catalyst.jl b/src/Catalyst.jl
index 41ebe97b57..4017f0e03d 100644
--- a/src/Catalyst.jl
+++ b/src/Catalyst.jl
@@ -56,7 +56,7 @@ export species, params, reactions, speciesmap, paramsmap, numspecies, numreactio
 export make_empty_network, addspecies!, addparam!, addreaction!
 export dependants, dependents, substoichmat, prodstoichmat, netstoichmat
 export conservationlaws, conservedquantities
-export reactioncomplexes, reactionrates, complexstoichmat, complexoutgoingmat
+export reactioncomplexmap, reactioncomplexes, reactionrates, complexstoichmat, complexoutgoingmat
     
 # for Latex printing of ReactionSystems
 include("latexify_recipes.jl")