up

Author: vyudu

src/networkapi.jl

@@ -1656,159 +1656,3 @@ end
 # Checks if a unit consist of exponents with base 1 (and is this unitless).
 unitless_exp(u) = istree(u) && (operation(u) == ^) && (arguments(u)[1] == 1)
 
-"""
-    iscomplexbalanced(rs::ReactionSystem, parametermap)
-
-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 ofvariable-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)
-    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)
-
-    # Compute ρ using the matrix-tree theorem
-    g = incidencematgraph(rs)
-    R = ratematrix(rs, rates)
-    ρ = matrixtree(g, R)
-
-    # Determine if 1) ρ is positive and 2) D^T Ln ρ lies in the image of S^T
-    if all(>(0), ρ)
-        img = D'*log.(ρ)
-        if rank(S') == rank(hcat(S', img)) return true else return false end 
-    else
-        return false
-    end
-end
-
-function iscomplexbalanced(rs::ReactionSystem, parametermap::Vector{Pair{Symbol, Float64}}) 
-    pdict = Dict(parametermap)
-    iscomplexbalanced(rs, pdict)
-end
-
-function iscomplexbalanced(rs::ReactionSystem, parametermap::Tuple{Pair{Symbol, Float64}}) 
-    pdict = Dict(parametermap)
-    iscomplexbalanced(rs, pdict)
-end
-
-iscomplexbalanced(rs::ReactionSystem, parametermap) = error("Parameter map must be a dictionary, tuple, or vector of symbol/value pairs.")
-
-
-"""
-    ratematrix(rs::ReactionSystem, parametermap)
-
-    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. 
-"""
-
-function ratematrix(rs::ReactionSystem, rates::Vector{Float64}) 
-    complexes, D = reactioncomplexes(rs)
-    n = length(complexes)
-    rxns = reactions(rs)
-    ratematrix = zeros(n, n)
-
-    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]
-    end
-    ratematrix
-end
-
-function ratematrix(rs::ReactionSystem, parametermap::Dict{Symbol, Float64}) 
-    if length(parametermap) != numparams(rs) 
-        error("The number of reaction rates must be equal to the number of parameters")
-    end
-    pmap = symmap_to_varmap(rs, parametermap)
-    rates = [substitute(rate, pmap) for rate in reactionrates(rs)]
-    ratematrix(rs, rates)
-end
-
-function ratematrix(rs::ReactionSystem, parametermap::Vector{Pair{Symbol, Float64}}) 
-    pdict = Dict(parametermap)
-    ratematrix(rs, pdict)
-end
-
-function ratematrix(rs::ReactionSystem, parametermap::Tuple{Pair{Symbol, Float64}}) 
-    pdict = Dict(parametermap)
-    ratematrix(rs, pdict)
-end
-
-ratematrix(rs::ReactionSystem, parametermap) = error("Parameter map must be a dictionary, tuple, or vector of symbol/value pairs.")
-
-### BELOW: Helper functions for iscomplexbalanced
-
-function matrixtree(g::SimpleDiGraph, distmx::Matrix)
-    n = nv(g)
-    if size(distmx) != (n, n)
-        error("Size of distance matrix is incorrect")
-    end
-
-    π = zeros(n) 
-
-    if !Graphs.is_connected(g)
-        ccs = Graphs.connected_components(g)
-        for cc in ccs
-            sg, vmap = Graphs.induced_subgraph(g, cc)
-            distmx_s = distmx[cc, cc]
-            π_j = matrixtree(sg, distmx_s)
-            π[cc] = π_j
-        end
-        return π
-    end
-
-    # generate all spanning trees
-    ug = SimpleGraph(SimpleDiGraph(g))
-    trees = collect(Combinatorics.combinations(collect(edges(ug)), n-1))
-    trees = SimpleGraph.(trees)
-    trees = filter!(t->isempty(Graphs.cycle_basis(t)), trees)
-
-    # constructed rooted trees for every vertex, compute sum
-    for v in 1:n
-        rootedTrees = [reverse(Graphs.bfs_tree(t, v, dir=:in)) for t in trees]
-        π[v] = sum([treeweight(t, g, distmx) for t in rootedTrees])
-    end
-    
-    # sum the contributions
-    return π 
-end
-
-function treeweight(t::SimpleDiGraph, g::SimpleDiGraph, distmx::Matrix) 
-    prod = 1
-    for e in edges(t)
-        s = Graphs.src(e); t = Graphs.dst(e)
-        prod *= distmx[s, t] 
-    end
-    prod 
-end