Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
1 change: 1 addition & 0 deletions src/SparseMatrixColorings.jl
Original file line number Diff line number Diff line change
Expand Up @@ -47,6 +47,7 @@ include("graph.jl")
include("forest.jl")
include("order.jl")
include("coloring.jl")
include("postprocessing.jl")
include("result.jl")
include("matrices.jl")
include("interface.jl")
Expand Down
162 changes: 10 additions & 152 deletions src/coloring.jl
Original file line number Diff line number Diff line change
Expand Up @@ -81,7 +81,7 @@ end
"""
star_coloring(
g::AdjacencyGraph, vertices_in_order::AbstractVector, postprocessing::Bool;
forced_colors::Union{AbstractVector,Nothing}=nothing
postprocessing_minimizes::Symbol=:all_colors, forced_colors::Union{AbstractVector,Nothing}=nothing
)

Compute a star coloring of all vertices in the adjacency graph `g` and return a tuple `(color, star_set)`, where
Expand Down Expand Up @@ -110,6 +110,7 @@ function star_coloring(
g::AdjacencyGraph{T},
vertices_in_order::AbstractVector{<:Integer},
postprocessing::Bool;
postprocessing_minimizes::Symbol=:all_colors,
forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing,
) where {T<:Integer}
# Initialize data structures
Expand Down Expand Up @@ -168,7 +169,7 @@ function star_coloring(
if postprocessing
# Reuse the vector forbidden_colors to compute offsets during post-processing
offsets = forbidden_colors
postprocess!(color, star_set, g, offsets)
postprocess!(color, star_set, g, offsets, postprocessing_minimizes)
end
return color, star_set
end
Expand Down Expand Up @@ -250,7 +251,8 @@ struct StarSet{T}
end

"""
acyclic_coloring(g::AdjacencyGraph, vertices_in_order::AbstractVector, postprocessing::Bool)
acyclic_coloring(g::AdjacencyGraph, vertices_in_order::AbstractVector, postprocessing::Bool;
postprocessing_minimizes::Symbol=:all_colors)

Compute an acyclic coloring of all vertices in the adjacency graph `g` and return a tuple `(color, tree_set)`, where

Expand All @@ -273,7 +275,10 @@ If `postprocessing=true`, some colors might be replaced with `0` (the "neutral"
> [_New Acyclic and Star Coloring Algorithms with Application to Computing Hessians_](https://epubs.siam.org/doi/abs/10.1137/050639879), Gebremedhin et al. (2007), Algorithm 3.1
"""
function acyclic_coloring(
g::AdjacencyGraph{T}, vertices_in_order::AbstractVector{<:Integer}, postprocessing::Bool
g::AdjacencyGraph{T},
vertices_in_order::AbstractVector{<:Integer},
postprocessing::Bool;
postprocessing_minimizes::Symbol=:all_colors,
) where {T<:Integer}
# Initialize data structures
nv = nb_vertices(g)
Expand Down Expand Up @@ -345,7 +350,7 @@ function acyclic_coloring(
if postprocessing
# Reuse the vector forbidden_colors to compute offsets during post-processing
offsets = forbidden_colors
postprocess!(color, tree_set, g, offsets)
postprocess!(color, tree_set, g, offsets, postprocessing_minimizes)
end
return color, tree_set
end
Expand Down Expand Up @@ -643,150 +648,3 @@ function TreeSet(

return TreeSet(reverse_bfs_orders, is_star, tree_edge_indices, nt)
end

## Postprocessing, mirrors decompression code

function postprocess!(
color::AbstractVector{<:Integer},
star_or_tree_set::Union{StarSet,TreeSet},
g::AdjacencyGraph,
offsets::AbstractVector{<:Integer},
)
S = pattern(g)
edge_to_index = edge_indices(g)
# flag which colors are actually used during decompression
nb_colors = maximum(color)
color_used = zeros(Bool, nb_colors)

# nonzero diagonal coefficients force the use of their respective color (there can be no neutral colors if the diagonal is fully nonzero)
if has_diagonal(g)
for i in axes(S, 1)
if !iszero(S[i, i])
color_used[color[i]] = true
end
end
end

if star_or_tree_set isa StarSet
# only the colors of the hubs are used
(; star, hub) = star_or_tree_set
nb_trivial_stars = 0

# Iterate through all non-trivial stars
for s in eachindex(hub)
h = hub[s]
if h > 0
color_used[color[h]] = true
else
nb_trivial_stars += 1
end
end

# Process the trivial stars (if any)
if nb_trivial_stars > 0
rvS = rowvals(S)
for j in axes(S, 2)
for k in nzrange(S, j)
i = rvS[k]
if i > j
index_ij = edge_to_index[k]
s = star[index_ij]
h = hub[s]
if h < 0
h = abs(h)
spoke = h == j ? i : j
if color_used[color[spoke]]
# Switch the hub and the spoke to possibly avoid adding one more used color
hub[s] = spoke
else
# Keep the current hub
color_used[color[h]] = true
end
end
end
end
end
end
else
# only the colors of non-leaf vertices are used
(; reverse_bfs_orders, is_star, tree_edge_indices, nt) = star_or_tree_set
nb_trivial_trees = 0

# Iterate through all non-trivial trees
for k in 1:nt
# Position of the first edge in the tree
first = tree_edge_indices[k]

# Total number of edges in the tree
ne_tree = tree_edge_indices[k + 1] - first

# Check if we have more than one edge in the tree (non-trivial tree)
if ne_tree > 1
# Determine if the tree is a star
if is_star[k]
# It is a non-trivial star and only the color of the hub is needed
(_, hub) = reverse_bfs_orders[first]
color_used[color[hub]] = true
else
# It is not a star and both colors are needed during the decompression
(i, j) = reverse_bfs_orders[first]
color_used[color[i]] = true
color_used[color[j]] = true
end
else
nb_trivial_trees += 1
end
end

# Process the trivial trees (if any)
if nb_trivial_trees > 0
for k in 1:nt
# Position of the first edge in the tree
first = tree_edge_indices[k]

# Total number of edges in the tree
ne_tree = tree_edge_indices[k + 1] - first

# Check if we have exactly one edge in the tree
if ne_tree == 1
(i, j) = reverse_bfs_orders[first]
if color_used[color[i]]
# Make i the root to avoid possibly adding one more used color
# Switch it with the (only) leaf
reverse_bfs_orders[first] = (j, i)
else
# Keep j as the root
color_used[color[j]] = true
end
end
end
end
end

# if at least one of the colors is useless, modify the color assignments of vertices
if any(!, color_used)
num_colors_useless = 0

# determine what are the useless colors and compute the offsets
for ci in 1:nb_colors
if color_used[ci]
offsets[ci] = num_colors_useless
else
num_colors_useless += 1
end
end

# assign the neutral color to every vertex with a useless color and remap the colors
for i in eachindex(color)
ci = color[i]
if !color_used[ci]
# assign the neutral color
color[i] = 0
else
# remap the color to not have any gap
color[i] -= offsets[ci]
end
end
end
return color
end
12 changes: 7 additions & 5 deletions src/graph.jl
Original file line number Diff line number Diff line change
Expand Up @@ -227,6 +227,7 @@ The adjacency graph of a symmetric matrix `A ∈ ℝ^{n × n}` is `G(A) = (V, E)
struct AdjacencyGraph{T<:Integer,has_diagonal}
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

has_diagonal should be renamed into augmented_graph or something to that effect

S::SparsityPatternCSC{T}
edge_to_index::Vector{T}
original_size::Tuple{Int,Int}
end

Base.eltype(::AdjacencyGraph{T}) where {T} = T
Expand All @@ -235,16 +236,17 @@ function AdjacencyGraph(
S::SparsityPatternCSC{T},
edge_to_index::Vector{T}=build_edge_to_index(S);
has_diagonal::Bool=true,
original_size::Tuple{Int,Int}=size(S),
) where {T}
return AdjacencyGraph{T,has_diagonal}(S, edge_to_index)
return AdjacencyGraph{T,has_diagonal}(S, edge_to_index, original_size)
end

function AdjacencyGraph(A::SparseMatrixCSC; has_diagonal::Bool=true)
return AdjacencyGraph(SparsityPatternCSC(A); has_diagonal)
function AdjacencyGraph(A::SparseMatrixCSC; has_diagonal::Bool=true, original_size::Tuple{Int,Int}=size(A))
return AdjacencyGraph(SparsityPatternCSC(A); has_diagonal, original_size)
end

function AdjacencyGraph(A::AbstractMatrix; has_diagonal::Bool=true)
return AdjacencyGraph(SparseMatrixCSC(A); has_diagonal)
function AdjacencyGraph(A::AbstractMatrix; has_diagonal::Bool=true, original_size::Tuple{Int,Int}=size(A))
return AdjacencyGraph(SparseMatrixCSC(A); has_diagonal, original_size)
end

pattern(g::AdjacencyGraph) = g.S
Expand Down
28 changes: 17 additions & 11 deletions src/interface.jl
Original file line number Diff line number Diff line change
Expand Up @@ -69,11 +69,12 @@ It is passed as an argument to the main function [`coloring`](@ref).

# Constructors

GreedyColoringAlgorithm{decompression}(order=NaturalOrder(); postprocessing=false)
GreedyColoringAlgorithm(order=NaturalOrder(); postprocessing=false, decompression=:direct)
GreedyColoringAlgorithm{decompression}(order=NaturalOrder(); postprocessing=false, postprocessing_minimizes=:all_colors)
GreedyColoringAlgorithm(order=NaturalOrder(); postprocessing=false, postprocessing_minimizes=:all_colors, decompression=:direct)

- `order::Union{AbstractOrder,Tuple}`: the order in which the columns or rows are colored, which can impact the number of colors. Can also be a tuple of different orders to try out, from which the best order (the one with the lowest total number of colors) will be used.
- `postprocessing::Bool`: whether or not the coloring will be refined by assigning the neutral color `0` to some vertices.
- `postprocessing::Bool`: whether or not the coloring will be refined by assigning the neutral color `0` to some vertices. This option does not affect row or column colorings.
- `postprocessing_minimizes::Symbol`: either `:all_colors`, `:row_colors` or `:column_colors`. This option only applies to star and acyclic bicolorings.
Copy link
Owner

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

Suggested change
- `postprocessing_minimizes::Symbol`: either `:all_colors`, `:row_colors` or `:column_colors`. This option only applies to star and acyclic bicolorings.
- `postprocessing_minimizes::Symbol`: which number of distinct colors is heuristically minimized by postprocessing, either `:all_colors`, `:row_colors` or `:column_colors`. This option only affects bidirectional colorings.

The user-facing documentation does not mention star/acyclic terminology

- `decompression::Symbol`: either `:direct` or `:substitution`. Usually `:substitution` leads to fewer colors, at the cost of a more expensive coloring (and decompression). When `:substitution` is not applicable, it falls back on `:direct` decompression.

!!! warning
Expand All @@ -98,27 +99,30 @@ struct GreedyColoringAlgorithm{decompression,N,O<:NTuple{N,AbstractOrder}} <:
ADTypes.AbstractColoringAlgorithm
orders::O
postprocessing::Bool
postprocessing_minimizes::Symbol

function GreedyColoringAlgorithm{decompression}(
order_or_orders::Union{AbstractOrder,Tuple}=NaturalOrder();
postprocessing::Bool=false,
postprocessing_minimizes::Symbol=:all_colors,
) where {decompression}
check_valid_algorithm(decompression)
if order_or_orders isa AbstractOrder
orders = (order_or_orders,)
else
orders = order_or_orders
end
return new{decompression,length(orders),typeof(orders)}(orders, postprocessing)
return new{decompression,length(orders),typeof(orders)}(orders, postprocessing, postprocessing_minimizes)
end
end

function GreedyColoringAlgorithm(
order_or_orders::Union{AbstractOrder,Tuple}=NaturalOrder();
postprocessing::Bool=false,
decompression::Symbol=:direct,
postprocessing_minimizes::Symbol=:all_colors,
)
return GreedyColoringAlgorithm{decompression}(order_or_orders; postprocessing)
return GreedyColoringAlgorithm{decompression}(order_or_orders; postprocessing, postprocessing_minimizes)
end

## Coloring
Expand Down Expand Up @@ -279,7 +283,7 @@ function _coloring(
symmetric_pattern::Bool;
forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing,
)
ag = AdjacencyGraph(A; has_diagonal=true)
ag = AdjacencyGraph(A; has_diagonal=true, original_size=size(A))
color_and_star_set_by_order = map(algo.orders) do order
vertices_in_order = vertices(ag, order)
return star_coloring(ag, vertices_in_order, algo.postprocessing; forced_colors)
Expand All @@ -300,7 +304,7 @@ function _coloring(
decompression_eltype::Type{R},
symmetric_pattern::Bool,
) where {R}
ag = AdjacencyGraph(A; has_diagonal=true)
ag = AdjacencyGraph(A; has_diagonal=true, original_size=size(A))
color_and_tree_set_by_order = map(algo.orders) do order
vertices_in_order = vertices(ag, order)
return acyclic_coloring(ag, vertices_in_order, algo.postprocessing)
Expand All @@ -323,11 +327,12 @@ function _coloring(
forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing,
) where {R}
A_and_Aᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern)
ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; has_diagonal=false)
ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; has_diagonal=false, original_size=size(A))
postprocessing_minimizes = algo.postprocessing_minimizes
outputs_by_order = map(algo.orders) do order
vertices_in_order = vertices(ag, order)
_color, _star_set = star_coloring(
ag, vertices_in_order, algo.postprocessing; forced_colors
ag, vertices_in_order, algo.postprocessing; postprocessing_minimizes, forced_colors
)
(_row_color, _column_color, _symmetric_to_row, _symmetric_to_column) = remap_colors(
eltype(ag), _color, maximum(_color), size(A)...
Expand Down Expand Up @@ -370,10 +375,11 @@ function _coloring(
symmetric_pattern::Bool,
) where {R}
A_and_Aᵀ, edge_to_index = bidirectional_pattern(A; symmetric_pattern)
ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; has_diagonal=false)
ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; has_diagonal=false, original_size=size(A))
postprocessing_minimizes = algo.postprocessing_minimizes
outputs_by_order = map(algo.orders) do order
vertices_in_order = vertices(ag, order)
_color, _tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing)
_color, _tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing; postprocessing_minimizes)
(_row_color, _column_color, _symmetric_to_row, _symmetric_to_column) = remap_colors(
eltype(ag), _color, maximum(_color), size(A)...
)
Expand Down
Loading
Loading