Update has_diagonal

Author: amontoison

src/coloring.jl

@@ -79,6 +79,7 @@ If `postprocessing=true`, some colors might be replaced with `0` (the "neutral"
 function star_coloring(g::AdjacencyGraph, order::AbstractOrder, postprocessing::Bool)
     # Initialize data structures
     S = pattern(g)
+    has_diag = has_diagonal(g)
     nv = nb_vertices(g)
     ne = nb_edges(g)
     color = zeros(Int, nv)
@@ -91,8 +92,9 @@ function star_coloring(g::AdjacencyGraph, order::AbstractOrder, postprocessing::
     vertices_in_order = vertices(g, order)
 
     for v in vertices_in_order
-        for (iw, w) in enumerate(neighbors2(g, v))
-            (v == w || iszero(color[w])) && continue
+        for (iw, w) in enumerate(neighbors(g, v))
+            !has_diag || (v == w && continue)
+            iszero(color[w]) && continue
             index_vw = edge_to_index[S.colptr[v] + iw - 1]
             forbidden_colors[color[w]] = v
             (p, q, _) = first_neighbor[color[w]]
@@ -105,8 +107,9 @@ function star_coloring(g::AdjacencyGraph, order::AbstractOrder, postprocessing::
                 _treat!(treated, forbidden_colors, g, v, w, color)
             else
                 first_neighbor[color[w]] = (v, w, index_vw)
-                for (ix, x) in enumerate(neighbors2(g, w))
-                    (w == x || x == v || iszero(color[x])) && continue
+                for (ix, x) in enumerate(neighbors(g, w))
+                    !has_diag || (w == x && continue)
+                    (x == v || iszero(color[x])) && continue
                     index_wx = edge_to_index[S.colptr[w] + ix - 1]
                     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
@@ -140,7 +143,9 @@ function _treat!(
     w::Integer,
     color::AbstractVector{<:Integer},
 )
+    has_diag = has_diagonal(g)
     for x in neighbors(g, w)
+        !has_diag || (w == x && continue)
         iszero(color[x]) && continue
         forbidden_colors[color[x]] = v
     end
@@ -160,12 +165,14 @@ function _update_stars!(
     edge_to_index::AbstractVector{<:Integer},
 )
     S = pattern(g)
-    for (iw, w) in enumerate(neighbors2(g, v))
-        (v == w || iszero(color[w])) && continue
+    has_diag = has_diagonal(g)
+    for (iw, w) in enumerate(neighbors(g, v))
+        !has_diag || (v == w && continue)
+        iszero(color[w]) && continue
         index_vw = edge_to_index[S.colptr[v] + iw - 1]
         x_exists = false
-        for (ix, x) in enumerate(neighbors2(g, w))
-            (w == x) && continue
+        for (ix, x) in enumerate(neighbors(g, w))
+            !has_diag || (w == x && continue)
             if x != v && color[x] == color[v]  # vw, wx ∈ E
                 index_wx = edge_to_index[S.colptr[w] + ix - 1]
                 star_wx = star[index_wx]
@@ -232,6 +239,7 @@ If `postprocessing=true`, some colors might be replaced with `0` (the "neutral"
 function acyclic_coloring(g::AdjacencyGraph, order::AbstractOrder, postprocessing::Bool)
     # Initialize data structures
     S = pattern(g)
+    has_diag = has_diagonal(g)
     nv = nb_vertices(g)
     ne = nb_edges(g)
     color = zeros(Int, nv)
@@ -244,13 +252,16 @@ function acyclic_coloring(g::AdjacencyGraph, order::AbstractOrder, postprocessin
 
     for v in vertices_in_order
         for w in neighbors(g, v)
+            !has_diag || (v == w && continue)
             iszero(color[w]) && continue
             forbidden_colors[color[w]] = v
         end
         for w in neighbors(g, v)
+            !has_diag || (v == w && continue)
             iszero(color[w]) && continue
-            for (ix, x) in enumerate(neighbors2(g, w))
-                (w == x || iszero(color[x])) && continue
+            for (ix, x) in enumerate(neighbors(g, w))
+                !has_diag || (w == x && continue)
+                iszero(color[x]) && continue
                 if forbidden_colors[color[x]] != v
                     index_wx = edge_to_index[S.colptr[w] + ix - 1]
                     _prevent_cycle!(
@@ -272,15 +283,18 @@ function acyclic_coloring(g::AdjacencyGraph, order::AbstractOrder, postprocessin
                 break
             end
         end
-        for (iw, w) in enumerate(neighbors2(g, v))  # grow two-colored stars around the vertex v
-            (v == w || iszero(color[w])) && continue
+        for (iw, w) in enumerate(neighbors(g, v))  # grow two-colored stars around the vertex v
+            !has_diag || (v == w && continue)
+            iszero(color[w]) && continue
             index_vw = edge_to_index[S.colptr[v] + iw - 1]
             _grow_star!(v, w, index_vw, color, first_neighbor, forest)
         end
-        for (iw, w) in enumerate(neighbors2(g, v))
-            (v == w || iszero(color[w])) && continue
-            for (ix, x) in enumerate(neighbors2(g, w))
-                (w == x || x == v || iszero(color[x])) && continue
+        for (iw, w) in enumerate(neighbors(g, v))
+            !has_diag || (v == w && continue)
+            iszero(color[w]) && continue
+            for (ix, x) in enumerate(neighbors(g, w))
+                !has_diag || (w == x && continue)
+                (x == v || iszero(color[x])) && continue
                 if color[x] == color[v]
                     index_vw = edge_to_index[S.colptr[v] + iw - 1]
                     index_wx = edge_to_index[S.colptr[w] + ix - 1]

src/graph.jl

@@ -98,10 +98,10 @@ end
 
 Return a [`SparsityPatternCSC`](@ref) corresponding to the matrix `[0 Aᵀ; A 0]`, with a minimum of allocations.
 """
-bidirectional_pattern(A::AbstractMatrix; symmetric_pattern) =
+bidirectional_pattern(A::AbstractMatrix; symmetric_pattern::Bool) =
     bidirectional_pattern(SparsityPatternCSC(SparseMatrixCSC(A)); symmetric_pattern)
 
-function bidirectional_pattern(S::SparsityPatternCSC; symmetric_pattern)
+function bidirectional_pattern(S::SparsityPatternCSC; symmetric_pattern::Bool)
     m, n = size(S)
     p = m + n
     nnzS = nnz(S)
@@ -194,7 +194,7 @@ end
 ## Adjacency graph
 
 """
-    AdjacencyGraph{T,has_diagonal}
+    AdjacencyGraph{T}
 
 Undirected graph without self-loops representing the nonzeros of a symmetric matrix (typically a Hessian matrix).
 
@@ -205,22 +205,24 @@ The adjacency graph of a symmetric matrix `A ∈ ℝ^{n × n}` is `G(A) = (V, E)
 
 # Constructors
 
-    AdjacencyGraph(A::SparseMatrixCSC)
+    AdjacencyGraph(A::SparseMatrixCSC; has_diagonal::Bool=true)
 
 # Fields
 
-- `S::SparsityPatternCSC{T}`: Underlying sparsity pattern, whose diagonal is empty whenever `has_diagonal` is `false`
+- `S::SparsityPatternCSC{T}`: Underlying sparsity pattern
+- `has_diagonal::Bool`: Diagonal is empty whenever `has_diagonal` is `false`
 
 # References
 
 > [_What Color Is Your Jacobian? SparsityPatternCSC Coloring for Computing Derivatives_](https://epubs.siam.org/doi/10.1137/S0036144504444711), Gebremedhin et al. (2005)
 """
-struct AdjacencyGraph{T,has_diagonal}
+struct AdjacencyGraph{T}
     S::SparsityPatternCSC{T}
+    has_diagonal::Bool
 end
 
 function AdjacencyGraph(S::SparsityPatternCSC{T}; has_diagonal::Bool=true) where {T}
-    return AdjacencyGraph{T,has_diagonal}(S)
+    return AdjacencyGraph{T}(S, has_diagonal)
 end
 
 function AdjacencyGraph(A::SparseMatrixCSC; has_diagonal::Bool=true)
@@ -234,22 +236,9 @@ end
 pattern(g::AdjacencyGraph) = g.S
 nb_vertices(g::AdjacencyGraph) = pattern(g).n
 vertices(g::AdjacencyGraph) = 1:nb_vertices(g)
+has_diagonal(g::AdjacencyGraph) = g.has_diagonal
 
-has_diagonal(::AdjacencyGraph{T,hl}) where {T,hl} = hl
-
-function neighbors(g::AdjacencyGraph{T,true}, v::Integer) where {T}
-    S = pattern(g)
-    neighbors_v = view(rowvals(S), nzrange(S, v))
-    return Iterators.filter(!=(v), neighbors_v)  # TODO: optimize
-end
-
-function neighbors(g::AdjacencyGraph{T,false}, v::Integer) where {T}
-    S = pattern(g)
-    neighbors_v = view(rowvals(S), nzrange(S, v))
-    return neighbors_v
-end
-
-function neighbors2(g::AdjacencyGraph, v::Integer)
+function neighbors(g::AdjacencyGraph, v::Integer)
     S = pattern(g)
     neighbors_v = view(rowvals(S), nzrange(S, v))
     return neighbors_v
@@ -266,20 +255,27 @@ function degree(g::AdjacencyGraph, v::Integer)
 end
 
 function nb_edges(g::AdjacencyGraph)
-    ne = 0
-    for v in vertices(g)
-        for u in neighbors(g, v)
-            ne += 1
+    has_diag = has_diagonal(g)
+    if !has_diag
+        S = pattern(g)
+        ne = nnz(S)
+        return ne ÷ 2
+    else
+        ne = 0
+        for v in vertices(g)
+            ne += degree(g, v)
         end
+        return ne ÷ 2
     end
-    return ne ÷ 2
 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)
+        # We don't need the next line if we know that u != v
+        v == w && continue
         if w == u
             return true
         end

src/order.jl

@@ -209,6 +209,7 @@ function vertices(
         u = pop_next_candidate!(db; direction)
         direction == :low2high ? push!(π, u) : pushfirst!(π, u)
         for v in neighbors(g, u)
+            u == v && continue
             already_ordered(db, v) && continue
             update_bucket!(db, v; degtype, direction)
         end