Use Flint via Nemo for integer matrix calculations

Project.toml

@@ -3,7 +3,6 @@ uuid = "479239e8-5488-4da2-87a7-35f2df7eef83"
 version = "9.0.0"
 
 [deps]
-AbstractAlgebra = "c3fe647b-3220-5bb0-a1ea-a7954cac585d"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqBase = "2b5f629d-d688-5b77-993f-72d75c75574e"
 DiffEqJump = "c894b116-72e5-5b58-be3c-e6d8d4ac2b12"
@@ -12,14 +11,14 @@ Latexify = "23fbe1c1-3f47-55db-b15f-69d7ec21a316"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
 MacroTools = "1914dd2f-81c6-5fcd-8719-6d5c9610ff09"
 ModelingToolkit = "961ee093-0014-501f-94e3-6117800e7a78"
+Nemo = "2edaba10-b0f1-5616-af89-8c11ac63239a"
 Parameters = "d96e819e-fc66-5662-9728-84c9c7592b0a"
 Reexport = "189a3867-3050-52da-a836-e630ba90ab69"
 Requires = "ae029012-a4dd-5104-9daa-d747884805df"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 Symbolics = "0c5d862f-8b57-4792-8d23-62f2024744c7"
 
 [compat]
-AbstractAlgebra = "0.19, 0.21"
 DataStructures = "0.18"
 DiffEqBase = "6.54.0"
 DiffEqJump = "7.0"
@@ -28,6 +27,7 @@ Latexify = "0.14, 0.15"
 LightGraphs = "1.3"
 MacroTools = "0.5.5"
 ModelingToolkit = "6.3"
+Nemo = "0.26"
 Parameters = "0.12"
 Reexport = "0.2, 1.0"
 Requires = "1.0"

src/Catalyst.jl

@@ -4,9 +4,11 @@ $(DocStringExtensions.README)
 module Catalyst
 
 using DocStringExtensions
-using SparseArrays, DiffEqBase, Reexport, ModelingToolkit, DiffEqJump, LightGraphs
+using SparseArrays, DiffEqBase, Reexport, ModelingToolkit, DiffEqJump
+using Latexify, Requires
 
 # ModelingToolkit imports and convenience functions we use
+using Symbolics
 using ModelingToolkit: Symbolic, value, istree, get_states, get_ps, get_iv, get_systems, 
                        get_eqs, get_defaults, toparam
 import ModelingToolkit: get_variables, namespace_expr, namespace_equation, get_variables!, 
@@ -17,12 +19,11 @@ import ModelingToolkit: check_variables, check_parameters, _iszero, _merge, chec
 
 const DEFAULT_IV = (@parameters t)[1]
 @reexport using ModelingToolkit
-import MacroTools
+
 import Base: (==), merge!, merge, hash, size, getindex, setindex, isless, Sort.defalg, length, show
-using Symbolics
-using Latexify, Requires
-import LightGraphs
-import AbstractAlgebra; const AA = AbstractAlgebra
+import MacroTools, LightGraphs
+const LG = LightGraphs
+import Nemo: FlintZZ, matrix, nullspace, rank
 
 # as used in Catlab
 const USE_GV_JLL = Ref(false)

src/networkapi.jl

@@ -493,32 +493,32 @@ incidencegraph   = incidencematgraph(incidencemat)
 function incidencematgraph(incidencemat::Matrix{Int})
    @assert all(∈([-1,0,1]) ,incidencemat)
    n = size(incidencemat,1)  # no. of nodes/complexes
-   graph = LightGraphs.DiGraph(n)
+   graph = LG.DiGraph(n)
    for col in eachcol(incidencemat)
        src = 0; dst = 0;
        for i in eachindex(col)
               (col[i] == -1) && (src = i)
               (col[i] == 1) && (dst = i)
               (src != 0) && (dst != 0) && break
        end
-       add_edge!(graph, src, dst)
+       LG.add_edge!(graph, src, dst)
     end
    return graph
 end
 function incidencematgraph(incidencemat::SparseMatrixCSC{Int,Int})
    @assert all(∈([-1,0,1]) ,incidencemat)
    m,n = size(incidencemat)  
-   graph = LightGraphs.DiGraph(m)
+   graph = LG.DiGraph(m)
    rows = rowvals(incidencemat)
    vals = nonzeros(incidencemat)
    for j = 1:n
       inds=nzrange(incidencemat, j)
       row = rows[inds];
       val = vals[inds];
       if val[1] == -1
-         add_edge!(graph, row[1], row[2])
+         LG.add_edge!(graph, row[1], row[2])
       else
-         add_edge!(graph, row[2], row[1])
+         LG.add_edge!(graph, row[2], row[1])
       end
    end
    return graph
@@ -540,8 +540,8 @@ julia> linkageclasses(incidencegraph)
  [3, 4]
 ```
 """
-function linkageclasses(incidencegraph::SimpleDiGraph{Int64})
-   LightGraphs.connected_components(incidencegraph)
+function linkageclasses(incidencegraph)
+   LG.connected_components(incidencegraph)
 end
 
 
@@ -570,8 +570,8 @@ netstoich_mat    = netstoichmat(sir)
 δ = deficiency(netstoich_mat, incidence_graph, linkage_classes)
 ```
 """
-function deficiency(ns::AbstractMatrix, ig::SimpleDiGraph, lc::AbstractArray)
-   LightGraphs.nv(ig) - length(lc) - AA.rank(AA.matrix(AA.zz,ns))
+function deficiency(ns, ig, lc)
+   LG.nv(ig) - length(lc) - rank(matrix(FlintZZ,ns))
 end
 
 
@@ -584,28 +584,20 @@ end
 Given the net stoichiometry matrix of a reaction system, computes a matrix of
 conservation laws, each represented as a row in the output. 
 """
-function conservationlaws(nsm::AbstractMatrix)::Matrix
-    n_spec,n_reac = size(nsm)
-
+function conservationlaws(nsm::AbstractMatrix)
+
     # We basically have to compute the left null space of the matrix
-    # over the integers; this is best done using its Smith Normal Form.
-    # note, we transpose as this was written when netstoichmat was reac by spec
-    nsm_conv = AA.matrix(AA.ZZ, nsm')
-    S, T, U = AA.snf_with_transform(nsm_conv)
-
-    # Zero columns of S (which occur after nonzero columns in SNF)
-    # correspond to conserved quantities
-    n = findfirst(i -> all(S[:,i] .== 0), 1:n_spec)
-    if n === nothing
-        return zeros(Int, 0, n_spec)
-    end
-
-    ret = Matrix(U[:,n:end]')
-
+    # over the integers. We do this using Nemo's Flint integer (ZZ) interface.
+    N = nullspace(matrix(FlintZZ, nsm'))[2]
+
+    # to save allocations we manually take the adjoint when converting back
+    # to a Julia integer matrix from the Nemo matrix. 
+    ret = [convert(Int,N[i,j]) for j=1:size(N,2), i=1:size(N,1)]  
+
     # If all coefficients for a conservation law are negative
     # we might as well flip them to become positive
-    for i in 1:size(ret,1)
-        all(ret[i,:] .<= 0) && (ret[i,:] .*= -1)
+    for retcol in eachcol(ret)
+        all(r -> r <= 0, retcol) && (retcol .*= -1)
     end
 
     ret