internalized CliqueTrees

Project.toml

@@ -5,7 +5,6 @@ repo = "https://github.com/lanl-ansi/PowerModels.jl"
 version = "0.21.3"
 
 [deps]
-CliqueTrees = "60701a23-6482-424a-84db-faee86b9b1f8"
 InfrastructureModels = "2030c09a-7f63-5d83-885d-db604e0e9cc0"
 JSON = "682c06a0-de6a-54ab-a142-c8b1cf79cde6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
@@ -15,7 +14,6 @@ NLsolve = "2774e3e8-f4cf-5e23-947b-6d7e65073b56"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
 
 [compat]
-CliqueTrees = "0.4.0"
 HiGHS = "~1"
 InfrastructureModels = "~0.6, ~0.7"
 Ipopt = "~1"

src/PowerModels.jl

@@ -1,7 +1,5 @@
 module PowerModels
 
-import CliqueTrees
-
 import LinearAlgebra, SparseArrays
 
 import JSON
@@ -62,6 +60,7 @@ include("io/json.jl")
 
 include("form/iv.jl")
 
+include("form/cliquetrees/CliqueTrees.jl")
 include("form/acp.jl")
 include("form/acr.jl")
 include("form/act.jl")

src/form/cliquetrees/CliqueTrees.jl

@@ -0,0 +1,16 @@
+module CliqueTrees
+
+using Base: oneto, @propagate_inbounds
+using SparseArrays
+
+const AbstractScalar{T} = AbstractArray{T,0}
+const Scalar{T} = Array{T,0}
+const MAX_ITEMS_PRINTED = 5
+
+export permutation, MCS
+
+include("./abstract_linked_lists.jl")
+include("./doubly_linked_lists.jl")
+include("./elimination_algorithms.jl")
+
+end

src/form/cliquetrees/LICENSE

@@ -0,0 +1,21 @@
+MIT License
+
+Copyright (c) 2023 AlgebraicJulia
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.

src/form/cliquetrees/abstract_linked_lists.jl

@@ -0,0 +1,40 @@
+abstract type AbstractLinkedList{I<:Integer} end
+
+function Base.show(io::IO, ::MIME"text/plain", list::L) where {L<:AbstractLinkedList}
+    println(io, "$L:")
+
+    for (i, v) in enumerate(take(list, MAX_ITEMS_PRINTED + 1))
+        if i <= MAX_ITEMS_PRINTED
+            println(io, " $v")
+        else
+            println(io, " ⋮")
+        end
+    end
+end
+
+function Base.empty!(list::AbstractLinkedList{I}) where {I}
+    list.head[] = zero(I)
+    return list
+end
+
+function Base.isempty(list::AbstractLinkedList)
+    return iszero(list.head[])
+end
+
+#######################
+# Iteration Interface #
+#######################
+
+function Base.iterate(list::AbstractLinkedList{I}, i::I=list.head[]) where {I}
+    if !iszero(i)
+        @inbounds (i, list.next[i])
+    end
+end
+
+function Base.IteratorSize(::Type{<:AbstractLinkedList})
+    return Base.SizeUnknown()
+end
+
+function Base.eltype(::Type{<:AbstractLinkedList{I}}) where {I}
+    return I
+end

src/form/cliquetrees/doubly_linked_lists.jl

@@ -0,0 +1,91 @@
+# A doubly linked list of distinct natural numbers.
+struct DoublyLinkedList{
+    I,Init<:AbstractScalar{I},Prev<:AbstractVector{I},Next<:AbstractVector{I}
+} <: AbstractLinkedList{I}
+    head::Init
+    prev::Prev
+    next::Next
+end
+
+function DoublyLinkedList{I}(n::Integer) where {I}
+    head = fill(zero(I))
+    prev = Vector{I}(undef, n)
+    next = Vector{I}(undef, n)
+    return DoublyLinkedList(head, prev, next)
+end
+
+function DoublyLinkedList{I}(vector::AbstractVector) where {I}
+    list = DoublyLinkedList{I}(length(vector))
+    return prepend!(list, vector)
+end
+
+function DoublyLinkedList(vector::AbstractVector{I}) where {I}
+    return DoublyLinkedList{I}(vector)
+end
+
+@propagate_inbounds function Base.pushfirst!(list::DoublyLinkedList, i::Integer)
+    @boundscheck checkbounds(list.prev, i)
+    @boundscheck checkbounds(list.next, i)
+
+    @inbounds begin
+        n = list.next[i] = list.head[]
+        list.head[] = i
+        list.prev[i] = 0
+
+        if !iszero(n)
+            list.prev[n] = i
+        end
+    end
+
+    return list
+end
+
+@propagate_inbounds function Base.popfirst!(list::DoublyLinkedList)
+    i = list.head[]
+    @boundscheck checkbounds(list.next, i)
+
+    @inbounds begin
+        n = list.head[] = list.next[i]
+
+        if !iszero(n)
+            list.prev[n] = 0
+        end
+    end
+
+    return i
+end
+
+@propagate_inbounds function Base.delete!(list::DoublyLinkedList, i::Integer)
+    @boundscheck checkbounds(list.prev, i)
+    @boundscheck checkbounds(list.next, i)
+
+    @inbounds begin
+        p = list.prev[i]
+        n = list.next[i]
+
+        if !iszero(p)
+            list.next[p] = n
+        else
+            list.head[] = n
+        end
+
+        if !iszero(n)
+            list.prev[n] = p
+        end
+    end
+
+    return list
+end
+
+function Base.prepend!(list::DoublyLinkedList, vector::AbstractVector)
+    @views list.next[vector[begin:(end - 1)]] = vector[(begin + 1):end]
+    @views list.prev[vector[(begin + 1):end]] = vector[begin:(end - 1)]
+
+    if !isempty(vector)
+        list.next[vector[end]] = list.head[]
+        list.prev[vector[begin]] = 0
+        list.head[] = vector[begin]
+    end
+
+    return list
+end

src/form/cliquetrees/elimination_algorithms.jl

@@ -0,0 +1,96 @@
+"""
+    EliminationAlgorithm
+
+A graph elimination algorithm. The options are
+
+| type             | name                                         | time     | space    |
+|:---------------- |:-------------------------------------------- |:-------- |:-------- |
+| [`MCS`](@ref)    | maximum cardinality search                   | O(m + n) | O(n)     |
+"""
+abstract type EliminationAlgorithm end
+
+"""
+    MCS <: EliminationAlgorithm
+
+    MCS()
+
+The maximum cardinality search algorithm.
+
+### Reference
+
+Tarjan, Robert E., and Mihalis Yannakakis. "Simple linear-time algorithms to test chordality of graphs, test acyclicity of hypergraphs, and selectively reduce acyclic hypergraphs." *SIAM Journal on Computing* 13.3 (1984): 566-579.
+"""
+struct MCS <: EliminationAlgorithm end
+
+function permutation(matrix, alg::MCS)
+    index, card = mcs(matrix)
+    return invperm(index), index
+end
+
+"""
+    mcs(matrix[, clique::AbstractVector])
+
+Perform a maximum cardinality search, optionally specifying a clique to be ordered last.
+Returns the inverse permutation.
+"""
+function mcs(matrix::SparseMatrixCSC{<:Any, V}) where V
+    mcs(matrix, oneto(zero(V)))
+end
+
+# Simple Linear-Time Algorithms to Test Chordality of BipartiteGraphs, Test Acyclicity of Hypergraphs, and Selectively Reduce Acyclic Hypergraphs
+# Tarjan and Yannakakis
+# Maximum Cardinality Search
+#
+# Construct a fill-reducing permutation of a graph.
+# The complexity is O(m + n), where m = |E| and n = |V|.
+function mcs(matrix::SparseMatrixCSC{<:Any, V}, clique::AbstractVector) where {V}
+    n = convert(V, size(matrix, 2))
+
+    # construct disjoint sets data structure
+    head = zeros(V, n + one(V))
+    prev = Vector{V}(undef, n + one(V))
+    next = Vector{V}(undef, n + one(V))
+
+    function set(i)
+        @inbounds DoublyLinkedList(view(head, i), prev, next)
+    end
+
+    # run algorithm
+    alpha = Vector{V}(undef, n)
+    card = ones(V, n)
+    prepend!(set(one(V)), oneto(n))
+
+    j = one(V)
+    k = convert(V, lastindex(clique))
+
+    @inbounds for i in reverse(oneto(n))
+        v = zero(V)
+
+        if k in eachindex(clique)
+            v = convert(V, clique[k])
+            k -= one(V)
+            delete!(set(j), v)
+        else
+            v = popfirst!(set(j))
+        end
+
+        alpha[v] = i
+        card[v] = one(V) - card[v]
+
+        for w in @view rowvals(matrix)[nzrange(matrix, v)]
+            if card[w] >= one(V)
+                delete!(set(card[w]), w)
+                card[w] += one(V)
+                pushfirst!(set(card[w]), w)
+            end
+        end
+
+        j += one(V)
+
+        while j >= one(V) && isempty(set(j))
+            j -= one(V)
+        end
+    end
+
+    return alpha, card
+end

src/form/wrm.jl

@@ -1,4 +1,5 @@
 ### sdp relaxations in the rectangular W-space
+import .CliqueTrees
 import LinearAlgebra: Hermitian, cholesky, Symmetric, diag, I
 import SparseArrays: SparseMatrixCSC, sparse, spdiagm, findnz, spzeros, nonzeros