Replace chordal completion implementation.

Author: samuelsonric

docs/src/constraints.md

@@ -344,7 +344,6 @@ SumOfSquares.PolyJuMP.bridges
 Chordal extension:
 ```@docs
 SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph.neighbors
-SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph.fill_in
 SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph.is_clique
 SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph.LabelledGraph
 SumOfSquares.Certificate.Sparsity.ChordalExtensionGraph.add_node!

src/Certificate/Sparsity/ChordalExtensionGraph.jl

@@ -1,9 +1,8 @@
 module ChordalExtensionGraph
 
-using CliqueTrees
-using CliqueTrees: EliminationAlgorithm
-using DataStructures
-using SparseArrays
+import CliqueTrees
+import DataStructures
+import SparseArrays
 
 export AbstractCompletion, ChordalCompletion, ClusterCompletion
 export LabelledGraph, add_node!, add_edge!, add_clique!, chordal_extension
@@ -12,35 +11,30 @@ abstract type AbstractGreedyAlgorithm end
 
 abstract type AbstractCompletion end
 
-struct ChordalCompletion{A<:AbstractGreedyAlgorithm} <: AbstractCompletion
-    algo::A
-end
-ChordalCompletion() = ChordalCompletion(GreedyFillIn())
-
 struct ClusterCompletion <: AbstractCompletion end
 
-struct CliqueTreesCompletion{A<:EliminationAlgorithm} <: AbstractCompletion
+struct ChordalCompletion{A<:CliqueTrees.EliminationAlgorithm} <: AbstractCompletion
     algo::A
 end
 
-function CliqueTreesCompletion()
-    # the default algorithm is the minimum fill heuristic
-    CliqueTreesCompletion(MF())
+# the default algorithm is the greedy minimum fill heuristic
+const GreedyFillIn = CliqueTrees.MF
+
+function ChordalCompletion()
+    ChordalCompletion(GreedyFillIn())
 end
 
 # With a `Vector{Vector{Int}}` with unsorted neighbors, computing fill-in
 # would be inefficient.
 # With a `Vector{Vector{Int}}` with sorted neighbors, adding an edge would be
 # inefficient.
 # With `Vector{Set{Int}}` both adding edges and computing fill-in is efficient.
-
 struct Graph
     neighbors::Vector{Set{Int}}
-    disabled::BitSet
 end
-Graph() = Graph(Set{Int}[], BitSet())
+Graph() = Graph(Set{Int}[])
 Base.broadcastable(g::Graph) = Ref(g)
-Base.copy(G::Graph) = Graph(deepcopy(G.neighbors), copy(G.disabled))
+Base.copy(G::Graph) = Graph(deepcopy(G.neighbors))
 function add_node!(G::Graph)
     push!(G.neighbors, Set{Int}())
     return length(G.neighbors)
@@ -54,15 +48,14 @@ function add_edge!(G::Graph, i::Int, j::Int)
     return (i, j)
 end
 function add_clique!(G::Graph, nodes::Vector{Int})
-    for i in 1:(length(nodes)-1)
-        for j in (i+1):length(nodes)
+    n = length(nodes)
+
+    for i in 1:(n-1)
+        for j in (i+1):n
             add_edge!(G, nodes[i], nodes[j])
         end
     end
 end
-disable_node!(G::Graph, node::Int) = push!(G.disabled, node)
-is_enabled(G::Graph, node::Int) = !(node in G.disabled)
-enable_all_nodes!(G::Graph) = empty!(G.disabled)
 
 """
     neighbors(G::Graph, node::Int}
@@ -73,94 +66,24 @@ function neighbors(G::Graph, node::Int)
     return G.neighbors[node]
 end
 
-function _num_edges_subgraph(
-    G::Graph,
-    nodes::Union{Vector{Int},Set{Int}},
-    node::Int,
-)
-    neighs = neighbors(G, node)
-    return count(nodes) do node
-        return is_enabled(G, node) && node in neighs
-    end
-end
-function num_edges_subgraph(G::Graph, nodes::Union{Vector{Int},Set{Int}})
-    return mapreduce(+, nodes; init = 0) do node
-        return is_enabled(G, node) ? _num_edges_subgraph(G, nodes, node) : 0
-    end
-end
-
-function num_missing_edges_subgraph(
-    G::Graph,
-    nodes::Union{Vector{Int},Set{Int}},
-)
-    n = count(node -> is_enabled(G, node), nodes)
-    # A clique is a completely connected graph. As such it has n*(n-1)/2 undirected
-    # or equivalently n*(n-1) directed edges.
-    return div(n * (n - 1) - num_edges_subgraph(G, nodes), 2)
-end
-
-"""
-    fill_in(G::Graph{T}, i::T}
-
-Return number of edges that need to be added to make the neighbors of `i` a clique.
-"""
-function fill_in(G::Graph, node::Int)
-    return num_missing_edges_subgraph(G, neighbors(G, node))
-end
-
 """
     is_clique(G::Graph{T}, x::Vector{T})
 
 Return a `Bool` indication whether `x` is a clique in `G`.
 """
 function is_clique(G::Graph, nodes::Vector{Int})
-    return iszero(num_missing_edges_subgraph(G, nodes))
-end
+    n = length(nodes)
 
-struct FillInCache
-    graph::Graph
-    fill_in::Vector{Int}
-end
-function FillInCache(graph::Graph)
-    return FillInCache(graph, [fill_in(graph, i) for i in 1:num_nodes(graph)])
-end
-Base.copy(G::FillInCache) = FillInCache(copy(G.graph), copy(G.fill_in))
-
-num_nodes(G::FillInCache) = num_nodes(G.graph)
-neighbors(G::FillInCache, node::Int) = neighbors(G.graph, node)
-function add_edge!(G::FillInCache, i::Int, j::Int)
-    ni = neighbors(G, i)
-    nj = neighbors(G, j)
-    if i in nj
-        @assert j in ni
-        return
-    end
-    @assert !(j in ni)
-    for node in ni
-        @assert node != i
-        @assert node != j
-        if node in nj
-            G.fill_in[node] -= 1
+    for i in 1:(n-1)
+        for j in (i+1):n
+            if nodes[j] ∉ neighbors(G, nodes[i])
+                return false
+            end
         end
     end
-    G.fill_in[i] +=
-        count(k -> is_enabled(G.graph, k), ni) -
-        _num_edges_subgraph(G.graph, ni, j)
-    G.fill_in[j] +=
-        count(k -> is_enabled(G.graph, k), nj) -
-        _num_edges_subgraph(G.graph, nj, i)
-    return add_edge!(G.graph, i, j)
-end
-fill_in(G::FillInCache, node::Int) = G.fill_in[node]
-function disable_node!(G::FillInCache, node::Int)
-    for neighbor in neighbors(G, node)
-        nodes = neighbors(G, neighbor)
-        G.fill_in[neighbor] -=
-            (length(nodes) - 1) - _num_edges_subgraph(G.graph, nodes, node)
-    end
-    return disable_node!(G.graph, node)
+
+    return true
 end
-is_enabled(G::FillInCache, node::Int) = is_enabled(G.graph, node)
 
 """
     struct LabelledGraph{T}
@@ -256,116 +179,24 @@ function add_clique!(G::LabelledGraph{T}, x::Vector{T}) where {T}
     return add_clique!(G.graph, add_node!.(G, x))
 end
 
-struct GreedyFillIn <: AbstractGreedyAlgorithm end
-cache(G::Graph, ::GreedyFillIn) = FillInCache(G)
-heuristic_value(G::FillInCache, node::Int, ::GreedyFillIn) = fill_in(G, node)
-
-function _greedy_triangulation!(G, algo::AbstractGreedyAlgorithm)
-    candidate_cliques = Vector{Int}[]
-    for i in 1:num_nodes(G)
-        node = argmin(map(1:num_nodes(G)) do node
-            if is_enabled(G, node)
-                heuristic_value(G, node, algo)
-            else
-                typemax(Int)
-            end
-        end)
-        #=
-                # look at all its neighbors. In G we want the neighbors to form a clique.
-                neighbor_nodes = collect(neighbors(G, node))
-
-                # add neighbors and node as a potentially maximal clique
-                candidate_clique = [neighbor for neighbor in neighbor_nodes if is_enabled(G, neighbor)]
-                push!(candidate_clique, node)
-                push!(candidate_cliques, candidate_clique)
-
-                # add edges to G to make the new candidate_clique at least potentially a clique
-                for i in eachindex(neighbor_nodes)
-                    if is_enabled(G, i)
-                        for j in (i + 1):length(neighbor_nodes)
-                            if is_enabled(G, j)
-                                add_edge!(G, neighbor_nodes[i], neighbor_nodes[j])
-                            end
-                        end
-                    end
-                end
-        =#
-        neighbor_nodes = [
-            neighbor for
-            neighbor in neighbors(G, node) if is_enabled(G, neighbor)
-        ]
-        push!(candidate_cliques, [node, neighbor_nodes...])
-        for i in eachindex(neighbor_nodes)
-            for j in (i+1):length(neighbor_nodes)
-                add_edge!(G, neighbor_nodes[i], neighbor_nodes[j])
-            end
-        end
-
-        disable_node!(G, node)
-    end
-    return candidate_cliques
-end
-
-"""
-    completion(G::Graph, comp::ChordalCompletion)
-
-Return a chordal extension of `G` and the corresponding maximal cliques.
-
-The algoritm is Algorithm 3 in [BA10] with the GreedyFillIn heuristic of Table I.
-
-[BA10] Bodlaender, Hans L., and Arie MCA Koster.
-Treewidth computations I. Upper bounds.
-Information and Computation 208, no. 3 (2010): 259-275.
-Utrecht University, Utrecht, The Netherlands www.cs.uu.nl
-"""
-function completion(G::Graph, comp::ChordalCompletion)
-    H = copy(G)
-
-    candidate_cliques = _greedy_triangulation!(cache(H, comp.algo), comp.algo)
-    enable_all_nodes!(H)
-
-    sort!.(candidate_cliques)
-    unique!(candidate_cliques)
-
-    # check whether candidate cliques are actually cliques
-    candidate_cliques = candidate_cliques[[
-        is_clique(H, clique) for clique in candidate_cliques
-    ]]
-    sort!(candidate_cliques, by = x -> length(x))
-    reverse!(candidate_cliques) # TODO use `rev=true` in `sort!`.
-
-    maximal_cliques = [first(candidate_cliques)]
-    for clique in Iterators.drop(candidate_cliques, 1)
-        if all(other_clique -> !(clique ⊆ other_clique), maximal_cliques)
-            push!(maximal_cliques, clique)
-        end
-    end
-
-    if length(Set(Iterators.flatten(maximal_cliques))) != num_nodes(H)
-        error("Maximal cliques do not cover all nodes.")
-    end
-
-    return H, maximal_cliques
-end
-
 """
-    completion(G::Graph, comp::ChordalCompletion)
+    completion(G::Graph, comp::ClusterCompletion)
 
 Return a cluster completion of `G` and the corresponding maximal cliques.
 """
 function completion(G::Graph, ::ClusterCompletion)
     H = copy(G)
-    union_find = IntDisjointSets(num_nodes(G))
+    union_find = DataStructures.IntDisjointSets(num_nodes(G))
     for from in 1:num_nodes(G)
         for to in neighbors(G, from)
-            union!(union_find, from, to)
+            DataStructures.union!(union_find, from, to)
         end
     end
-    cliques = [Int[] for i in 1:num_groups(union_find)]
+    cliques = [Int[] for i in 1:DataStructures.num_groups(union_find)]
     ids = zeros(Int, num_nodes(G))
     k = 0
     for node in 1:num_nodes(G)
-        root = find_root!(union_find, node)
+        root = DataStructures.find_root!(union_find, node)
         if iszero(ids[root])
             k += 1
             ids[root] = k
@@ -379,38 +210,46 @@ function completion(G::Graph, ::ClusterCompletion)
     return H, cliques
 end
 
-function completion(G::Graph, comp::CliqueTreesCompletion)
+"""
+    completion(G::Graph, comp::ChordalCompletion)
+
+Return a chordal extension of `G` and the corresponding maximal cliques.
+
+The algoritm is Algorithm 3 in [BA10] with the GreedyFillIn heuristic of Table I.
+
+[BA10] Bodlaender, Hans L., and Arie MCA Koster.
+Treewidth computations I. Upper bounds.
+Information and Computation 208, no. 3 (2010): 259-275.
+Utrecht University, Utrecht, The Netherlands www.cs.uu.nl
+"""
+function completion(G::Graph, comp::ChordalCompletion)
     # construct a copy H of the graph G
     H = copy(G)
 
     # construct adjacency matrix of H
-    matrix = spzeros(Bool, num_nodes(H), num_nodes(H))
+    matrix = SparseArrays.spzeros(Bool, num_nodes(H), num_nodes(H))
 
     for j in axes(matrix, 2)
-        matrix.colptr[j + 1] = matrix.colptr[j] + length(neighbors(H, j))
-        append!(rowvals(matrix), neighbors(H, j))
+        list = neighbors(H, j)
+        SparseArrays.getcolptr(matrix)[j + 1] = SparseArrays.getcolptr(matrix)[j] + length(list)
+        append!(SparseArrays.rowvals(matrix), list)
+        append!(SparseArrays.nonzeros(matrix), zeros(Bool, length(list)))
     end
 
     # sort row indices of adjacency matrix
     matrix = copy(transpose(matrix))
 
     # construct tree decomposition of H
-    label, tree = cliquetree(matrix; alg=comp.algo)
+    label, tree = CliqueTrees.cliquetree(matrix; alg=comp.algo)
 
-    # construct chordal extension of H
-    F = FilledGraph(tree)
+    # construct vector of cliques and append fill edges to H
+    cliques = Vector{Vector{Int}}(undef, length(tree))
 
-    # append fill edges to H
-    for v in CliqueTrees.vertices(F)
-        for w in CliqueTrees.neighbors(F, v)
-            push!(neighbors(H, label[v]), label[w])
-            push!(neighbors(H, label[w]), label[v])
-        end
+    for (i, v) in enumerate(tree)
+        clique = cliques[i] = label[v]
+        add_clique!(H, clique)
     end
 
-    # compute maximal cliques of H
-    cliques = [label[clique] for clique in tree]
-
     return H, cliques
 end