Revamp graph structures

README.md

@@ -17,11 +17,13 @@ pkg> add https://github.com/gdalle/SparseMatrixColorings.jl
 
 ## Background
 
-The algorithms implemented in this package are all taken from the following survey:
+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)
 
-Some parts of the survey (like definitions and theorems) are also copied verbatim or referred to by their number in the documentation.
+> [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.
 
 ## Alternatives
 

docs/Project.toml

@@ -1,3 +1,4 @@
 [deps]
 Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4"
+DocumenterInterLinks = "d12716ef-a0f6-4df4-a9f1-a5a34e75c656"
 SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35"

docs/make.jl

@@ -1,6 +1,9 @@
 using Documenter
+using DocumenterInterLinks
 using SparseMatrixColorings
 
+links = InterLinks("ADTypes" => "https://sciml.github.io/ADTypes.jl/stable/")
+
 cp(joinpath(@__DIR__, "..", "README.md"), joinpath(@__DIR__, "src", "index.md"); force=true)
 
 makedocs(;
@@ -9,6 +12,7 @@ makedocs(;
     sitename="SparseMatrixColorings.jl",
     format=Documenter.HTML(),
     pages=["Home" => "index.md", "API reference" => "api.md"],
+    plugins=[links],
 )
 
 deploydocs(; repo="github.com/gdalle/SparseMatrixColorings.jl", devbranch="main")

src/SparseMatrixColorings.jl

@@ -3,27 +3,24 @@
 
 Coloring algorithms for sparse Jacobian and Hessian matrices.
 
-The algorithms implemented in this package are all taken from the following survey:
+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)
 
-Some parts of the survey (like definitions and theorems) are also copied verbatim or referred to by their number in the documentation.
+> [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.
 """
 module SparseMatrixColorings
 
 using ADTypes: ADTypes, AbstractColoringAlgorithm
-using LinearAlgebra: Transpose, parent, transpose
-using Random: AbstractRNG
-using SparseArrays: SparseMatrixCSC, nzrange, rowvals
-
-include("utils.jl")
-
-include("bipartite_graph.jl")
-include("adjacency_graph.jl")
-
-include("distance2_coloring.jl")
-include("star_coloring.jl")
+using LinearAlgebra: Diagonal, Transpose, checksquare, parent, transpose
+using Random: AbstractRNG, default_rng, randperm
+using SparseArrays: SparseArrays, SparseMatrixCSC, nzrange, rowvals, spzeros
 
+include("graph.jl")
+include("order.jl")
+include("coloring.jl")
 include("adtypes.jl")
 include("check.jl")
 

src/adtypes.jl

@@ -5,25 +5,32 @@ Matrix coloring algorithm for sparse Jacobians and Hessians.
 
 Compatible with the [ADTypes.jl coloring framework](https://sciml.github.io/ADTypes.jl/stable/#Coloring-algorithm).
 
+# Constructor
+
+    GreedyColoringAlgorithm(order::AbstractOrder=NaturalOrder())
+
 # Implements
 
-- `ADTypes.column_coloring` with a partial distance-2 coloring of the bipartite graph
-- `ADTypes.row_coloring` with a partial distance-2 coloring of the bipartite graph
-- `ADTypes.symmetric_coloring` with a star coloring of the adjacency graph
+- [`ADTypes.column_coloring`](@extref ADTypes) and [`ADTypes.row_coloring`](@extref ADTypes) with a partial distance-2 coloring of the bipartite graph
+- [`ADTypes.symmetric_coloring`](@extref ADTypes) with a star coloring of the adjacency graph
 """
-struct GreedyColoringAlgorithm <: ADTypes.AbstractColoringAlgorithm end
+struct GreedyColoringAlgorithm{O<:AbstractOrder} <: ADTypes.AbstractColoringAlgorithm
+    order::O
+end
+
+GreedyColoringAlgorithm() = GreedyColoringAlgorithm(NaturalOrder())
 
-function ADTypes.column_coloring(A::AbstractMatrix, ::GreedyColoringAlgorithm)
-    g = BipartiteGraph(A)
-    return distance2_column_coloring(g)
+function ADTypes.column_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
+    bg = BipartiteGraph(A)
+    return partial_distance2_coloring(bg, Val(2), algo.order)
 end
 
-function ADTypes.row_coloring(A::AbstractMatrix, ::GreedyColoringAlgorithm)
-    g = BipartiteGraph(A)
-    return distance2_row_coloring(g)
+function ADTypes.row_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
+    bg = BipartiteGraph(A)
+    return partial_distance2_coloring(bg, Val(1), algo.order)
 end
 
-function ADTypes.symmetric_coloring(A::AbstractMatrix, ::GreedyColoringAlgorithm)
-    g = AdjacencyGraph(A)
-    return star_coloring(g)
+function ADTypes.symmetric_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm)
+    ag = AdjacencyGraph(A)
+    return star_coloring(ag, algo.order)
 end

src/check.jl

@@ -1,14 +1,18 @@
 """
-    check_structurally_orthogonal_columns(A, colors; verbose=false)
+    check_structurally_orthogonal_columns(
+        A::AbstractMatrix, colors::AbstractVector{<:Integer}
+        verbose=false
+    )
 
 Return `true` if coloring the columns of the matrix `A` with the vector `colors` results in a partition that is structurally orthogonal, and `false` otherwise.
 
-Def 3.2: A partition of the columns of a matrix `A` is _structurally orthogonal_ if, for every nonzero element `A[i, j]`, the group containing column `A[:, j]` has no other column with a nonzero in row `i`.
+A partition of the columns of a matrix `A` is _structurally orthogonal_ if, for every nonzero element `A[i, j]`, the group containing column `A[:, j]` has no other column with a nonzero in row `i`.
 
-Thm 3.5: The function [`distance2_column_coloring`](@ref) applied to the [`BipartiteGraph`](@ref) of `A` should return a suitable coloring.
+!!! warning
+    This function is not coded with efficiency in mind, it is designed for small-scale tests.
 """
 function check_structurally_orthogonal_columns(
-    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose=false
+    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=false
 )
     for c in unique(colors)
         js = filter(j -> colors[j] == c, axes(A, 2))
@@ -23,16 +27,20 @@ function check_structurally_orthogonal_columns(
 end
 
 """
-    check_structurally_orthogonal_rows(A, colors; verbose=false)
+    check_structurally_orthogonal_rows(
+        A::AbstractMatrix, colors::AbstractVector{<:Integer};
+        verbose=false
+    )
 
 Return `true` if coloring the rows of the matrix `A` with the vector `colors` results in a partition that is structurally orthogonal, and `false` otherwise.
 
-Def 3.2: A partition of the rows of a matrix `A` is _structurally orthogonal_ if, for every nonzero element `A[i, j]`, the group containing row `A[i, :]` has no other row with a nonzero in column `j`.
+A partition of the rows of a matrix `A` is _structurally orthogonal_ if, for every nonzero element `A[i, j]`, the group containing row `A[i, :]` has no other row with a nonzero in column `j`.
 
-Thm 3.5: The function [`distance2_row_coloring`](@ref) applied to the [`BipartiteGraph`](@ref) of `A` should return a suitable coloring.
+!!! warning
+    This function is not coded with efficiency in mind, it is designed for small-scale tests.
 """
 function check_structurally_orthogonal_rows(
-    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose=false
+    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=false
 )
     for c in unique(colors)
         is = filter(i -> colors[i] == c, axes(A, 1))
@@ -47,17 +55,23 @@ function check_structurally_orthogonal_rows(
 end
 
 """
-    check_symmetrically_orthogonal(A, colors; verbose=false)
+    check_symmetrically_orthogonal(
+        A::AbstractMatrix, colors::AbstractVector{<:Integer};
+        verbose=false
+    )
 
 Return `true` if coloring the columns of the symmetric matrix `A` with the vector `colors` results in a partition that is symmetrically orthogonal, and `false` otherwise.
 
-Def 4.2: A partition of the columns of a symmetrix matrix `A` is _symmetrically orthogonal_ if, for every nonzero element `A[i, j]`, either
+A partition of the columns of a symmetrix matrix `A` is _symmetrically orthogonal_ if, for every nonzero element `A[i, j]`, either
 
 1. the group containing the column `A[:, j]` has no other column with a nonzero in row `i`
 2. the group containing the column `A[:, i]` has no other column with a nonzero in row `j`
+
+!!! warning
+    This function is not coded with efficiency in mind, it is designed for small-scale tests.
 """
 function check_symmetrically_orthogonal(
-    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose=false
+    A::AbstractMatrix, colors::AbstractVector{<:Integer}; verbose::Bool=false
 )
     for i in axes(A, 2), j in axes(A, 2)
         if !iszero(A[i, j])

src/coloring.jl

@@ -0,0 +1,85 @@
+"""
+    partial_distance2_coloring(bg::BipartiteGraph, ::Val{side}, order::AbstractOrder)
+
+Compute a distance-2 coloring of the given `side` (`1` or `2`) in the bipartite graph `bg` and return a vector of integer colors.
+
+A _distance-2 coloring_ is such that two vertices have different colors if they are at distance at most 2.
+
+The vertices are colored in a greedy fashion, following the `order` supplied.
+
+# See also
+
+- [`BipartiteGraph`](@ref)
+- [`AbstractOrder`](@ref)
+"""
+function partial_distance2_coloring(
+    bg::BipartiteGraph, ::Val{side}, order::AbstractOrder
+) where {side}
+    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 w in neighbors(bg, Val(side), v)
+            for x in neighbors(bg, Val(other_side), w)
+                if !iszero(colors[x])
+                    forbidden_colors[colors[x]] = v
+                end
+            end
+        end
+        for c in eachindex(forbidden_colors)
+            if forbidden_colors[c] != v
+                colors[v] = c
+                break
+            end
+        end
+    end
+    return colors
+end
+
+"""
+    star_coloring(ag::AdjacencyGraph, order::AbstractOrder)
+
+Compute a star coloring of all vertices in the adjacency graph `ag` 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)
+- [`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)
+            if !iszero(colors[w])  # w is colored
+                forbidden_colors[colors[w]] = v
+            end
+            for x in neighbors(ag, w)
+                if !iszero(colors[x]) && iszero(colors[w])  # w is not colored
+                    forbidden_colors[colors[x]] = v
+                else
+                    for y in neighbors(ag, x)
+                        if !iszero(colors[y]) && y != w
+                            if colors[y] == colors[w]
+                                forbidden_colors[colors[x]] = v
+                                break
+                            end
+                        end
+                    end
+                end
+            end
+        end
+        for c in eachindex(forbidden_colors)
+            if forbidden_colors[c] != v
+                colors[v] = c
+                break
+            end
+        end
+    end
+    return colors
+end