diff --git a/Project.toml b/Project.toml index 1ebf50d..44a0a2d 100644 --- a/Project.toml +++ b/Project.toml @@ -1,7 +1,7 @@ name = "SparseMatrixColorings" uuid = "0a514795-09f3-496d-8182-132a7b665d35" authors = ["Guillaume Dalle", "Alexis Montoison"] -version = "0.4.12" +version = "0.4.13" [deps] ADTypes = "47edcb42-4c32-4615-8424-f2b9edc5f35b" diff --git a/docs/src/vis.md b/docs/src/vis.md index 408e439..72c76f2 100644 --- a/docs/src/vis.md +++ b/docs/src/vis.md @@ -61,19 +61,17 @@ Finally, a background color can be passed via the `background` keyword argument. We demonstrate this on a bidirectional coloring. ```@example img - S = sparse([ 1 1 1 1 1 1 1 1 1 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 1 0 0 0 0 0 0 0 1 - 1 0 0 0 0 0 0 0 1 1 1 1 1 1 1 1 1 1 ]) problem_bi = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional) -algo_bi = GreedyColoringAlgorithm(RandomOrder(StableRNG(0)); decompression=:direct) +algo_bi = GreedyColoringAlgorithm(RandomOrder(StableRNG(0)); postprocessing=true, decompression=:direct) result_bi = coloring(S, problem_bi, algo_bi) A_img, Br_img, Bc_img = show_colors( diff --git a/ext/SparseMatrixColoringsColorsExt.jl b/ext/SparseMatrixColoringsColorsExt.jl index 99e1eff..8733e04 100644 --- a/ext/SparseMatrixColoringsColorsExt.jl +++ b/ext/SparseMatrixColoringsColorsExt.jl @@ -131,7 +131,9 @@ function show_colors!( if !iszero(A[I]) r, c = Tuple(I) area = matrix_entry_area(I, scale, pad) - A_img[area] .= A_colors[c] + if column_colors(res)[c] > 0 + A_img[area] .= A_colors[c] + end end end for I in CartesianIndices(B) @@ -163,7 +165,9 @@ function show_colors!( if !iszero(A[I]) r, c = Tuple(I) area = matrix_entry_area(I, scale, pad) - A_img[area] .= A_colors[r] + if row_colors(res)[r] > 0 + A_img[area] .= A_colors[r] + end end end for I in CartesianIndices(B) @@ -205,9 +209,13 @@ function show_colors!( area = matrix_entry_area(I, scale, pad) for i in axes(area, 1), j in axes(area, 2) if j > i - A_img[area[i, j]] = A_ccolors[c] + if column_colors(res)[c] > 0 + A_img[area[i, j]] = A_ccolors[c] + end elseif i > j - A_img[area[i, j]] = A_rcolors[r] + if row_colors(res)[r] > 0 + A_img[area[i, j]] = A_rcolors[r] + end end end end diff --git a/src/check.jl b/src/check.jl index 5210a3f..ad7c552 100644 --- a/src/check.jl +++ b/src/check.jl @@ -105,23 +105,10 @@ function symmetrically_orthogonal_columns( for i in axes(A, 2), j in axes(A, 2) iszero(A[i, j]) && continue ci, cj = color[i], color[j] - gi, gj = group[ci], group[cj] - A_gj_rowi = view(A, i, gj) - A_gi_rowj = view(A, j, gi) - nonzeros_gj_rowi = count(!iszero, A_gj_rowi) - nonzeros_gi_rowj = count(!iszero, A_gi_rowj) - if nonzeros_gj_rowi > 1 && nonzeros_gi_rowj > 1 - if verbose - gj_incompatible_columns = gj[findall(!iszero, A_gj_rowi)] - gi_incompatible_columns = gi[findall(!iszero, A_gi_rowj)] - @warn """ - For coefficient (i=$i, j=$j) with column colors (ci=$ci, cj=$cj): - - In color ci=$ci, columns $gi_incompatible_columns all have nonzeros in row j=$j. - - In color cj=$cj, columns $gj_incompatible_columns all have nonzeros in row i=$i. - """ - end - return false - end + check = _bilateral_check( + A; i, j, ci, cj, row_group=group, column_group=group, verbose + ) + !check && return false end return true end @@ -156,6 +143,62 @@ function structurally_biorthogonal( for i in axes(A, 1), j in axes(A, 2) iszero(A[i, j]) && continue ci, cj = row_color[i], column_color[j] + check = _bilateral_check(A; i, j, ci, cj, row_group, column_group, verbose) + !check && return false + end + return true +end + +function _bilateral_check( + A::AbstractMatrix; + i::Integer, + j::Integer, + ci::Integer, + cj::Integer, + row_group::AbstractVector, + column_group::AbstractVector, + verbose::Bool, +) + if ci == 0 && cj == 0 + if verbose + @warn """ + For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj): + - Row color ci=$ci is neutral. + - Column color cj=$cj is neutral. + """ + end + return false + elseif ci == 0 && cj != 0 + gj = column_group[cj] + A_gj_rowi = view(A, i, gj) + nonzeros_gj_rowi = count(!iszero, A_gj_rowi) + if nonzeros_gj_rowi > 1 + if verbose + gj_incompatible_columns = gj[findall(!iszero, A_gj_rowi)] + @warn """ + For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj): + - Row color ci=$ci is neutral. + - In column color cj=$cj, columns $gj_incompatible_columns all have nonzeros in row i=$i. + """ + end + return false + end + elseif ci != 0 && cj == 0 + gi = row_group[ci] + A_gi_columnj = view(A, gi, j) + nonzeros_gi_columnj = count(!iszero, A_gi_columnj) + if nonzeros_gi_columnj > 1 + if verbose + gi_incompatible_rows = gi[findall(!iszero, A_gi_columnj)] + @warn """ + For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj): + - In row color ci=$ci, rows $gi_incompatible_rows all have nonzeros in column j=$j. + - Column color cj=$cj is neutral. + """ + end + return false + end + else gi, gj = row_group[ci], column_group[cj] A_gj_rowi = view(A, i, gj) A_gi_columnj = view(A, gi, j) @@ -166,7 +209,7 @@ function structurally_biorthogonal( gj_incompatible_columns = gj[findall(!iszero, A_gj_rowi)] gi_incompatible_rows = gi[findall(!iszero, A_gi_columnj)] @warn """ - For coefficient (i=$i, j=$j) with row color ci=$ci and column color cj=$cj: + For coefficient (i=$i, j=$j) with colors (ci=$ci, cj=$cj): - In row color ci=$ci, rows $gi_incompatible_rows all have nonzeros in column j=$j. - In column color cj=$cj, columns $gj_incompatible_columns all have nonzeros in row i=$i. """ @@ -199,12 +242,16 @@ function directly_recoverable_columns( return false end group = group_by_color(color) - B = stack(group; dims=2) do g - dropdims(sum(A[:, g]; dims=2); dims=2) + B = if isempty(group) + similar(A, size(A, 1), 0) + else + stack(group; dims=2) do g + dropdims(sum(A[:, g]; dims=2); dims=2) + end end A_unique = Set(unique(A)) B_unique = Set(unique(B)) - if !issubset(A_unique, B_unique) + if !issubset(A_unique, push!(B_unique, zero(eltype(B)))) if verbose @warn "Coefficients $(sort(collect(setdiff(A_unique, B_unique)))) are not directly recoverable." return false diff --git a/src/coloring.jl b/src/coloring.jl index c406eef..1a53202 100644 --- a/src/coloring.jl +++ b/src/coloring.jl @@ -54,7 +54,7 @@ function partial_distance2_coloring!( end """ - star_coloring(g::AdjacencyGraph, order::AbstractOrder) + star_coloring(g::AdjacencyGraph, order::AbstractOrder; postprocessing::Bool) Compute a star coloring of all vertices in the adjacency graph `g` and return a tuple `(color, star_set)`, where @@ -65,6 +65,8 @@ A _star coloring_ is a distance-1 coloring such that every path on 4 vertices us The vertices are colored in a greedy fashion, following the `order` supplied. +If `postprocessing=true`, some colors might be replaced with `0` (the "neutral" color) as long as they are not needed during decompression. + # See also - [`AdjacencyGraph`](@ref) @@ -74,7 +76,7 @@ The vertices are colored in a greedy fashion, following the `order` supplied. > [_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 4.1 """ -function star_coloring(g::AdjacencyGraph, order::AbstractOrder) +function star_coloring(g::AdjacencyGraph, order::AbstractOrder; postprocessing::Bool) # Initialize data structures nv = nb_vertices(g) color = zeros(Int, nv) @@ -82,7 +84,7 @@ function star_coloring(g::AdjacencyGraph, order::AbstractOrder) first_neighbor = fill((0, 0), nv) # at first no neighbors have been encountered treated = zeros(Int, nv) star = Dict{Tuple{Int,Int},Int}() - hub = Int[] + hub = Int[] # one hub for each star, including the trivial ones vertices_in_order = vertices(g, order) for v in vertices_in_order @@ -102,7 +104,7 @@ function star_coloring(g::AdjacencyGraph, order::AbstractOrder) for x in neighbors(g, w) (x == v || iszero(color[x])) && continue wx = _sort(w, x) - if x == hub[star[wx]] # potential Case 2 + if x == hub[star[wx]] # potential Case 2 (which is always false for trivial stars with two vertices, since the associated hub is negative) forbidden_colors[color[x]] = v end end @@ -116,7 +118,12 @@ function star_coloring(g::AdjacencyGraph, order::AbstractOrder) end _update_stars!(star, hub, g, v, color, first_neighbor) end - return color, StarSet(star, hub) + hub .= abs.(hub) # make all hubs actual (positive) hubs before creating StarSet + star_set = StarSet(star, hub) + if postprocessing + postprocess!(color, star_set, g) + end + return color, star_set end """ @@ -131,7 +138,7 @@ $TYPEDFIELDS struct StarSet "a mapping from edges (pair of vertices) to their star index" star::Dict{Tuple{Int,Int},Int} - "a mapping from star indices to their hub (undefined hubs for single-edge stars are the negative value of one of the vertices, picked arbitrarily)" + "a mapping from star indices to their hub (hubs of single-edge stars are picked arbitrarily)" hub::Vector{Int} "a mapping from star indices to the vector of their spokes" spokes::Vector{Vector{Int}} @@ -141,9 +148,9 @@ function StarSet(star, hub) spokes = [Int[] for s in eachindex(hub)] for ((i, j), s) in pairs(star) h = hub[s] - if i == abs(h) + if i == h push!(spokes[s], j) - elseif j == abs(h) + elseif j == h push!(spokes[s], i) end end @@ -200,7 +207,7 @@ function _update_stars!( hub[star[vq]] = v # this may already be true star[vw] = star[vq] else # vw forms a new star - push!(hub, -max(v, w)) # hub is undefined so we set it to a negative value, but it allows us to remember one of the two vertices + push!(hub, -max(v, w)) # star is trivial (composed only of two vertices) so we set the hub to a negative value, but it allows us to choose one of the two vertices star[vw] = length(hub) end end @@ -247,7 +254,7 @@ function symmetric_coefficient( end """ - acyclic_coloring(g::AdjacencyGraph, order::AbstractOrder) + acyclic_coloring(g::AdjacencyGraph, order::AbstractOrder; postprocessing::Bool) Compute an acyclic coloring of all vertices in the adjacency graph `g` and return a tuple `(color, tree_set)`, where @@ -258,6 +265,8 @@ An _acyclic coloring_ is a distance-1 coloring with the further restriction that The vertices are colored in a greedy fashion, following the `order` supplied. +If `postprocessing=true`, some colors might be replaced with `0` (the "neutral" color) as long as they are not needed during decompression. + # See also - [`AdjacencyGraph`](@ref) @@ -267,7 +276,7 @@ The vertices are colored in a greedy fashion, following the `order` supplied. > [_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, order::AbstractOrder) +function acyclic_coloring(g::AdjacencyGraph, order::AbstractOrder; postprocessing::Bool) # Initialize data structures nv = nb_vertices(g) ne = nb_edges(g) @@ -319,8 +328,11 @@ function acyclic_coloring(g::AdjacencyGraph, order::AbstractOrder) for edge in forest.revmap find_root!(forest, edge) end - - return color, TreeSet(forest) + tree_set = TreeSet(forest, nb_vertices(g)) + if postprocessing + postprocess!(color, tree_set, g) + end + return color, tree_set end function _prevent_cycle!( @@ -401,4 +413,163 @@ $TYPEDFIELDS struct TreeSet "a forest of two-colored trees" forest::DisjointSets{Tuple{Int,Int}} + vertices_by_tree::Vector{Vector{Int}} + reverse_bfs_orders::Vector{Vector{Tuple{Int,Int}}} +end + +function TreeSet(forest::DisjointSets{Tuple{Int,Int}}, nvertices::Int) + # forest is a structure DisjointSets from DataStructures.jl + # - forest.intmap: a dictionary that maps an edge (i, j) to an integer k + # - forest.revmap: a dictionary that does the reverse of intmap, mapping an integer k to an edge (i, j) + # - forest.internal.ngroups: the number of trees in the forest + ntrees = forest.internal.ngroups + + # dictionary that maps a tree's root to the index of the tree + roots = Dict{Int,Int}() + + # vector of dictionaries where each dictionary stores the neighbors of each vertex in a tree + trees = [Dict{Int,Vector{Int}}() for i in 1:ntrees] + + # counter of the number of roots found + k = 0 + for edge in forest.revmap + i, j = edge + # forest has already been compressed so this doesn't change its state + # I wanted to use find_root but it is deprecated + root_edge = find_root!(forest, edge) + root = forest.intmap[root_edge] + + # Update roots + if !haskey(roots, root) + k += 1 + roots[root] = k + end + + # index of the tree T that contains this edge + index_tree = roots[root] + + # Update the neighbors of i in the tree T + if !haskey(trees[index_tree], i) + trees[index_tree][i] = [j] + else + push!(trees[index_tree][i], j) + end + + # Update the neighbors of j in the tree T + if !haskey(trees[index_tree], j) + trees[index_tree][j] = [i] + else + push!(trees[index_tree][j], i) + end + end + + # degrees is a vector of integers that stores the degree of each vertex in a tree + degrees = Vector{Int}(undef, nvertices) + + # list of vertices for each tree in the forest + vertices_by_tree = [collect(keys(trees[i])) for i in 1:ntrees] + + # reverse breadth first (BFS) traversal order for each tree in the forest + reverse_bfs_orders = [Tuple{Int,Int}[] for i in 1:ntrees] + + # nvmax is the number of vertices of the biggest tree in the forest + nvmax = mapreduce(length, max, vertices_by_tree; init=0) + + # Create a queue with a fixed size nvmax + queue = Vector{Int}(undef, nvmax) + + for k in 1:ntrees + tree = trees[k] + + # Initialize the queue to store the leaves + queue_start = 1 + queue_end = 0 + + # compute the degree of each vertex in the tree + for (vertex, neighbors) in tree + degree = length(neighbors) + degrees[vertex] = degree + + # the vertex is a leaf + if degree == 1 + queue_end += 1 + queue[queue_end] = vertex + end + end + + # continue until all leaves are treated + while queue_start <= queue_end + leaf = queue[queue_start] + queue_start += 1 + + # Mark the vertex as removed + degrees[leaf] = 0 + + for neighbor in tree[leaf] + if degrees[neighbor] != 0 + # (leaf, neighbor) represents the next edge to visit during decompression + push!(reverse_bfs_orders[k], (leaf, neighbor)) + + # reduce the degree of the neighbor + degrees[neighbor] -= 1 + + # check if the neighbor is now a leaf + if degrees[neighbor] == 1 + queue_end += 1 + queue[queue_end] = neighbor + end + end + end + end + end + + return TreeSet(forest, vertices_by_tree, reverse_bfs_orders) +end + +## Postprocessing, mirrors decompression code + +function postprocess!( + color::AbstractVector{<:Integer}, + star_or_tree_set::Union{StarSet,TreeSet}, + g::AdjacencyGraph, +) + (; S) = g + # flag which colors are actually used during decompression + color_used = zeros(Bool, maximum(color)) + + # nonzero diagonal coefficients force the use of their respective color (there can be no neutral colors if the diagonal is fully nonzero) + for i in axes(S, 1) + if !iszero(S[i, i]) + color_used[color[i]] = true + end + end + if star_or_tree_set isa StarSet + # only the colors of the hubs are used + (; hub) = star_or_tree_set + for s in eachindex(hub) + j = hub[s] + color_used[color[j]] = true + end + else + # only the colors of non-leaf vertices are used + (; reverse_bfs_orders) = star_or_tree_set + for k in eachindex(reverse_bfs_orders) + for (i, j) in reverse_bfs_orders[k] + color_used[color[j]] = true + end + end + end + # if at least one of the colors is useless, modify the color assignments of vertices + if any(!, color_used) + # assign the neutral color to every vertex with a useless color + for i in eachindex(color) + ci = color[i] + if !color_used[ci] + color[i] = 0 + end + end + # remap colors to decrease the highest one by filling gaps + color .= remap_colors(color)[1] + end + return color end diff --git a/src/decompression.jl b/src/decompression.jl index 6cba735..9b8d397 100644 --- a/src/decompression.jl +++ b/src/decompression.jl @@ -46,16 +46,24 @@ function compress end function compress(A, result::AbstractColoringResult{structure,:column}) where {structure} group = column_groups(result) - B = stack(group; dims=2) do g - dropdims(sum(A[:, g]; dims=2); dims=2) + if isempty(group) + B = similar(A, size(A, 1), 0) + else + B = stack(group; dims=2) do g + dropdims(sum(A[:, g]; dims=2); dims=2) + end end return B end function compress(A, result::AbstractColoringResult{structure,:row}) where {structure} group = row_groups(result) - B = stack(group; dims=1) do g - dropdims(sum(A[g, :]; dims=1); dims=1) + if isempty(group) + B = similar(A, 0, size(A, 2)) + else + B = stack(group; dims=1) do g + dropdims(sum(A[g, :]; dims=1); dims=1) + end end return B end @@ -65,11 +73,19 @@ function compress( ) where {structure} row_group = row_groups(result) column_group = column_groups(result) - Br = stack(row_group; dims=1) do g - dropdims(sum(A[g, :]; dims=1); dims=1) + if isempty(row_group) + Br = similar(A, 0, size(A, 2)) + else + Br = stack(row_group; dims=1) do g + dropdims(sum(A[g, :]; dims=1); dims=1) + end end - Bc = stack(column_group; dims=2) do g - dropdims(sum(A[:, g]; dims=2); dims=2) + if isempty(column_group) + Bc = similar(A, size(A, 1), 0) + else + Bc = stack(column_group; dims=2) do g + dropdims(sum(A[:, g]; dims=2); dims=2) + end end return Br, Bc end @@ -410,7 +426,7 @@ function decompress!( end end for s in eachindex(hub, spokes) - j = abs(hub[s]) + j = hub[s] cj = color[j] for i in spokes[s] if in_triangle(i, j, uplo) diff --git a/src/interface.jl b/src/interface.jl index f2b601f..2bcb488 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -69,10 +69,11 @@ It is passed as an argument to the main function [`coloring`](@ref). # Constructors - GreedyColoringAlgorithm{decompression}(order=NaturalOrder()) - GreedyColoringAlgorithm(order=NaturalOrder(); decompression=:direct) + GreedyColoringAlgorithm{decompression}(order=NaturalOrder(); postprocessing=false) + GreedyColoringAlgorithm(order=NaturalOrder(); postprocessing=false, decompression=:direct) - `order::AbstractOrder`: the order in which the columns or rows are colored, which can impact the number of colors. +- `postprocessing::Bool`: whether or not the coloring will be refined by assigning the neutral color `0` to some vertices. - `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 @@ -96,20 +97,23 @@ See their respective docstrings for details. struct GreedyColoringAlgorithm{decompression,O<:AbstractOrder} <: ADTypes.AbstractColoringAlgorithm order::O + postprocessing::Bool end function GreedyColoringAlgorithm{decompression}( - order::AbstractOrder=NaturalOrder() + order::AbstractOrder=NaturalOrder(); postprocessing::Bool=false ) where {decompression} check_valid_algorithm(decompression) - return GreedyColoringAlgorithm{decompression,typeof(order)}(order) + return GreedyColoringAlgorithm{decompression,typeof(order)}(order, postprocessing) end function GreedyColoringAlgorithm( - order::AbstractOrder=NaturalOrder(); decompression::Symbol=:direct + order::AbstractOrder=NaturalOrder(); + postprocessing::Bool=false, + decompression::Symbol=:direct, ) check_valid_algorithm(decompression) - return GreedyColoringAlgorithm{decompression,typeof(order)}(order) + return GreedyColoringAlgorithm{decompression,typeof(order)}(order, postprocessing) end """ @@ -208,7 +212,7 @@ function coloring( decompression_eltype::Type=Float64, ) ag = AdjacencyGraph(A) - color, star_set = star_coloring(ag, algo.order) + color, star_set = star_coloring(ag, algo.order; postprocessing=algo.postprocessing) return StarSetColoringResult(A, ag, color, star_set) end @@ -219,7 +223,7 @@ function coloring( decompression_eltype::Type=Float64, ) ag = AdjacencyGraph(A) - color, tree_set = acyclic_coloring(ag, algo.order) + color, tree_set = acyclic_coloring(ag, algo.order; postprocessing=algo.postprocessing) return TreeSetColoringResult(A, ag, color, tree_set, decompression_eltype) end @@ -243,10 +247,12 @@ function coloring( ] # TODO: slow ag = AdjacencyGraph(A_and_Aᵀ) if decompression == :direct - color, star_set = star_coloring(ag, algo.order) + color, star_set = star_coloring(ag, algo.order; postprocessing=algo.postprocessing) symmetric_result = StarSetColoringResult(A_and_Aᵀ, ag, color, star_set) else - color, tree_set = acyclic_coloring(ag, algo.order) + color, tree_set = acyclic_coloring( + ag, algo.order; postprocessing=algo.postprocessing + ) symmetric_result = TreeSetColoringResult( A_and_Aᵀ, ag, color, tree_set, decompression_eltype ) @@ -270,6 +276,7 @@ end function ADTypes.symmetric_coloring(A::AbstractMatrix, algo::GreedyColoringAlgorithm) ag = AdjacencyGraph(A) - color, star_set = star_coloring(ag, algo.order) + # never postprocess because end users do not expect zeros + color, star_set = star_coloring(ag, algo.order; postprocessing=false) return color end diff --git a/src/result.jl b/src/result.jl index 13a867d..e3716cb 100644 --- a/src/result.jl +++ b/src/result.jl @@ -44,23 +44,23 @@ function row_colors end """ column_groups(result::AbstractColoringResult) -Return a vector `group` such that for every color `c`, `group[c]` contains the indices of all columns that are colored with `c`. +Return a vector `group` such that for every non-zero color `c`, `group[c]` contains the indices of all columns that are colored with `c`. """ function column_groups end """ row_groups(result::AbstractColoringResult) -Return a vector `group` such that for every color `c`, `group[c]` contains the indices of all rows that are colored with `c`. +Return a vector `group` such that for every non-zero color `c`, `group[c]` contains the indices of all rows that are colored with `c`. """ function row_groups end """ ncolors(result::AbstractColoringResult) -Return the number of different colors used to color the matrix. +Return the number of different non-zero colors used to color the matrix. -For bidirectional partitions, this number is the sum of the number of row colors and the number of column colors. +For bidirectional partitions, this number is the sum of the number of non-zero row colors and the number of non-zero column colors. """ function ncolors(res::AbstractColoringResult{structure,:column}) where {structure} return length(column_groups(res)) @@ -77,32 +77,35 @@ end """ group_by_color(color::Vector{Int}) -Create a color-indexed vector `group` such that `i ∈ group[c]` iff `color[i] == c`. +Create a color-indexed vector `group` such that `i ∈ group[c]` iff `color[i] == c` for all `c > 0`. -Assumes the colors are contiguously numbered from `1` to some `cmax`. +Assumes the colors are contiguously numbered from `0` to some `cmax`. """ function group_by_color(color::AbstractVector{<:Integer}) cmin, cmax = extrema(color) - @assert cmin >= 1 + @assert cmin >= 0 # Compute group sizes and offsets for a joint storage group_sizes = zeros(Int, cmax) # allocation 1, size cmax for c in color - group_sizes[c] += 1 + if c > 0 + group_sizes[c] += 1 + end end group_offsets = cumsum(group_sizes) # allocation 2, size cmax # Concatenate all groups inside a single vector - group_flat = similar(color) # allocation 3, size n + group_flat = similar(color, sum(group_sizes)) # allocation 3, size <= n for (k, c) in enumerate(color) - i = group_offsets[c] - group_sizes[c] + 1 - group_flat[i] = k - group_sizes[c] -= 1 + if c > 0 + i = group_offsets[c] - group_sizes[c] + 1 + group_flat[i] = k + group_sizes[c] -= 1 + end end # Create views into contiguous blocks of the group vector - group = Vector{typeof(view(group_flat, 1:1))}(undef, cmax) # allocation 4, size cmax - for c in 1:cmax + group = map(1:cmax) do c i = 1 + (c == 1 ? 0 : group_offsets[c - 1]) j = group_offsets[c] - group[c] = view(group_flat, i:j) + view(group_flat, i:j) end return group end @@ -291,6 +294,7 @@ function TreeSetColoringResult( tree_set::TreeSet, decompression_eltype::Type{R}, ) where {R} + (; vertices_by_tree, reverse_bfs_orders) = tree_set S = ag.S nvertices = length(color) group = group_by_color(color) @@ -300,7 +304,6 @@ function TreeSetColoringResult( diagonal_nzind = Int[] ndiag = 0 - n = size(S, 1) rv = rowvals(S) for j in axes(S, 2) for k in nzrange(S, j) @@ -318,142 +321,40 @@ function TreeSetColoringResult( lower_triangle_offsets = Vector{Int}(undef, nedges) upper_triangle_offsets = Vector{Int}(undef, nedges) - # forest is a structure DisjointSets from DataStructures.jl - # - forest.intmap: a dictionary that maps an edge (i, j) to an integer k - # - forest.revmap: a dictionary that does the reverse of intmap, mapping an integer k to an edge (i, j) - # - forest.internal.ngroups: the number of trees in the forest - forest = tree_set.forest - ntrees = forest.internal.ngroups - - # dictionary that maps a tree's root to the index of the tree - roots = Dict{Int,Int}() - - # vector of dictionaries where each dictionary stores the neighbors of each vertex in a tree - trees = [Dict{Int,Vector{Int}}() for i in 1:ntrees] - - # counter of the number of roots found - k = 0 - for edge in forest.revmap - i, j = edge - # forest has already been compressed so this doesn't change its state - # I wanted to use find_root but it is deprecated - root_edge = find_root!(forest, edge) - root = forest.intmap[root_edge] - - # Update roots - if !haskey(roots, root) - k += 1 - roots[root] = k - end - - # index of the tree T that contains this edge - index_tree = roots[root] - - # Update the neighbors of i in the tree T - if !haskey(trees[index_tree], i) - trees[index_tree][i] = [j] - else - push!(trees[index_tree][i], j) - end - - # Update the neighbors of j in the tree T - if !haskey(trees[index_tree], j) - trees[index_tree][j] = [i] - else - push!(trees[index_tree][j], i) - end - end - - # degrees is a vector of integers that stores the degree of each vertex in a tree - degrees = Vector{Int}(undef, nvertices) - - # list of vertices for each tree in the forest - vertices_by_tree = [collect(keys(trees[i])) for i in 1:ntrees] - - # reverse breadth first (BFS) traversal order for each tree in the forest - reverse_bfs_orders = [Tuple{Int,Int}[] for i in 1:ntrees] - - # nvmax is the number of vertices of the biggest tree in the forest - nvmax = mapreduce(length, max, vertices_by_tree; init=0) - - # Create a queue with a fixed size nvmax - queue = Vector{Int}(undef, nvmax) - # Index in lower_triangle_offsets and upper_triangle_offsets index_offsets = 0 - for k in 1:ntrees - tree = trees[k] - - # Initialize the queue to store the leaves - queue_start = 1 - queue_end = 0 - - # compute the degree of each vertex in the tree - for (vertex, neighbors) in tree - degree = length(neighbors) - degrees[vertex] = degree - - # the vertex is a leaf - if degree == 1 - queue_end += 1 - queue[queue_end] = vertex - end - end - - # continue until all leaves are treated - while queue_start <= queue_end - leaf = queue[queue_start] - queue_start += 1 - - # Mark the vertex as removed - degrees[leaf] = 0 - - for neighbor in tree[leaf] - if degrees[neighbor] != 0 - # (leaf, neighbor) represents the next edge to visit during decompression - push!(reverse_bfs_orders[k], (leaf, neighbor)) - - # reduce the degree of the neighbor - degrees[neighbor] -= 1 - - # check if the neighbor is now a leaf - if degrees[neighbor] == 1 - queue_end += 1 - queue[queue_end] = neighbor - end - - # Update lower_triangle_offsets and upper_triangle_offsets - i = leaf - j = neighbor - col_i = view(rv, nzrange(S, i)) - col_j = view(rv, nzrange(S, j)) - index_offsets += 1 - - #! format: off - # S[i,j] is in the lower triangular part of S - if in_triangle(i, j, :L) - # uplo = :L or uplo = :F - # S[i,j] is stored at index_ij = (S.colptr[j+1] - offset_L) in S.nzval - lower_triangle_offsets[index_offsets] = length(col_j) - searchsortedfirst(col_j, i) + 1 - - # uplo = :U or uplo = :F - # S[j,i] is stored at index_ji = (S.colptr[i] + offset_U) in S.nzval - upper_triangle_offsets[index_offsets] = searchsortedfirst(col_i, j)::Int - 1 - - # S[i,j] is in the upper triangular part of S - else - # uplo = :U or uplo = :F - # S[i,j] is stored at index_ij = (S.colptr[j] + offset_U) in S.nzval - upper_triangle_offsets[index_offsets] = searchsortedfirst(col_j, i)::Int - 1 - - # uplo = :L or uplo = :F - # S[j,i] is stored at index_ji = (S.colptr[i+1] - offset_L) in S.nzval - lower_triangle_offsets[index_offsets] = length(col_i) - searchsortedfirst(col_i, j) + 1 - end - #! format: on - end + for k in eachindex(reverse_bfs_orders) + for (leaf, neighbor) in reverse_bfs_orders[k] + # Update lower_triangle_offsets and upper_triangle_offsets + i = leaf + j = neighbor + col_i = view(rv, nzrange(S, i)) + col_j = view(rv, nzrange(S, j)) + index_offsets += 1 + + #! format: off + # S[i,j] is in the lower triangular part of S + if in_triangle(i, j, :L) + # uplo = :L or uplo = :F + # S[i,j] is stored at index_ij = (S.colptr[j+1] - offset_L) in S.nzval + lower_triangle_offsets[index_offsets] = length(col_j) - searchsortedfirst(col_j, i) + 1 + + # uplo = :U or uplo = :F + # S[j,i] is stored at index_ji = (S.colptr[i] + offset_U) in S.nzval + upper_triangle_offsets[index_offsets] = searchsortedfirst(col_i, j)::Int - 1 + + # S[i,j] is in the upper triangular part of S + else + # uplo = :U or uplo = :F + # S[i,j] is stored at index_ij = (S.colptr[j] + offset_U) in S.nzval + upper_triangle_offsets[index_offsets] = searchsortedfirst(col_j, i)::Int - 1 + + # uplo = :L or uplo = :F + # S[j,i] is stored at index_ji = (S.colptr[i+1] - offset_L) in S.nzval + lower_triangle_offsets[index_offsets] = length(col_i) - searchsortedfirst(col_i, j) + 1 end + #! format: on end end @@ -528,10 +429,14 @@ function LinearSystemColoringResult( for (l, (i, j)) in enumerate(strict_upper_nonzero_inds) ci = color[i] cj = color[j] - ki = (ci - 1) * n + j # A[i, j] appears in B[j, ci] - kj = (cj - 1) * n + i # A[i, j] appears in B[i, cj] - T[ki, l] = 1 - T[kj, l] = 1 + if ci > 0 + ki = (ci - 1) * n + j # A[i, j] appears in B[j, ci] + T[ki, l] = 1 + end + if cj > 0 + kj = (cj - 1) * n + i # A[i, j] appears in B[i, cj] + T[kj, l] = 1 + end end T_factorization = factorize(T) @@ -553,7 +458,7 @@ end """ remap_colors(color::Vector{Int}) -Renumber the colors in `color` using their index in the vector `sort(unique(color))`, so that they are forced to go from `1` to some `cmax` contiguously. +Renumber the colors in `color` using their index in the vector `sort(unique(color))`, so that the non-zero colors are forced to go from `1` to some `cmax` contiguously. Return a tuple `(remapped_colors, color_to_ind)` such that `remapped_colors` is a vector containing the renumbered colors and `color_to_ind` is a dictionary giving the translation between old and new color numberings. @@ -562,7 +467,8 @@ For all vertex indices `i` we have: remapped_color[i] = color_to_ind[color[i]] """ function remap_colors(color::Vector{Int}) - color_to_ind = Dict(c => i for (i, c) in enumerate(sort(unique(color)))) + color_to_ind = Dict(c => i for (i, c) in enumerate(sort(unique(color[color .> 0])))) + color_to_ind[0] = 0 remapped_colors = [color_to_ind[c] for c in color] return remapped_colors, color_to_ind end diff --git a/test/check.jl b/test/check.jl index 60e7ecc..4e9e5a6 100644 --- a/test/check.jl +++ b/test/check.jl @@ -100,9 +100,9 @@ end @test_logs ( :warn, """ -For coefficient (i=2, j=3) with column colors (ci=3, cj=1): -- In color ci=3, columns [2, 4] all have nonzeros in row j=3. -- In color cj=1, columns [1, 3, 5, 6] all have nonzeros in row i=2. +For coefficient (i=2, j=3) with colors (ci=3, cj=1): +- In row color ci=3, rows [2, 4] all have nonzeros in column j=3. +- In column color cj=1, columns [1, 3, 5, 6] all have nonzeros in row i=2. """, ) symmetrically_orthogonal_columns(A, [1, 3, 1, 3, 1, 1]; verbose=true) @@ -151,9 +151,36 @@ end @test_logs ( :warn, """ -For coefficient (i=1, j=1) with row color ci=1 and column color cj=1: +For coefficient (i=1, j=1) with colors (ci=1, cj=1): - In row color ci=1, rows [1, 2, 3] all have nonzeros in column j=1. - In column color cj=1, columns [1, 2, 3, 4] all have nonzeros in row i=1. """, ) !structurally_biorthogonal(A, [1, 1, 1, 2], [1, 1, 1, 1, 2]; verbose=true) + + @test_logs ( + :warn, + """ +For coefficient (i=1, j=2) with colors (ci=0, cj=2): +- Row color ci=0 is neutral. +- In column color cj=2, columns [2, 3, 4] all have nonzeros in row i=1. +""", + ) structurally_biorthogonal(A, [0, 2, 2, 3], [1, 2, 2, 2, 3], verbose=true) + + @test_logs ( + :warn, + """ +For coefficient (i=2, j=1) with colors (ci=2, cj=0): +- In row color ci=2, rows [2, 3] all have nonzeros in column j=1. +- Column color cj=0 is neutral. +""", + ) structurally_biorthogonal(A, [1, 2, 2, 3], [0, 2, 2, 2, 3], verbose=true) + + @test_logs ( + :warn, + """ +For coefficient (i=1, j=1) with colors (ci=0, cj=0): +- Row color ci=0 is neutral. +- Column color cj=0 is neutral. +""", + ) structurally_biorthogonal(A, [0, 2, 2, 3], [0, 2, 2, 2, 3], verbose=true) end diff --git a/test/random.jl b/test/random.jl index 06e8565..02c1a3a 100644 --- a/test/random.jl +++ b/test/random.jl @@ -51,46 +51,70 @@ end; @testset "Symmetric coloring & direct decompression" begin problem = ColoringProblem(; structure=:symmetric, partition=:column) - algo = GreedyColoringAlgorithm(; decompression=:direct) - @testset "$((; n, p))" for (n, p) in symmetric_params - A0 = sparse(Symmetric(sprand(rng, n, n, p))) - color0 = symmetric_coloring(A0, algo) - test_coloring_decompression(A0, problem, algo; color0) + @testset for algo in ( + GreedyColoringAlgorithm(; postprocessing=false, decompression=:direct), + GreedyColoringAlgorithm(; postprocessing=true, decompression=:direct), + ) + @testset "$((; n, p))" for (n, p) in symmetric_params + A0 = sparse(Symmetric(sprand(rng, n, n, p))) + color0 = algo.postprocessing ? nothing : symmetric_coloring(A0, algo) + test_coloring_decompression(A0, problem, algo; color0) + end end end; @testset "Symmetric coloring & substitution decompression" begin problem = ColoringProblem(; structure=:symmetric, partition=:column) - algo = GreedyColoringAlgorithm(; decompression=:substitution) - @testset "$((; n, p))" for (n, p) in symmetric_params - A0 = sparse(Symmetric(sprand(rng, n, n, p))) - # TODO: find tests for recoverability - test_coloring_decompression(A0, problem, algo) + @testset for algo in ( + GreedyColoringAlgorithm(; postprocessing=false, decompression=:substitution), + GreedyColoringAlgorithm(; postprocessing=true, decompression=:substitution), + ) + @testset "$((; n, p))" for (n, p) in symmetric_params + A0 = sparse(Symmetric(sprand(rng, n, n, p))) + # TODO: find tests for recoverability + test_coloring_decompression(A0, problem, algo) + end end end; @testset "Bicoloring & direct decompression" begin problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional) - algo = GreedyColoringAlgorithm(RandomOrder(rng); decompression=:direct) - @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params - A0 = sprand(rng, m, n, p) - test_bicoloring_decompression(A0, problem, algo) - end - @testset "$((; n, p))" for (n, p) in symmetric_params - A0 = sparse(Symmetric(sprand(rng, n, n, p))) - test_bicoloring_decompression(A0, problem, algo) + @testset for algo in ( + GreedyColoringAlgorithm( + RandomOrder(rng); postprocessing=false, decompression=:direct + ), + GreedyColoringAlgorithm( + RandomOrder(rng); postprocessing=true, decompression=:direct + ), + ) + @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params + A0 = sprand(rng, m, n, p) + test_bicoloring_decompression(A0, problem, algo) + end + @testset "$((; n, p))" for (n, p) in symmetric_params + A0 = sparse(Symmetric(sprand(rng, n, n, p))) + test_bicoloring_decompression(A0, problem, algo) + end end end; @testset "Bicoloring & substitution decompression" begin problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional) - algo = GreedyColoringAlgorithm(RandomOrder(rng); decompression=:substitution) - @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params - A0 = sprand(rng, m, n, p) - test_bicoloring_decompression(A0, problem, algo) - end - @testset "$((; n, p))" for (n, p) in symmetric_params - A0 = sparse(Symmetric(sprand(rng, n, n, p))) - test_bicoloring_decompression(A0, problem, algo) + @testset for algo in ( + GreedyColoringAlgorithm( + RandomOrder(rng); postprocessing=false, decompression=:substitution + ), + GreedyColoringAlgorithm( + RandomOrder(rng); postprocessing=true, decompression=:substitution + ), + ) + @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params + A0 = sprand(rng, m, n, p) + test_bicoloring_decompression(A0, problem, algo) + end + @testset "$((; n, p))" for (n, p) in symmetric_params + A0 = sparse(Symmetric(sprand(rng, n, n, p))) + test_bicoloring_decompression(A0, problem, algo) + end end end; diff --git a/test/result.jl b/test/result.jl index 288c129..661297f 100644 --- a/test/result.jl +++ b/test/result.jl @@ -9,5 +9,25 @@ using Test @test all(1:maximum(color)) do c all(color[group[c]] .== c) && issorted(group[c]) end + + color = rand(0:cmax, n) + group = group_by_color(color) + @test length(group) == maximum(color) + @test all(1:maximum(color)) do c + all(color[group[c]] .== c) && issorted(group[c]) + end end end + +@testset "Empty compression" begin + A = rand(10, 10) + color = zeros(Int, 10) + problem = ColoringProblem{:nonsymmetric,:column}() + algo = ConstantColoringAlgorithm(A, color; partition=:column) + B = compress(A, coloring(A, problem, algo)) + @test size(B, 2) == 0 + problem = ColoringProblem{:nonsymmetric,:row}() + algo = ConstantColoringAlgorithm(A, color; partition=:row) + B = compress(A, coloring(A, problem, algo)) + @test size(B, 1) == 0 +end diff --git a/test/suitesparse.jl b/test/suitesparse.jl index 3950629..4f67995 100644 --- a/test/suitesparse.jl +++ b/test/suitesparse.jl @@ -95,7 +95,7 @@ what_table_41_42 = CSV.read( @test nb_edges(ag) == row[:E] @test maximum_degree(ag) == row[:Δ] @test minimum_degree(ag) == row[:δ] - color_N, _ = star_coloring(ag, NaturalOrder()) + color_N, _ = star_coloring(ag, NaturalOrder(); postprocessing=false) @test_skip row[:KS1] <= length(unique(color_N)) <= row[:KS2] # TODO: find better yield() end diff --git a/test/utils.jl b/test/utils.jl index d58461c..b367226 100644 --- a/test/utils.jl +++ b/test/utils.jl @@ -86,6 +86,7 @@ function test_coloring_decompression( A2 = respectful_similar(A) A2 .= zero(eltype(A2)) for c in unique(color) + c == 0 && continue if partition == :column decompress_single_color!(A2, B[:, c], c, result) elseif partition == :row @@ -125,6 +126,7 @@ function test_coloring_decompression( A4both .= zero(eltype(A)) for c in unique(color) + c == 0 && continue decompress_single_color!(A4upper, B[:, c], c, result, :U) decompress_single_color!(A4lower, B[:, c], c, result, :L) decompress_single_color!(A4both, B[:, c], c, result, :F) @@ -149,6 +151,9 @@ function test_coloring_decompression( @testset "Coherence between all colorings" begin @test all(color_vec .== Ref(color_vec[1])) + if !all(color_vec .== Ref(color_vec[1])) + @show color_vec + end end end @@ -168,8 +173,8 @@ function test_bicoloring_decompression( end Br, Bc = compress(A, result) row_color, column_color = row_colors(result), column_colors(result) - @test size(Br, 1) == length(unique(row_color)) - @test size(Bc, 2) == length(unique(column_color)) + @test size(Br, 1) == length(unique(row_color[row_color .> 0])) + @test size(Bc, 2) == length(unique(column_color[column_color .> 0])) @test ncolors(result) == size(Br, 1) + size(Bc, 2) if decompression == :direct