diff --git a/Project.toml b/Project.toml index 24ac9ff9..b239314b 100644 --- a/Project.toml +++ b/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" diff --git a/src/SparseDiffTools.jl b/src/SparseDiffTools.jl index c72b3053..e99ace8b 100644 --- a/src/SparseDiffTools.jl +++ b/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") diff --git a/src/coloring/acyclic_coloring.jl b/src/coloring/acyclic_coloring.jl new file mode 100644 index 00000000..894d9fdb --- /dev/null +++ b/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 diff --git a/src/coloring/backtracking_coloring.jl b/src/coloring/backtracking_coloring.jl index f1d60f6c..510a3d72 100644 --- a/src/coloring/backtracking_coloring.jl +++ b/src/coloring/backtracking_coloring.jl @@ -1,5 +1,3 @@ -using LightGraphs - """ color_graph(g::LightGraphs, ::BacktrackingColor) diff --git a/src/coloring/high_level.jl b/src/coloring/high_level.jl index 25077205..70ff9fa0 100644 --- a/src/coloring/high_level.jl +++ b/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 diff --git a/test/runtests.jl b/test/runtests.jl index 4b1aaab5..2d124e68 100644 --- a/test/runtests.jl +++ b/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 diff --git a/test/test_acyclic.jl b/test/test_acyclic.jl new file mode 100644 index 00000000..3ee1ce8f --- /dev/null +++ b/test/test_acyclic.jl @@ -0,0 +1,128 @@ +using SparseDiffTools +using LightGraphs +using Test + +using Random +Random.seed!(123) + +#= Test data =# +test_graphs = Vector{SimpleGraph}(undef, 0) +test_graphs_dir = Vector{SimpleDiGraph}(undef, 0) + +for _ in 1:5 + nv = rand(5:20) + ne = rand(1:100) + graph = SimpleGraph(nv) + for e in 1:ne + v1 = rand(1:nv) + v2 = rand(1:nv) + while v1 == v2 + v2 = rand(1:nv) + end + add_edge!(graph, v1, v2) + end + push!(test_graphs, copy(graph)) +end + +#= + Coloring needs to satisfy two conditions: + +1. every pair of adjacent vertices receives distinct colors +(a distance-1 coloring) + +2. For any vertex v, any color that leads to a two-colored path +involving v and three other vertices is impermissible for v. +In other words, every path on four vertices uses at least three +colors. +=# + + +#Sample graph from Gebremedhin AH, Manne F, Pothen A. **New Acyclic and Star Coloring Algorithms with Application to Computing Hessians** + +#= + +(1) ----- (2) ----- (3) ---- (4) + | \ / | | + | \ / | | + | \/ | | + (5)---(6) -------- (7) ----(8) + | | / | | + | | (11) | | + | | / | | + (9) ---(10) ----- (12)-----(13) + \ / + \ / + ------------ +=# + +gx = SimpleGraph(13) + +add_edge!(gx,1,2) +add_edge!(gx,1,5) +add_edge!(gx,1,6) +add_edge!(gx,2,3) +add_edge!(gx,2,6) +add_edge!(gx,3,4) +add_edge!(gx,3,7) +add_edge!(gx,4,8) +add_edge!(gx,5,6) +add_edge!(gx,5,9) +add_edge!(gx,6,7) +add_edge!(gx,6,10) +add_edge!(gx,7,8) +add_edge!(gx,7,11) +add_edge!(gx,7,12) +add_edge!(gx,8,13) +add_edge!(gx,9,10) +add_edge!(gx,9,12) +add_edge!(gx,10,11) +add_edge!(gx,10,12) +add_edge!(gx,12,13) + +push!(test_graphs, gx) + +#create directed copies +for g in test_graphs + dg = SimpleDiGraph(nv(g)) + for e in edges(g) + src1 = src(e) + dst1 = dst(e) + add_edge!(dg, src1, dst1) + add_edge!(dg, dst1, src1) + end + push!(test_graphs_dir, dg) +end + + +for i in 1:6 + g = test_graphs[i] + dg = test_graphs_dir[i] + out_colors = SparseDiffTools.color_graph(g, SparseDiffTools.AcyclicColoring()) + + #test condition 1 + for v in vertices(g) + color = out_colors[v] + for j in inneighbors(g, v) + @test out_colors[j] != color + end + end +end + +for i in 3:6 + g = test_graphs[i] + dg = test_graphs_dir[i] + out_colors = SparseDiffTools.color_graph(g, SparseDiffTools.AcyclicColoring()) + + #test condition 2 + cycles = simplecycles(dg) + for c in cycles + colors = zeros(Int, 0) + if length(c) > 2 + for v in c + push!(colors, out_colors[v]) + end + @test length(unique(colors)) >= 3 + end + end + +end