From 30cf870d621bf3cd595927c60dcb655531ffba62 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 10 Jan 2025 16:49:42 -0500 Subject: [PATCH 1/3] simplify iscomplexbalanced, rename adjacencymat, adjust implementation to be more similar --- src/network_analysis.jl | 125 +++++++++++++--------------- test/network_analysis/crn_theory.jl | 12 +-- 2 files changed, 66 insertions(+), 71 deletions(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 18bbcdb6f9..86869de253 100644 --- a/src/network_analysis.jl +++ b/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 diff --git a/test/network_analysis/crn_theory.jl b/test/network_analysis/crn_theory.jl index deb53b21b5..4b52171a5f 100644 --- a/test/network_analysis/crn_theory.jl +++ b/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 From 12e7b2b4f7bd040f4fb407817427f67a1cd7bb88 Mon Sep 17 00:00:00 2001 From: vyudu Date: Fri, 17 Jan 2025 14:00:59 -0500 Subject: [PATCH 2/3] add validation --- src/network_analysis.jl | 19 ++++++++++++++++++- 1 file changed, 18 insertions(+), 1 deletion(-) diff --git a/src/network_analysis.jl b/src/network_analysis.jl index 86869de253..8cc558d3e9 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -215,6 +215,13 @@ end **Warning**: Unlike other Catalyst functions, the `fluxmat` function will return a `Matrix{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs, and to ensure that multiplying the sparse form of the matrix will work. """ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) + deps = Set() + for (i, rx) in enumerate(reactions(rn)) + empty!(deps) + get_variables!(deps, rx.rate, species(rn)) + (!isempty(deps)) && (error("Reaction $rx's rate constant depends on species $(join(deps, ", ")). `adjacencymat` cannot support rate constants of this form.")) + end + rates = if isempty(pmap) reactionrates(rn) else @@ -290,6 +297,10 @@ function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric rxs = reactions(rn) sm = speciesmap(rn) + for rx in rxs + !ismassaction(rx, rn) && error("The supplied ReactionSystem has non-mass action reaction $rx. The `massactionvector` can only be constructed for mass action networks.") + end + specs = if isempty(scmap) species(rn) else @@ -1001,7 +1012,6 @@ function iscomplexbalanced(rs::ReactionSystem, parametermap::Dict; sparse = fals 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 - L = laplacianmat(rs, parametermap; sparse) D = incidencemat(rs; sparse) S = netstoichmat(rs; sparse) @@ -1053,6 +1063,13 @@ end - 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 adjacencymat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false) + deps = Set() + for (i, rx) in enumerate(reactions(rn)) + empty!(deps) + get_variables!(deps, rx.rate, species(rn)) + (!isempty(deps)) && (error("Reaction $rx's rate constant depends on species $(join(deps, ", ")). `adjacencymat` cannot support rate constants of this form.")) + end + rates = if isempty(pmap) reactionrates(rn) else From 49a56fae839275476a049c960470d55f04123be2 Mon Sep 17 00:00:00 2001 From: Vincent Du <54586336+vyudu@users.noreply.github.com> Date: Fri, 17 Jan 2025 14:05:16 -0500 Subject: [PATCH 3/3] Update network_analysis.jl --- 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 8cc558d3e9..2800f36ddc 100644 --- a/src/network_analysis.jl +++ b/src/network_analysis.jl @@ -219,7 +219,7 @@ function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false) for (i, rx) in enumerate(reactions(rn)) empty!(deps) get_variables!(deps, rx.rate, species(rn)) - (!isempty(deps)) && (error("Reaction $rx's rate constant depends on species $(join(deps, ", ")). `adjacencymat` cannot support rate constants of this form.")) + (!isempty(deps)) && (error("Reaction $rx's rate constant depends on species $(join(deps, ", ")). `fluxmat` cannot support rate constants of this form.")) end rates = if isempty(pmap)