From bedcd88344a60678ab1099bb9a8d414a9b606632 Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Fri, 10 Oct 2025 18:52:45 +0200
Subject: [PATCH 01/10] Allow picking the best of several orders for
 GreedyColoringAlgorithm

---
 src/interface.jl | 73 +++++++++++++++++++++++++++++++-----------------
 test/utils.jl    | 26 +++++++++++++++++
 2 files changed, 74 insertions(+), 25 deletions(-)

diff --git a/src/interface.jl b/src/interface.jl
index 79fe71e9..435da353 100644
--- a/src/interface.jl
+++ b/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,NTuple}`: 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,NTuple{N,AbstractOrder}}=NaturalOrder();
+        postprocessing::Bool=false,
+    ) where {decompression,N}
+        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,NTuple{N,AbstractOrder}}=NaturalOrder();
     postprocessing::Bool=false,
     decompression::Symbol=:direct,
-)
-    check_valid_algorithm(decompression)
-    return GreedyColoringAlgorithm{decompression,typeof(order)}(order, postprocessing)
+) where {N}
+    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,8 +320,11 @@ 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)
+    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
         symmetric_result = StarSetColoringResult(A_and_Aᵀ, ag, color, star_set)
         return BicoloringResult(A, ag, symmetric_result, R)
@@ -326,8 +346,11 @@ 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)
+    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
         symmetric_result = TreeSetColoringResult(A_and_Aᵀ, ag, color, tree_set, R)
         return BicoloringResult(A, ag, symmetric_result, R)
diff --git a/test/utils.jl b/test/utils.jl
index bb80f95f..32f94c4d 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -15,6 +15,10 @@ using SparseMatrixColorings:
     structurally_biorthogonal
 using Test
 
+const _ALL_ORDERS = (
+    NaturalOrder(), LargestFirst(), SmallestLast(), IncidenceDegree(), DynamicLargestFirst()
+)
+
 function test_coloring_decompression(
     A0::AbstractMatrix,
     problem::ColoringProblem{structure,partition},
@@ -169,6 +173,17 @@ function test_coloring_decompression(
             @show color_vec
         end
     end
+
+    @testset "More orders is better" begin
+        if algo.orders ⊆ _ALL_ORDERS
+            better_algo = GreedyColoringAlgorithm{decompression}(
+                _ALL_ORDERS; algo.postprocessing
+            )
+            result = coloring(A0, problem, algo)
+            better_result = coloring(A0, problem, better_algo)
+            @test ncolors(result) >= ncolors(better_result)
+        end
+    end
 end
 
 function test_bicoloring_decompression(
@@ -216,6 +231,17 @@ function test_bicoloring_decompression(
             end
         end
     end
+
+    @testset "More orders is better" begin
+        if algo.orders ⊆ _ALL_ORDERS
+            better_algo = GreedyColoringAlgorithm{decompression}(
+                _ALL_ORDERS; algo.postprocessing
+            )
+            result = coloring(A0, problem, algo)
+            better_result = coloring(A0, problem, better_algo)
+            @test ncolors(result) >= ncolors(better_result)
+        end
+    end
 end
 
 function test_structured_coloring_decompression(A::AbstractMatrix)

From 9d7403215dcb8723a579c5498c70291c85f9dc7f Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Fri, 10 Oct 2025 19:11:17 +0200
Subject: [PATCH 02/10] Fix type params

---
 src/interface.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/src/interface.jl b/src/interface.jl
index 435da353..2d749595 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -100,9 +100,9 @@ struct GreedyColoringAlgorithm{decompression,N,O<:NTuple{N,AbstractOrder}} <:
     postprocessing::Bool
 
     function GreedyColoringAlgorithm{decompression}(
-        order_or_orders::Union{AbstractOrder,NTuple{N,AbstractOrder}}=NaturalOrder();
+        order_or_orders::Union{AbstractOrder,NTuple}=NaturalOrder();
         postprocessing::Bool=false,
-    ) where {decompression,N}
+    ) where {decompression}
         check_valid_algorithm(decompression)
         if order_or_orders isa AbstractOrder
             orders = (order_or_orders,)
@@ -114,10 +114,10 @@ struct GreedyColoringAlgorithm{decompression,N,O<:NTuple{N,AbstractOrder}} <:
 end
 
 function GreedyColoringAlgorithm(
-    order_or_orders::Union{AbstractOrder,NTuple{N,AbstractOrder}}=NaturalOrder();
+    order_or_orders::Union{AbstractOrder,NTuple}=NaturalOrder();
     postprocessing::Bool=false,
     decompression::Symbol=:direct,
-) where {N}
+)
     return GreedyColoringAlgorithm{decompression}(order_or_orders; postprocessing)
 end
 

From c104a5344250e246509ea00c80f256dc46cb438b Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Fri, 10 Oct 2025 19:20:22 +0200
Subject: [PATCH 03/10] Tuple

---
 src/interface.jl | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/src/interface.jl b/src/interface.jl
index 2d749595..44700f5b 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -100,7 +100,7 @@ struct GreedyColoringAlgorithm{decompression,N,O<:NTuple{N,AbstractOrder}} <:
     postprocessing::Bool
 
     function GreedyColoringAlgorithm{decompression}(
-        order_or_orders::Union{AbstractOrder,NTuple}=NaturalOrder();
+        order_or_orders::Union{AbstractOrder,Tuple}=NaturalOrder();
         postprocessing::Bool=false,
     ) where {decompression}
         check_valid_algorithm(decompression)
@@ -114,7 +114,7 @@ struct GreedyColoringAlgorithm{decompression,N,O<:NTuple{N,AbstractOrder}} <:
 end
 
 function GreedyColoringAlgorithm(
-    order_or_orders::Union{AbstractOrder,NTuple}=NaturalOrder();
+    order_or_orders::Union{AbstractOrder,Tuple}=NaturalOrder();
     postprocessing::Bool=false,
     decompression::Symbol=:direct,
 )

From 77b920aad8f338cea00ea6141f7e4be1fa62ad65 Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Mon, 13 Oct 2025 10:34:02 +0200
Subject: [PATCH 04/10] Better tests

---
 docs/src/tutorial.md   |  2 +-
 test/type_stability.jl | 10 ++++++++++
 test/utils.jl          | 42 ++++++++++++++++++++++++++----------------
 3 files changed, 37 insertions(+), 17 deletions(-)

diff --git a/docs/src/tutorial.md b/docs/src/tutorial.md
index fec54f13..6e803511 100644
--- a/docs/src/tutorial.md
+++ b/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
diff --git a/test/type_stability.jl b/test/type_stability.jl
index 5a488596..3c6f9f02 100644
--- a/test/type_stability.jl
+++ b/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;
diff --git a/test/utils.jl b/test/utils.jl
index 32f94c4d..2462f8c1 100644
--- a/test/utils.jl
+++ b/test/utils.jl
@@ -175,14 +175,19 @@ function test_coloring_decompression(
     end
 
     @testset "More orders is better" begin
-        if algo.orders ⊆ _ALL_ORDERS
-            better_algo = GreedyColoringAlgorithm{decompression}(
-                _ALL_ORDERS; algo.postprocessing
-            )
-            result = coloring(A0, problem, algo)
-            better_result = coloring(A0, problem, better_algo)
-            @test ncolors(result) >= ncolors(better_result)
-        end
+        more_orders = (algo.orders..., _ALL_ORDERS...)
+        better_algo = GreedyColoringAlgorithm{decompression}(
+            more_orders; algo.postprocessing
+        )
+        all_algos = [
+            GreedyColoringAlgorithm{decompression}(order; algo.postprocessing) for
+            order in more_orders
+        ]
+        result = coloring(A0, problem, algo)
+        better_result = coloring(A0, problem, better_algo)
+        all_results = [coloring(A0, problem, _algo) for _algo in all_algos]
+        @test ncolors(better_result) <= ncolors(result)
+        @test ncolors(better_result) == minimum(ncolors, all_results)
     end
 end
 
@@ -233,14 +238,19 @@ function test_bicoloring_decompression(
     end
 
     @testset "More orders is better" begin
-        if algo.orders ⊆ _ALL_ORDERS
-            better_algo = GreedyColoringAlgorithm{decompression}(
-                _ALL_ORDERS; algo.postprocessing
-            )
-            result = coloring(A0, problem, algo)
-            better_result = coloring(A0, problem, better_algo)
-            @test ncolors(result) >= ncolors(better_result)
-        end
+        more_orders = (algo.orders..., _ALL_ORDERS...)
+        better_algo = GreedyColoringAlgorithm{decompression}(
+            more_orders; algo.postprocessing
+        )
+        all_algos = [
+            GreedyColoringAlgorithm{decompression}(order; algo.postprocessing) for
+            order in more_orders
+        ]
+        result = coloring(A0, problem, algo)
+        better_result = coloring(A0, problem, better_algo)
+        all_results = [coloring(A0, problem, _algo) for _algo in all_algos]
+        @test ncolors(better_result) <= ncolors(result)
+        @test ncolors(better_result) == minimum(ncolors, all_results)
     end
 end
 

From 7a5d6fe31fa47eaaba5c72964865624cbeba0482 Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Mon, 13 Oct 2025 11:42:27 +0200
Subject: [PATCH 05/10] Better test

---
 src/interface.jl | 30 +++++++++------
 test/order.jl    | 96 ++++++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 114 insertions(+), 12 deletions(-)

diff --git a/src/interface.jl b/src/interface.jl
index 44700f5b..76874237 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -320,18 +320,21 @@ 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)
-    color_and_star_set_by_order = map(algo.orders) do order
+    outputs_by_order = map(algo.orders) do order
         vertices_in_order = vertices(ag, order)
-        return star_coloring(ag, vertices_in_order, algo.postprocessing)
+        color, star_set = star_coloring(ag, vertices_in_order, algo.postprocessing)
+        row_color, column_color, _ = remap_colors(
+            eltype(ag), color, maximum(color), size(A)...
+        )
+        return (; color, star_set, row_color, column_color)
     end
-    color, star_set = argmin(maximum ∘ first, color_and_star_set_by_order)
+    (; color, star_set, row_color, column_color) = argmin(
+        t -> maximum(t.row_color) + maximum(t.column_color), 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 row_color, column_color
     end
 end
@@ -346,18 +349,21 @@ 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)
-    color_and_tree_set_by_order = map(algo.orders) do order
+    outputs_by_order = map(algo.orders) do order
         vertices_in_order = vertices(ag, order)
-        return acyclic_coloring(ag, vertices_in_order, algo.postprocessing)
+        color, tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing)
+        row_color, column_color, _ = remap_colors(
+            eltype(ag), color, maximum(color), size(A)...
+        )
+        return (; color, tree_set, row_color, column_color)
     end
-    color, tree_set = argmin(maximum ∘ first, color_and_tree_set_by_order)
+    (; color, tree_set, row_color, column_color) = argmin(
+        t -> maximum(t.row_color) + maximum(t.column_color), 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 row_color, column_color
     end
 end
diff --git a/test/order.jl b/test/order.jl
index 9bb9e6dc..838c32a9 100644
--- a/test/order.jl
+++ b/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

From 54622dd0a291f623ad4142fb47930625be2a7a50 Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Mon, 13 Oct 2025 11:45:17 +0200
Subject: [PATCH 06/10] Fix foc

---
 src/interface.jl | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/src/interface.jl b/src/interface.jl
index 76874237..47b00d14 100644
--- a/src/interface.jl
+++ b/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::Union{AbstractOrder,NTuple}`: 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.
+- `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.
 

From be5a06a5d977b567721e25ae40c740aa04a02369 Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Mon, 13 Oct 2025 12:03:45 +0200
Subject: [PATCH 07/10] Fix type inference inside closure

---
 src/interface.jl | 24 ++++++++++++------------
 1 file changed, 12 insertions(+), 12 deletions(-)

diff --git a/src/interface.jl b/src/interface.jl
index 47b00d14..ca8fee1f 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -322,14 +322,14 @@ function _coloring(
     ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; has_diagonal=false)
     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, _ = remap_colors(
-            eltype(ag), color, maximum(color), size(A)...
+        _color, _star_set = star_coloring(ag, vertices_in_order, algo.postprocessing)
+        _row_color, _column_color, _ = remap_colors(
+            eltype(ag), _color, maximum(_color), size(A)...
         )
-        return (; color, star_set, row_color, column_color)
+        return (_color, _star_set, _row_color, _column_color)
     end
-    (; color, star_set, row_color, column_color) = argmin(
-        t -> maximum(t.row_color) + maximum(t.column_color), outputs_by_order
+    (color, star_set, row_color, column_color) = 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)
@@ -351,14 +351,14 @@ function _coloring(
     ag = AdjacencyGraph(A_and_Aᵀ, edge_to_index; has_diagonal=false)
     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, _ = remap_colors(
-            eltype(ag), color, maximum(color), size(A)...
+        _color, _tree_set = acyclic_coloring(ag, vertices_in_order, algo.postprocessing)
+        _row_color, _column_color, _ = remap_colors(
+            eltype(ag), _color, maximum(_color), size(A)...
         )
-        return (; color, tree_set, row_color, column_color)
+        return (; _color, _tree_set, _row_color, _column_color)
     end
-    (; color, tree_set, row_color, column_color) = argmin(
-        t -> maximum(t.row_color) + maximum(t.column_color), outputs_by_order
+    (color, tree_set, row_color, column_color) = 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)

From 8e132dc29cd2597390a4ebc21337fd3b905bf0b3 Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Mon, 13 Oct 2025 13:50:52 +0200
Subject: [PATCH 08/10] Fix seed

---
 test/random.jl | 8 ++++----
 1 file changed, 4 insertions(+), 4 deletions(-)

diff --git a/test/random.jl b/test/random.jl
index 02c1a3a1..406dfea9 100644
--- a/test/random.jl
+++ b/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

From ace3de0913a52f1ca463bd3b7478d5c3614cd21f Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Mon, 13 Oct 2025 17:31:34 +0200
Subject: [PATCH 09/10] Test on 1.11

---
 .github/workflows/Test.yml | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

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

From a1b5b5eea18a3f3e7b8801c8f797509f1f1cf63f Mon Sep 17 00:00:00 2001
From: Guillaume Dalle <22795598+gdalle@users.noreply.github.com>
Date: Tue, 14 Oct 2025 07:01:51 +0200
Subject: [PATCH 10/10] Avoid duplicate remap_colors

---
 src/interface.jl | 48 ++++++++++++++++++++++++++++++++++++++++--------
 src/result.jl    |  7 ++++---
 2 files changed, 44 insertions(+), 11 deletions(-)

diff --git a/src/interface.jl b/src/interface.jl
index ca8fee1f..c1a56c29 100644
--- a/src/interface.jl
+++ b/src/interface.jl
@@ -323,17 +323,33 @@ function _coloring(
     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, _ = remap_colors(
+        (_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)
+        return (
+            _color,
+            _star_set,
+            _row_color,
+            _column_color,
+            _symmetric_to_row,
+            _symmetric_to_column,
+        )
     end
-    (color, star_set, row_color, column_color) = argmin(
+    (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)
+        return BicoloringResult(
+            A,
+            ag,
+            symmetric_result,
+            row_color,
+            column_color,
+            symmetric_to_row,
+            symmetric_to_column,
+            R,
+        )
     else
         return row_color, column_color
     end
@@ -352,17 +368,33 @@ function _coloring(
     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, _ = remap_colors(
+        (_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)
+        return (
+            _color,
+            _tree_set,
+            _row_color,
+            _column_color,
+            _symmetric_to_row,
+            _symmetric_to_column,
+        )
     end
-    (color, tree_set, row_color, column_color) = argmin(
+    (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)
+        return BicoloringResult(
+            A,
+            ag,
+            symmetric_result,
+            row_color,
+            column_color,
+            symmetric_to_row,
+            symmetric_to_column,
+            R,
+        )
     else
         return row_color, column_color
     end
diff --git a/src/result.jl b/src/result.jl
index b5b9cdd2..7a3e3b68 100644
--- a/src/result.jl
+++ b/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)