Merge pull request #4 from SciML/docs_test_build

Docs test build

Author: vyudu

Project.toml

@@ -21,6 +21,7 @@ Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 Oscar = "f1435218-dba5-11e9-1e4d-f1a5fab5fc13"
 Polyhedra = "67491407-f73d-577b-9b50-8179a7c68029"
 PolynomialRoots = "3a141323-8675-5d76-9d11-e1df1406c778"
+PrecompileTools = "aea7be01-6a6a-4083-8856-8a6e6704d82a"
 ReactionNetworkImporters = "b4db0fb7-de2a-5028-82bf-5021f5cfa881"
 SBMLToolkit = "86080e66-c8ac-44c2-a1a0-9adaadfe4a4e"
 Satisfiability = "160ab843-0bc6-4ba4-9585-b7478b70f443"
@@ -45,6 +46,7 @@ Nemo = "0.47"
 Oscar = "1.1.1"
 Polyhedra = "0.7.8"
 PolynomialRoots = "1.0.0"
+PrecompileTools = "1.2.1"
 ReactionNetworkImporters = "0.15.1"
 SBMLToolkit = "0.1.29"
 SafeTestsets = "0.1.0"

ext/CNACDDLib.jl

@@ -0,0 +1,7 @@
+module CNACDDLibExt
+
+using ModelingToolkit, Setfield
+#import BifurcationKit
+#
+#### Observable Plotting Handling ###
+end

src/CatalystNetworkAnalysis.jl

@@ -1,28 +1,28 @@
 module CatalystNetworkAnalysis
 
+using PrecompileTools: @setup_workload, @compile_workload
+
 using Catalyst
 using Satisfiability # For siphon detection
 
 # Algebraic functionality
-using Oscar, Nemo
-import Hecke: n_positive_roots
+using Oscar
+
+# Linear programming (for concordance + deficiency)
+using JuMP, HiGHS 
+const M::Float64 = 1E4
+const ϵ::Float64 = 1E-4
 
-using JuMP, HiGHS # For concordance and deficiency algorithms
-using MixedSubdivisions, DynamicPolynomials # For polytope analysis
 using LinearAlgebra
 using Graphs
-
-using SparseArrays, StaticArrays
 using IterTools, Combinatorics
+using SparseArrays
 
+# Polytope analysis (EFMs)
+using MixedSubdivisions, DynamicPolynomials
 using Polyhedra
 import CDDLib
 
-import ModelingToolkit as MT
-
-const M::Float64 = 1E4
-const ϵ::Float64 = 1E-1
-
 include("persistence.jl")
 export ispersistent, minimalsiphons, iscritical, isconservative, isconsistent
 include("concordance.jl")

src/concentrationrobustness.jl

@@ -6,13 +6,13 @@
     - :GLOBAL_ACR - this species is absolutely concentration-robust for every choice of rate constants
     - :INCONCLUSIVE - the algorithm currently cannot decide whether this network has ACR. One could try calling this function with rate constants provided. 
     - :NO_ACR - the reaction network does not have ACR. 
-    - :INEXACTPARAMS - the algorithm cannot conclude concentration-robustness due to inexact parameters (Floats that are too small)
+    - :INEXACTPARAMS - the algorithm cannot conclude concentration-robustness due to inexact parameters (floats that are too small)
 
     Follows the approach outlined in [Puente et al. 2023](https://arxiv.org/abs/2401.00078).
 """
 function isconcentrationrobust(rn::ReactionSystem; p::VarMapType = Dict()) 
     nps = Catalyst.get_networkproperties(rn)
-    Catalyst.deficiency(rn) == 1 && (isempty(robustspecies(rn)) : return :GLOBAL_ACR)
+    Catalyst.deficiency(rn) == 1 && !isempty(robustspecies_δ1(rn)) && return :GLOBAL_ACR
 
     # Check necessary condition
     possibly_ACR = false
@@ -28,7 +28,6 @@ function isconcentrationrobust(rn::ReactionSystem; p::VarMapType = Dict())
     end
 
     # Convert parameter values to rational values. 
-
     try
         ftype, stype = eltype(p).parameters
         stype <: Rational || (p = Dict{ftype, Rational}(p))
@@ -67,8 +66,10 @@ function isconcentrationrobust(rn::ReactionSystem; p::VarMapType = Dict())
     r_ξ, ξ = polynomial_ring(QQ, :ξ)
 
     !isempty(p) && for i in 1:numspecies(rn)
-        IQ = eliminate(I, vcat(specs[1:i-1], specs[i+1:end]))
+        IQ = Oscar.eliminate(I, vcat(specs[1:i-1], specs[i+1:end]))
         for g in gens(IQ)
+            iszero(g) && continue
+
             # Generate the univariate polynomial corresponding to the generator
             coeff_pos = [exponent_vector(g, j)[i]+1 for j in 1:length(Oscar.terms(g))]
             coeffs = zeros(QQFieldElem, first(coeff_pos))
@@ -83,7 +84,6 @@ function isconcentrationrobust(rn::ReactionSystem; p::VarMapType = Dict())
         end
     end
 
-    
     return :INCONCLUSIVE
 end
 
@@ -118,7 +118,7 @@ end
 
     For a network of deficiency one, return a vector of indices corresponding to species that are concentration robust, i.e. for every positive equilbrium, the concentration of species s will be the same. 
 """
-function robustspecies(rn::ReactionSystem)
+function robustspecies_δ1(rn::ReactionSystem)
     complexes, D = reactioncomplexes(rn)
     nps = Catalyst.get_networkproperties(rn)
 

src/concordance.jl

@@ -49,7 +49,7 @@ end
 """
     isconcordant(rn::ReactionSystem, atol=1e-12)
 
-    Given a reaction network (and an absolute tolerance for the nullspace matrix below which entries should be zero), test whether the reaction network's graph has a property called concordance. A concordant network will not admit multiple equilibria in any stoichiometric compatibility class. The algorithm for this check follows Haixia Ji's PhD thesis, (Ji, 2011).  
+Given a reaction network (and an absolute tolerance for the nullspace matrix below which entries should be zero), test whether the reaction network's graph has a property called concordance. A concordant network will not admit multiple equilibria in any stoichiometric compatibility class. The algorithm for this check follows Haixia Ji's PhD thesis, (Ji, 2011).  
 """
 function isconcordant(rn::ReactionSystem) 
     S = netstoichmat(rn)
@@ -66,8 +66,8 @@ function isconcordant(rn::ReactionSystem)
     #   2. If α[r] == 0 for some reaction r, either σ[s] == 0 for all s in the reactant complex, 
     #   3. or else there are two species s1, s2 in the reactant complex for which sign(σ[s1]) != sign(σ[s2])
 
-    model = add_sign_constraints(S; varName = "σ") 
-    add_subspace_constraints(kerS; model, varName = "α")
+    model = add_sign_constraints(S; var_name = "σ") 
+    add_subspace_constraints(kerS; model, var_name = "α")
     add_concordance_constraints(model, rn)
 
     optimize!(model)

src/cycles.jl

@@ -4,7 +4,7 @@
     Checks if a reaction network is consistent, i.e. admits a positive equilibrium for some choice of rate constants. Equivalent to [`ispositivelydependent`](@ref).
 """
 function isconsistent(rs::ReactionSystem)
-    cyclemat = cycles(rs)
+    cyclemat = Catalyst.cycles(rs)
     has_positive_solution(cyclemat)
 end
 

src/deficiencytheory.jl

@@ -17,17 +17,17 @@ function deficiencyonealgorithm(rn::ReactionSystem)
     s, c = size(Y); r = size(S, 2)
 
     # Initialize sign compatibility model. 
-    model = C.add_sign_constraints(S, var = "μ")
+    model = C.add_sign_constraints(S, var_name = "μ")
     @variable(model, n)
 
     # Iterate over partitions. For each pair of UML partition and confluence vector, 
     # check if there exists a μ that satisfies the resulting constraints. 
     for partition in partitions
-        feasible = solveconstraints(rn, model, g, partition, cutdict)
+        sol, feasible = solveconstraints(rn, model, g, partition, cutdict)
         feasible && return true
 
         if reversible
-            feasible = solveconstraints(rn, model, -g, partition, cutdict)
+            sol, feasible = solveconstraints(rn, model, -g, partition, cutdict)
             feasible && return true
         end
     end
@@ -36,7 +36,8 @@ function deficiencyonealgorithm(rn::ReactionSystem)
 end
 
 function isregular(rn::ReactionSystem) 
-    lcs = linkageclasses(rn); tlcs = terminallinkageclasses(rn)
+    lcs = Catalyst.linkageclasses(rn)
+    tlcs = terminallinkageclasses(rn)
 
     img = incidencematgraph(rn)
     tlc_graphs = [Graphs.induced_subgraph(img, tlc)[1] for tlc in tlcs] 
@@ -57,7 +58,8 @@ function confluencevector(rn::ReactionSystem)
 
     idx = findfirst(i -> @view(g[:, i]) != zeros(nc), 1:nc) 
     g = g[:, idx]
-    tlcs = terminallinkageclasses(rn); lcs = linkageclasses(rn)
+    tlcs = terminallinkageclasses(rn)
+    lcs = Catalyst.linkageclasses(rn)
 
     absorptive = true
     for tlc in tlcs
@@ -127,7 +129,7 @@ end
 # created by removing each cut-link connecting two terminal complexes. 
 function cutlinkpartitions(rn::ReactionSystem) 
     tlcs = terminallinkageclasses(rn)
-    lcs = linkageclasses(rn)
+    lcs = Catalyst.linkageclasses(rn)
     inclusions = [findfirst(lc -> issubset(tlc, lc), lcs) for tlc in tlcs]
 
     img = Graphs.SimpleGraph(incidencematgraph(rn))

src/lp_utils.jl

@@ -5,11 +5,11 @@
 #################################
 
 """
-    add_subspace_constraints(S; model = nothing, varName = "")
+    add_subspace_constraints(S; model = nothing, var_name = "")
 
     Given a matrix S, add the constraints that the solution vector μ will lie in the image space of S. If the model does not already exists, initialize one.
 """
-function add_subspace_constraints(S::M; model = nothing, varName::String = "") where M <: AbstractMatrix
+function add_subspace_constraints(S::T; model = nothing, var_name::String = "") where T <: AbstractMatrix
     (s, r) = size(S)
 
     # Initialize model if none provided. 
@@ -20,20 +20,20 @@ function add_subspace_constraints(S::M; model = nothing, varName::String = "") w
         @objective(model, Min, 0)
     end
 
-    coeffs = varName*"_coeffs"
+    coeffs = var_name*"_coeffs"
     model[Symbol(coeffs)] = @variable(model, [i = 1:r], base_name = coeffs) 
-    model[Symbol(varName)] = @variable(model, [i = 1:s], base_name = varName)
+    model[Symbol(var_name)] = @variable(model, [i = 1:s], base_name = var_name)
 
-    @constraint(model, S*model[Symbol(coeffs)] == model[Symbol(var)])
+    @constraint(model, S*model[Symbol(coeffs)] == model[Symbol(var_name)])
     return model
 end
 
 """
-    add_sign_constraints(S; model = nothing, varName = "")
+    add_sign_constraints(S; model = nothing, var_name = "")
 
-    Given a matrix S, add the constraints that the solution vector var will be sign-compatible with the image space of S. If the model does not already exists, initialize one.
+Given a matrix S, add the constraints that the solution vector var will be sign-compatible with the image space of S. If the model does not already exists, initialize one.
 """
-function add_sign_constraints(S::M; model = nothing, varName::String = "") where M <: AbstractMatrix
+function add_sign_constraints(S::T; model = nothing, var_name::String = "") where T <: AbstractMatrix
     (s, r) = size(S)
 
     # Initialize model if none provided. 
@@ -44,13 +44,13 @@ function add_sign_constraints(S::M; model = nothing, varName::String = "") where
         @objective(model, Min, 0)
     end
 
-    coeffs = varName*"_coeffs"
+    coeffs = var_name*"_coeffs"
     model[Symbol(coeffs)] = @variable(model, [i = 1:r], base_name = coeffs) 
-    model[Symbol(varName)] = @variable(model, [i = 1:s], base_name = varName)
+    model[Symbol(var_name)] = @variable(model, [i = 1:s], base_name = var_name)
 
-    ispos = var*"_ispos"
-    isneg = var*"_isneg"
-    iszer = var*"_iszero"
+    ispos = var_name*"_ispos"
+    isneg = var_name*"_isneg"
+    iszer = var_name*"_iszero"
     model[Symbol(ispos)] = @variable(model, [i = 1:s], Bin, base_name = ispos)
     model[Symbol(isneg)] = @variable(model, [i = 1:s], Bin, base_name = isneg)
     model[Symbol(iszer)] = @variable(model, [i = 1:s], Bin, base_name = iszer) 
@@ -60,17 +60,17 @@ function add_sign_constraints(S::M; model = nothing, varName::String = "") where
        sum(model[Symbol(iszer)]) <= s - 1 # Ensure that var is not the zero vector.
 
        # iszero = 1 --> var[i] == 0 <--> (S * coeffs)[i] == 0
-       model[Symbol(var)] + M * (ones(s) - model[Symbol(iszer)]) ≥ zeros(s)
-       model[Symbol(var)] - M * (ones(s) - model[Symbol(iszer)]) ≤ zeros(s)
+       model[Symbol(var_name)] + M * (ones(s) - model[Symbol(iszer)]) ≥ zeros(s)
+       model[Symbol(var_name)] - M * (ones(s) - model[Symbol(iszer)]) ≤ zeros(s)
        (S*model[Symbol(coeffs)]) + M * (ones(s) - model[Symbol(iszer)]) ≥ zeros(s)
        (S*model[Symbol(coeffs)]) - M * (ones(s) - model[Symbol(iszer)]) ≤ zeros(s)
 
        # isnegative = 1 --> var[i] < 0 <--> (S * coeffs)[i] < 0
-       model[Symbol(var)] - M * (ones(s) - model[Symbol(isneg)]) ≤ -ones(s) * ϵ
+       model[Symbol(var_name)] - M * (ones(s) - model[Symbol(isneg)]) ≤ -ones(s) * ϵ
        (S*model[Symbol(coeffs)]) - M * (ones(s) - model[Symbol(isneg)]) ≤ -ones(s) * ϵ
 
        # ispositive = 1 --> var[i] > 0 <--> (S * coeffs)[i] > 0
-       model[Symbol(var)] + M * (ones(s) - model[Symbol(ispos)]) ≥ ones(s) * ϵ
+       model[Symbol(var_name)] + M * (ones(s) - model[Symbol(ispos)]) ≥ ones(s) * ϵ
        (S*model[Symbol(coeffs)]) + M * (ones(s) - model[Symbol(ispos)]) ≥ ones(s) * ϵ
     end)
 
@@ -101,15 +101,15 @@ function has_positive_solution(M::T; nonneg = false) where {T <: AbstractMatrix}
 end
 
 # Suppose we have a cone defined by the intersection of the nullspace of S and the positive orthant. Then is_extreme_ray tells whether a given vector in this cone is an extreme ray.
-function is_extreme_ray(S::M, x::V; atol = 1e-9) where {M <: AbstractMatrix, V <: AbstractVector}
+function is_extreme_ray(S::T, x::V; atol = 1e-9) where {T <: AbstractMatrix, V <: AbstractVector}
     m, n = size(S)
     length(x) != n && error("The length of x is not correct, expected $n and received $(length(x)).")
     !isapprox(S*x, zeros(m); atol) && error("The provided vector $x is not a solution to Sx = 0.")
     idxset = findall(!=(0), x)
     is_extreme_idxset(S, idxset)
 end
 
-function is_extreme_idxset(S::M, idxs::Vector{Int}) where {M <: AbstractMatrix}
+function is_extreme_idxset(S::T, idxs::Vector{Int}) where {T <: AbstractMatrix}
     m, n = size(S)
     cone_mat = [I; S; -S] 
     cone_mat_eq  = Matrix{eltype(S)}(undef, n+2*m - length(idxs), n)