Compare colorings against values from papers (#3)

* Test against reference papers

* Rng

* Using

* Rename

* CSV

* Add performance tests

* Correct column intersection graph

* Missing using

Author: gdalle

Project.toml

@@ -5,12 +5,14 @@ version = "0.1.0"
 
 [deps]
 ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b"
+DocStringExtensions = "ffbed154-4ef7-542d-bbb7-c09d3a79fcae"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
 ADTypes = "1.2.1"
+DocStringExtensions = "0.9"
 LinearAlgebra = "1"
 Random = "1"
 SparseArrays = "1"

README.md

@@ -21,7 +21,7 @@ The algorithms implemented in this package are mainly taken from the following a
 
 > [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005)
 
-> [ColPack: Software for graph coloring and related problems in scientific computing](https://dl.acm.org/doi/10.1145/2513109.2513110), Gebremedhin et al. (2013)
+> [_ColPack: Software for graph coloring and related problems in scientific computing_](https://dl.acm.org/doi/10.1145/2513109.2513110), Gebremedhin et al. (2013)
 
 Some parts of the articles (like definitions) are thus copied verbatim in the documentation.
 

src/SparseMatrixColorings.jl

@@ -1,22 +1,24 @@
 """
     SparseMatrixColorings
 
-Coloring algorithms for sparse Jacobian and Hessian matrices.
-
-The algorithms implemented in this package are mainly taken from the following articles:
-
-> [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005)
-
-> [ColPack: Software for graph coloring and related problems in scientific computing](https://dl.acm.org/doi/10.1145/2513109.2513110), Gebremedhin et al. (2013)
-
-Some parts of the articles (like definitions) are thus copied verbatim in the documentation.
+$README
 """
 module SparseMatrixColorings
 
 using ADTypes: ADTypes, AbstractColoringAlgorithm
+using DocStringExtensions
 using LinearAlgebra: Diagonal, Transpose, checksquare, parent, transpose
 using Random: AbstractRNG, default_rng, randperm
-using SparseArrays: SparseArrays, SparseMatrixCSC, nzrange, rowvals, spzeros
+using SparseArrays:
+    SparseArrays,
+    SparseMatrixCSC,
+    dropzeros,
+    dropzeros!,
+    nnz,
+    nzrange,
+    rowvals,
+    sparse,
+    spzeros
 
 include("graph.jl")
 include("order.jl")

src/adtypes.jl

@@ -20,17 +20,21 @@ end
 
 GreedyColoringAlgorithm() = GreedyColoringAlgorithm(NaturalOrder())
 
+function Base.show(io::IO, algo::GreedyColoringAlgorithm)
+    return print(io, "GreedyColoringAlgorithm($(algo.order))")
+end
+
 function ADTypes.column_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
-    bg = BipartiteGraph(A)
+    bg = bipartite_graph(A)
     return partial_distance2_coloring(bg, Val(2), algo.order)
 end
 
 function ADTypes.row_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
-    bg = BipartiteGraph(A)
+    bg = bipartite_graph(A)
     return partial_distance2_coloring(bg, Val(1), algo.order)
 end
 
 function ADTypes.symmetric_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
-    ag = AdjacencyGraph(A)
-    return star_coloring(ag, algo.order)
+    ag = adjacency_graph(A)
+    return star_coloring1(ag, algo.order)
 end

src/coloring.jl

@@ -15,10 +15,24 @@ The vertices are colored in a greedy fashion, following the `order` supplied.
 function partial_distance2_coloring(
     bg::BipartiteGraph, ::Val{side}, order::AbstractOrder
 ) where {side}
+    colors = Vector{Int}(undef, length(bg, Val(side)))
+    forbidden_colors = Vector{Int}(undef, length(bg, Val(side)))
+    vertices_in_order = vertices(bg, Val(side), order)
+    partial_distance2_coloring!(colors, forbidden_colors, bg, Val(side), vertices_in_order)
+    return colors
+end
+
+function partial_distance2_coloring!(
+    colors::Vector{Int},
+    forbidden_colors::Vector{Int},
+    bg::BipartiteGraph,
+    ::Val{side},
+    vertices_in_order::AbstractVector{<:Integer},
+) where {side}
+    colors .= 0
+    forbidden_colors .= 0
     other_side = 3 - side
-    colors = zeros(Int, length(bg, Val(side)))
-    forbidden_colors = zeros(Int, length(bg, Val(side)))
-    for v in vertices(bg, Val(side), order)
+    for v in vertices_in_order
         for w in neighbors(bg, Val(side), v)
             for x in neighbors(bg, Val(other_side), w)
                 if !iszero(colors[x])
@@ -33,37 +47,48 @@ function partial_distance2_coloring(
             end
         end
     end
-    return colors
 end
 
 """
-    star_coloring(ag::AdjacencyGraph, order::AbstractOrder)
+    star_coloring1(g::Graph, order::AbstractOrder)
 
-Compute a star coloring of all vertices in the adjacency graph `ag` and return a vector of integer colors.
+Compute a star coloring of all vertices in the adjacency graph `g` and return a vector of integer colors.
 
 A _star coloring_ is a distance-1 coloring such that every path on 4 vertices uses at least 3 colors.
 
 The vertices are colored in a greedy fashion, following the `order` supplied.
 
 # See also
 
-- [`AdjacencyGraph`](@ref)
+- [`Graph`](@ref)
 - [`AbstractOrder`](@ref)
 """
-function star_coloring(ag::AdjacencyGraph, order::AbstractOrder)
-    n = length(ag)
-    colors = zeros(Int, n)
-    forbidden_colors = zeros(Int, n)
-    for v in vertices(ag, order)
-        for w in neighbors(ag, v)
+function star_coloring1(g::Graph, order::AbstractOrder)
+    colors = Vector{Int}(undef, length(g))
+    forbidden_colors = Vector{Int}(undef, length(g))
+    vertices_in_order = vertices(g, order)
+    star_coloring1!(colors, forbidden_colors, g, vertices_in_order)
+    return colors
+end
+
+function star_coloring1!(
+    colors::Vector{Int},
+    forbidden_colors::Vector{Int},
+    g::Graph,
+    vertices_in_order::AbstractVector{<:Integer},
+)
+    colors .= 0
+    forbidden_colors .= 0
+    for v in vertices_in_order
+        for w in neighbors(g, v)
             if !iszero(colors[w])  # w is colored
                 forbidden_colors[colors[w]] = v
             end
-            for x in neighbors(ag, w)
+            for x in neighbors(g, w)
                 if !iszero(colors[x]) && iszero(colors[w])  # w is not colored
                     forbidden_colors[colors[x]] = v
                 else
-                    for y in neighbors(ag, x)
+                    for y in neighbors(g, x)
                         if !iszero(colors[y]) && y != w
                             if colors[y] == colors[w]
                                 forbidden_colors[colors[x]] = v
@@ -81,5 +106,4 @@ function star_coloring(ag::AdjacencyGraph, order::AbstractOrder)
             end
         end
     end
-    return colors
 end

src/graph.jl

@@ -1,81 +1,129 @@
+## Standard graph
+
+"""
+    Graph{T}
+
+Undirected graph structure stored in Compressed Sparse Column (CSC) format.
+
+# Fields
+
+- `colptr::Vector{T}`: same as for `SparseMatrixCSC`
+- `rowval::Vector{T}`: same as for `SparseMatrixCSC`
+"""
 struct Graph{T<:Integer}
     colptr::Vector{T}
     rowval::Vector{T}
 end
 
 Graph(A::SparseMatrixCSC) = Graph(A.colptr, A.rowval)
+Graph(A::AbstractMatrix) = Graph(sparse(A))
 
 Base.length(g::Graph) = length(g.colptr) - 1
+SparseArrays.nnz(g::Graph) = length(g.rowval)
 
+vertices(g::Graph) = 1:length(g)
 neighbors(g::Graph, v::Integer) = view(g.rowval, g.colptr[v]:(g.colptr[v + 1] - 1))
+degree(g::Graph, v::Integer) = length(g.colptr[v]:(g.colptr[v + 1] - 1))
 
-## Adjacency graph
+maximum_degree(g::Graph) = maximum(Base.Fix1(degree, g), vertices(g))
+minimum_degree(g::Graph) = minimum(Base.Fix1(degree, g), vertices(g))
+
+## Bipartite graph
 
 """
-    AdjacencyGraph
+    BipartiteGraph{T}
 
-Undirected graph representing the nonzeros of a symmetrix matrix (typically a Hessian matrix).
+Undirected bipartite graph structure stored in bidirectional Compressed Sparse Column format (redundancy allows for faster access).
 
-The adjacency graph of a symmetrix matrix `A ∈ ℝ^{n × n}` is `G(A) = (V, E)` where
+A bipartite graph has two "sides", which we number `1` and `2`.
 
-- `V = 1:n` is the set of rows or columns
-- `(i, j) ∈ E` whenever `A[i, j] ≠ 0` and `i ≠ j`
+# Fields
+
+- `g1::Graph{T}`: contains the neighbors for vertices on side `1`
+- `g2::Graph{T}`: contains the neighbors for vertices on side `2`
 """
-struct AdjacencyGraph{T}
-    g::Graph{T}
+struct BipartiteGraph{T<:Integer}
+    g1::Graph{T}
+    g2::Graph{T}
 end
 
-function AdjacencyGraph(H::SparseMatrixCSC)
-    g = Graph(H - Diagonal(H))
-    return AdjacencyGraph(g)
-end
+Base.length(bg::BipartiteGraph, ::Val{1}) = length(bg.g1)
+Base.length(bg::BipartiteGraph, ::Val{2}) = length(bg.g2)
+SparseArrays.nnz(bg::BipartiteGraph) = nnz(bg.g1)
 
-Base.length(ag::AdjacencyGraph) = length(ag.g)
+"""
+    vertices(bg::BipartiteGraph, Val(side))
 
+Return the list of vertices of `bg` from the specified `side` as a range `1:n`.
 """
-    neighbors(ag::AdjacencyGraph, v::Integer)
+vertices(bg::BipartiteGraph, ::Val{side}) where {side} = 1:length(bg, Val(side))
 
-Return the neighbors of `v` in the graph `ag`.
 """
-neighbors(ag::AdjacencyGraph, v::Integer) = neighbors(ag.g, v)
+    neighbors(bg::BipartiteGraph, Val(side), v::Integer)
 
-degree(ag::AdjacencyGraph, v::Integer) = length(neighbors(ag, v))
+Return the neighbors of `v` (a vertex from the specified `side`, `1` or `2`), in the graph `bg`.
+"""
+neighbors(bg::BipartiteGraph, ::Val{1}, v::Integer) = neighbors(bg.g1, v)
+neighbors(bg::BipartiteGraph, ::Val{2}, v::Integer) = neighbors(bg.g2, v)
 
-## Bipartite graph
+degree(bg::BipartiteGraph, ::Val{1}, v::Integer) = degree(bg.g1, v)
+degree(bg::BipartiteGraph, ::Val{2}, v::Integer) = degree(bg.g2, v)
+
+function maximum_degree(bg::BipartiteGraph, ::Val{side}) where {side}
+    return maximum(v -> degree(bg, Val(side), v), vertices(bg, Val(side)))
+end
+
+function minimum_degree(bg::BipartiteGraph, ::Val{side}) where {side}
+    return minimum(v -> degree(bg, Val(side), v), vertices(bg, Val(side)))
+end
+
+## Construct from matrices
+
+"""
+    adjacency_graph(H::AbstractMatrix)
+
+Return a [`Graph`](@ref) representing the nonzeros of a symmetric matrix (typically a Hessian matrix).
+
+The adjacency graph of a symmetrix matric `A ∈ ℝ^{n × n}` is `G(A) = (V, E)` where
+
+- `V = 1:n` is the set of rows or columns `i`/`j`
+- `(i, j) ∈ E` whenever `A[i, j] ≠ 0` and `i ≠ j`
+"""
+adjacency_graph(H::SparseMatrixCSC) = Graph(H - Diagonal(H))
+adjacency_graph(H::AbstractMatrix) = adjacency_graph(sparse(H))
 
 """
-    BipartiteGraph
+    bipartite_graph(J::AbstractMatrix)
 
-Undirected bipartite graph representing the nonzeros of a non-symmetric matrix (typically a Jacobian matrix).
+Return a [`BipartiteGraph`](@ref) representing the nonzeros of a non-symmetric matrix (typically a Jacobian matrix).
 
 The bipartite graph of a matrix `A ∈ ℝ^{m × n}` is `Gb(A) = (V₁, V₂, E)` where
 
 - `V₁ = 1:m` is the set of rows `i`
 - `V₂ = 1:n` is the set of columns `j`
 - `(i, j) ∈ E` whenever `A[i, j] ≠ 0`
 """
-struct BipartiteGraph{T}
-    g1::Graph{T}
-    g2::Graph{T}
-end
-
-function BipartiteGraph(J::SparseMatrixCSC)
-    g1 = Graph(SparseMatrixCSC(transpose(J)))  # rows to columns
+function bipartite_graph(J::SparseMatrixCSC)
+    g1 = Graph(transpose(J))  # rows to columns
     g2 = Graph(J)  # columns to rows
     return BipartiteGraph(g1, g2)
 end
 
-Base.length(bg::BipartiteGraph, ::Val{1}) = length(bg.g1)
-Base.length(bg::BipartiteGraph, ::Val{2}) = length(bg.g2)
+bipartite_graph(J::AbstractMatrix) = bipartite_graph(sparse(J))
 
 """
-    neighbors(bg::BipartiteGraph, Val(side), v::Integer)
+    column_intersection_graph(J::AbstractMatrix)
 
-Return the neighbors of `v`, which is a vertex from the specified `side` (`1` or `2`), in the graph `bg`.
-"""
-neighbors(bg::BipartiteGraph, ::Val{1}, v::Integer) = neighbors(bg.g1, v)
-neighbors(bg::BipartiteGraph, ::Val{2}, v::Integer) = neighbors(bg.g2, v)
+Return a [`Graph`](@ref) representing the column intersections of a non-symmetric matrix (typically a Jacobian matrix).
+
+The column intersection graph of a matrix `A ∈ ℝ^{m × n}` is `Gc(A) = (V, E)` where
 
-function degree(bg::BipartiteGraph, ::Val{side}, v::Integer) where {side}
-    return length(neighbors(bg, Val(side), v))
+- `V = 1:n` is the set of columns `j`
+- `(j1, j2) ∈ E` whenever `A[:, j1] ∩ A[:, j2] ≠ ∅`
+"""
+function column_intersection_graph(J::SparseMatrixCSC)
+    A = transpose(J) * J
+    return adjacency_graph(A - Diagonal(A))
 end
+
+column_intersection_graph(J::AbstractMatrix) = column_intersection_graph(sparse(J))

src/order.jl

@@ -18,12 +18,12 @@ Order vertices as they come in the graph.
 """
 struct NaturalOrder <: AbstractOrder end
 
-function vertices(ag::AdjacencyGraph, ::NaturalOrder)
-    return 1:length(ag)
+function vertices(g::Graph, ::NaturalOrder)
+    return vertices(g)
 end
 
 function vertices(bg::BipartiteGraph, ::Val{side}, ::NaturalOrder) where {side}
-    return 1:length(bg, Val(side))
+    return vertices(bg, Val(side))
 end
 
 """
@@ -37,8 +37,8 @@ end
 
 RandomOrder() = RandomOrder(default_rng())
 
-function vertices(ag::AdjacencyGraph, order::RandomOrder)
-    return randperm(order.rng, length(ag))
+function vertices(g::Graph, order::RandomOrder)
+    return randperm(order.rng, length(g))
 end
 
 function vertices(bg::BipartiteGraph, ::Val{side}, order::RandomOrder) where {side}
@@ -52,12 +52,12 @@ Order vertices by decreasing degree.
 """
 struct LargestFirst <: AbstractOrder end
 
-function vertices(ag::AdjacencyGraph, ::LargestFirst)
-    criterion(v) = degree(ag, v)
-    return sort(1:length(ag); by=criterion, rev=true)
+function vertices(g::Graph, ::LargestFirst)
+    criterion(v) = degree(g, v)
+    return sort(vertices(g); by=criterion, rev=true)
 end
 
 function vertices(bg::BipartiteGraph, ::Val{side}, ::LargestFirst) where {side}
     criterion(v) = degree(bg, Val(side), v)
-    return sort(1:length(bg, Val(side)); by=criterion, rev=true)
+    return sort(vertices(bg, Val(side)); by=criterion, rev=true)
 end