Allow picking the best of several orders for GreedyColoringAlgorithm

.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
+          fail_ci_if_error: false

docs/src/tutorial.md

@@ -31,7 +31,7 @@ problem = ColoringProblem()
 
 The algorithm defines how you want to solve it. It can be either a [`GreedyColoringAlgorithm`](@ref) or a [`ConstantColoringAlgorithm`](@ref). For `GreedyColoringAlgorithm`, you can select options such as
 
-- the order in which vertices are processed (a subtype of [`AbstractOrder`](@ref SparseMatrixColorings.AbstractOrder))
+- the order in which vertices are processed (a subtype of [`AbstractOrder`](@ref SparseMatrixColorings.AbstractOrder) , or a tuple of such objects)
 - the type of decompression you want (`:direct` or `:substitution`)
 
 ```@example tutorial

src/interface.jl

@@ -72,7 +72,7 @@ It is passed as an argument to the main function [`coloring`](@ref).
     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.
+- `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.
 - `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.
 
@@ -94,26 +94,31 @@ See their respective docstrings for details.
 - [`AbstractOrder`](@ref)
 - [`decompress`](@ref)
 """
-struct GreedyColoringAlgorithm{decompression,O<:AbstractOrder} <:
+struct GreedyColoringAlgorithm{decompression,N,O<:NTuple{N,AbstractOrder}} <:
        ADTypes.AbstractColoringAlgorithm
-    order::O
+    orders::O
     postprocessing::Bool
-end
 
-function GreedyColoringAlgorithm{decompression}(
-    order::AbstractOrder=NaturalOrder(); postprocessing::Bool=false
-) where {decompression}
-    check_valid_algorithm(decompression)
-    return GreedyColoringAlgorithm{decompression,typeof(order)}(order, postprocessing)
+    function GreedyColoringAlgorithm{decompression}(
+        order_or_orders::Union{AbstractOrder,Tuple}=NaturalOrder();
+        postprocessing::Bool=false,
+    ) 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)
+    end
 end
 
 function GreedyColoringAlgorithm(
-    order::AbstractOrder=NaturalOrder();
+    order_or_orders::Union{AbstractOrder,Tuple}=NaturalOrder();
     postprocessing::Bool=false,
     decompression::Symbol=:direct,
 )
-    check_valid_algorithm(decompression)
-    return GreedyColoringAlgorithm{decompression,typeof(order)}(order, postprocessing)
+    return GreedyColoringAlgorithm{decompression}(order_or_orders; postprocessing)
 end
 
 ## Coloring
@@ -229,8 +234,11 @@ function _coloring(
 )
     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_by_order = map(algo.orders) do order
+        vertices_in_order = vertices(bg, Val(2), order)
+        return partial_distance2_coloring(bg, Val(2), vertices_in_order)
+    end
+    color = argmin(maximum, color_by_order)
     if speed_setting isa WithResult
         return ColumnColoringResult(A, bg, color)
     else
@@ -248,8 +256,11 @@ function _coloring(
 )
     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_by_order = map(algo.orders) do order
+        vertices_in_order = vertices(bg, Val(1), order)
+        return partial_distance2_coloring(bg, Val(1), vertices_in_order)
+    end
+    color = argmin(maximum, color_by_order)
     if speed_setting isa WithResult
         return RowColoringResult(A, bg, color)
     else
@@ -266,8 +277,11 @@ function _coloring(
     symmetric_pattern::Bool,
 )
     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_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)
+    end
+    color, star_set = argmin(maximum ∘ first, color_and_star_set_by_order)
     if speed_setting isa WithResult
         return StarSetColoringResult(A, ag, color, star_set)
     else
@@ -284,8 +298,11 @@ function _coloring(
     symmetric_pattern::Bool,
 ) 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_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)
+    end
+    color, tree_set = argmin(maximum ∘ first, color_and_tree_set_by_order)
     if speed_setting isa WithResult
         return TreeSetColoringResult(A, ag, color, tree_set, R)
     else
@@ -303,15 +320,37 @@ function _coloring(
 ) 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)
+    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)
+        (_row_color, _column_color, _symmetric_to_row, _symmetric_to_column) = remap_colors(
+            eltype(ag), _color, maximum(_color), size(A)...
+        )
+        return (
+            _color,
+            _star_set,
+            _row_color,
+            _column_color,
+            _symmetric_to_row,
+            _symmetric_to_column,
+        )
+    end
+    (color, star_set, row_color, column_color, symmetric_to_row, symmetric_to_column) = argmin(
+        t -> maximum(t[3]) + maximum(t[4]), outputs_by_order
+    )  # can't use ncolors without computing the full result
     if speed_setting isa WithResult
         symmetric_result = StarSetColoringResult(A_and_Aᵀ, ag, color, star_set)
-        return BicoloringResult(A, ag, symmetric_result, R)
-    else
-        row_color, column_color, _ = remap_colors(
-            eltype(ag), color, maximum(color), size(A)...
+        return BicoloringResult(
+            A,
+            ag,
+            symmetric_result,
+            row_color,
+            column_color,
+            symmetric_to_row,
+            symmetric_to_column,
+            R,
         )
+    else
         return row_color, column_color
     end
 end
@@ -326,15 +365,37 @@ function _coloring(
 ) 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)
+    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)
+        (_row_color, _column_color, _symmetric_to_row, _symmetric_to_column) = remap_colors(
+            eltype(ag), _color, maximum(_color), size(A)...
+        )
+        return (
+            _color,
+            _tree_set,
+            _row_color,
+            _column_color,
+            _symmetric_to_row,
+            _symmetric_to_column,
+        )
+    end
+    (color, tree_set, row_color, column_color, symmetric_to_row, symmetric_to_column) = argmin(
+        t -> maximum(t[3]) + maximum(t[4]), outputs_by_order
+    )  # can't use ncolors without computing the full result
     if speed_setting isa WithResult
         symmetric_result = TreeSetColoringResult(A_and_Aᵀ, ag, color, tree_set, R)
-        return BicoloringResult(A, ag, symmetric_result, R)
-    else
-        row_color, column_color, _ = remap_colors(
-            eltype(ag), color, maximum(color), size(A)...
+        return BicoloringResult(
+            A,
+            ag,
+            symmetric_result,
+            row_color,
+            column_color,
+            symmetric_to_row,
+            symmetric_to_column,
+            R,
         )
+    else
         return row_color, column_color
     end
 end

src/result.jl

@@ -686,14 +686,15 @@ function BicoloringResult(
     A::AbstractMatrix,
     ag::AdjacencyGraph{T},
     symmetric_result::AbstractColoringResult{:symmetric,:column},
+    row_color::Vector{T},
+    column_color::Vector{T},
+    symmetric_to_row::Vector{T},
+    symmetric_to_column::Vector{T},
     decompression_eltype::Type{R},
 ) where {T,R}
     m, n = size(A)
     symmetric_color = column_colors(symmetric_result)
     num_sym_colors = maximum(symmetric_color)
-    row_color, column_color, symmetric_to_row, symmetric_to_column = remap_colors(
-        T, symmetric_color, num_sym_colors, m, n
-    )
     column_group = group_by_color(T, column_color)
     row_group = group_by_color(T, row_color)
     Br_and_Bc = Matrix{R}(undef, n + m, num_sym_colors)

test/order.jl

@@ -146,3 +146,99 @@ end;
         @test isperm(π)
     end
 end
+
+@testset "Multiple orders" begin
+    # I used brute force to find examples where LargestFirst is *strictly* better than NaturalOrder, just to check that the best order is indeed selected when multiple orders are provided
+    @testset "Column coloring" begin
+        A = [
+            0 0 1 1
+            0 1 0 1
+            0 0 1 1
+            1 1 0 0
+        ]
+        problem = ColoringProblem{:nonsymmetric,:column}()
+        algo = GreedyColoringAlgorithm(NaturalOrder())
+        better_algo = GreedyColoringAlgorithm((NaturalOrder(), LargestFirst()))
+        @test ncolors(coloring(A, problem, better_algo)) <
+            ncolors(coloring(A, problem, algo))
+    end
+    @testset "Row coloring" begin
+        A = [
+            1 0 0 0
+            0 0 1 0
+            0 1 1 1
+            1 0 0 1
+        ]
+        problem = ColoringProblem{:nonsymmetric,:row}()
+        algo = GreedyColoringAlgorithm(NaturalOrder())
+        better_algo = GreedyColoringAlgorithm((NaturalOrder(), LargestFirst()))
+        @test ncolors(coloring(A, problem, better_algo)) <
+            ncolors(coloring(A, problem, algo))
+    end
+    @testset "Star coloring" begin
+        A = [
+            0 1 0 1 1
+            1 1 0 1 0
+            0 0 1 0 1
+            1 1 0 1 0
+            1 0 1 0 0
+        ]
+        problem = ColoringProblem{:symmetric,:column}()
+        algo = GreedyColoringAlgorithm(NaturalOrder())
+        better_algo = GreedyColoringAlgorithm((NaturalOrder(), LargestFirst()))
+        @test ncolors(coloring(A, problem, better_algo)) <
+            ncolors(coloring(A, problem, algo))
+    end
+    @testset "Acyclic coloring" begin
+        A = [
+            1 0 0 0 0 1 0
+            0 0 0 1 0 0 0
+            0 0 0 1 0 0 0
+            0 1 1 1 0 1 1
+            0 0 0 0 0 0 1
+            1 0 0 1 0 0 1
+            0 0 0 1 1 1 1
+        ]
+        problem = ColoringProblem{:symmetric,:column}()
+        algo = GreedyColoringAlgorithm{:substitution}(NaturalOrder())
+        better_algo = GreedyColoringAlgorithm{:substitution}((
+            NaturalOrder(), LargestFirst()
+        ))
+        @test ncolors(coloring(A, problem, better_algo)) <
+            ncolors(coloring(A, problem, algo))
+    end
+    @testset "Star bicoloring" begin
+        A = [
+            0 1 0 0 0
+            1 0 1 0 0
+            0 1 0 0 1
+            0 0 0 0 0
+            0 0 1 0 1
+        ]
+        problem = ColoringProblem{:nonsymmetric,:bidirectional}()
+        algo = GreedyColoringAlgorithm(NaturalOrder())
+        better_algo = GreedyColoringAlgorithm((NaturalOrder(), LargestFirst()))
+        @test ncolors(coloring(A, problem, better_algo)) <
+            ncolors(coloring(A, problem, algo))
+    end
+    @testset "Acyclic bicoloring" begin
+        A = [
+            0 1 0 1 1 0 1 0 1
+            1 0 0 0 0 0 0 0 1
+            0 0 0 0 0 0 0 0 0
+            1 0 0 1 1 0 1 0 0
+            1 0 0 1 0 0 0 0 0
+            0 0 0 0 0 0 0 0 0
+            1 0 0 1 0 0 0 0 0
+            0 0 0 0 0 0 0 0 0
+            1 1 0 0 0 0 0 0 0
+        ]
+        problem = ColoringProblem{:nonsymmetric,:bidirectional}()
+        algo = GreedyColoringAlgorithm{:substitution}(NaturalOrder())
+        better_algo = GreedyColoringAlgorithm{:substitution}((
+            NaturalOrder(), LargestFirst()
+        ))
+        @test ncolors(coloring(A, problem, better_algo)) <
+            ncolors(coloring(A, problem, algo))
+    end
+end

test/random.jl

@@ -81,10 +81,10 @@ end;
     problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional)
     @testset for algo in (
         GreedyColoringAlgorithm(
-            RandomOrder(rng); postprocessing=false, decompression=:direct
+            RandomOrder(StableRNG(0), 0); postprocessing=false, decompression=:direct
         ),
         GreedyColoringAlgorithm(
-            RandomOrder(rng); postprocessing=true, decompression=:direct
+            RandomOrder(StableRNG(0), 0); postprocessing=true, decompression=:direct
         ),
     )
         @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params
@@ -102,10 +102,10 @@ end;
     problem = ColoringProblem(; structure=:nonsymmetric, partition=:bidirectional)
     @testset for algo in (
         GreedyColoringAlgorithm(
-            RandomOrder(rng); postprocessing=false, decompression=:substitution
+            RandomOrder(StableRNG(0), 0); postprocessing=false, decompression=:substitution
         ),
         GreedyColoringAlgorithm(
-            RandomOrder(rng); postprocessing=true, decompression=:substitution
+            RandomOrder(StableRNG(0), 0); postprocessing=true, decompression=:substitution
         ),
     )
         @testset "$((; m, n, p))" for (m, n, p) in asymmetric_params

test/type_stability.jl

@@ -40,11 +40,21 @@ rng = StableRNG(63)
                 ColoringProblem(; structure, partition),
                 GreedyColoringAlgorithm(order; decompression),
             )
+            @test_opt coloring(
+                A,
+                ColoringProblem(; structure, partition),
+                GreedyColoringAlgorithm((NaturalOrder(), order); decompression),
+            )
             @inferred coloring(
                 A,
                 ColoringProblem(; structure, partition),
                 GreedyColoringAlgorithm(order; decompression),
             )
+            @inferred coloring(
+                A,
+                ColoringProblem(; structure, partition),
+                GreedyColoringAlgorithm((NaturalOrder(), order); decompression),
+            )
         end
     end
 end;