cycle wip

Author: vyudu

src/network_analysis.jl

@@ -718,3 +718,43 @@ function conservationlaw_errorcheck(rs, pre_varmap)
     isempty(conservedequations(Catalyst.flatten(rs))) ||
         error("The system has conservation laws but initial conditions were not provided for some species.")
 end
+
+"""
+    cycles(rs::ReactionSystem)
+
+    Returns the matrix of cycles, or right eigenvectors of the stoichiometric matrix. 
+"""
+
+function cycles(rs::ReactionSystem) 
+    # nps = get_networkproperties(rs)
+    nsm = netstoichmat(rs)
+    cycles(nsm)
+end
+
+function cycles(nsm::T; col_order = nothing) where {T <: AbstractMatrix}
+
+    # compute the left nullspace over the integers
+    N = MT.nullspace(nsm; col_order)
+
+    # if all coefficients for a conservation law are negative, make positive
+    for Nrow in eachcol(N)
+        all(r -> r <= 0, Nrow) && (Nrow .*= -1)
+    end
+
+    # check we haven't overflowed
+    iszero(nsm * N) || error("Calculation of the cycle matrix was inaccurate, "
+          * "likely due to numerical overflow. Please use a larger integer "
+          * "type like Int128 or BigInt for the net stoichiometry matrix.")
+
+    T(N)
+end
+
+"""
+    fluxmodes(rs::ReactionSystem)
+
+    See documentation for [`cycles`](@ref). 
+"""
+
+function fluxmodes(rs::ReactionSystem) 
+    cycles(rs)
+end

src/reactionsystem.jl

@@ -80,6 +80,7 @@ Base.@kwdef mutable struct NetworkProperties{I <: Integer, V <: BasicSymbolic{Re
     isempty::Bool = true
     netstoichmat::Union{Matrix{Int}, SparseMatrixCSC{Int, Int}} = Matrix{Int}(undef, 0, 0)
     conservationmat::Matrix{I} = Matrix{I}(undef, 0, 0)
+    cyclemat::Matrix{I} = Matrix{I}(undef, 0, 0)
     col_order::Vector{Int} = Int[]
     rank::Int = 0
     nullity::Int = 0
@@ -1456,4 +1457,4 @@ function validate(rs::ReactionSystem, info::String = "")
 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)
+unitless_exp(u) = istree(u) && (operation(u) == ^) && (arguments(u)[1] == 1)