Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Add API functions for the other kinds of matrixes that a CRN ODE system can be factored into #1134

Open
wants to merge 28 commits into
base: master
Choose a base branch
from
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
4 changes: 2 additions & 2 deletions src/Catalyst.jl
Original file line number Diff line number Diff line change
Expand Up @@ -128,8 +128,8 @@ export @reaction_network, @network_component, @reaction, @species
# Network analysis functionality.
include("network_analysis.jl")
export reactioncomplexmap, reactioncomplexes, incidencemat
export complexstoichmat
export complexoutgoingmat, incidencematgraph, linkageclasses, stronglinkageclasses,
export complexstoichmat, laplacianmat, fluxmat, massactionvector, complexoutgoingmat
export incidencematgraph, linkageclasses, stronglinkageclasses,
terminallinkageclasses, deficiency, subnetworks
export linkagedeficiencies, isreversible, isweaklyreversible
export conservationlaws, conservedquantities, conservedequations, conservationlaw_constants
Expand Down
130 changes: 129 additions & 1 deletion src/network_analysis.jl
Original file line number Diff line number Diff line change
Expand Up @@ -192,6 +192,134 @@ function complexstoichmat(::Type{Matrix{Int}}, rn::ReactionSystem, rcs)
Z
end

@doc raw"""
laplacianmat(rn::ReactionSystem, pmap=Dict(); sparse=false)

Return the negative of the graph Laplacian of the reaction network. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the graph Laplacian, and ``Φ`` is the [`massactionvector`](@ref). ``A_k`` is an n-by-n matrix, where n is the number of complexes, where ``A_{ij} = k_{ij}`` if a reaction exists between the two complexes and 0 otherwise.
Returns a symbolic matrix by default, but will return a numerical matrix if parameter values are specified via pmap.
"""
function laplacianmat(rn::ReactionSystem, pmap = Dict(); sparse = false)
D = incidencemat(rn; sparse)
K = fluxmat(rn, pmap; sparse)
D*K
end

@doc raw"""
fluxmat(rn::ReactionSystem, pmap = Dict(); sparse=false)

Return an r×c matrix ``K`` such that, if complex ``j`` is the substrate complex of reaction ``i``, then ``K_{ij} = k``, the rate constant for this reaction. Mostly a helper function for the network Laplacian, [`laplacianmat`](@ref). Has the useful property that ``\frac{dx}{dt} = S*K*Φ(x)``, where S is the [`netstoichmat`](@ref) or net stoichiometry matrix and ``Φ(x)`` is the [`massactionvector`](@ref).
Returns a symbolic matrix by default, but will return a numerical matrix if rate constants are specified as a `Tuple`, `Vector`, or `Dict` of symbol-value pairs via `pmap`.
"""
function fluxmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse=false)
rates = if isempty(pmap)
reactionrates(rn)
else
substitutevals(rn, pmap, parameters(rn), reactionrates(rn))
end

rcmap = reactioncomplexmap(rn)
nc = length(rcmap)
nr = length(rates)
mtype = eltype(rates) <: Symbolics.BasicSymbolic ? Num : eltype(rates)
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Should this be Num and not the unwrapped type?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

I think it has to, doesn't seem like it's actually possible to have zeros in this matrix if the eltype is BasicSymbolic

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Can't it by type Any?

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Is that better than being Num? I assumed more specific was better

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

We don't want to return mixtures of Nums and non-Nums across different methods, so we should not re-wrap internal symbolics.

Copy link
Collaborator Author

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Could I just wrap it in an Any when I return it? It's useful to have it be Num at least temporarily because I use zero(T) and one(T) in various places (and product over an empty Num[] returns 1, which is useful).

Copy link
Collaborator Author

@vyudu vyudu Dec 13, 2024

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Actually this creates some issues for the sparse method (can't really do any arithmetic with a SparseMatrixCSC{Any, T} since zero(Any) is undefined). Would Union{Float64, BasicSymbolic} be okay? If not it might just be better to define two functions, one symbolic and one numeric

mat = if sparse
fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates)
else
fluxmat(Matrix{mtype}, rcmap, rates)
end
mtype == Num ? Matrix{Any}(mat) : mat
end

function fluxmat(::Type{SparseMatrixCSC{T, Int}}, rcmap, rates) where T
Is = Int[]
Js = Int[]
Vs = T[]
for (i, (complex, rxs)) in enumerate(rcmap)
for (rx, dir) in rxs
dir == -1 && begin
push!(Is, rx)
push!(Js, i)
push!(Vs, rates[rx])
end
end
end
Z = sparse(Is, Js, Vs, length(rates), length(rcmap))
end

function fluxmat(::Type{Matrix{T}}, rcmap, rates) where T
nr = length(rates)
nc = length(rcmap)
K = zeros(T, nr, nc)
for (i, (complex, rxs)) in enumerate(rcmap)
for (rx, dir) in rxs dir == -1 && (K[rx, i] = rates[rx])
end
end
K
end

function fluxmat(rn::ReactionSystem, pmap::Vector; sparse = false)
pdict = Dict(pmap)
fluxmat(rn, pdict; sparse)
end

function fluxmat(rn::ReactionSystem, pmap::Tuple; sparse = false)
pdict = Dict(pmap)
fluxmat(rn, pdict; sparse)
end

# Helper to substitute values into a (vector of) symbolic expressions. The syms are the symbols to substitute and the symexprs are the expressions to substitute into.
function substitutevals(rn::ReactionSystem, map::Dict, syms, symexprs)
length(map) != length(syms) && error("Incorrect number of parameter-value pairs were specified.")
map = symmap_to_varmap(rn, map)
map = Dict(ModelingToolkit.value(k) => v for (k, v) in map)
vals = [substitute(expr, map) for expr in symexprs]
end

"""
massactionvector(rn::ReactionSystem, scmap = Dict(); combinatoric_ratelaws = true)

Return the vector whose entries correspond to the "mass action products" of each complex. For example, given the complex A + B, the corresponding entry of the vector would be ``A*B``, and for the complex 2X + Y, the corresponding entry would be ``X^2*Y``. The ODE system of a chemical reaction network can be factorized as ``\frac{dx}{dt} = Y A_k Φ(x)``, where ``Y`` is the [`complexstoichmat`](@ref) and ``A_k`` is the negative of the [`laplacianmat`](@ref). This utility returns ``Φ(x)``.
Returns a symbolic vector by default, but will return a numerical vector if species concentrations are specified as a tuple, vector, or dictionary via scmap.
If the `combinatoric_ratelaws` option is set, will include prefactors for that (see [introduction to Catalyst's rate laws](@ref introduction_to_catalyst_ratelaws). Will default to the default for the system.
"""
function massactionvector(rn::ReactionSystem, scmap::Dict = Dict(); combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn))
r = numreactions(rn)
rxs = reactions(rn)
sm = speciesmap(rn)

specs = if isempty(scmap)
species(rn)
else
substitutevals(rn, scmap, species(rn), species(rn))
end

if !all(r -> ismassaction(r, rn), rxs)
error("The supplied ReactionSystem has reactions that are not ismassaction. The mass action vector is only defined for pure mass action networks.")
end

vtype = eltype(specs) <: Symbolics.BasicSymbolic ? Num : eltype(specs)
Φ = Vector{vtype}()
rcmap = reactioncomplexmap(rn)
for comp in keys(reactioncomplexmap(rn))
subs = map(ce -> getfield(ce, :speciesid), comp)
stoich = map(ce -> getfield(ce, :speciesstoich), comp)
maprod = prod(vtype[specs[s]^α for (s, α) in zip(subs, stoich)])
combinatoric_ratelaws && (maprod /= prod(map(factorial, stoich)))
push!(Φ, maprod)
end

vtype == Num ? Vector{Any}(Φ) : Φ
end

function massactionvector(rn::ReactionSystem, scmap::Tuple; combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn))
sdict = Dict(scmap)
massactionvector(rn, sdict; combinatoric_ratelaws)
end

function massactionvector(rn::ReactionSystem, scmap::Vector; combinatoric_ratelaws = Catalyst.get_combinatoric_ratelaws(rn))
sdict = Dict(scmap)
massactionvector(rn, sdict; combinatoric_ratelaws)
end

@doc raw"""
complexoutgoingmat(network::ReactionSystem; sparse=false)

Expand Down Expand Up @@ -787,7 +915,7 @@ function isdetailedbalanced(rs::ReactionSystem, parametermap::Dict; abstol=0, re
elseif !isreversible(rs)
return false
elseif !all(r -> ismassaction(r, rs), reactions(rs))
error("The supplied ReactionSystem has reactions that are not ismassaction. Testing for being complex balanced is currently only supported for pure mass action networks.")
error("The supplied ReactionSystem has reactions that are not ismassaction. Testing for being detailed balanced is currently only supported for pure mass action networks.")
end

isforestlike(rs) && deficiency(rs) == 0 && return true
Expand Down
Loading
Loading