Added acyclic coloring

Project.toml

@@ -8,6 +8,7 @@ Adapt = "79e6a3ab-5dfb-504d-930d-738a2a938a0e"
 ArrayInterface = "4fba245c-0d91-5ea0-9b3e-6abc04ee57a9"
 BandedMatrices = "aae01518-5342-5314-be14-df237901396f"
 BlockBandedMatrices = "ffab5731-97b5-5995-9138-79e8c1846df0"
+DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
 DiffEqDiffTools = "01453d9d-ee7c-5054-8395-0335cb756afa"
 ForwardDiff = "f6369f11-7733-5829-9624-2563aa707210"
 LightGraphs = "093fc24a-ae57-5d10-9952-331d41423f4d"
@@ -19,6 +20,7 @@ VertexSafeGraphs = "19fa3120-7c27-5ec5-8db8-b0b0aa330d6f"
 [compat]
 ArrayInterface = "1.1"
 julia = "1"
+DataStructures = "0.17"
 
 [extras]
 IterativeSolvers = "42fd0dbc-a981-5370-80f2-aaf504508153"

src/SparseDiffTools.jl

@@ -14,6 +14,7 @@ using SparseArrays, ArrayInterface
 
 using BlockBandedMatrices: blocksize, nblocks
 using ForwardDiff: Dual, jacobian, partials, DEFAULT_CHUNK_THRESHOLD
+using DataStructures: DisjointSets, find_root, union!
 
 using ArrayInterface: matrix_colors
 
@@ -39,6 +40,7 @@ include("coloring/high_level.jl")
 include("coloring/backtracking_coloring.jl")
 include("coloring/contraction_coloring.jl")
 include("coloring/greedy_d1_coloring.jl")
+include("coloring/acyclic_coloring.jl")
 include("coloring/greedy_star1_coloring.jl")
 include("coloring/greedy_star2_coloring.jl")
 include("coloring/matrix2graph.jl")

src/coloring/acyclic_coloring.jl

@@ -0,0 +1,201 @@
+"""
+        color_graph(g::LightGraphs.AbstractGraphs, ::AcyclicColoring)
+
+Returns a coloring vector following the acyclic coloring rules (1) the coloring
+corresponds to a distance-1 coloring, and (2) vertices in every cycle of the
+graph are assigned at least three distinct colors. This variant of coloring is
+called acyclic since every subgraph induced by vertices assigned any two colors
+is a collection of trees—and hence is acyclic.
+
+Reference: Gebremedhin AH, Manne F, Pothen A. **New Acyclic and Star Coloring Algorithms with Application to Computing Hessians**
+"""
+function color_graph(g::LightGraphs.AbstractGraph, ::AcyclicColoring)
+
+    color = zeros(Int, nv(g))
+    forbidden_colors = zeros(Int, nv(g))
+
+    set = DisjointSets{LightGraphs.Edge}([])
+
+    first_visit_to_tree = Array{Tuple{Int, Int}, 1}(undef, ne(g))
+    first_neighbor = Array{Tuple{Int, Int}, 1}(undef, nv(g))
+
+    for v in vertices(g)
+        #enforces the first condition of acyclic coloring
+        for w in outneighbors(g, v)
+            if color[w] != 0
+                forbidden_colors[color[w]] = v
+            end
+        end
+        #enforces the second condition of acyclic coloring
+        for w in outneighbors(g, v)
+            if color[w] != 0 #colored neighbor
+                for x in outneighbors(g, w)
+                    if color[x] != 0 #colored x
+                        if forbidden_colors[color[x]] != v
+                            prevent_cycle(v, w, x, g, color, forbidden_colors, first_visit_to_tree, set)
+                        end
+                    end
+                end
+            end
+        end
+
+        color[v] = min_index(forbidden_colors, v)
+
+        #grow star for every edge connecting colored vertices v and w
+        for w in outneighbors(g, v)
+            if color[w] != 0
+                grow_star!(set, v, w, g, first_neighbor, color)
+            end
+        end
+
+        #merge the newly formed stars into existing trees if possible
+        for w in outneighbors(g, v)
+            if color[w] != 0
+                for x in outneighbors(g, w)
+                    if color[x] != 0 && x != v
+                        if color[x] == color[v]
+                            merge_trees!(set, v, w, x, g)
+                        end
+                    end
+                end
+            end
+        end
+    end
+
+    return color
+end
+
+"""
+        prevent_cycle(v::Integer,
+                    w::Integer,
+                    x::Integer,
+                    g::LightGraphs.AbstractGraph,
+                    color::AbstractVector{<:Integer},
+                    forbidden_colors::AbstractVector{<:Integer},
+                    first_visit_to_tree::Array{Tuple{Integer, Integer}, 1},
+                    set::DisjointSets{LightGraphs.Edge})
+
+Subroutine to avoid generation of 2-colored cycle due to coloring of vertex v,
+which is adjacent to vertices w and x in graph g. Disjoint set is used to store
+the induced 2-colored subgraphs/trees where the id of set is a key edge of g
+"""
+function prevent_cycle(v::Integer,
+                        w::Integer,
+                        x::Integer,
+                        g::LightGraphs.AbstractGraph,
+                        color::AbstractVector{<:Integer},
+                        forbidden_colors::AbstractVector{<:Integer},
+                        first_visit_to_tree::AbstractVector{<: Tuple{Integer, Integer}},
+                        set::DisjointSets{LightGraphs.Edge})
+
+    edge = find_edge(g, w, x)
+    e = find_root(set, edge)
+    p, q = first_visit_to_tree[edge_index(g, e)]
+    if p != v
+        first_visit_to_tree[edge_index(g, e)] = (v, w)
+    elseif q != w
+        forbidden_colors[color[x]] = v
+    end
+end
+
+"""
+        min_index(forbidden_colors::AbstractVector{<:Integer}, v::Integer)
+
+Returns min{i > 0 such that forbidden_colors[i] != v}
+"""
+function min_index(forbidden_colors::AbstractVector{<:Integer}, v::Integer)
+    return findfirst(!isequal(v), forbidden_colors)
+end
+
+"""
+        grow_star!(set::DisjointSets{LightGraphs.Edge},
+                v::Integer,
+                w::Integer,
+                g::LightGraphs.AbstractGraph
+                first_neighbor::Array{Tuple{Integer, Integer}, 1})
+
+Subroutine to grow a 2-colored star after assigning a new color to the
+previously uncolored vertex v, by comparing it with the adjacent vertex w.
+Disjoint set is used to store stars in sets, which are identified through key
+edges present in g.
+"""
+function grow_star!(set::DisjointSets{LightGraphs.Edge},
+                   v::Integer,
+                   w::Integer,
+                   g::LightGraphs.AbstractGraph,
+                   first_neighbor::AbstractArray{<: Tuple{Integer, Integer}, 1},
+                   color::AbstractVector{<: Integer})
+    edge = find_edge(g, v, w)
+    push!(set, edge)
+    p, q = first_neighbor[color[w]]
+    if p != v
+        first_neighbor[color[w]] = (v, w)
+    else
+        edge1 = find_edge(g, v, w)
+        edge2 = find_edge(g, p, q)
+        e1 = find_root(set, edge1)
+        e2 = find_root(set, edge2)
+        union!(set, e1, e2)
+    end
+    return nothing
+end
+
+
+"""
+        merge_trees!(v::Integer,
+                w::Integer,
+                x::Integer,
+                g::LightGraphs.AbstractGraph,
+                set::DisjointSets{LightGraphs.Edge})
+
+Subroutine to merge trees present in the disjoint set which have a
+common edge.
+"""
+function merge_trees!(set::DisjointSets{LightGraphs.Edge},
+                    v::Integer,
+                    w::Integer,
+                    x::Integer,
+                    g::LightGraphs.AbstractGraph)
+    edge1 = find_edge(g, v, w)
+    edge2 = find_edge(g, w, x)
+    e1 = find_root(set, edge1)
+    e2 = find_root(set, edge2)
+    if (e1 != e2)
+        union!(set, e1, e2)
+    end
+end
+
+
+"""
+        find_edge(g::LightGraphs.AbstractGraph, v::Integer, w::Integer)
+
+Returns an edge object of the type LightGraphs.Edge which represents the
+edge connecting vertices v and w of the undirected graph g
+"""
+function find_edge(g::LightGraphs.AbstractGraph,
+                   v::Integer,
+                   w::Integer)
+    for e in edges(g)
+        if (src(e) == v && dst(e) == w) || (src(e) == w && dst(e) == v)
+            return e
+        end
+    end
+    throw(ArgumentError("$v and $w are not connected in graph g"))
+end
+
+"""
+        edge_index(g::LightGraphs.AbstractGraph, e::LightGraphs.Edge)
+
+Returns an Integer value which uniquely identifies the edge e in graph
+g. Used as an index in main function to avoid custom arrays with non-
+numerical indices.
+"""
+function edge_index(g::LightGraphs.AbstractGraph,
+                    e::LightGraphs.Edge)
+    for (i, edge) in enumerate(edges(g))
+        if edge == e
+            return i
+        end
+    end
+    throw(ArgumentError("Edge $e is not present in graph g"))
+end

src/coloring/backtracking_coloring.jl

@@ -1,5 +1,3 @@
-using LightGraphs
-
 """
     color_graph(g::LightGraphs, ::BacktrackingColor)
 

src/coloring/high_level.jl

@@ -1,20 +1,21 @@
-abstract type SparseDiffToolsColoringAlgorithm <: ArrayInterface.ColoringAlgorithm end
-struct GreedyD1Color <: SparseDiffToolsColoringAlgorithm end
-struct BacktrackingColor <: SparseDiffToolsColoringAlgorithm end
-struct ContractionColor <: SparseDiffToolsColoringAlgorithm end
-struct GreedyStar1Color <: SparseDiffToolsColoringAlgorithm end
-struct GreedyStar2Color <: SparseDiffToolsColoringAlgorithm end
-
-"""
-    matrix_colors(A,alg::ColoringAlgorithm = GreedyD1Color())
-
-    Returns the colorvec vector for the matrix A using the chosen coloring
-    algorithm. If a known analytical solution exists, that is used instead.
-    The coloring defaults to a greedy distance-1 coloring.
-
-"""
-function ArrayInterface.matrix_colors(A::AbstractMatrix,alg::SparseDiffToolsColoringAlgorithm = GreedyD1Color(); partition_by_rows::Bool = false)
-    _A = A isa SparseMatrixCSC ? A : sparse(A) # Avoid the copy
-    A_graph = matrix2graph(_A, partition_by_rows)
-    color_graph(A_graph,alg)
-end
+abstract type SparseDiffToolsColoringAlgorithm <: ArrayInterface.ColoringAlgorithm end
+struct GreedyD1Color <: SparseDiffToolsColoringAlgorithm end
+struct BacktrackingColor <: SparseDiffToolsColoringAlgorithm end
+struct ContractionColor <: SparseDiffToolsColoringAlgorithm end
+struct GreedyStar1Color <: SparseDiffToolsColoringAlgorithm end
+struct GreedyStar2Color <: SparseDiffToolsColoringAlgorithm end
+struct AcyclicColoring <: SparseDiffToolsColoringAlgorithm end
+
+"""
+    matrix_colors(A,alg::ColoringAlgorithm = GreedyD1Color())
+
+    Returns the colorvec vector for the matrix A using the chosen coloring
+    algorithm. If a known analytical solution exists, that is used instead.
+    The coloring defaults to a greedy distance-1 coloring.
+
+"""
+function ArrayInterface.matrix_colors(A::AbstractMatrix, alg::SparseDiffToolsColoringAlgorithm = GreedyD1Color(); partition_by_rows::Bool = false)
+    _A = A isa SparseMatrixCSC ? A : sparse(A) # Avoid the copy
+    A_graph = matrix2graph(_A, partition_by_rows)
+    return color_graph(A_graph, alg)
+end

test/runtests.jl

@@ -8,6 +8,7 @@ if GROUP == "All"
     @time @safetestset "Exact coloring via contraction" begin include("test_contraction.jl") end
     @time @safetestset "Greedy distance-1 coloring" begin include("test_greedy_d1.jl") end
     @time @safetestset "Greedy star coloring" begin include("test_greedy_star.jl") end
+    @time @safetestset "Acyclic coloring" begin include("test_acyclic.jl") end
     @time @safetestset "Matrix to graph conversion" begin include("test_matrix2graph.jl") end
     @time @safetestset "AD using colorvec vector" begin include("test_ad.jl") end
     @time @safetestset "Integration test" begin include("test_integration.jl") end