Degree-based orderings, take 2 (#125)

* Degree-based orderings

* Fix docs

* Clearly separate `AdjacencyGraph` and `BipartiteGraph` from the underlying sparsity pattern

* Typos

* Fix

* Fixes

* Implement for BipartiteGraph

* No abstract

* Fix

* Codecov

* Fix compat

* Tested dynamic orders

* Fixes

* Add experimental warning

Author: gdalle

docs/src/api.md

@@ -47,4 +47,8 @@ AbstractOrder
 NaturalOrder
 RandomOrder
 LargestFirst
+SmallestLast
+IncidenceDegree
+DynamicLargestFirst
+DynamicDegreeBasedOrder
 ```

docs/src/dev.md

@@ -46,6 +46,7 @@ SparseMatrixColorings.LinearSystemColoringResult
 SparseMatrixColorings.directly_recoverable_columns
 SparseMatrixColorings.symmetrically_orthogonal_columns
 SparseMatrixColorings.structurally_orthogonal_columns
+SparseMatrixColorings.valid_dynamic_order
 ```
 
 ## Matrix handling

src/SparseMatrixColorings.jl

@@ -56,6 +56,7 @@ include("examples.jl")
 include("show_colors.jl")
 
 export NaturalOrder, RandomOrder, LargestFirst
+export DynamicDegreeBasedOrder, SmallestLast, IncidenceDegree, DynamicLargestFirst
 export ColoringProblem, GreedyColoringAlgorithm, AbstractColoringResult
 export ConstantColoringAlgorithm
 export coloring

src/check.jl

@@ -138,3 +138,67 @@ function directly_recoverable_columns(
     end
     return true
 end
+
+"""
+    valid_dynamic_order(g::AdjacencyGraph, π::AbstractVector{Int}, order::DynamicDegreeBasedOrder)
+    valid_dynamic_order(bg::AdjacencyGraph, ::Val{side}, π::AbstractVector{Int}, order::DynamicDegreeBasedOrder)
+
+Check that a permutation `π` corresponds to a valid application of a [`DynamicDegreeBasedOrder`](@ref).
+
+This is done by checking, for each ordered vertex, that its back- or forward-degree was the smallest or largest among the remaining vertices (the specifics depend on the order parameters).
+
+!!! warning
+    This function is not coded with efficiency in mind, it is designed for small-scale tests.
+"""
+function valid_dynamic_order(
+    g::AdjacencyGraph, π::AbstractVector{Int}, ::DynamicDegreeBasedOrder{degtype,direction}
+) where {degtype,direction}
+    length(π) != nb_vertices(g) && return false
+    length(unique(π)) != nb_vertices(g) && return false
+    for i in eachindex(π)
+        vi = π[i]
+        yet_to_be_ordered = direction == :low2high ? π[i:end] : π[begin:i]
+        considered_for_degree = degtype == :back ? π[begin:(i - 1)] : π[(i + 1):end]
+        di = degree_in_subset(g, vi, considered_for_degree)
+        considered_for_degree_switched = copy(considered_for_degree)
+        for vj in yet_to_be_ordered
+            replace!(considered_for_degree_switched, vj => vi)
+            dj = degree_in_subset(g, vj, considered_for_degree_switched)
+            replace!(considered_for_degree_switched, vi => vj)
+            if direction == :low2high
+                dj > di && return false
+            else
+                dj < di && return false
+            end
+        end
+    end
+    return true
+end
+
+function valid_dynamic_order(
+    g::BipartiteGraph,
+    ::Val{side},
+    π::AbstractVector{Int},
+    ::DynamicDegreeBasedOrder{degtype,direction},
+) where {side,degtype,direction}
+    length(π) != nb_vertices(g, Val(side)) && return false
+    length(unique(π)) != nb_vertices(g, Val(side)) && return false
+    for i in eachindex(π)
+        vi = π[i]
+        yet_to_be_ordered = direction == :low2high ? π[i:end] : π[begin:i]
+        considered_for_degree = degtype == :back ? π[begin:(i - 1)] : π[(i + 1):end]
+        di = degree_dist2_in_subset(g, Val(side), vi, considered_for_degree)
+        considered_for_degree_switched = copy(considered_for_degree)
+        for vj in yet_to_be_ordered
+            replace!(considered_for_degree_switched, vj => vi)
+            dj = degree_dist2_in_subset(g, Val(side), vj, considered_for_degree_switched)
+            replace!(considered_for_degree_switched, vi => vj)
+            if direction == :low2high
+                dj > di && return false
+            else
+                dj < di && return false
+            end
+        end
+    end
+    return true
+end

src/graph.jl

@@ -161,6 +161,25 @@ end
 maximum_degree(g::AdjacencyGraph) = maximum(Base.Fix1(degree, g), vertices(g))
 minimum_degree(g::AdjacencyGraph) = minimum(Base.Fix1(degree, g), vertices(g))
 
+function has_neighbor(g::AdjacencyGraph, v::Integer, u::Integer)
+    for w in neighbors(g, v)
+        if w == u
+            return true
+        end
+    end
+    return false
+end
+
+function degree_in_subset(g::AdjacencyGraph, v::Integer, subset::AbstractVector{Int})
+    d = 0
+    for u in subset
+        if has_neighbor(g, v, u)
+            d += 1
+        end
+    end
+    return d
+end
+
 ## Bipartite graph
 
 """
@@ -235,6 +254,18 @@ function neighbors(bg::BipartiteGraph, ::Val{side}, v::Integer) where {side}
     return view(rowvals(S), nzrange(S, v))
 end
 
+function neighbors_dist2(bg::BipartiteGraph{T}, ::Val{side}, v::Integer) where {T,side}
+    # TODO: make more efficient
+    other_side = 3 - side
+    neigh = Set{T}()
+    for u in neighbors(bg, Val(side), v)
+        for w in neighbors(bg, Val(other_side), u)
+            w != v && push!(neigh, w)
+        end
+    end
+    return neigh
+end
+
 function degree(bg::BipartiteGraph, ::Val{side}, v::Integer) where {side}
     return length(neighbors(bg, Val(side), v))
 end
@@ -248,13 +279,31 @@ function minimum_degree(bg::BipartiteGraph, ::Val{side}) where {side}
 end
 
 function degree_dist2(bg::BipartiteGraph{T}, ::Val{side}, v::Integer) where {T,side}
-    # not efficient, for testing purposes only
+    return length(neighbors_dist2(bg, Val(side), v))
+end
+
+function has_neighbor_dist2(
+    bg::BipartiteGraph, ::Val{side}, v::Integer, u::Integer
+) where {side}
     other_side = 3 - side
-    neighbors_dist2 = Set{T}()
-    for u in neighbors(bg, Val(side), v)
-        for w in neighbors(bg, Val(other_side), u)
-            w != v && push!(neighbors_dist2, w)
+    for w1 in neighbors(bg, Val(side), v)
+        for w2 in neighbors(bg, Val(other_side), w1)
+            if w2 == u
+                return true
+            end
         end
     end
-    return length(neighbors_dist2)
+    return false
+end
+
+function degree_dist2_in_subset(
+    bg::BipartiteGraph, ::Val{side}, v::Integer, subset::AbstractVector{Int}
+) where {side}
+    d = 0
+    for u in subset
+        if has_neighbor_dist2(bg, Val(side), v, u)
+            d += 1
+        end
+    end
+    return d
 end

src/order.jl

@@ -6,11 +6,14 @@ Abstract supertype for the vertex order used inside [`GreedyColoringAlgorithm`](
 In this algorithm, the rows and columns of a matrix form a graph, and the vertices are colored one after the other in a greedy fashion.
 Depending on how the vertices are ordered, the number of colors necessary may vary.
 
-# Subtypes
+# Options
 
 - [`NaturalOrder`](@ref)
 - [`RandomOrder`](@ref)
 - [`LargestFirst`](@ref)
+- [`IncidenceDegree`](@ref) (experimental)
+- [`SmallestLast`](@ref) (experimental)
+- [`DynamicLargestFirst`](@ref) (experimental)
 """
 abstract type AbstractOrder end
 
@@ -79,3 +82,189 @@ function vertices(bg::BipartiteGraph, ::Val{side}, ::LargestFirst) where {side}
     criterion(v) = degrees_dist2[v]
     return sort(vertices(bg, Val(side)); by=criterion, rev=true)
 end
+
+"""
+    DynamicDegreeBasedOrder{degtype,direction}
+
+Instance of [`AbstractOrder`](@ref) which sorts vertices using a dynamically computed degree.
+
+!!! danger
+    This order is still experimental and needs more tests, correctness is not yet guaranteed.
+
+# Type parameters
+
+- `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`)
+
+# Concrete variants
+
+- [`IncidenceDegree`](@ref)
+- [`SmallestLast`](@ref)
+- [`DynamicLargestFirst`](@ref)
+
+# References
+
+- [_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 DegreeBuckets{B}
+    degrees::Vector{Int}
+    buckets::B
+    positions::Vector{Int}
+end
+
+function DegreeBuckets(degrees::Vector{Int}, dmax)
+    buckets = Dict(d => Int[] for d in 0:dmax)
+    positions = similar(degrees, Int)
+    for v in eachindex(degrees, positions)
+        d = degrees[v]
+        push!(buckets[d], v)  # TODO: optimize
+        positions[v] = length(buckets[d])
+    end
+    return DegreeBuckets(degrees, buckets, positions)
+end
+
+function degree_increasing(; degtype, direction)
+    increasing =
+        (degtype == :back && direction == :low2high) ||
+        (degtype == :forward && direction == :high2low)
+    return increasing
+end
+
+function mark_ordered!(db::DegreeBuckets, v::Integer)
+    db.degrees[v] = -1
+    db.positions[v] = -1
+    return nothing
+end
+
+already_ordered(db::DegreeBuckets, v::Integer) = db.degrees[v] == -1
+
+function pop_next_candidate!(db::DegreeBuckets; direction::Symbol)
+    (; buckets) = db
+    if direction == :low2high
+        candidate_degree = maximum(d for (d, bucket) in pairs(buckets) if !isempty(bucket))
+    else
+        candidate_degree = minimum(d for (d, bucket) in pairs(buckets) if !isempty(bucket))
+    end
+    candidate_bucket = buckets[candidate_degree]
+    candidate = pop!(candidate_bucket)
+    mark_ordered!(db, candidate)
+    return candidate
+end
+
+function update_bucket!(db::DegreeBuckets, v::Integer; degtype, direction)
+    (; degrees, buckets, positions) = db
+    d, p = degrees[v], positions[v]
+    bucket = buckets[d]
+    # select previous or next bucket for the move
+    d_new = degree_increasing(; degtype, direction) ? d + 1 : d - 1
+    bucket_new = buckets[d_new]
+    # put v at the end of its bucket by swapping
+    w = bucket[end]
+    bucket[p] = w
+    positions[w] = p
+    bucket[end] = v
+    positions[v] = length(bucket)
+    # move v from the old bucket to the new one
+    @assert pop!(bucket) == v
+    push!(bucket_new, v)
+    degrees[v] = d_new
+    positions[v] = length(bucket_new)
+    return nothing
+end
+
+function vertices(
+    g::AdjacencyGraph, ::DynamicDegreeBasedOrder{degtype,direction}
+) where {degtype,direction}
+    if degree_increasing(; degtype, direction)
+        degrees = zeros(Int, nb_vertices(g))
+    else
+        degrees = [degree(g, v) for v in vertices(g)]
+    end
+    db = DegreeBuckets(degrees, maximum_degree(g))
+    π = Int[]
+    for _ in 1:nb_vertices(g)
+        u = pop_next_candidate!(db; direction)
+        direction == :low2high ? push!(π, u) : pushfirst!(π, u)
+        for v in neighbors(g, u)
+            already_ordered(db, v) && continue
+            update_bucket!(db, v; degtype, direction)
+        end
+    end
+    return π
+end
+
+function vertices(
+    g::BipartiteGraph, ::Val{side}, ::DynamicDegreeBasedOrder{degtype,direction}
+) where {side,degtype,direction}
+    other_side = 3 - side
+    if degree_increasing(; degtype, direction)
+        degrees = zeros(Int, nb_vertices(g, Val(side)))
+    else
+        degrees = [degree_dist2(g, Val(side), v) for v in vertices(g, Val(side))]  # TODO: optimize
+    end
+    maxd2 = maximum(v -> degree_dist2(g, Val(side), v), vertices(g, Val(side)))  # TODO: optimize
+    db = DegreeBuckets(degrees, maxd2)
+    π = Int[]
+    visited = falses(nb_vertices(g, Val(side)))
+    for _ in 1:nb_vertices(g, Val(side))
+        u = pop_next_candidate!(db; direction)
+        direction == :low2high ? push!(π, u) : pushfirst!(π, u)
+        for w in neighbors(g, Val(side), u)
+            for v in neighbors(g, Val(other_side), w)
+                if v == u || visited[v]
+                    continue
+                else
+                    visited[v] = true
+                end
+                already_ordered(db, v) && continue
+                update_bucket!(db, v; degtype, direction)
+            end
+        end
+        fill!(visited, false)
+    end
+    return π
+end
+
+"""
+    IncidenceDegree()
+
+Instance of [`AbstractOrder`](@ref) which sorts vertices from lowest to highest using the dynamic back degree.
+
+!!! danger
+    This order is still experimental and needs more tests, correctness is not yet guaranteed.
+
+# See also
+
+- [`DynamicDegreeBasedOrder`](@ref)
+"""
+const IncidenceDegree = DynamicDegreeBasedOrder{:back,:low2high}
+
+"""
+    SmallestLast()
+
+Instance of [`AbstractOrder`](@ref) which sorts vertices from highest to lowest using the dynamic back degree.
+
+!!! danger
+    This order is still experimental and needs more tests, correctness is not yet guaranteed.
+
+# See also
+
+- [`DynamicDegreeBasedOrder`](@ref)
+"""
+const SmallestLast = DynamicDegreeBasedOrder{:back,:high2low}
+
+"""
+    DynamicLargestFirst()
+
+Instance of [`AbstractOrder`](@ref) which sorts vertices from lowest to highest using the dynamic forward degree.
+
+!!! danger
+    This order is still experimental and needs more tests, correctness is not yet guaranteed.
+
+# See also
+    
+- [`DynamicDegreeBasedOrder`](@ref)
+"""
+const DynamicLargestFirst = DynamicDegreeBasedOrder{:forward,:low2high}