Allow exact reproduction of ColPack's dynamic orders

src/order.jl

@@ -105,7 +105,7 @@ function vertices(bg::BipartiteGraph{T}, ::Val{side}, ::LargestFirst) where {T,s
 end
 
 """
-    DynamicDegreeBasedOrder{degtype,direction}
+    DynamicDegreeBasedOrder{degtype,direction}(; reproduce_colpack=false)
 
 Instance of [`AbstractOrder`](@ref) which sorts vertices using a dynamically computed degree.
 
@@ -117,6 +117,10 @@ Instance of [`AbstractOrder`](@ref) which sorts vertices using a dynamically com
 - `degtype::Symbol`: can be `:forward` (for the forward degree) or `:back` (for the back degree)
 - `direction::Symbol`: can be `:low2high` (if the order is defined from lowest to highest, i.e. `1` to `n`) or `:high2low` (if the order is defined from highest to lowest, i.e. `n` to `1`)
 
+# Settings
+
+- `reproduce_colpack::Bool`: whether to manage the buckets in the same way as the original ColPack implementation. When `reproduce_colpack=true`, we always append and remove vertices from the end of a bucket, which incurs a large performance penalty because every modification requires a circular permutation of the corresponding bucket. This setting is mostly for the purpose of reproducing past research results which rely on implementation details.
+
 # Concrete variants
 
 - [`IncidenceDegree`](@ref)
@@ -127,17 +131,28 @@ Instance of [`AbstractOrder`](@ref) which sorts vertices using a dynamically com
 
 - [_ColPack: Software for graph coloring and related problems in scientific computing_](https://dl.acm.org/doi/10.1145/2513109.2513110), Gebremedhin et al. (2013), Section 5
 """
-struct DynamicDegreeBasedOrder{degtype,direction} <: AbstractOrder end
+struct DynamicDegreeBasedOrder{degtype,direction} <: AbstractOrder
+    reproduce_colpack::Bool
+end
+
+function DynamicDegreeBasedOrder{degtype,direction}(;
+    reproduce_colpack::Bool=false
+) where {degtype,direction}
+    return DynamicDegreeBasedOrder{degtype,direction}(reproduce_colpack)
+end
 
 struct DegreeBuckets{T}
     degrees::Vector{T}
     bucket_storage::Vector{T}
     bucket_low::Vector{T}
     bucket_high::Vector{T}
     positions::Vector{T}
+    reproduce_colpack::Bool
 end
 
-function DegreeBuckets(::Type{T}, degrees::Vector{T}, dmax::Integer) where {T}
+function DegreeBuckets(
+    ::Type{T}, degrees::Vector{T}, dmax::Integer; reproduce_colpack::Bool
+) where {T}
     # number of vertices per degree class
     deg_count = zeros(T, dmax + 1)
     for d in degrees
@@ -157,7 +172,9 @@ function DegreeBuckets(::Type{T}, degrees::Vector{T}, dmax::Integer) where {T}
         bucket_storage[positions[v]] = v
         deg_count[d + 1] -= 1
     end
-    return DegreeBuckets(degrees, bucket_storage, bucket_low, bucket_high, positions)
+    return DegreeBuckets(
+        degrees, bucket_storage, bucket_low, bucket_high, positions, reproduce_colpack
+    )
 end
 
 maxdeg(db::DegreeBuckets) = length(db.bucket_low) - 1
@@ -205,8 +222,42 @@ function pop_next_candidate!(db::DegreeBuckets; direction::Symbol)
     return candidate
 end
 
+function rotate_bucket_left!(db::DegreeBuckets, d::Integer)
+    (; bucket_storage, bucket_high, bucket_low, positions) = db
+    low, high = bucket_low[d + 1], bucket_high[d + 1]
+    # remember first element v
+    v = bucket_storage[low]
+    # shift everyone else one index down
+    for i in (low + 1):high
+        w = bucket_storage[i]
+        bucket_storage[i - 1] = w
+        positions[w] = i - 1
+    end
+    # put v back at the end
+    bucket_storage[high] = v
+    positions[v] = high
+    return nothing
+end
+
+function rotate_bucket_right!(db::DegreeBuckets, d::Integer)
+    (; bucket_storage, bucket_high, bucket_low, positions) = db
+    low, high = bucket_low[d + 1], bucket_high[d + 1]
+    # remember last element v
+    v = bucket_storage[high]
+    # shift everyone else one index up
+    for i in (high - 1):-1:low
+        w = bucket_storage[i]
+        bucket_storage[i + 1] = w
+        positions[w] = i + 1
+    end
+    # put v back at the start
+    bucket_storage[low] = v
+    positions[v] = low
+    return nothing
+end
+
 function update_bucket!(db::DegreeBuckets, v::Integer; degtype, direction)
-    (; degrees, bucket_storage, bucket_low, bucket_high, positions) = db
+    (; degrees, bucket_storage, bucket_low, bucket_high, positions, reproduce_colpack) = db
     d, p = degrees[v], positions[v]
     low, high = bucket_low[d + 1], bucket_high[d + 1]
     # select previous or next bucket for the move
@@ -227,11 +278,27 @@ function update_bucket!(db::DegreeBuckets, v::Integer; degtype, direction)
         # update v's stats
         degrees[v] = d_new
         positions[v] = low_new - 1
+        if reproduce_colpack
+            # move v from start to end of the next bucket, preserving order
+            rotate_bucket_left!(db, d_new)  # expensive
+        end
     else
-        # move the vertex w located at the start of the current bucket to v's position (!= ColPack)
-        w = bucket_storage[low]
-        bucket_storage[p] = w
-        positions[w] = p
+        if reproduce_colpack
+            # move the vertex w located at the end of the current bucket to v's position
+            w = bucket_storage[high]
+            bucket_storage[p] = w
+            positions[w] = p
+            # explicitly put v at the end
+            bucket_storage[high] = v
+            positions[v] = high
+            # move v from end to start of the current bucket, preserving order
+            rotate_bucket_right!(db, d)  # expensive
+        else
+            # move the vertex w located at the start of the current bucket to v's position (!= ColPack)
+            w = bucket_storage[low]
+            bucket_storage[p] = w
+            positions[w] = p
+        end
         # shrink current bucket from the left
         # morally we put v at the start and then ignore it
         bucket_low[d + 1] += 1
@@ -249,14 +316,16 @@ function update_bucket!(db::DegreeBuckets, v::Integer; degtype, direction)
 end
 
 function vertices(
-    g::AdjacencyGraph{T}, ::DynamicDegreeBasedOrder{degtype,direction}
+    g::AdjacencyGraph{T}, order::DynamicDegreeBasedOrder{degtype,direction}
 ) where {T<:Integer,degtype,direction}
     if degree_increasing(; degtype, direction)
         degrees = zeros(T, nb_vertices(g))
     else
         degrees = T[degree(g, v) for v in vertices(g)]
     end
-    db = DegreeBuckets(T, degrees, maximum_degree(g))
+    db = DegreeBuckets(
+        T, degrees, maximum_degree(g); reproduce_colpack=order.reproduce_colpack
+    )
     π = T[]
     sizehint!(π, nb_vertices(g))
     for _ in 1:nb_vertices(g)
@@ -272,7 +341,7 @@ function vertices(
 end
 
 function vertices(
-    g::BipartiteGraph{T}, ::Val{side}, ::DynamicDegreeBasedOrder{degtype,direction}
+    g::BipartiteGraph{T}, ::Val{side}, order::DynamicDegreeBasedOrder{degtype,direction}
 ) where {T<:Integer,side,degtype,direction}
     other_side = 3 - side
     # compute dist-2 degrees in an optimized way
@@ -294,7 +363,7 @@ function vertices(
         degrees = degrees_dist2
     end
     maxd2 = maximum(degrees_dist2)
-    db = DegreeBuckets(T, degrees, maxd2)
+    db = DegreeBuckets(T, degrees, maxd2; reproduce_colpack=order.reproduce_colpack)
     π = T[]
     sizehint!(π, n)
     visited = falses(n)
@@ -318,7 +387,7 @@ function vertices(
 end
 
 """
-    IncidenceDegree()
+    IncidenceDegree(; reproduce_colpack=false)
 
 Instance of [`AbstractOrder`](@ref) which sorts vertices from lowest to highest using the dynamic back degree.
 
@@ -332,7 +401,7 @@ Instance of [`AbstractOrder`](@ref) which sorts vertices from lowest to highest
 const IncidenceDegree = DynamicDegreeBasedOrder{:back,:low2high}
 
 """
-    SmallestLast()
+    SmallestLast(; reproduce_colpack=false)
 
 Instance of [`AbstractOrder`](@ref) which sorts vertices from highest to lowest using the dynamic back degree.
 
@@ -346,7 +415,7 @@ Instance of [`AbstractOrder`](@ref) which sorts vertices from highest to lowest
 const SmallestLast = DynamicDegreeBasedOrder{:back,:high2low}
 
 """
-    DynamicLargestFirst()
+    DynamicLargestFirst(; reproduce_colpack=false)
 
 Instance of [`AbstractOrder`](@ref) which sorts vertices from lowest to highest using the dynamic forward degree.
 

test/order.jl

@@ -89,8 +89,14 @@ end;
 end;
 
 @testset "Dynamic degree-based orders" begin
-    @testset "$order" for order in
-                          [SmallestLast(), IncidenceDegree(), DynamicLargestFirst()]
+    @testset "$order" for order in [
+        SmallestLast(),
+        SmallestLast(; reproduce_colpack=true),
+        IncidenceDegree(),
+        IncidenceDegree(; reproduce_colpack=true),
+        DynamicLargestFirst(),
+        DynamicLargestFirst(; reproduce_colpack=true),
+    ]
         @testset "AdjacencyGraph" begin
             for (n, p) in Iterators.product(20:20:100, 0.0:0.1:0.2)
                 yield()

test/reference/colpack_table_6_7.csv

@@ -1,11 +1,11 @@
-row,modified,group,name,V1,V2,E,Δ1,Δ2,LF2,N2,LF1,N1
-1,false,LPnetlib,lp_cre_a,3516,7248,18168,360,14,360,360,14,16
-2,false,LPnetlib,lp_ken_11,14694,21349,49058,122,3,128,130,4,5
-3,false,LPnetlib,lp_ken_13,28632,42659,97246,170,3,174,176,5,4
-4,false,LPnetlib,lp_maros_r7,3136,9408,144848,48,46,70,74,100,72
-5,false,LPnetlib,lp_cre_d,8926,73948,246614,808,13,808,813,13,15
-6,false,LPnetlib,lp_ken_18,105127,154699,358171,325,3,325,330,5,5
-7,false,Bai,af23560,23560,23560,460598,21,21,59,44,43,43
-8,false,Shen,e40r0100,17281,17281,553562,62,62,85,95,85,87
-9,false,vanHeukelum,cage11,39082,39082,559722,31,31,70,81,67,81
-10,false,vanHeukelum,cage12,130228,130228,2032536,33,33,79,96,73,96
+row,modified,group,name,V1,V2,E,Δ1,Δ2,SL2,ID2,LF2,DLF2,N2,SL1,ID1,LF1,DLF1,N1
+1,false,LPnetlib,lp_cre_a,3516,7248,18168,360,14,360,360,360,360,360,14,14,14,14,16
+2,false,LPnetlib,lp_ken_11,14694,21349,49058,122,3,125,124,128,122,130,4,4,4,5,5
+3,false,LPnetlib,lp_ken_13,28632,42659,97246,170,3,171,171,174,170,176,4,5,5,5,4
+4,false,LPnetlib,lp_maros_r7,3136,9408,144848,48,46,83,90,70,114,74,80,88,100,113,72
+5,false,LPnetlib,lp_cre_d,8926,73948,246614,808,13,808,808,808,808,813,15,15,13,14,15
+6,false,LPnetlib,lp_ken_18,105127,154699,358171,325,3,325,326,328,325,330,4,5,5,5,5
+7,false,Bai,af23560,23560,23560,460598,21,21,42,42,43,59,44,42,43,43,60,43
+8,false,Shen,e40r0100,17281,17281,553562,62,62,70,71,87,85,95,70,71,85,86,87
+9,false,vanHeukelum,cage11,39082,39082,559722,31,31,64,67,67,70,81,64,67,67,70,81
+10,false,vanHeukelum,cage12,130228,130228,2032536,33,33,67,72,73,79,96,67,72,73,79,96

test/suitesparse.jl

@@ -3,11 +3,10 @@ using DataFrames
 using LinearAlgebra
 using MatrixDepot
 using SparseArrays
+using SparseMatrixColorings
 using SparseMatrixColorings:
     AdjacencyGraph,
     BipartiteGraph,
-    LargestFirst,
-    NaturalOrder,
     degree,
     minimum_degree,
     maximum_degree,
@@ -19,6 +18,14 @@ using SparseMatrixColorings:
     vertices
 using Test
 
+nbunique(x) = length(unique(x))
+
+_N() = NaturalOrder()
+_LF() = LargestFirst()
+_SL() = SmallestLast(; reproduce_colpack=true)
+_ID() = IncidenceDegree(; reproduce_colpack=true)
+_DLF() = DynamicLargestFirst(; reproduce_colpack=true)
+
 ## Distance-2 coloring
 
 #=
@@ -29,20 +36,50 @@ colpack_table_6_7 = CSV.read(
     joinpath(@__DIR__, "reference", "colpack_table_6_7.csv"), DataFrame
 )
 
-@testset "Distance-2 coloring (ColPack paper)" begin
+@testset verbose = true "Distance-2 coloring (ColPack paper)" begin
     @testset "$(row[:name])" for row in eachrow(colpack_table_6_7)
         original_mat = matrixdepot("$(row[:group])/$(row[:name])")
         mat = dropzeros(original_mat)
         bg = BipartiteGraph(mat)
-        @test nb_vertices(bg, Val(1)) == row[:V1]
-        @test nb_vertices(bg, Val(2)) == row[:V2]
-        @test nb_edges(bg) == row[:E]
-        @test maximum_degree(bg, Val(1)) == row[:Δ1]
-        @test maximum_degree(bg, Val(2)) == row[:Δ2]
-        color_N1 = partial_distance2_coloring(bg, Val(1), NaturalOrder())
-        color_N2 = partial_distance2_coloring(bg, Val(2), NaturalOrder())
-        @test length(unique(color_N1)) == row[:N1]
-        @test length(unique(color_N2)) == row[:N2]
+        @testset "Graph features" begin
+            @test nb_vertices(bg, Val(1)) == row[:V1]
+            @test nb_vertices(bg, Val(2)) == row[:V2]
+            @test nb_edges(bg) == row[:E]
+            @test maximum_degree(bg, Val(1)) == row[:Δ1]
+            @test maximum_degree(bg, Val(2)) == row[:Δ2]
+        end
+        @testset "Natural" begin
+            @test nbunique(partial_distance2_coloring(bg, Val(1), _N())) == row[:N1]
+            @test nbunique(partial_distance2_coloring(bg, Val(2), _N())) == row[:N2]
+        end
+        yield()
+        @testset "LargestFirst" begin
+            @test nbunique(partial_distance2_coloring(bg, Val(1), _LF())) == row[:LF1]
+            @test nbunique(partial_distance2_coloring(bg, Val(2), _LF())) == row[:LF2]
+        end
+        yield()
+        if row[:name] == "af23560"
+            # orders differ for this one, not sure why
+            continue
+        end
+        if row[:E] > 200_000
+            # just to spare computational resources, but the larger tests pass too
+            continue
+        end
+        @testset "SmallestLast" begin
+            @test nbunique(partial_distance2_coloring(bg, Val(1), _SL())) == row[:SL1]
+            @test nbunique(partial_distance2_coloring(bg, Val(2), _SL())) == row[:SL2]
+        end
+        yield()
+        @testset "IncidenceDegree" begin
+            @test nbunique(partial_distance2_coloring(bg, Val(1), _ID())) == row[:ID1]
+            @test nbunique(partial_distance2_coloring(bg, Val(2), _ID())) == row[:ID2]
+        end
+        yield()
+        @testset "DynamicLargestFirst" begin
+            @test nbunique(partial_distance2_coloring(bg, Val(1), _DLF())) == row[:DLF1]
+            @test nbunique(partial_distance2_coloring(bg, Val(2), _DLF())) == row[:DLF2]
+        end
         yield()
     end
 end;