simplify iscomplexbalanced, rename adjacencymat, adjust implementation to be more similar

Author: vyudu

src/network_analysis.jl

@@ -206,9 +206,6 @@ function laplacianmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false)
     D * K
 end
 
-Base.zero(::Type{Union{R, Symbolics.BasicSymbolic{Real}}}) where R <: Real = zero(R)
-Base.one(::Type{Union{R, Symbolics.BasicSymbolic{Real}}}) where R <: Real = one(R)
-
 @doc raw"""
     fluxmat(rn::ReactionSystem, pmap = Dict(); sparse=false)
 
@@ -225,8 +222,6 @@ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false)
     end
 
     rcmap = reactioncomplexmap(rn)
-    nc = length(rcmap)
-    nr = length(rates)
     mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates)
     if sparse
         return fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates)
@@ -936,7 +931,7 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict; abstol=0, re
     complexes, D = reactioncomplexes(rs)
     img = incidencematgraph(rs)
     undir_img = SimpleGraph(incidencematgraph(rs))
-    K = ratematrix(rs, pmap)
+    K = adjacencymat(rs, pmap)
 
     spanning_forest = Graphs.kruskal_mst(undir_img)
     outofforest_edges = setdiff(collect(edges(undir_img)), spanning_forest)
@@ -993,52 +988,26 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap)
 end
 
 """
-    iscomplexbalanced(rs::ReactionSystem, parametermap)
+    iscomplexbalanced(rs::ReactionSystem, parametermap; sparse = false)
 
 Constructively compute whether a network will have complex-balanced equilibrium
 solutions, following the method in van der Schaft et al., [2015](https://link.springer.com/article/10.1007/s10910-015-0498-2#Sec3). 
-Accepts a dictionary, vector, or tuple of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
-"""
-function iscomplexbalanced(rs::ReactionSystem, parametermap::Dict)
-    if length(parametermap) != numparams(rs)
-        error("Incorrect number of parameters specified.")
-    end
 
-    pmap = symmap_to_varmap(rs, parametermap)
-    pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap)
-
-    sm = speciesmap(rs)
-    cm = reactioncomplexmap(rs)
-    complexes, D = reactioncomplexes(rs)
+Requires the specification of values for the parameters/rate constants. Accepts a dictionary, vector, or tuple of parameter-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
+"""
+function iscomplexbalanced(rs::ReactionSystem, parametermap::Dict; sparse = false)
     rxns = reactions(rs)
-    nc = length(complexes)
-    nr = numreactions(rs)
-    nm = numspecies(rs)
-
     if !all(r -> ismassaction(r, rs), rxns)
         error("The supplied ReactionSystem has reactions that are not ismassaction. Testing for being complex balanced is currently only supported for pure mass action networks.")
     end
 
-    rates = [substitute(rate, pmap) for rate in reactionrates(rs)]
-
-    # Construct kinetic matrix, K
-    K = zeros(nr, nc)
-    for c in 1:nc
-        complex = complexes[c]
-        for (r, dir) in cm[complex]
-            rxn = rxns[r]
-            if dir == -1
-                K[r, c] = rates[r]
-            end
-        end
-    end
-
-    L = -D * K
-    S = netstoichmat(rs)
+    L = laplacianmat(rs, parametermap; sparse)
+    D = incidencemat(rs; sparse)
+    S = netstoichmat(rs; sparse)
 
     # Compute ρ using the matrix-tree theorem
     g = incidencematgraph(rs)
-    R = ratematrix(rs, rates)
+    R = adjacencymat(rs, parametermap; sparse)
     ρ = matrixtree(g, R)
 
     # Determine if 1) ρ is positive and 2) D^T Ln ρ lies in the image of S^T
@@ -1069,50 +1038,74 @@ function iscomplexbalanced(rs::ReactionSystem, parametermap)
 end
 
 """
-    ratematrix(rs::ReactionSystem, parametermap)
+    adjacencymat(rs::ReactionSystem, pmap = Dict(); sparse = false)
 
     Given a reaction system with n complexes, outputs an n-by-n matrix where R_{ij} is the rate 
     constant of the reaction between complex i and complex j. Accepts a dictionary, vector, or tuple 
     of variable-to-value mappings, e.g. [k1 => 1.0, k2 => 2.0,...]. 
+
+    Equivalent to the adjacency matrix of the weighted graph whose weights are the 
+    reaction rates. 
+    Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap. 
+
+    The difference between this matrix and [`fluxmat`](@ref) is quite subtle. The non-zero entries to both matrices are the rate constants. However:
+        - In `fluxmat`, the rows represent complexes and columns represent reactions, and an entry (c, r) is non-zero if reaction r's source complex is c. 
+        - In `adjacencymat`, the rows and columns both represent complexes, and an entry (c1, c2) is non-zero if there is a reaction c1 --> c2.
 """
-function ratematrix(rs::ReactionSystem, rates::Vector{Float64})
-    complexes, D = reactioncomplexes(rs)
-    n = length(complexes)
-    rxns = reactions(rs)
-    ratematrix = zeros(n, n)
+function adjacencymat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false)
+    rates = if isempty(pmap)
+        reactionrates(rn)
+    else
+        substitutevals(rn, pmap, parameters(rn), reactionrates(rn))
+    end
+    mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates)
 
-    for r in 1:length(rxns)
-        rxn = rxns[r]
-        s = findfirst(==(-1), @view D[:, r])
-        p = findfirst(==(1), @view D[:, r])
-        ratematrix[s, p] = rates[r]
+    if sparse
+        return adjacencymat(SparseMatrixCSC{mtype, Int}, incidencemat(rn), rates)
+    else
+        return adjacencymat(Matrix{mtype}, incidencemat(rn), rates)
     end
-    ratematrix
 end
 
-function ratematrix(rs::ReactionSystem, parametermap::Dict)
-    if length(parametermap) != numparams(rs)
-        error("Incorrect number of parameters specified.")
+function adjacencymat(::Type{SparseMatrixCSC{T, Int}}, D, rates) where T
+    Is = Int[]
+    Js = Int[]
+    Vs = T[]
+    nc = size(D, 1)
+
+    for r in 1:length(rates)
+        s = findfirst(==(-1), @view D[:, r])
+        p = findfirst(==(1), @view D[:, r])
+        push!(Is, s)
+        push!(Js, p)
+        push!(Vs, rates[r])
     end
+    A = sparse(Is, Js, Vs, nc, nc)
+end
 
-    pmap = symmap_to_varmap(rs, parametermap)
-    pmap = Dict(ModelingToolkit.value(k) => v for (k, v) in pmap)
+function adjacencymat(::Type{Matrix{T}}, D, rates) where T
+    nc = size(D, 1)
+    A = zeros(T, nc, nc)
 
-    rates = [substitute(rate, pmap) for rate in reactionrates(rs)]
-    ratematrix(rs, rates)
+    for r in 1:length(rates)
+        s = findfirst(==(-1), @view D[:, r])
+        p = findfirst(==(1), @view D[:, r])
+        A[s, p] = rates[r]
+    end
+    A
 end
 
-function ratematrix(rs::ReactionSystem, parametermap::Vector{<:Pair})
-    pdict = Dict(parametermap)
-    ratematrix(rs, pdict)
+function adjacencymat(rn::ReactionSystem, pmap::Vector{<:Pair}; sparse=false)
+    pdict = Dict(pmap)
+    adjacencymat(rn, pdict; sparse)
 end
 
-function ratematrix(rs::ReactionSystem, parametermap::Tuple)
-    pdict = Dict(parametermap)
-    ratematrix(rs, pdict)
+function adjacencymat(rn::ReactionSystem, pmap::Tuple; sparse=false)
+    pdict = Dict(pmap)
+    adjacencymat(rn, pdict; sparse)
 end
 
-function ratematrix(rs::ReactionSystem, parametermap)
+function adjacencymat(rn::ReactionSystem, pmap)
     error("Parameter map must be a dictionary, tuple, or vector of symbol/value pairs.")
 end
 

test/network_analysis/crn_theory.jl

@@ -1,7 +1,6 @@
 # Tests for properties from chemical reaction network theory: deficiency theorems, complex/detailed balance, etc.
 using Catalyst, StableRNGs, LinearAlgebra, Test
 rng = StableRNG(514)
-
 # Tests that `iscomplexbalanced` works for different rate inputs.
 # Tests that non-valid rate input yields and error
 let
@@ -57,11 +56,14 @@ let
     rates_invalid = reshape(rate_vals, 1, 8)
 
     # Tests that all input types generates the correct rate matrix.
-    Catalyst.ratematrix(rn, rate_vals) == rate_mat
-    Catalyst.ratematrix(rn, rates_vec) == rate_mat
-    Catalyst.ratematrix(rn, rates_tup) == rate_mat
-    Catalyst.ratematrix(rn, rates_dict) == rate_mat
+    @test Catalyst.adjacencymat(rn, rates_vec) == rate_mat
+    @test Catalyst.adjacencymat(rn, rates_tup) == rate_mat
+    @test Catalyst.adjacencymat(rn, rates_dict) == rate_mat
+    @test_throws Exception Catalyst.adjacencymat(rn, rate_vals)
     @test_throws Exception Catalyst.iscomplexbalanced(rn, rates_invalid)
+
+    # Test sparse matrix
+    @test Catalyst.adjacencymat(rn, rates_vec; sparse = true) == rate_mat
 end
 
 ### CONCENTRATION ROBUSTNESS TESTS