Support for Sparse Matrices

Project.toml

@@ -9,12 +9,14 @@ AbstractTrees = "1520ce14-60c1-5f80-bbc7-55ef81b5835c"
 Catlab = "134e5e36-593f-5add-ad60-77f754baafbe"
 CuthillMcKee = "17f17636-5e38-52e3-a803-7ae3aaaf3da9"
 DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8"
+Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 LinearAlgebra = "37e2e46d-f89d-539d-b4ee-838fcccc9c8e"
 LinkedLists = "70f5e60a-1556-5f34-a19e-a48b3e4aaee9"
 MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078"
 Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b"
 PartialFunctions = "570af359-4316-4cb7-8c74-252c00c2016b"
 SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf"
+SymRCM = "286e6d88-80af-4590-acc9-0001b223b9bd"
 TreeWidthSolver = "7d267fc5-9ace-409f-a54c-cd2374872a55"
 
 [compat]

src/Decompositions.jl

@@ -214,7 +214,7 @@ end
 #    jtree    junction tree
 # ----------------------------------------
 function StrDecomp(graph::AbstractSymmetricGraph, jtree::JunctionTree)
-    n = length(jtree)
+    n = treesize(jtree)
     tree = Graph(n)
 
     for i in 1:n - 1
@@ -293,7 +293,7 @@ end
 
 
 function homomorphisms(graph::AbstractSymmetricGraph, jtree::JunctionTree)
-    n = length(jtree)
+    n = treesize(jtree)
     subgraph = Vector{Any}(undef, 2n - 1)
     homomorphism = Vector{Any}(undef, 2n - 2)
 
@@ -310,20 +310,20 @@ function homomorphisms(graph::AbstractSymmetricGraph, jtree::JunctionTree)
     for i in 1:n - 1
         # seperator(i) → clique(parent(i))
         j = parentindex(jtree, i)
-        homomorphism[i] = induced_homomorphism(subgraph[n + i], subgraph[j], seperator_to_parent(jtree, i))
+        homomorphism[i] = induced_homomorphism(subgraph[n + i], subgraph[j], lift_par(jtree, i))
     end
 
     for i in 1:n - 1
         # seperator(i) → clique(i)
-        homomorphism[n + i - 1] = induced_homomorphism(subgraph[n + i], subgraph[i], seperator_to_self(jtree, i))
+        homomorphism[n + i - 1] = induced_homomorphism(subgraph[n + i], subgraph[i], lift_sep(jtree, i))
     end
 
     subgraph, homomorphism
 end
 
 
 function induced_order(order::Order, elements::AbstractVector)
-    Order(sortperm(inverse(order, elements)))
+    Order(sortperm(view(inv(order), elements)))
 end
 
 

src/JunctionTrees.jl

@@ -2,41 +2,51 @@ module JunctionTrees
 
 
 import AMD
+import Catlab.BasicGraphs
 import CuthillMcKee
-import LinkedLists
+import LinkedLists: ListNode, LinkedList
 import Metis
+import SymRCM
 import TreeWidthSolver
 
 using AbstractTrees
-using Catlab.BasicGraphs
+using Base.Order: Ordering
 using DataStructures
+using Graphs
+using Graphs.SimpleGraphs
+using LinearAlgebra
 using SparseArrays
+using SparseArrays: AbstractSparseMatrixCSC
 
 
 # Orders
-export Order, inverse
+export Order
 
 
 # Elimination Algorithms
-export AMDJL_AMD, CuthillMcKeeJL_RCM, MetisJL_ND, TreeWidthSolverJL_BT, MCS
+export AMDJL_AMD, CuthillMcKeeJL_RCM, SymRCMJL_RCM, MetisJL_ND, TreeWidthSolverJL_BT, MCS
+
+
+# Ordered Graphs
+export OrderedGraph
 
 
 # Supernode Types
 export Node, Maximal, Fundamental 
 
 
 # Junction Trees
-export JunctionTree, width, height, seperator, residual, clique, seperator_to_parent, seperator_to_self
+export JunctionTree, treewidth, seperator, residual, clique, find_clique, lift_par, lift_sep, lift
 
 
+include("junction_trees/fixed_stacks.jl")
+include("junction_trees/disjoint_sets.jl")
 include("junction_trees/orders.jl")
 include("junction_trees/elimination_algorithms.jl")
-include("junction_trees/ordered_graphs.jl")
 include("junction_trees/trees.jl")
-include("junction_trees/postorder_trees.jl")
-include("junction_trees/elimination_trees.jl")
+include("junction_trees/child_indices.jl")
+include("junction_trees/ordered_graphs.jl")
 include("junction_trees/supernode_types.jl")
-include("junction_trees/supernode_trees.jl")
 include("junction_trees/junction_trees.jl")
 
 

src/junction_trees/child_indices.jl

@@ -0,0 +1,33 @@
+struct ChildIndices
+    tree::Tree
+    index::Int
+end
+
+
+function Base.iterate(iterator::ChildIndices)
+    iterate(iterator, iterator.tree.child[iterator.index])
+end
+
+
+function Base.iterate(iterator::ChildIndices, i::Integer)
+    if iszero(i)
+        nothing
+    else
+        i, iterator.tree.brother[i]
+    end
+end
+
+
+function Base.IteratorSize(::Type{ChildIndices})
+    Base.SizeUnknown()
+end
+
+
+function Base.eltype(::Type{ChildIndices})
+    Int
+end
+
+
+function AbstractTrees.childindices(tree::Tree, i::Integer)
+    ChildIndices(tree, i)
+end

src/junction_trees/disjoint_sets.jl

@@ -0,0 +1,21 @@
+struct DisjointSets
+    sets::IntDisjointSets{Int}
+    index::Vector{Int}
+    root::Vector{Int}
+
+    function DisjointSets(n::Integer)
+        new(IntDisjointSets(n), collect(1:n), collect(1:n))
+    end
+end
+
+
+function find!(sets::DisjointSets, u::Integer)
+    sets.index[find_root!(sets.sets, u)]
+end
+
+
+function Base.union!(sets::DisjointSets, u::Integer, v::Integer)
+    w = max(u, v)
+    sets.root[w] = root_union!(sets.sets, sets.root[u], sets.root[v])
+    sets.index[sets.root[w]] = w
+end

src/junction_trees/elimination_algorithms.jl

@@ -3,6 +3,7 @@
 
 A graph elimination algorithm. The options are
 - [`CuthillMcKeeJL_RCM`](@ref)
+- [`SymRCMJL_RCM`](@ref)
 - [`AMDJL_AMD`](@ref)
 - [`MetisJL_ND`](@ref)
 - [`TreeWidthSolverJL_BT`](@ref)
@@ -19,6 +20,14 @@
 struct CuthillMcKeeJL_RCM <: EliminationAlgorithm end
 
 
+"""
+    SymRCMJL_RCM <: EliminationAlgorithm
+
+The reverse Cuthill-McKee algorithm. Uses SymRCM.jl.
+"""
+struct SymRCMJL_RCM <: EliminationAlgorithm end
+
+
 """
     AMDJL_AMD <: EliminationAlgorithm
 
@@ -52,96 +61,109 @@
 
 
 """
-    Order(graph::AbstractSymmetricGraph[, ealg::EliminationAlgorithm])
+    Order(graph[, ealg::EliminationAlgorithm])
 
 Construct an elimination order for a simple graph, optionally specifying an elimination algorithm.
 """
-function Order(graph::AbstractSymmetricGraph)
+function Order(graph, ealg::EliminationAlgorithm=DEFAULT_ELIMINATION_ALGORITHM)
+    Order(adjacencymatrix(graph), ealg)
+end
+
+
+# Construct an elimination order.
+function Order(graph::AbstractMatrix)
     Order(graph, DEFAULT_ELIMINATION_ALGORITHM)
 end
 
 
 # Construct an order using the reverse Cuthill-McKee algorithm. Uses CuthillMcKee.jl.
-function Order(graph::AbstractSymmetricGraph, ealg::CuthillMcKeeJL_RCM)
-    order = CuthillMcKee.symrcm(adjacencymatrix(graph))
+function Order(graph::AbstractMatrix, ealg::CuthillMcKeeJL_RCM)
+    order = CuthillMcKee.symrcm(graph)
+    Order(order)
+end
+
+
+function Order(graph::AbstractMatrix, ealg::SymRCMJL_RCM)
+    order = SymRCM.symrcm(graph)
     Order(order)
 end
 
 
 # Construct an order using the approximate minimum degree algorithm. Uses AMD.jl.
-function Order(graph::AbstractSymmetricGraph, ealg::AMDJL_AMD)
-    order = AMD.symamd(adjacencymatrix(graph))
+function Order(graph::AbstractMatrix, ealg::AMDJL_AMD)
+    order = AMD.symamd(graph)
     Order(order)
 end
 
 
 # Construct an order using the nested dissection heuristic. Uses Metis.jl.
-function Order(graph::AbstractSymmetricGraph, ealg::MetisJL_ND)
-    order, index = Metis.permutation(adjacencymatrix(graph))
+function Order(graph::AbstractMatrix, ealg::MetisJL_ND)
+    order, index = Metis.permutation(graph)
     Order(order, index)
 end
 
 
 # Construct an order using the Bouchitte-Todinca algorithm. Uses TreeWidthSolver.jl.
-function Order(graph::AbstractSymmetricGraph, ealg::TreeWidthSolverJL_BT)
-    n = nv(graph)
+function Order(graph::AbstractSparseMatrixCSC, ealg::TreeWidthSolverJL_BT)
+    n = size(graph, 1)
     T = TreeWidthSolver.LongLongUInt{n ÷ 64 + 1}
     fadjlist = Vector{Vector{Int}}(undef, n)
     bitfadjlist = Vector{T}(undef, n)
 
     for i in 1:n
-        fadjlist[i] = sort(collect(outneighbors(graph, i)))
+        fadjlist[i] = rowvals(graph)[nzrange(graph, i)]
         bitfadjlist[i] = TreeWidthSolver.bmask(T, fadjlist[i])
     end
 
     bitgraph = TreeWidthSolver.MaskedBitGraph(bitfadjlist, fadjlist, TreeWidthSolver.bmask(T, 1:n))
     decomposition = TreeWidthSolver.bt_algorithm(bitgraph, TreeWidthSolver.all_pmc_enmu(bitgraph, false), ones(n), false, true)
-    order = reverse(vcat(TreeWidthSolver.EliminationOrder(decomposition.tree).order...))
+    order = reverse(reduce(vcat, TreeWidthSolver.EliminationOrder(decomposition.tree).order))
     Order(order)
 end
 
 
 # Construct an order using the maximum cardinality search algorithm.
-function Order(graph::AbstractSymmetricGraph, ealg::MCS)
-    order, index = mcs(graph)
-    Order(order, index)
+function Order(graph::AbstractMatrix, ealg::MCS)
+    mcs(graph)
 end
 
 
 # Construct the adjacency matrix of a graph.
-function adjacencymatrix(graph::AbstractSymmetricGraph)
-    m = ne(graph)
-    n = nv(graph)
-
-    colptr = ones(Int, n + 1)
-    rowval = sizehint!(Vector{Int}(), 2m)
-
-    for j in 1:n
-        ns = collect(neighbors(graph, j))
-        sort!(ns)
-        colptr[j + 1] = colptr[j] + length(ns)
-        append!(rowval, ns)
+function adjacencymatrix(graph::BasicGraphs.AbstractSymmetricGraph)
+    rowval = Vector{Int}(undef, BasicGraphs.ne(graph))
+    colptr = Vector{Int}(undef, BasicGraphs.nv(graph) + 1)
+    colptr[1] = 1
+
+    for i in BasicGraphs.vertices(graph)
+        column = BasicGraphs.all_neighbors(graph, i)
+        colptr[i + 1] = colptr[i] + length(column)
+        rowval[colptr[i]:colptr[i + 1] - 1] = sort!(collect(column))
     end
 
-    nzval = ones(Int, length(rowval))
-    SparseMatrixCSC(n, n, colptr, rowval, nzval)
+    nzval = ones(Bool, BasicGraphs.ne(graph))
+    SparseMatrixCSC(BasicGraphs.nv(graph), BasicGraphs.nv(graph), colptr, rowval, nzval)
+end
+
+
+# Construct the adjacency matrix of a graph.
+function adjacencymatrix(graph::AbstractGraph)
+    adjacency_matrix(graph; dir=:both)
 end
 
 
 # Simple Linear-Time Algorithms to Test Chordality of Graphs, Test Acyclicity of Hypergraphs, and Selectively Reduce Acyclic Hypergraphs
 # Tarjan and Yannakakis
 # Maximum Cardinality Search
-function mcs(graph::AbstractSymmetricGraph)
-    n = nv(graph)
-    α = Vector{Int}(undef, n)
-    β = Vector{Int}(undef, n)
-    size = Vector{Int}(undef, n)
-    set = Vector{LinkedLists.LinkedList{Int}}(undef, n)
-    pointer = Vector{LinkedLists.ListNode{Int}}(undef, n)
+function mcs(graph::AbstractSparseMatrixCSC)
+    n = size(graph, 1)
+    α = Order(undef, n)
+    len = Vector{Int}(undef, n)
+    set = Vector{LinkedList{Int}}(undef, n)
+    pointer = Vector{ListNode{Int}}(undef, n)
 
     for i in 1:n
-        size[i] = 1
-        set[i] = LinkedLists.LinkedList{Int}()
+        len[i] = 1
+        set[i] = LinkedList{Int}()
         pointer[i] = push!(set[1], i)
     end
 
@@ -151,15 +173,14 @@
     while i >= 1
         v = first(set[j])
         deleteat!(set[j], pointer[v])        
-        α[v] = i
-        β[i] = v
-        size[v] = 0
-
-        for w in neighbors(graph, v)
-            if size[w] >= 1
-                deleteat!(set[size[w]], pointer[w])
-                size[w] += 1
-                pointer[w] = push!(set[size[w]], w)
+        α[i] = v
+        len[v] = 0
+
+        for w in view(rowvals(graph), nzrange(graph, v))
+            if len[w] >= 1
+                deleteat!(set[len[w]], pointer[w])
+                len[w] += 1
+                pointer[w] = push!(set[len[w]], w)
             end
         end
 
@@ -171,7 +192,7 @@
         end
     end
 
-    β, α
+    α
 end