Star coloring without dictionaries

Author: amontoison

src/SparseMatrixColorings.jl

@@ -40,7 +40,8 @@ using SparseArrays:
     nzrange,
     rowvals,
     sparse,
-    spzeros
+    spzeros,
+    triu
 
 include("graph.jl")
 include("order.jl")

src/coloring.jl

@@ -84,7 +84,7 @@ function star_coloring(g::AdjacencyGraph, order::AbstractOrder; postprocessing::
     forbidden_colors = zeros(Int, nv)
     first_neighbor = fill((0, 0), nv)  # at first no neighbors have been encountered
     treated = zeros(Int, nv)
-    star = Dict{Tuple{Int,Int},Int}()
+    star = Vector{Int}(undef, ne)
     sizehint!(star, ne)
     hub = Int[]  # one hub for each star, including the trivial ones
     nb_spokes = Int[]  # number of spokes for each star
@@ -107,7 +107,8 @@ function star_coloring(g::AdjacencyGraph, order::AbstractOrder; postprocessing::
                 for x in neighbors(g, w)
                     (x == v || iszero(color[x])) && continue
                     wx = _sort(w, x)
-                    if x == hub[star[wx]]  # potential Case 2 (which is always false for trivial stars with two vertices, since the associated hub is negative)
+                    index_wx = g.M[wx...]
+                    if x == hub[star[index_wx]]  # potential Case 2 (which is always false for trivial stars with two vertices, since the associated hub is negative)
                         forbidden_colors[color[x]] = v
                     end
                 end
@@ -121,7 +122,7 @@ function star_coloring(g::AdjacencyGraph, order::AbstractOrder; postprocessing::
         end
         _update_stars!(star, hub, nb_spokes, g, v, color, first_neighbor)
     end
-    star_set = StarSet(star, hub, nb_spokes)
+    star_set = StarSet(g, star, hub, nb_spokes)
     if postprocessing
         # Reuse the vector forbidden_colors to compute offsets during post-processing
         postprocess!(color, star_set, g, forbidden_colors)
@@ -140,33 +141,44 @@ $TYPEDFIELDS
 """
 struct StarSet
     "a mapping from edges (pair of vertices) to their star index"
-    star::Dict{Tuple{Int,Int},Int}
+    star::Vector{Int}
     "a mapping from star indices to their hub (undefined hubs for single-edge stars are the negative value of one of the vertices, picked arbitrarily)"
     hub::Vector{Int}
     "a mapping from star indices to the vector of their spokes"
     spokes::Vector{Vector{Int}}
+    M::SparseMatrixCSC{Int,Int}
 end
 
-function StarSet(star::Dict{Tuple{Int,Int},Int}, hub::Vector{Int}, nb_spokes::Vector{Int})
+function StarSet(
+    g::AdjacencyGraph{Int}, star::Vector{Int}, hub::Vector{Int}, nb_spokes::Vector{Int}
+)
     # Create a list of spokes for each star, preallocating their sizes based on nb_spokes
     spokes = [Vector{Int}(undef, ns) for ns in nb_spokes]
 
     # Reuse nb_spokes as counters to track the current index while filling the spokes
     fill!(nb_spokes, 0)
 
-    for ((i, j), s) in pairs(star)
-        h = abs(hub[s])
-        nb_spokes[s] += 1
-        index = nb_spokes[s]
-
-        # Assign the non-hub vertex (spoke) to the correct position in spokes
-        if i == h
-            spokes[s][index] = j
-        elseif j == h
-            spokes[s][index] = i
+    nv = nb_vertices(g)
+    k = 0
+    rows = rowvals(g.M)
+    for j in 1:nv
+        for p in nzrange(g.M, j)
+            i = rows[p]
+            k += 1
+            s = star[k]
+            h = abs(hub[s])
+            nb_spokes[s] += 1
+            index = nb_spokes[s]
+
+            # Assign the non-hub vertex (spoke) to the correct position in spokes
+            if i == h
+                spokes[s][index] = j
+            elseif j == h
+                spokes[s][index] = i
+            end
         end
     end
-    return StarSet(star, hub, spokes)
+    return StarSet(star, hub, spokes, g.M)
 end
 
 _sort(u, v) = (min(u, v), max(u, v))
@@ -191,7 +203,7 @@ end
 
 function _update_stars!(
     # modified
-    star::Dict{<:Tuple,<:Integer},
+    star::AbstractVector{<:Integer},
     hub::AbstractVector{<:Integer},
     nb_spokes::AbstractVector{<:Integer},
     # not modified
@@ -203,14 +215,16 @@ function _update_stars!(
     for w in neighbors(g, v)
         iszero(color[w]) && continue
         vw = _sort(v, w)
+        index_vw = g.M[vw...]
         x_exists = false
         for x in neighbors(g, w)
             if x != v && color[x] == color[v]  # vw, wx ∈ E
                 wx = _sort(w, x)
-                star_wx = star[wx]
+                index_wx = g.M[wx...]
+                star_wx = star[index_wx]
                 hub[star_wx] = w  # this may already be true
                 nb_spokes[star_wx] += 1
-                star[vw] = star_wx
+                star[index_vw] = star_wx
                 x_exists = true
                 break
             end
@@ -219,14 +233,15 @@ function _update_stars!(
             (p, q) = first_neighbor[color[w]]
             if p == v && q != w  # vw, vq ∈ E and color[w] = color[q]
                 vq = _sort(v, q)
-                star_vq = star[vq]
+                index_vq = g.M[vq...]
+                star_vq = star[index_vq]
                 hub[star_vq] = v  # this may already be true
                 nb_spokes[star_vq] += 1
-                star[vw] = star_vq
+                star[index_vw] = star_vq
             else  # vw forms a new star
                 push!(hub, -max(v, w))  # star is trivial (composed only of two vertices) so we set the hub to a negative value, but it allows us to choose one of the two vertices
                 push!(nb_spokes, 1)
-                star[vw] = length(hub)
+                star[index_vw] = length(hub)
             end
         end
     end
@@ -251,7 +266,7 @@ This function corresponds to algorithm `DirectRecover2` in the paper.
 function symmetric_coefficient(
     i::Integer, j::Integer, color::AbstractVector{<:Integer}, star_set::StarSet
 )
-    (; star, hub) = star_set
+    (; M, star, hub) = star_set
     if i == j
         # diagonal
         return i, color[j]
@@ -260,7 +275,8 @@ function symmetric_coefficient(
         # star only contains one triangle
         i, j = j, i
     end
-    star_id = star[i, j]
+    index_ij = M[i, j]
+    star_id = star[index_ij]
     h = abs(hub[star_id])
     if h == j
         # i is the spoke

src/graph.jl

@@ -189,14 +189,23 @@ The adjacency graph of a symmetric matrix `A ∈ ℝ^{n × n}` is `G(A) = (V, E)
 """
 struct AdjacencyGraph{T,has_diagonal}
     S::SparsityPatternCSC{T}
+    M::SparseMatrixCSC{T,T}
 end
 
-function AdjacencyGraph(S::SparsityPatternCSC{T}; has_diagonal::Bool=true) where {T}
-    return AdjacencyGraph{T,has_diagonal}(S)
+function AdjacencyGraph(
+    S::SparsityPatternCSC{T}, M::SparseMatrixCSC{T,T}; has_diagonal::Bool=true
+) where {T}
+    return AdjacencyGraph{T,has_diagonal}(S, M)
 end
 
 function AdjacencyGraph(A::SparseMatrixCSC; has_diagonal::Bool=true)
-    return AdjacencyGraph(SparsityPatternCSC(A); has_diagonal)
+    S = SparsityPatternCSC(A)
+    M = triu(A, 1)
+    nnzM = Int(nnz(M))
+    nzval = collect(Base.OneTo(nnzM))
+    m, n = size(A)
+    M = SparseMatrixCSC(m, n, M.colptr, M.rowval, nzval)
+    return AdjacencyGraph(S, M; has_diagonal)
 end
 
 function AdjacencyGraph(A::AbstractMatrix; has_diagonal::Bool=true)