Merge branch 'master' of https://github.com/SciML/Catalyst.jl into graph-docs

comment out HC

Author: vyudu

Project.toml

@@ -67,7 +67,7 @@ RuntimeGeneratedFunctions = "0.5.12"
 SciMLBase = "2.57.2"
 Setfield = "1"
 # StructuralIdentifiability = "0.5.8"
-SymbolicUtils = "< 3.8"
+SymbolicUtils = "3.8.1"
 Symbolics = "6.22"
 Unitful = "1.12.4"
 julia = "1.10"

README.md

@@ -72,10 +72,10 @@ be found in its corresponding research paper, [Catalyst: Fast and flexible model
 - Model steady states can be [computed through homotopy continuation](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/homotopy_continuation/) using [HomotopyContinuation.jl](https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl) (which can find *all* steady states of systems with multiple ones), by [forward ODE simulations](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/nonlinear_solve/#steady_state_solving_simulation) using [SteadyStateDiffEq.jl](https://github.com/SciML/SteadyStateDiffEq.jl), or by [numerically solving steady-state nonlinear equations](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/nonlinear_solve/#steady_state_solving_nonlinear) using [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl).
 - [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) can be used to [compute bifurcation diagrams](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/bifurcation_diagrams/) of model steady states (including finding periodic orbits).
 - [DynamicalSystems.jl](https://github.com/JuliaDynamics/DynamicalSystems.jl) can be used to compute model [basins of attraction](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/dynamical_systems/#dynamical_systems_basins_of_attraction), [Lyapunov spectrums](https://docs.sciml.ai/Catalyst/stable/steady_state_functionality/dynamical_systems/#dynamical_systems_lyapunov_exponents), and other dynamical system properties.
-- [StructuralIdentifiability.jl](https://github.com/SciML/StructuralIdentifiability.jl) can be used to [perform structural identifiability analysis](https://docs.sciml.ai/Catalyst/stable/inverse_problems/structural_identifiability/).
 - [Optimization.jl](https://github.com/SciML/Optimization.jl), [DiffEqParamEstim.jl](https://github.com/SciML/DiffEqParamEstim.jl), and [PEtab.jl](https://github.com/sebapersson/PEtab.jl) can all be used to [fit model parameters to data](https://sebapersson.github.io/PEtab.jl/stable/Define_in_julia/).
 - [GlobalSensitivity.jl](https://github.com/SciML/GlobalSensitivity.jl) can be used to perform [global sensitivity analysis](https://docs.sciml.ai/Catalyst/stable/inverse_problems/global_sensitivity_analysis/) of model behaviors.
 - [SciMLSensitivity.jl](https://github.com/SciML/SciMLSensitivity.jl) can be used to compute local sensitivities of functions containing forward model simulations.
+<!-- - [StructuralIdentifiability.jl](https://github.com/SciML/StructuralIdentifiability.jl) can be used to [perform structural identifiability analysis](https://docs.sciml.ai/Catalyst/stable/inverse_problems/structural_identifiability/). -->
 
 #### Features of packages built upon Catalyst
 - Catalyst [`ReactionSystem`](@ref)s can be [imported from SBML files](https://docs.sciml.ai/Catalyst/stable/model_creation/model_file_loading_and_export/#Loading-SBML-files-using-SBMLImporter.jl-and-SBMLToolkit.jl) via [SBMLImporter.jl](https://github.com/sebapersson/SBMLImporter.jl) and [SBMLToolkit.jl](https://github.com/SciML/SBMLToolkit.jl), and [from BioNetGen .net files](https://docs.sciml.ai/Catalyst/stable/model_creation/model_file_loading_and_export/#file_loading_rni_net) and various stoichiometric matrix network representations using [ReactionNetworkImporters.jl](https://github.com/SciML/ReactionNetworkImporters.jl).

docs/src/index.md

@@ -41,10 +41,10 @@ etc).
 - Model steady states can be [computed through homotopy continuation](@ref homotopy_continuation) using [HomotopyContinuation.jl](https://github.com/JuliaHomotopyContinuation/HomotopyContinuation.jl) (which can find *all* steady states of systems with multiple ones), by [forward ODE simulations](@ref steady_state_solving_simulation) using [SteadyStateDiffEq.jl)](https://github.com/SciML/SteadyStateDiffEq.jl), or by [numerically solving steady-state nonlinear equations](@ref steady_state_solving_nonlinear) using [NonlinearSolve.jl](https://github.com/SciML/NonlinearSolve.jl).
 [BifurcationKit.jl](https://github.com/bifurcationkit/BifurcationKit.jl) can be used to compute bifurcation diagrams of model steady states (including finding periodic orbits).
 - [DynamicalSystems.jl](https://github.com/JuliaDynamics/DynamicalSystems.jl) can be used to compute model [basins of attraction](@ref dynamical_systems_basins_of_attraction), [Lyapunov spectrums](@ref dynamical_systems_lyapunov_exponents), and other dynamical system properties.
-<!--- [StructuralIdentifiability.jl](https://github.com/SciML/StructuralIdentifiability.jl) can be used to perform structural identifiability analysis.-->
 - [Optimization.jl](https://github.com/SciML/Optimization.jl), [DiffEqParamEstim.jl](https://github.com/SciML/DiffEqParamEstim.jl), and [PEtab.jl](https://github.com/sebapersson/PEtab.jl) can all be used to [fit model parameters to data](https://sebapersson.github.io/PEtab.jl/stable/Define_in_julia/).
 - [GlobalSensitivity.jl](https://github.com/SciML/GlobalSensitivity.jl) can be used to perform [global sensitivity analysis](@ref global_sensitivity_analysis) of model behaviors.
 - [SciMLSensitivity.jl](https://github.com/SciML/SciMLSensitivity.jl) can be used to compute local sensitivities of functions containing forward model simulations.
+<!--- [StructuralIdentifiability.jl](https://github.com/SciML/StructuralIdentifiability.jl) can be used to perform structural identifiability analysis.-->
 
 #### [Features of packages built upon Catalyst](@id doc_index_features_other_packages)
 - Catalyst [`ReactionSystem`](@ref)s can be [imported from SBML files](@ref model_file_import_export_sbml) via [SBMLImporter.jl](https://github.com/sebapersson/SBMLImporter.jl) and [SBMLToolkit.jl](https://github.com/SciML/SBMLToolkit.jl), and [from BioNetGen .net files](@ref model_file_import_export_sbml_rni_net) and various stoichiometric matrix network representations using [ReactionNetworkImporters.jl](https://github.com/SciML/ReactionNetworkImporters.jl).

src/Catalyst.jl

@@ -115,8 +115,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

src/network_analysis.jl

@@ -192,6 +192,143 @@ 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. 
+
+    **Warning**: Unlike other Catalyst functions, the `laplacianmat` 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 laplacianmat(rn::ReactionSystem, pmap::Dict = Dict(); sparse = false)
+    D = incidencemat(rn; sparse)
+    K = fluxmat(rn, pmap; sparse)
+    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)
+
+    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`.
+
+    **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)
+    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)
+    if sparse
+        return fluxmat(SparseMatrixCSC{mtype, Int}, rcmap, rates)
+    else
+        return fluxmat(Matrix{mtype}, rcmap, rates)
+    end
+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) 
+    pdict = Dict(pmap)
+    fluxmat(rn, pdict)
+end
+
+function fluxmat(rn::ReactionSystem, pmap::Tuple) 
+    pdict = Dict(pmap)
+    fluxmat(rn, pdict)
+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.
+
+    **Warning**: Unlike other Catalyst functions, the `massactionvector` function will return a `Vector{Num}` in the symbolic case. This is to allow easier computation of the matrix decomposition of the ODEs.
+"""
+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
+
+    Φ
+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)
 
@@ -787,7 +924,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