Skip to content

Commit

Permalink
Compare colorings against values from papers (#3)
Browse files Browse the repository at this point in the history
* Test against reference papers

* Rng

* Using

* Rename

* CSV

* Add performance tests

* Correct column intersection graph

* Missing using
  • Loading branch information
gdalle authored May 21, 2024
1 parent 947a213 commit 35c22cd
Show file tree
Hide file tree
Showing 16 changed files with 412 additions and 114 deletions.
2 changes: 2 additions & 0 deletions Project.toml
Original file line number Diff line number Diff line change
Expand Up @@ -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"
Expand Down
2 changes: 1 addition & 1 deletion README.md
Original file line number Diff line number Diff line change
Expand Up @@ -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.

Expand Down
22 changes: 12 additions & 10 deletions src/SparseMatrixColorings.jl
Original file line number Diff line number Diff line change
@@ -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")
Expand Down
12 changes: 8 additions & 4 deletions src/adtypes.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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
56 changes: 40 additions & 16 deletions src/coloring.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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])
Expand All @@ -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
Expand All @@ -81,5 +106,4 @@ function star_coloring(ag::AdjacencyGraph, order::AbstractOrder)
end
end
end
return colors
end
120 changes: 84 additions & 36 deletions src/graph.jl
Original file line number Diff line number Diff line change
@@ -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))
18 changes: 9 additions & 9 deletions src/order.jl
Original file line number Diff line number Diff line change
Expand Up @@ -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

"""
Expand All @@ -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}
Expand All @@ -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
Loading

0 comments on commit 35c22cd

Please sign in to comment.