diff --git a/.github/workflows/Test.yml b/.github/workflows/Test.yml index ef161b41..3e17ae97 100644 --- a/.github/workflows/Test.yml +++ b/.github/workflows/Test.yml @@ -20,7 +20,7 @@ jobs: runs-on: ubuntu-latest strategy: matrix: - julia-version: ['1.10', '1'] + julia-version: ['1.10', '1.11'] steps: - uses: actions/checkout@v5 @@ -40,4 +40,4 @@ jobs: with: files: lcov.info token: ${{ secrets.CODECOV_TOKEN }} - fail_ci_if_error: false \ No newline at end of file + fail_ci_if_error: false diff --git a/Project.toml b/Project.toml index 6656e125..612ff5a1 100644 --- a/Project.toml +++ b/Project.toml @@ -15,11 +15,14 @@ SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" CUDA = "052768ef-5323-5732-b1bb-66c8b64840ba" CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" +MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee" [extensions] SparseMatrixColoringsCUDAExt = "CUDA" SparseMatrixColoringsCliqueTreesExt = "CliqueTrees" SparseMatrixColoringsColorsExt = "Colors" +SparseMatrixColoringsJuMPExt = ["JuMP", "MathOptInterface"] [compat] ADTypes = "1.2.1" @@ -27,7 +30,9 @@ CUDA = "5.8.2" CliqueTrees = "1" Colors = "0.12.11, 0.13" DocStringExtensions = "0.8,0.9" +JuMP = "1.29.1" LinearAlgebra = "<0.0.1, 1" +MathOptInterface = "1.45.0" PrecompileTools = "1.2.1" Random = "<0.0.1, 1" SparseArrays = "<0.0.1, 1" diff --git a/docs/src/api.md b/docs/src/api.md index c121d1e3..f90c888f 100644 --- a/docs/src/api.md +++ b/docs/src/api.md @@ -17,8 +17,14 @@ SparseMatrixColorings coloring fast_coloring ColoringProblem +``` + +## Coloring algorithms + +```@docs GreedyColoringAlgorithm ConstantColoringAlgorithm +OptimalColoringAlgorithm ``` ## Result analysis diff --git a/ext/SparseMatrixColoringsJuMPExt.jl b/ext/SparseMatrixColoringsJuMPExt.jl new file mode 100644 index 00000000..ca4f5e7b --- /dev/null +++ b/ext/SparseMatrixColoringsJuMPExt.jl @@ -0,0 +1,80 @@ +module SparseMatrixColoringsJuMPExt + +using ADTypes: ADTypes +using JuMP: + Model, + is_solved_and_feasible, + optimize!, + primal_status, + set_silent, + set_start_value, + value, + @variable, + @constraint, + @objective +using JuMP +import MathOptInterface as MOI +using SparseMatrixColorings: + BipartiteGraph, OptimalColoringAlgorithm, nb_vertices, neighbors, pattern, vertices + +function optimal_distance2_coloring( + bg::BipartiteGraph, + ::Val{side}, + optimizer::O; + silent::Bool=true, + assert_solved::Bool=true, +) where {side,O} + other_side = 3 - side + n = nb_vertices(bg, Val(side)) + model = Model(optimizer) + silent && set_silent(model) + # one variable per vertex to color, removing some renumbering symmetries + @variable(model, 1 <= color[i=1:n] <= i, Int) + # one variable to count the number of distinct colors + @variable(model, ncolors, Int) + @constraint(model, [ncolors; color] in MOI.CountDistinct(n + 1)) + # distance-2 coloring: neighbors of the same vertex must have distinct colors + for i in vertices(bg, Val(other_side)) + neigh = neighbors(bg, Val(other_side), i) + @constraint(model, color[neigh] in MOI.AllDifferent(length(neigh))) + end + # minimize the number of distinct colors (can't use maximum because they are not necessarily numbered contiguously) + @objective(model, Min, ncolors) + # actual solving step where time is spent + optimize!(model) + if assert_solved + # assert feasibility and optimality + @assert is_solved_and_feasible(model) + else + # only assert feasibility + @assert primal_status(model) == MOI.FEASIBLE_POINT + end + # native solver solutions are floating point numbers + color_int = round.(Int, value.(color)) + # remap to 1:cmax in case they are not contiguous + true_ncolors = 0 + remap = fill(0, maximum(color_int)) + for c in color_int + if remap[c] == 0 + true_ncolors += 1 + remap[c] = true_ncolors + end + end + return remap[color_int] +end + +function ADTypes.column_coloring(A::AbstractMatrix, algo::OptimalColoringAlgorithm) + bg = BipartiteGraph(A) + return optimal_distance2_coloring( + bg, Val(2), algo.optimizer; algo.silent, algo.assert_solved + ) +end + +function ADTypes.row_coloring(A::AbstractMatrix, algo::OptimalColoringAlgorithm) + bg = BipartiteGraph(A) + return optimal_distance2_coloring( + bg, Val(1), algo.optimizer; algo.silent, algo.assert_solved + ) +end + +end diff --git a/src/SparseMatrixColorings.jl b/src/SparseMatrixColorings.jl index 7c25bac4..fca70229 100644 --- a/src/SparseMatrixColorings.jl +++ b/src/SparseMatrixColorings.jl @@ -56,6 +56,7 @@ include("decompression.jl") include("check.jl") include("examples.jl") include("show_colors.jl") +include("optimal.jl") include("precompile.jl") @@ -64,6 +65,7 @@ export DynamicDegreeBasedOrder, SmallestLast, IncidenceDegree, DynamicLargestFir export PerfectEliminationOrder export ColoringProblem, GreedyColoringAlgorithm, AbstractColoringResult export ConstantColoringAlgorithm +export OptimalColoringAlgorithm export coloring, fast_coloring export column_colors, row_colors, ncolors export column_groups, row_groups diff --git a/src/adtypes.jl b/src/adtypes.jl index 7500b1e2..7c464433 100644 --- a/src/adtypes.jl +++ b/src/adtypes.jl @@ -1,21 +1,33 @@ -function coloring( - A::AbstractMatrix, - ::ColoringProblem{:nonsymmetric,:column}, - algo::ADTypes.NoColoringAlgorithm; - kwargs..., -) - bg = BipartiteGraph(A) - color = convert(Vector{eltype(bg)}, ADTypes.column_coloring(A, algo)) - return ColumnColoringResult(A, bg, color) -end +## From ADTypes to SMC function coloring( A::AbstractMatrix, - ::ColoringProblem{:nonsymmetric,:row}, - algo::ADTypes.NoColoringAlgorithm; - kwargs..., -) - bg = BipartiteGraph(A) - color = convert(Vector{eltype(bg)}, ADTypes.row_coloring(A, algo)) - return RowColoringResult(A, bg, color) + problem::ColoringProblem{structure,partition}, + algo::ADTypes.AbstractColoringAlgorithm; + decompression_eltype::Type{R}=Float64, + symmetric_pattern::Bool=false, +) where {structure,partition,R} + symmetric_pattern = symmetric_pattern || A isa Union{Symmetric,Hermitian} + if structure == :nonsymmetric + if partition == :column + forced_colors = ADTypes.column_coloring(A, algo) + elseif partition == :row + forced_colors = ADTypes.row_coloring(A, algo) + else + # TODO: improve once https://github.com/SciML/ADTypes.jl/issues/69 is done + A_and_Aᵀ, _ = bidirectional_pattern(A; symmetric_pattern) + forced_colors = ADTypes.symmetric_coloring(A_and_Aᵀ, algo) + end + else + forced_colors = ADTypes.symmetric_coloring(A, algo) + end + return _coloring( + WithResult(), + A, + problem, + GreedyColoringAlgorithm(), + R, + symmetric_pattern; + forced_colors, + ) end diff --git a/src/coloring.jl b/src/coloring.jl index 7eec55c0..df204915 100644 --- a/src/coloring.jl +++ b/src/coloring.jl @@ -1,5 +1,10 @@ +struct InvalidColoringError <: Exception end + """ - partial_distance2_coloring(bg::BipartiteGraph, ::Val{side}, vertices_in_order::AbstractVector) + partial_distance2_coloring( + bg::BipartiteGraph, ::Val{side}, vertices_in_order::AbstractVector; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing + ) Compute a distance-2 coloring of the given `side` (`1` or `2`) in the bipartite graph `bg` and return a vector of integer colors. @@ -7,6 +12,8 @@ A _distance-2 coloring_ is such that two vertices have different colors if they The vertices are colored in a greedy fashion, following the order supplied. +The optional `forced_colors` keyword argument is used to enforce predefined vertex colors (e.g. coming from another optimization algorithm) but still run the distance-2 coloring procedure to verify correctness. + # See also - [`BipartiteGraph`](@ref) @@ -17,11 +24,16 @@ The vertices are colored in a greedy fashion, following the order supplied. > [_What Color Is Your Jacobian? Graph Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005), Algorithm 3.2 """ function partial_distance2_coloring( - bg::BipartiteGraph{T}, ::Val{side}, vertices_in_order::AbstractVector{<:Integer} + bg::BipartiteGraph{T}, + ::Val{side}, + vertices_in_order::AbstractVector{<:Integer}; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) where {T,side} color = Vector{T}(undef, nb_vertices(bg, Val(side))) forbidden_colors = Vector{T}(undef, nb_vertices(bg, Val(side))) - partial_distance2_coloring!(color, forbidden_colors, bg, Val(side), vertices_in_order) + partial_distance2_coloring!( + color, forbidden_colors, bg, Val(side), vertices_in_order; forced_colors + ) return color end @@ -30,7 +42,8 @@ function partial_distance2_coloring!( forbidden_colors::AbstractVector{<:Integer}, bg::BipartiteGraph, ::Val{side}, - vertices_in_order::AbstractVector{<:Integer}, + vertices_in_order::AbstractVector{<:Integer}; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) where {side} color .= 0 forbidden_colors .= 0 @@ -44,17 +57,32 @@ function partial_distance2_coloring!( end end end - for i in eachindex(forbidden_colors) - if forbidden_colors[i] != v - color[v] = i - break + if isnothing(forced_colors) + for i in eachindex(forbidden_colors) + if forbidden_colors[i] != v + color[v] = i + break + end + end + else + f = forced_colors[v] + if ( + (f == 0 && length(neighbors(bg, Val(side), v)) > 0) || + (f > 0 && forbidden_colors[f] == v) + ) + throw(InvalidColoringError()) + else + color[v] = f end end end end """ - star_coloring(g::AdjacencyGraph, vertices_in_order::AbstractVector, postprocessing::Bool) + star_coloring( + g::AdjacencyGraph, vertices_in_order::AbstractVector, postprocessing::Bool; + 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 @@ -67,6 +95,8 @@ 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. +The optional `forced_colors` keyword argument is used to enforce predefined vertex colors (e.g. coming from another optimization algorithm) but still run the star coloring procedure to verify correctness and build auxiliary data structures, useful during decompression. + # See also - [`AdjacencyGraph`](@ref) @@ -77,7 +107,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 4.1 """ function star_coloring( - g::AdjacencyGraph{T}, vertices_in_order::AbstractVector{<:Integer}, postprocessing::Bool + g::AdjacencyGraph{T}, + vertices_in_order::AbstractVector{<:Integer}, + postprocessing::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) where {T<:Integer} # Initialize data structures nv = nb_vertices(g) @@ -115,10 +148,18 @@ function star_coloring( end end end - for i in eachindex(forbidden_colors) - if forbidden_colors[i] != v - color[v] = i - break + if isnothing(forced_colors) + for i in eachindex(forbidden_colors) + if forbidden_colors[i] != v + color[v] = i + break + end + end + else + if forbidden_colors[forced_colors[v]] == v # TODO: handle forced_colors[v] == 0 + throw(InvalidColoringError()) + else + color[v] = forced_colors[v] end end _update_stars!(star, hub, g, v, color, first_neighbor) @@ -209,7 +250,10 @@ struct StarSet{T} end """ - acyclic_coloring(g::AdjacencyGraph, vertices_in_order::AbstractVector, postprocessing::Bool) + acyclic_coloring( + g::AdjacencyGraph, vertices_in_order::AbstractVector, postprocessing::Bool; + forced_colors::Union{AbstractVector,Nothing}=nothing + ) Compute an acyclic coloring of all vertices in the adjacency graph `g` and return a tuple `(color, tree_set)`, where @@ -222,6 +266,8 @@ 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. +The optional `forced_colors` keyword argument is used to enforce predefined vertex colors (e.g. coming from another optimization algorithm) but still run the acyclic coloring procedure to verify correctness and build auxiliary data structures, useful during decompression. + # See also - [`AdjacencyGraph`](@ref) @@ -232,7 +278,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; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) where {T<:Integer} # Initialize data structures nv = nb_vertices(g) @@ -271,10 +320,18 @@ function acyclic_coloring( end end end - for i in eachindex(forbidden_colors) - if forbidden_colors[i] != v - color[v] = i - break + if isnothing(forced_colors) + for i in eachindex(forbidden_colors) + if forbidden_colors[i] != v + color[v] = i + break + end + end + else + if forbidden_colors[forced_colors[v]] == v # TODO: handle forced_colors[v] == 0 + throw(InvalidColoringError()) + else + color[v] = forced_colors[v] end end for (w, index_vw) in neighbors_with_edge_indices(g, v) # grow two-colored stars around the vertex v diff --git a/src/constant.jl b/src/constant.jl index d79caab0..56bf4101 100644 --- a/src/constant.jl +++ b/src/constant.jl @@ -4,20 +4,24 @@ Coloring algorithm which always returns the same precomputed vector of colors. Useful when the optimal coloring of a matrix can be determined a priori due to its specific structure (e.g. banded). -It is passed as an argument to the main function [`coloring`](@ref), but will only work if the associated `problem` has `:nonsymmetric` structure. -Indeed, for symmetric coloring problems, we need more than just the vector of colors to allow fast decompression. +It is passed as an argument to the main function [`coloring`](@ref), but will only work if the associated `problem` has a `:column` or `:row` partition. # Constructors ConstantColoringAlgorithm{partition}(matrix_template, color) - ConstantColoringAlgorithm(matrix_template, color; partition=:column) + ConstantColoringAlgorithm{partition,structure}(matrix_template, color) + ConstantColoringAlgorithm( + matrix_template, color; + structure=:nonsymmetric, partition=:column + ) - `partition::Symbol`: either `:row` or `:column`. +- `structure::Symbol`: either `:nonsymmetric` or `:symmetric`. - `matrix_template::AbstractMatrix`: matrix for which the vector of colors was precomputed (the algorithm will only accept matrices of the exact same size). - `color::Vector{<:Integer}`: vector of integer colors, one for each row or column (depending on `partition`). !!! warning - The second constructor (based on keyword arguments) is type-unstable. + The constructor based on keyword arguments is type-unstable if these arguments are not compile-time constants. We do not necessarily verify consistency between the matrix template and the vector of colors, this is the responsibility of the user. @@ -63,71 +67,68 @@ julia> column_colors(result) - [`ADTypes.column_coloring`](@extref ADTypes.column_coloring) - [`ADTypes.row_coloring`](@extref ADTypes.row_coloring) +- [`ADTypes.symmetric_coloring`](@extref ADTypes.symmetric_coloring) """ -struct ConstantColoringAlgorithm{ - partition, - M<:AbstractMatrix, - T<:Integer, - R<:AbstractColoringResult{:nonsymmetric,partition,:direct}, -} <: ADTypes.AbstractColoringAlgorithm +struct ConstantColoringAlgorithm{partition,structure,M<:AbstractMatrix,T<:Integer} <: + ADTypes.AbstractColoringAlgorithm matrix_template::M color::Vector{T} - result::R -end -function ConstantColoringAlgorithm{:column}( - matrix_template::AbstractMatrix, color::Vector{<:Integer} -) - bg = BipartiteGraph(matrix_template) - result = ColumnColoringResult(matrix_template, bg, color) - T, M, R = eltype(bg), typeof(matrix_template), typeof(result) - return ConstantColoringAlgorithm{:column,M,T,R}(matrix_template, color, result) + function ConstantColoringAlgorithm{partition,structure}( + matrix_template::AbstractMatrix, color::Vector{<:Integer} + ) where {partition,structure} + check_valid_problem(structure, partition) + return new{partition,structure,typeof(matrix_template),eltype(color)}( + matrix_template, color + ) + end end -function ConstantColoringAlgorithm{:row}( +function ConstantColoringAlgorithm{partition}( matrix_template::AbstractMatrix, color::Vector{<:Integer} -) - bg = BipartiteGraph(matrix_template) - result = RowColoringResult(matrix_template, bg, color) - T, M, R = eltype(bg), typeof(matrix_template), typeof(result) - return ConstantColoringAlgorithm{:row,M,T,R}(matrix_template, color, result) +) where {partition} + return ConstantColoringAlgorithm{partition,:nonsymmetric}(matrix_template, color) end function ConstantColoringAlgorithm( - matrix_template::AbstractMatrix, color::Vector{<:Integer}; partition::Symbol=:column + matrix_template::AbstractMatrix, + color::Vector{<:Integer}; + structure::Symbol=:nonsymmetric, + partition::Symbol=:column, ) - return ConstantColoringAlgorithm{partition}(matrix_template, color) + return ConstantColoringAlgorithm{partition,structure}(matrix_template, color) end -function coloring( - A::AbstractMatrix, - ::ColoringProblem{:nonsymmetric,partition}, - algo::ConstantColoringAlgorithm{partition}; - decompression_eltype::Type=Float64, - symmetric_pattern::Bool=false, -) where {partition} - (; matrix_template, result) = algo +function check_template(algo::ConstantColoringAlgorithm, A::AbstractMatrix) + (; matrix_template) = algo if size(A) != size(matrix_template) throw( DimensionMismatch( "`ConstantColoringAlgorithm` expected matrix of size $(size(matrix_template)) but got matrix of size $(size(A))", ), ) - else - return result end end function ADTypes.column_coloring( - A::AbstractMatrix, algo::ConstantColoringAlgorithm{:column} + A::AbstractMatrix, algo::ConstantColoringAlgorithm{:column,:nonsymmetric} +) + check_template(algo, A) + return algo.color +end + +function ADTypes.row_coloring( + A::AbstractMatrix, algo::ConstantColoringAlgorithm{:row,:nonsymmetric} ) - problem = ColoringProblem{:nonsymmetric,:column}() - result = coloring(A, problem, algo) - return column_colors(result) + check_template(algo, A) + return algo.color end -function ADTypes.row_coloring(A::AbstractMatrix, algo::ConstantColoringAlgorithm) - problem = ColoringProblem{:nonsymmetric,:row}() - result = coloring(A, problem, algo) - return row_colors(result) +function ADTypes.symmetric_coloring( + A::AbstractMatrix, algo::ConstantColoringAlgorithm{:column,:symmetric} +) + check_template(algo, A) + return algo.color end + +# TODO: handle bidirectional once https://github.com/SciML/ADTypes.jl/issues/69 is done diff --git a/src/interface.jl b/src/interface.jl index 79fe71e9..882b70bc 100644 --- a/src/interface.jl +++ b/src/interface.jl @@ -225,12 +225,13 @@ function _coloring( ::ColoringProblem{:nonsymmetric,:column}, algo::GreedyColoringAlgorithm, decompression_eltype::Type, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) symmetric_pattern = symmetric_pattern || A isa Union{Symmetric,Hermitian} bg = BipartiteGraph(A; symmetric_pattern) vertices_in_order = vertices(bg, Val(2), algo.order) - color = partial_distance2_coloring(bg, Val(2), vertices_in_order) + color = partial_distance2_coloring(bg, Val(2), vertices_in_order; forced_colors) if speed_setting isa WithResult return ColumnColoringResult(A, bg, color) else @@ -244,12 +245,13 @@ function _coloring( ::ColoringProblem{:nonsymmetric,:row}, algo::GreedyColoringAlgorithm, decompression_eltype::Type, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) symmetric_pattern = symmetric_pattern || A isa Union{Symmetric,Hermitian} bg = BipartiteGraph(A; symmetric_pattern) vertices_in_order = vertices(bg, Val(1), algo.order) - color = partial_distance2_coloring(bg, Val(1), vertices_in_order) + color = partial_distance2_coloring(bg, Val(1), vertices_in_order; forced_colors) if speed_setting isa WithResult return RowColoringResult(A, bg, color) else @@ -263,11 +265,14 @@ function _coloring( ::ColoringProblem{:symmetric,:column}, algo::GreedyColoringAlgorithm{:direct}, decompression_eltype::Type, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) ag = AdjacencyGraph(A; has_diagonal=true) vertices_in_order = vertices(ag, algo.order) - color, star_set = star_coloring(ag, vertices_in_order, algo.postprocessing) + color, star_set = star_coloring( + ag, vertices_in_order, algo.postprocessing; forced_colors + ) if speed_setting isa WithResult return StarSetColoringResult(A, ag, color, star_set) else @@ -281,11 +286,14 @@ function _coloring( ::ColoringProblem{:symmetric,:column}, algo::GreedyColoringAlgorithm{:substitution}, decompression_eltype::Type{R}, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + forced_colors::Union{AbstractVector{<:Integer},Nothing}=nothing, ) where {R} ag = AdjacencyGraph(A; has_diagonal=true) vertices_in_order = vertices(ag, algo.order) - color, tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing) + color, tree_set = acyclic_coloring( + ag, vertices_in_order, algo.postprocessing; forced_colors + ) if speed_setting isa WithResult return TreeSetColoringResult(A, ag, color, tree_set, R) else @@ -299,12 +307,15 @@ function _coloring( ::ColoringProblem{:nonsymmetric,:bidirectional}, algo::GreedyColoringAlgorithm{:direct}, decompression_eltype::Type{R}, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + 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) vertices_in_order = vertices(ag, algo.order) - color, star_set = star_coloring(ag, vertices_in_order, algo.postprocessing) + color, star_set = star_coloring( + ag, vertices_in_order, algo.postprocessing; forced_colors + ) if speed_setting isa WithResult symmetric_result = StarSetColoringResult(A_and_Aᵀ, ag, color, star_set) return BicoloringResult(A, ag, symmetric_result, R) @@ -322,12 +333,15 @@ function _coloring( ::ColoringProblem{:nonsymmetric,:bidirectional}, algo::GreedyColoringAlgorithm{:substitution}, decompression_eltype::Type{R}, - symmetric_pattern::Bool, + symmetric_pattern::Bool; + 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) vertices_in_order = vertices(ag, algo.order) - color, tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing) + color, tree_set = acyclic_coloring( + ag, vertices_in_order, algo.postprocessing; forced_colors + ) if speed_setting isa WithResult symmetric_result = TreeSetColoringResult(A_and_Aᵀ, ag, color, tree_set, R) return BicoloringResult(A, ag, symmetric_result, R) diff --git a/src/optimal.jl b/src/optimal.jl new file mode 100644 index 00000000..58e92c9c --- /dev/null +++ b/src/optimal.jl @@ -0,0 +1,33 @@ +""" + OptimalColoringAlgorithm + +Coloring algorithm that relies on mathematical programming with [JuMP](https://jump.dev/) to find an optimal coloring. + +!!! warning + This algorithm is only available when JuMP is loaded. If you encounter a method error, run `import JuMP` in your REPL and try again. + +!!! danger + The coloring problem is NP-hard, so it is unreasonable to expect an optimal solution in reasonable time for large instances. + +# Constructor + + OptimalColoringAlgorithm(optimizer; silent::Bool=true, assert_solved::Bool=true) + +The `optimizer` argument can be any JuMP-compatible optimizer. +However, the problem formulation is best suited to CP-SAT optimizers like [MiniZinc](https://github.com/jump-dev/MiniZinc.jl). +You can use [`optimizer_with_attributes`](https://jump.dev/JuMP.jl/stable/api/JuMP/#optimizer_with_attributes) to set solver-specific parameters. + +# Keyword arguments + +- `silent`: whether to suppress solver output +- `assert_solved`: whether to check that the solver found an optimal solution (as opposed to running out of time for example) +""" +struct OptimalColoringAlgorithm{O} <: ADTypes.AbstractColoringAlgorithm + optimizer::O + silent::Bool + assert_solved::Bool +end + +function OptimalColoringAlgorithm(optimizer; silent::Bool=true, assert_solved::Bool=true) + return OptimalColoringAlgorithm(optimizer, silent, assert_solved) +end diff --git a/test/Project.toml b/test/Project.toml index ab595aba..217d62c9 100644 --- a/test/Project.toml +++ b/test/Project.toml @@ -12,10 +12,13 @@ CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8" Colors = "5ae59095-9a9b-59fe-a467-6f913c188581" DataFrames = "a93c6f00-e57d-5684-b7b6-d8193f3e46c0" Documenter = "e30172f5-a6a5-5a46-863b-614d45cd2de4" +HiGHS = "87dc4568-4c63-4d18-b0c0-bb2238e4078b" JET = "c3a54625-cd67-489e-a8e7-0a5a0ff4e31b" +JuMP = "4076af6c-e467-56ae-b986-b466b2749572" JuliaFormatter = "98e50ef6-434e-11e9-1051-2b60c6c9e899" LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e" MatrixDepot = "b51810bb-c9f3-55da-ae3c-350fc1fbce05" +MiniZinc = "a7f392d2-6c35-496e-b8cc-0974fbfcbf91" Random = "9a3f8284-a2c9-5f02-9a11-845980a1fd5c" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" SparseMatrixColorings = "0a514795-09f3-496d-8182-132a7b665d35" diff --git a/test/adtypes.jl b/test/adtypes.jl index bbb124a6..da0807d3 100644 --- a/test/adtypes.jl +++ b/test/adtypes.jl @@ -1,26 +1,49 @@ using ADTypes: ADTypes using SparseArrays +using LinearAlgebra using SparseMatrixColorings using Test -@testset "Column coloring" begin - problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) - algo = ADTypes.NoColoringAlgorithm() - A = sprand(10, 20, 0.1) - result = coloring(A, problem, algo) - B = compress(A, result) - @test size(B) == size(A) - @test column_colors(result) == ADTypes.column_coloring(A, algo) - @test decompress(B, result) == A -end +@testset "NoColoringAlgorithm" begin + @testset "Column coloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) + algo = ADTypes.NoColoringAlgorithm() + A = sprand(10, 20, 0.3) + result = coloring(A, problem, algo) + B = compress(A, result) + @test size(B) == size(A) + @test column_colors(result) == ADTypes.column_coloring(A, algo) + @test decompress(B, result) == A + end + + @testset "Row coloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:row) + algo = ADTypes.NoColoringAlgorithm() + A = sprand(10, 20, 0.3) + result = coloring(A, problem, algo) + B = compress(A, result) + @test size(B) == size(A) + @test row_colors(result) == ADTypes.row_coloring(A, algo) + @test decompress(B, result) == A + end + + @testset "Symmetric coloring" begin + problem = ColoringProblem(; structure=:symmetric, partition=:column) + algo = ADTypes.NoColoringAlgorithm() + A = Symmetric(sprand(20, 20, 0.3)) + result = coloring(A, problem, algo) + B = compress(A, result) + @test size(B) == size(A) + @test column_colors(result) == ADTypes.column_coloring(A, algo) + @test decompress(B, result) == A + end -@testset "Row coloring" begin - problem = ColoringProblem(; structure=:nonsymmetric, partition=:row) - algo = ADTypes.NoColoringAlgorithm() - A = sprand(10, 20, 0.1) - result = coloring(A, problem, algo) - B = compress(A, result) - @test size(B) == size(A) - @test row_colors(result) == ADTypes.row_coloring(A, algo) - @test decompress(B, result) == A + @testset "Bicoloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional) + algo = ADTypes.NoColoringAlgorithm() + A = sprand(10, 20, 0.3) + result = coloring(A, problem, algo) + Br, Bc = compress(A, result) + @test decompress(Br, Bc, result) == A + end end diff --git a/test/constant.jl b/test/constant.jl index 5de881b8..97c7f238 100644 --- a/test/constant.jl +++ b/test/constant.jl @@ -1,16 +1,20 @@ using ADTypes: ADTypes using SparseMatrixColorings +using SparseMatrixColorings: InvalidColoringError using Test -matrix_template = ones(100, 200) +matrix_template = ones(Bool, 10, 20) +sym_matrix_template = ones(Bool, 10, 10) @testset "Column coloring" begin problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) - color = rand(1:5, size(matrix_template, 2)) + color = collect(1:20) algo = ConstantColoringAlgorithm(matrix_template, color; partition=:column) - wrong_algo = ConstantColoringAlgorithm(matrix_template, color; partition=:row) + wrong_algo = ConstantColoringAlgorithm{:row}(matrix_template, color) + wrong_color = ConstantColoringAlgorithm{:column}(matrix_template, ones(Int, 20)) @test_throws DimensionMismatch coloring(transpose(matrix_template), problem, algo) @test_throws MethodError coloring(matrix_template, problem, wrong_algo) + @test_throws InvalidColoringError coloring(matrix_template, problem, wrong_color) result = coloring(matrix_template, problem, algo) @test column_colors(result) == color @test ADTypes.column_coloring(matrix_template, algo) == color @@ -19,9 +23,13 @@ end @testset "Row coloring" begin problem = ColoringProblem(; structure=:nonsymmetric, partition=:row) - color = rand(1:5, size(matrix_template, 1)) + color = collect(1:10) algo = ConstantColoringAlgorithm(matrix_template, color; partition=:row) + wrong_algo = ConstantColoringAlgorithm{:column}(matrix_template, color) + wrong_color = ConstantColoringAlgorithm{:row}(matrix_template, ones(Int, 10)) @test_throws DimensionMismatch coloring(transpose(matrix_template), problem, algo) + @test_throws MethodError coloring(matrix_template, problem, wrong_algo) + @test_throws InvalidColoringError coloring(matrix_template, problem, wrong_color) result = coloring(matrix_template, problem, algo) @test row_colors(result) == color @test ADTypes.row_coloring(matrix_template, algo) == color @@ -29,9 +37,22 @@ end end @testset "Symmetric coloring" begin - wrong_problem = ColoringProblem(; structure=:symmetric, partition=:column) - color = rand(1:5, size(matrix_template, 2)) - algo = ConstantColoringAlgorithm(matrix_template, color; partition=:column) - @test_throws MethodError coloring(matrix_template, wrong_problem, algo) - @test_throws MethodError ADTypes.symmetric_coloring(matrix_template, algo) + problem = ColoringProblem(; structure=:symmetric, partition=:column) + color = collect(1:10) + algo = ConstantColoringAlgorithm( + sym_matrix_template, color; partition=:column, structure=:symmetric + ) + wrong_algo = ConstantColoringAlgorithm{:column,:nonsymmetric}( + sym_matrix_template, color + ) + wrong_color = ConstantColoringAlgorithm{:column,:symmetric}( + sym_matrix_template, ones(Int, 20) + ) + @test_throws DimensionMismatch coloring(matrix_template, problem, algo) + @test_throws MethodError coloring(sym_matrix_template, problem, wrong_algo) + @test_throws InvalidColoringError coloring(sym_matrix_template, problem, wrong_color) + result = coloring(sym_matrix_template, problem, algo) + @test column_colors(result) == color + @test ADTypes.symmetric_coloring(sym_matrix_template, algo) == color + @test_throws MethodError ADTypes.column_coloring(sym_matrix_template, algo) end diff --git a/test/optimal.jl b/test/optimal.jl new file mode 100644 index 00000000..10fc32ef --- /dev/null +++ b/test/optimal.jl @@ -0,0 +1,46 @@ +using SparseArrays +using SparseMatrixColorings +using StableRNGs +using Test +using JuMP +using MiniZinc +using HiGHS + +rng = StableRNG(0) + +asymmetric_params = vcat( + [(10, 20, p) for p in (0.0:0.1:0.5)], [(20, 10, p) for p in (0.0:0.1:0.5)] +) + +algo = GreedyColoringAlgorithm() +optalgo = OptimalColoringAlgorithm(() -> MiniZinc.Optimizer{Float64}("highs"); silent=false) + +@testset "Column coloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:column) + for (m, n, p) in asymmetric_params + A = sprand(rng, m, n, p) + result = coloring(A, problem, algo) + optresult = coloring(A, problem, optalgo) + @test ncolors(result) >= ncolors(optresult) + end +end + +@testset "Row coloring" begin + problem = ColoringProblem(; structure=:nonsymmetric, partition=:row) + for (m, n, p) in asymmetric_params + A = sprand(rng, m, n, p) + result = coloring(A, problem, algo) + optresult = coloring(A, problem, optalgo) + @test ncolors(result) >= ncolors(optresult) + end +end + +@testset "Too big" begin + A = sprand(rng, Bool, 100, 100, 0.1) + optalgo_timelimit = OptimalColoringAlgorithm( + optimizer_with_attributes(HiGHS.Optimizer, "time_limit" => 10.0); # 1 second + silent=false, + assert_solved=false, + ) + @test_throws AssertionError coloring(A, ColoringProblem(), optalgo_timelimit) +end diff --git a/test/result.jl b/test/result.jl index 0bea3353..49344aab 100644 --- a/test/result.jl +++ b/test/result.jl @@ -1,3 +1,4 @@ +using SparseMatrixColorings using SparseMatrixColorings: group_by_color, UnsupportedDecompressionError using Test @@ -20,7 +21,7 @@ using Test end @testset "Empty compression" begin - A = rand(10, 10) + A = zeros(Bool, 10, 10) color = zeros(Int, 10) problem = ColoringProblem{:nonsymmetric,:column}() algo = ConstantColoringAlgorithm(A, color; partition=:column) diff --git a/test/runtests.jl b/test/runtests.jl index 81c4ecc9..87b12f57 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -58,6 +58,9 @@ include("utils.jl") @testset "Constant coloring" begin include("constant.jl") end + @testset "Optimal coloring" begin + include("optimal.jl") + end @testset "ADTypes coloring algorithms" begin include("adtypes.jl") end