From fad62b297d3f846a933b086fb2ff7ed6ecc900ee Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sat, 16 Nov 2024 21:14:40 -0500 Subject: [PATCH 01/48] temp commit... --- src/StructuredDecompositions.jl | 1 - src/junction_trees/JunctionTrees.jl | 66 +--- src/junction_trees/abstract_trees.jl | 25 ++ src/junction_trees/elimination_algorithms.jl | 47 ++- src/junction_trees/elimination_trees.jl | 388 +++++++------------ src/junction_trees/junction_trees.jl | 189 --------- src/junction_trees/ordered_graphs.jl | 143 +++++++ src/junction_trees/orders.jl | 156 ++------ src/junction_trees/postorder_trees.jl | 104 +++++ src/junction_trees/supernode_trees.jl | 205 ++++++++++ src/junction_trees/supernode_types.jl | 70 ++++ src/junction_trees/supernodes.jl | 36 -- src/junction_trees/trees.jl | 205 +--------- src/nested_uwds/NestedUWDs.jl | 34 -- src/nested_uwds/nested_uwds.jl | 312 --------------- test/JunctionTrees.jl | 27 +- test/NestedUWDs.jl | 122 ------ test/runtests.jl | 4 - 18 files changed, 815 insertions(+), 1319 deletions(-) create mode 100644 src/junction_trees/abstract_trees.jl delete mode 100644 src/junction_trees/junction_trees.jl create mode 100644 src/junction_trees/ordered_graphs.jl create mode 100644 src/junction_trees/postorder_trees.jl create mode 100644 src/junction_trees/supernode_trees.jl create mode 100644 src/junction_trees/supernode_types.jl delete mode 100644 src/junction_trees/supernodes.jl delete mode 100644 src/nested_uwds/NestedUWDs.jl delete mode 100644 src/nested_uwds/nested_uwds.jl delete mode 100644 test/NestedUWDs.jl diff --git a/src/StructuredDecompositions.jl b/src/StructuredDecompositions.jl index 8c90cc6..2e82fef 100644 --- a/src/StructuredDecompositions.jl +++ b/src/StructuredDecompositions.jl @@ -5,7 +5,6 @@ include("Decompositions.jl") include("FunctorUtils.jl") include("DecidingSheaves.jl") include("junction_trees/JunctionTrees.jl") -include("nested_uwds/NestedUWDs.jl") end diff --git a/src/junction_trees/JunctionTrees.jl b/src/junction_trees/JunctionTrees.jl index 87ae3d3..8fbfe33 100644 --- a/src/junction_trees/JunctionTrees.jl +++ b/src/junction_trees/JunctionTrees.jl @@ -10,74 +10,32 @@ using Catlab.BasicGraphs using DataStructures using SparseArrays -# Elimination Algorithms -export EliminationAlgorithm, AMDJL_AMD, CuthillMcKeeJL_RCM, MetisJL_ND, MCS - -# Supernodes -export Supernode, Node, MaximalSupernode, FundamentalSupernode # Orders export Order -# Elimination Trees -export EliminationTree -export getwidth, getsupernode, getsubtree, getlevel - -# Junction Trees -export JunctionTree -export getseperator, getresidual - - -# Add an element x to a sorted set v. -# Returns true if x ∉ v. -# Returns false if x ∈ v. -function insertsorted!(v::Vector, x::Integer) - i = searchsortedfirst(v, x) - - if i > length(v) || v[i] != x - insert!(v, i, x) - true - else - false - end -end - -# Delete an element x from a sorted set v. -# Returns true if x ∈ v. -# Returns false if x ∉ v. -function deletesorted!(v::Vector, x::Integer) - i = searchsortedfirst(v, x) +# Elimination Algorithms +export AMDJL_AMD, CuthillMcKeeJL_RCM, MetisJL_ND - if i <= length(v) && v[i] == x - deleteat!(v, i) - true - else - false - end -end +# Supernode Types +export Node, Maximal, Fundamental -# Delete the elements xs from a sorted set v. -# Returns true if xs and v intersect. -# Returns false if xs and v are disjoint. -function deletesorted!(v::Vector, xs::AbstractVector) - isintersecting = true - for x in xs - isintersecting = deletesorted!(v, x) || isintersecting - end +# Supernode Trees +export SupernodeTree, width, height - isintersecting -end - -include("elimination_algorithms.jl") -include("supernodes.jl") include("orders.jl") +include("elimination_algorithms.jl") +include("ordered_graphs.jl") +include("abstract_trees.jl") include("trees.jl") +include("postorder_trees.jl") include("elimination_trees.jl") -include("junction_trees.jl") +include("supernode_types.jl") +include("supernode_trees.jl") end diff --git a/src/junction_trees/abstract_trees.jl b/src/junction_trees/abstract_trees.jl new file mode 100644 index 0000000..618e23d --- /dev/null +++ b/src/junction_trees/abstract_trees.jl @@ -0,0 +1,25 @@ +# A rooted tree. +abstract type AbstractTree end + + +# Compute a postordering of tree's vertices. +function postorder(tree::AbstractTree) + n = length(tree) + order = Vector{Int}(undef, n) + index = Vector{Int}(undef, n) + + for node in PreOrderDFS(IndexNode(tree)) + order[n] = node.index + index[node.index] = n + n -= 1 + end + + Order(order, index) +end + + +# Get the height of a tree. +function height(tree::AbstractTree) + n = length(tree) + maximum(map(i -> level(tree, i), 1:n)) +end diff --git a/src/junction_trees/elimination_algorithms.jl b/src/junction_trees/elimination_algorithms.jl index d302a74..f87dc5b 100644 --- a/src/junction_trees/elimination_algorithms.jl +++ b/src/junction_trees/elimination_algorithms.jl @@ -5,7 +5,7 @@ A graph elimination algorithm. The options are - [`CuthillMcKeeJL_RCM`](@ref) - [`AMDJL_AMD`](@ref) - [`MetisJL_ND`](@ref) -- [`MCS`](@ref) +- [`Order`](@ref) """ abstract type EliminationAlgorithm end @@ -34,12 +34,47 @@ The nested dissection heuristic. Uses Metis.jl. struct MetisJL_ND <: EliminationAlgorithm end -""" - MCS <: EliminationAlgorithm +# Construct an order using the reverse Cuthill-McKee algorithm. Uses +# CuthillMcKee.jl. +function Order(graph::AbstractSymmetricGraph, ::CuthillMcKeeJL_RCM) + order = CuthillMcKee.symrcm(adjacencymatrix(graph)) + Order(order) +end -The maximum cardinality search algorithm. -""" -struct MCS <: EliminationAlgorithm end + +# Construct an order using the approximate minimum degree algorithm. Uses +# AMD.jl. +function Order(graph::AbstractSymmetricGraph, ::AMDJL_AMD) + order = AMD.symamd(adjacencymatrix(graph)) + Order(order) +end + + +# Construct an order using the nested dissection heuristic. Uses Metis.jl. +function Order(graph::AbstractSymmetricGraph, ::MetisJL_ND) + order, index = Metis.permutation(adjacencymatrix(graph)) + Order(order, index) +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) + end + + nzval = ones(Int, length(rowval)) + SparseMatrixCSC(n, n, colptr, rowval, nzval) +end const DEFAULT_ELIMINATION_ALGORITHM = AMDJL_AMD() diff --git a/src/junction_trees/elimination_trees.jl b/src/junction_trees/elimination_trees.jl index c5d57d9..0102eff 100644 --- a/src/junction_trees/elimination_trees.jl +++ b/src/junction_trees/elimination_trees.jl @@ -1,268 +1,117 @@ -# A supernodal elimination tree. -struct EliminationTree - order::Order - tree::Tree - firstsupernodelist::Vector{Int} - lastsupernodelist::Vector{Int} - subtreelist::Vector{Int} - width::Int +# An ordered graph (G, σ) equipped with the elimination tree T of its elimination graph. +# Nodes i in T correspond to vertices σ(i) in G. +struct EliminationTree{T <: AbstractTree} <: AbstractTree + tree::T + ograph::OrderedGraph end +# Construct an elimination tree using an elimination algorithm. function EliminationTree( - order::Order, - tree::Tree, - supernodelist::AbstractVector, - subtreelist::AbstractVector, - width::Integer) - - n = length(order) - m = length(tree) - postorder = Order(n) - firstsupernodelist = Vector{Int}(undef, m) - lastsupernodelist = Vector{Int}(undef, m) - - i₂ = 0 - - for j in 1:m - supernode = supernodelist[j] - i₁ = i₂ + 1 - i₂ = i₂ + length(supernode) - firstsupernodelist[j] = i₁ - lastsupernodelist[j] = i₂ - postorder[i₁:i₂] .= supernode - end - - order = compose(postorder, order) - subtreelist = subtreelist[postorder] - - EliminationTree( - order, - tree, - firstsupernodelist, - lastsupernodelist, - subtreelist, - width) -end - - -# Construct a supernodal elimination tree. -# -# The complexity is -# 𝒪(m α(m, n) + n) -# where m = |E|, n = |V|, and α is the inverse Ackermann function. -function EliminationTree( - graph::AbstractSymmetricGraph, - order::Order, - supernode::Supernode=DEFAULT_SUPERNODE) - - etree = Tree(graph, order) - _, outdegreelist = getdegrees(graph, order, etree) - - supernodelist, subtreelist, parentlist = makestree( - etree, - outdegreelist, - supernode) - - n = nv(graph) - tree = Tree(subtreelist[n], parentlist) - postorder, tree = makepostorder(tree) - - supernodelist = supernodelist[postorder] - subtreelist = postorder.index[subtreelist] - width = maximum(outdegreelist) - - EliminationTree( - order, - tree, - supernodelist, - subtreelist, - width) -end - - -# Construct a supernodal elimination tree, first computing an elimination order. -function EliminationTree( - graph::AbstractSymmetricGraph, - algorithm::EliminationAlgorithm=DEFAULT_ELIMINATION_ALGORITHM, - supernode::Supernode=DEFAULT_SUPERNODE) - - order = Order(graph, algorithm) - EliminationTree(graph, order, supernode) -end - - -# Get the number of nodes in a supernodal elimination tree. -function Base.length(stree::EliminationTree) - length(stree.tree) -end - - -# Get the width of a supernodal elimination tree. -function getwidth(stree::EliminationTree) - stree.width + graph::AbstractSymmetricGraph, + ealg::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM) + EliminationTree(OrderedGraph(graph, ealg)) end -# Get the supernode at node i. -function getsupernode(stree::EliminationTree, i::Integer) - i₁ = stree.firstsupernodelist[i] - i₂ = stree.lastsupernodelist[i] - stree.order[i₁:i₂] +# Construct the elimination tree of an ordered graph. +function EliminationTree(ograph::OrderedGraph) + EliminationTree(Tree(etree(ograph)), ograph) end -# Get the highest node containing a vertex v. -function getsubtree(stree::EliminationTree, v::Integer) - stree.subtreelist[stree.order.index[v]] +# Postorder an elimination tree. +function EliminationTree{PostorderTree}(etree::EliminationTree, order::Order) + EliminationTree(PostorderTree(etree.tree, order), OrderedGraph(etree.ograph, order)) end -# Get the highest node containing vertices vs. -function getsubtree(stree::EliminationTree, vs::AbstractVector) - init = length(stree.order) - stree.subtreelist[minimum(stree.order.index[vs]; init)] -end - - -# Get the level of node i. -function getlevel(stree::EliminationTree, i::Integer) - getlevel(stree.tree, i) -end - - -# Evaluate whether node i₁ is a descendant of node i₂. -function AbstractTrees.isdescendant(stree::EliminationTree, i₁::Integer, i₂::Integer) - isdescendant(stree.tree, i₁, i₂) -end +# A Compact Row Storage Scheme for Cholesky Factors Using Elimination Trees +# Liu +# Algorithm 4.2: Elimination Tree by Path Compression. +function etree(ograph::OrderedGraph) + n = nv(graph) + parent = collect(1:n) + ancestor = collect(1:n) + for i in 1:n + for k in inneighbors(ograph, i) + r = k -# Compute the supernodes, parent function, and first ancestor of a -# supernodal elimination tree. -# -# The complexity is -# 𝒪(n) -# where n = |V|. -# -# doi:10.1561/2400000006 -# Algorithm 4.1: Maximal supernodes and supernodal elimination tree. -function makestree(etree::Tree, outdegrees::AbstractVector, supernode::Supernode) - n = length(etree) - sbt = Vector{Int}(undef, n) - snd = Vector{Int}[] - q = Int[] - a = Int[] - - for v in 1:n - w′ = findchild(etree, outdegrees, v, supernode) - - if isnothing(w′) - i = length(snd) + 1 - sbt[v] = i - push!(snd, [v]) - push!(q, 0) - push!(a, 0) - else - i = sbt[w′] - sbt[v] = i - push!(snd[i], v) - end + while ancestor[r] != r && ancestor[r] != i + t = ancestor[r] + ancestor[r] = i + r = t + end - for w in childindices(etree, v) - if w !== w′ - j = sbt[w] - q[j] = i - a[j] = v + if ancestor[r] == r + ancestor[r] = i + parent[r] = i end end end - snd, sbt, q, a + parent end -# Find a child w of v such that -# v ∈ snd(w). -# If no such child exists, return nothing. -function findchild(etree::Tree, outdegrees::AbstractVector, v::Integer, ::Supernode) end - - -function findchild(etree::Tree, outdegrees::AbstractVector, v::Integer, ::MaximalSupernode) - for w in childindices(etree, v) - if outdegrees[w] == outdegrees[v] + 1 - return w - end - end -end - - -function findchild(etree::Tree, outdegrees::AbstractVector, v::Integer, ::FundamentalSupernode) - ws = childindices(etree, v) - - if length(ws) == 1 - w = only(ws) - - if outdegrees[w] == outdegrees[v] + 1 - return w - end - end +# An Efficient Algorithm to Compute Row and Column Counts for Sparse Cholesky Factorization +# Gilbert, Ng, and Peyton +# Figure 3: Implementation of algorithm to compute row and column counts. +function supcnt(etree::EliminationTree) + order = postorder(etree) + index = inverse(order) + rc, cc = supcnt(EliminationTree{PostorderTree}(etree, order)) + rc[index], cc[index] end -# Compute the row and column counts of a graph's elimination graph. -# -# The complexity is -# 𝒪(m α(m, n)) -# where m = |E|, n = |V|, and α is the inverse Ackermann function. -# -# doi:10.1137/S089547989223692 -# Figure 3: Implementation of algorithm to compute row and column counts -function getdegrees(graph::AbstractSymmetricGraph, order::Order, etree::Tree) - n = nv(graph) - forest = IntDisjointSets(n) - rvert = Vector{Int}(undef, n) - index = Vector{Int}(undef, n) - rvert .= index .= 1:n - - function FIND(p) - index[find_root!(forest, p)] +# An Efficient Algorithm to Compute Row and Column Counts for Sparse Cholesky Factorization +# Gilbert, Ng, and Peyton +# Figure 3: Implementation of algorithm to compute row and column counts. +function supcnt(etree::EliminationTree{PostorderTree}) + n = length(etree) + + #### Disjoint Set Union #### + + rvert = collect(1:n) + index = collect(1:n) + forest = IntDisjointSets(n) + + function find(u) + index[find_root!(forest, u)] end - - function UNION(u, v) + + function union(u, v) w = max(u, v) rvert[w] = root_union!(forest, rvert[u], rvert[v]) index[rvert[w]] = w end - postorder, etree = makepostorder(etree) - graph = Graph(graph, compose(postorder, order)) - prev_p = Vector{Int}(undef, n) - prev_nbr = Vector{Int}(undef, n) - rc = Vector{Int}(undef, n) - wt = Vector{Int}(undef, n) - - for u in 1:n - prev_p[u] = 0 - prev_nbr[u] = 0 - rc[u] = 1 - wt[u] = isempty(childindices(etree, u)) - end + ############################ + + prev_p = zeros(Int, n) + prev_nbr = zeros(Int, n) + rc = ones(Int, n) + wt = ones(Int, n) - for p in 1:n - if p != n - wt[parentindex(etree, p)] -= 1 - end + for u in 1:n - 1 + wt[parentindex(etree, u)] = 0 + end + + for p in 1:n - 1 + wt[parentindex(etree, p)] -= 1 - for u in neighbors(graph, p) - if getfirstdescendant(etree, p) > prev_nbr[u] + for u in outneighbors(etree, p) + if firstdescendant(etree, p) > prev_nbr[u] wt[p] += 1 - p′ = prev_p[u] + pp = prev_p[u] - if p′ == 0 - rc[u] += getlevel(etree, p) - getlevel(etree, u) + if iszero(pp) + rc[u] += level(etree, p) - level(etree, u) else - q = FIND(p′) - rc[u] += getlevel(etree, p) - getlevel(etree, q) + q = find(pp) + rc[u] += level(etree, p) - level(etree, q) wt[q] -= 1 end @@ -272,9 +121,7 @@ function getdegrees(graph::AbstractSymmetricGraph, order::Order, etree::Tree) prev_nbr[u] = p end - if p != n - UNION(p, parentindex(etree, p)) - end + union(p, parentindex(etree, p)) end cc = wt @@ -283,9 +130,51 @@ function getdegrees(graph::AbstractSymmetricGraph, order::Order, etree::Tree) cc[parentindex(etree, v)] += cc[v] end - indegrees = rc[postorder.index] .- 1 - outdegrees = cc[postorder.index] .- 1 - indegrees, outdegrees + rc, cc +end + + +# Compute higher degree of every vertex in the elimination graph of +# (G, σ). +function outdegrees(etree::EliminationTree) + rc, cc = supcnt(etree) + cc .- 1 +end + + +# Get the number of nodes in T. +function Base.length(etree::EliminationTree) + length(etree.tree) +end + + +# Get the level of a node i. +function level(etree::EliminationTree, i::Integer) + level(etree.tree, i) +end + + +# Get the first descendant of a node i. +function firstdescendant(etree::EliminationTree, i::Integer) + firstdescendant(etree.tree, i) +end + + +# Determine whether a vertex i is a descendant of a node j. +function isdescendant(etree::EliminationTree, i::Integer, j::Integer) + isdescendant(etree.tree, i, j) +end + + +# Get the vertex σ(i). +function order(etree::EliminationTree, i) + order(etree.ograph, i) +end + + +# Get the index σ⁻¹(v), +function inverse(etree::EliminationTree, i) + inverse(etree.ograph, i) end @@ -294,18 +183,18 @@ end ########################## -function AbstractTrees.rootindex(stree::EliminationTree) - rootindex(stree.tree) +function AbstractTrees.rootindex(etree::EliminationTree) + rootindex(etree.tree) end -function AbstractTrees.parentindex(stree::EliminationTree, i::Integer) - parentindex(stree.tree, i) +function AbstractTrees.parentindex(etree::EliminationTree, i::Integer) + parentindex(etree.tree, i) end -function AbstractTrees.childindices(stree::EliminationTree, i::Integer) - childindices(stree.tree, i) +function AbstractTrees.childindices(etree::EliminationTree, i::Integer) + childindices(etree.tree, i) end @@ -314,6 +203,21 @@ function AbstractTrees.NodeType(::Type{IndexNode{EliminationTree, Int}}) end -function AbstractTrees.nodetype(::Type{IndexNode{EliminationTree, Int}}) - IndexNode{EliminationTree, Int} +function AbstractTrees.nodetype(::Type{IndexNode{EliminationTree{T}, Int}}) where T + IndexNode{EliminationTree{T}, Int} +end + + +############################ +# Abstract Graph Interface # +############################ + + +function BasicGraphs.inneighbors(etree::EliminationTree, i::Integer) + inneighbors(etree.ograph, i) +end + + +function BasicGraphs.outneighbors(etree::EliminationTree, i::Integer) + outneighbors(etree.ograph, i) end diff --git a/src/junction_trees/junction_trees.jl b/src/junction_trees/junction_trees.jl deleted file mode 100644 index 87e22e0..0000000 --- a/src/junction_trees/junction_trees.jl +++ /dev/null @@ -1,189 +0,0 @@ -# A junction tree. -struct JunctionTree - stree::EliminationTree - seperatorlist::Vector{Vector{Int}} -end - - -# Construct a tree decomposition. -function JunctionTree(graph::AbstractSymmetricGraph, stree::EliminationTree) - graph = makeeliminationgraph(graph, stree) - - n = length(stree) - seperatorlist = Vector{Vector{Int}}(undef, n) - seperatorlist[n] = [] - - for i in 1:n - 1 - v₁ = stree.firstsupernodelist[i] - v₂ = stree.lastsupernodelist[i] - bag = collect(neighbors(graph, v₁)) - sort!(bag) - seperatorlist[i] = bag[v₂ - v₁ + 1:end] - end - - JunctionTree(stree, seperatorlist) -end - - -# Reorient a juncton tree towards the given root. -function JunctionTree(root::Integer, jtree::JunctionTree) - m = length(jtree.stree.order) - n = length(jtree) - seperatorlist = Vector{Vector{Int}}(undef, n) - supernodelist = Vector{Vector{Int}}(undef, n) - subtreelist = Vector{Int}(undef, m) - - v₁ = jtree.stree.firstsupernodelist[root] - v₂ = jtree.stree.lastsupernodelist[root] - seperatorlist[n] = [] - supernodelist[n] = [v₁:v₂; jtree.seperatorlist[root]] - subtreelist[supernodelist[n]] .= n - - tree = Tree(root, jtree.stree.tree) - postorder, tree = makepostorder(tree) - - for i in 1:n - 1 - j = postorder[i] - v₁ = jtree.stree.firstsupernodelist[j] - v₂ = jtree.stree.lastsupernodelist[j] - - if isdescendant(jtree, root, j) - seperatorlist[i] = jtree.seperatorlist[postorder[parentindex(tree, i)]] - supernodelist[i] = [v₁:v₂; jtree.seperatorlist[j]] - deletesorted!(supernodelist[i], seperatorlist[i]) - else - seperatorlist[i] = jtree.seperatorlist[j] - supernodelist[i] = v₁:v₂ - end - - subtreelist[supernodelist[i]] .= i - end - - order = jtree.stree.order - width = jtree.stree.width - stree = EliminationTree(order, tree, supernodelist, subtreelist, width) - - for i in 1:n - seperatorlist[i] = stree.order.index[order[seperatorlist[i]]] - sort!(seperatorlist[i]) - end - - JunctionTree(stree, seperatorlist) -end - - -# Construct a tree decomposition, first computing an elimination order and a supernodal -# elimination tree. -function JunctionTree( - graph::AbstractSymmetricGraph, - algorithm::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, - supernode::Supernode=DEFAULT_SUPERNODE) - - stree = EliminationTree(graph, algorithm, supernode) - JunctionTree(graph, stree) -end - - -# Get the number of nodes in a junction tree. -function Base.length(jtree::JunctionTree) - length(jtree.stree) -end - - -# Get the width of a junction tree. -function getwidth(jtree::JunctionTree) - getwidth(jtree.stree) -end - - -# Get the seperator at node i. -function getseperator(jtree::JunctionTree, i::Integer) - jtree.stree.order[jtree.seperatorlist[i]] -end - - -# Get the residual at node i. -function getresidual(jtree::JunctionTree, i::Integer) - getsupernode(jtree.stree, i) -end - - -# Get the highest node containing the vertex v. -function getsubtree(jtree::JunctionTree, v::Union{Integer, AbstractVector}) - getsubtree(jtree.stree, v) -end - - -# Get the level of node i. -function getlevel(jtree::JunctionTree, i::Integer) - getlevel(jtree.stree, i) -end - - -# Evaluate whether node i₁ is a descendant of node i₂. -function AbstractTrees.isdescendant(jtree::JunctionTree, i₁::Integer, i₂::Integer) - isdescendant(jtree.stree, i₁, i₂) -end - - -# Construct an elimination graph. -function makeeliminationgraph(graph::AbstractSymmetricGraph, stree::EliminationTree) - n = length(stree) - graph = Graph(graph, stree.order) - - for i in 1:n - 1 - u₁ = stree.firstsupernodelist[i] - u₂ = stree.lastsupernodelist[i] - - for u in u₁:u₂ - 1 - v = u + 1 - - for w in neighbors(graph, u) - if v != w && !has_edge(graph, v, w) - add_edge!(graph, v, w) - end - end - end - - u = u₂ - v = stree.firstsupernodelist[parentindex(stree, i)] - - for w in neighbors(graph, u) - if v != w && !has_edge(graph, v, w) - add_edge!(graph, v, w) - end - end - end - - graph -end - - -########################## -# Indexed Tree Interface # -########################## - - -function AbstractTrees.rootindex(jtree::JunctionTree) - rootindex(jtree.stree) -end - - -function AbstractTrees.parentindex(jtree::JunctionTree, i::Integer) - parentindex(jtree.stree, i) -end - - -function AbstractTrees.childindices(jtree::JunctionTree, i::Integer) - childindices(jtree.stree, i) -end - - -function AbstractTrees.NodeType(::Type{IndexNode{JunctionTree, Int}}) - HasNodeType() -end - - -function AbstractTrees.nodetype(::Type{IndexNode{JunctionTree, Int}}) - IndexNode{JunctionTree, Int} -end diff --git a/src/junction_trees/ordered_graphs.jl b/src/junction_trees/ordered_graphs.jl new file mode 100644 index 0000000..1915e69 --- /dev/null +++ b/src/junction_trees/ordered_graphs.jl @@ -0,0 +1,143 @@ +# An ordered graph (G, σ). +struct OrderedGraph + graph::Graph + order::Order +end + + +# Given a graph G and permutation σ, construct the ordered graph +# (G, σ). +function OrderedGraph(sgraph::AbstractSymmetricGraph, order::Order) + n = nv(sgraph) + graph = Graph(n) + + for e in edges(sgraph) + u = src(sgraph, e) + v = tgt(sgraph, e) + + if order(u, v) + add_edge!(graph, inverse(order, u), inverse(order, v)) + end + end + + OrderedGraph(graph, order) +end + + +# Given an ordered graph (G, σ) and permutation μ, construct the ordered graph +# (G, σ ∘ μ). +function OrderedGraph(ograph::OrderedGraph, order::Order) + n = nv(ograph) + graph = Graph(n) + + for e in edges(ograph) + u = src(ograph, e) + v = tgt(ograph, e) + + if order(u, v) + add_edge!(graph, inverse(order, u), inverse(order, v)) + else + add_edge!(graph, inverse(order, v), inverse(order, u)) + end + end + + OrderedGraph(graph, compose(order, ograph.order)) +end + + +# Construct the elimination graph of an ordered graph. +function eliminationgraph(ograph::OrderedGraph) + ograph = deepcopy(ograph) + + for i in vertices(ograph) + ns = collect(outneighbors(ograph, i)) + n = length(ns) + + for j in 1:n + for k in j + 1:n + add_edge!(ograph, ns[j], ns[k]) + end + end + end + + ograph +end + + +function Base.deepcopy(ograph::OrderedGraph) + order = deepcopy(ograph.order) + graph = deepcopy(ograph.graph) + OrderedGraph(graph, order) +end + + +# Get the vertex σ(i). +function order(ograph::OrderedGraph, i) + ograph.order[i] +end + + +# Get the index σ⁻¹(v), +function inverse(ograph::OrderedGraph, v) + inverse(ograph.order, v) +end + + +############################ +# Abstract Graph Interface # +############################ + + +function BasicGraphs.nv(ograph::OrderedGraph) + nv(ograph.graph) +end + + +function BasicGraphs.inneighbors(ograph::OrderedGraph, i) + inneighbors(ograph.graph, i) +end + + +function BasicGraphs.outneighbors(ograph::OrderedGraph, i) + outneighbors(ograph.graph, i) +end + + +function BasicGraphs.edges(ograph::OrderedGraph) + edges(ograph.graph) +end + + +function BasicGraphs.has_edge(ograph::OrderedGraph, i, j) + u = min(i, j) + v = max(i, j) + has_edge(ograph.graph, u, v) +end + + +function BasicGraphs.add_edge!(ograph::OrderedGraph, i, j) + u = min(i, j) + v = max(i, j) + + if !has_edge(ograph.graph, u, v) + add_edge!(ograph.graph, u, v) + true + else + false + end +end + + +function BasicGraphs.vertices(ograph::OrderedGraph) + vertices(ograph.graph) +end + + +function BasicGraphs.src(ograph::OrderedGraph, i) + src(ograph.graph, i) +end + + +function BasicGraphs.tgt(ograph::OrderedGraph, i) + tgt(ograph.graph, i) +end diff --git a/src/junction_trees/orders.jl b/src/junction_trees/orders.jl index 48f2a67..121ed64 100644 --- a/src/junction_trees/orders.jl +++ b/src/junction_trees/orders.jl @@ -1,170 +1,70 @@ -# A total ordering of the numbers {1, ..., n}. +# A permutation σ of the set {1, ..., n}. struct Order <: AbstractVector{Int} order::Vector{Int} index::Vector{Int} end -# Given a vector σ, construct the order ≺, where -# σ(i₁) ≺ σ(i₂) -# if -# i₁ < i₂. +# Construct an permutation σ from a vector +# (σ(1), ..., σ(n)). function Order(order::AbstractVector) n = length(order) index = Vector{Int}(undef, n) for i in 1:n index[order[i]] = i - end - + end + Order(order, index) end -# Construct an empty order of length n. -function Order(n::Integer) - order = Vector{Int}(undef, n) - index = Vector{Int}(undef, n) - Order(order, index) +# Determine if i < j, where +# u = σ(i) +# v = σ(j) +function (order::Order)(u, v) + inverse(order, u) < inverse(order, v) end -# Construct an elimination order using the reverse Cuthill-McKee algorithm. Uses -# CuthillMcKee.jl. -function Order(graph::AbstractSymmetricGraph, ::CuthillMcKeeJL_RCM) - order = CuthillMcKee.symrcm(adjacencymatrix(graph)) - Order(order) +# Compose two permutations. +function compose(left::Order, right::Order) + Order(right.order[left.order], left.index[right.index]) end -# Construct an elimination order using the approximate minimum degree algorithm. Uses -# AMD.jl. -function Order(graph::AbstractSymmetricGraph, ::AMDJL_AMD) - order = AMD.symamd(adjacencymatrix(graph)) - Order(order) +# Construct the inverse permutation. +function inverse(order::Order) + Order(order.index, order.order) end -# Construct an elimination order using the nested dissection heuristic. Uses Metis.jl. -function Order(graph::AbstractSymmetricGraph, ::MetisJL_ND) - order, index = Metis.permutation(adjacencymatrix(graph)) - Order(order, index) +# Get the index σ⁻¹(v), +function inverse(order::Order, v) + order.index[v] end -# Construct an elimination order using the maximum cardinality search algorithm. -function Order(graph::AbstractSymmetricGraph, ::MCS) - order, index = mcs(graph) - Order(order, index) -end +############################# +# Abstract Vector Interface # +############################# -# Compose as permutations. -function compose(order₁::Order, order₂::Order) - order = order₂.order[order₁.order] - index = order₁.index[order₂.index] - Order(order, index) -end - - -# Evaluate whether -# n₁ < n₂ -# in the given order. -function Base.isless(order::Order, n₁::Integer, n₂::Integer) - order.index[n₁] < order.index[n₂] -end - - -# Compute a vertex elimination order using the maximum cardinality search algorithm. -# -# The complexity is -# 𝒪(m + n), -# where m = |E| and n = |V|. -# -# https://doi.org/10.1137/0213035 -# 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{Vector{Int}}(undef, n) - - set .= [[]] - size .= 1 - append!(set[1], vertices(graph)) - - i = n - j = 1 - - while i >= 1 - v = pop!(set[j]) - α[v] = i - α⁻¹[i] = v - size[v] = 0 - - for w in neighbors(graph, v) - if size[w] >= 1 - deletesorted!(set[size[w]], w) - size[w] += 1 - insertsorted!(set[size[w]], w) - end - end - - i -= 1 - j += 1 - - while j >= 1 && isempty(set[j]) - j -= 1 - end - end - - α⁻¹, α +function Base.getindex(order::Order, i) + order.order[i] 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) - end - - nzval = ones(Int, length(rowval)) - SparseMatrixCSC(n, n, colptr, rowval, nzval) +function Base.IndexStyle(::Type{Order}) + IndexLinear() end -############################ -# AbstractVector Interface # -############################ - - function Base.size(order::Order) (length(order.order),) end -function Base.getindex(order::Order, i::Integer) - order.order[i] -end - - -function Base.setindex!(order::Order, v::Integer, i::Integer) - order.order[i] = v - order.index[v] = i - v -end - - -function Base.IndexStyle(::Type{Order}) - IndexLinear() +function Base.deepcopy(order::Order) + Order(copy(order.order), copy(order.index)) end diff --git a/src/junction_trees/postorder_trees.jl b/src/junction_trees/postorder_trees.jl new file mode 100644 index 0000000..450282c --- /dev/null +++ b/src/junction_trees/postorder_trees.jl @@ -0,0 +1,104 @@ +# A postordered rooted tree. +struct PostorderTree <: AbstractTree + parent::Vector{Int} # parent + children::Vector{Vector{Int}} # children + level::Vector{Int} # level + descendant::Vector{Int} # first descendant +end + + +# Construct a tree from a postordered list of parents. +function PostorderTree(parent::AbstractVector) + n = length(parent) + children = Vector{Vector{Int}}(undef, n) + level = Vector{Int}(undef, n) + descendant = Vector{Int}(undef, n) + + for i in 1:n + children[i] = [] + level[i] = 0 + descendant[i] = i + end + + for i in 1:n - 1 + j = parent[i] + push!(children[j], i) + descendant[j] = min(descendant[i], descendant[j]) + end + + for i in n - 1:-1:1 + j = parent[i] + level[i] = level[j] + 1 + end + + PostorderTree(parent, children, level, descendant) +end + + +# Postorder a tree. +function PostorderTree(tree::AbstractTree, order::Order) + n = length(tree) + parent = collect(1:n) + + for i in 1:n - 1 + parent[i] = inverse(order, parentindex(tree, order[i])) + end + + PostorderTree(parent) +end + + +# The number of node in a tree. +function Base.length(tree::PostorderTree) + length(tree.parent) +end + + +# Get the level of a node i. +function level(tree::PostorderTree, i::Integer) + tree.level[i] +end + + +# Get the first descendant of a node i. +function firstdescendant(tree::PostorderTree, i::Integer) + tree.descendant[i] +end + + +# Determine whether the node i is a descendant of the node j. +function isdescendant(tree::PostorderTree, i::Integer, j::Integer) + getdescendant(tree, j) <= i < j +end + + +########################## +# Indexed Tree Interface # +########################## + + +function AbstractTrees.parentindex(tree::PostorderTree, i::Integer) + if i != rootindex(tree) + tree.parent[i] + end +end + + +function AbstractTrees.childindices(tree::PostorderTree, i::Integer) + tree.children[i] +end + + +function AbstractTrees.rootindex(tree::PostorderTree) + length(tree) +end + + +function AbstractTrees.NodeType(::Type{IndexNode{PostorderTree, Int}}) + HasNodeType() +end + + +function AbstractTrees.nodetype(::Type{IndexNode{PostorderTree, Int}}) + IndexNode{PostorderTree, Int} +end diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl new file mode 100644 index 0000000..5bb5ba6 --- /dev/null +++ b/src/junction_trees/supernode_trees.jl @@ -0,0 +1,205 @@ +# An ordered graph (G, σ) equipped with a supernodal elimination tree T. +struct SupernodeTree <: AbstractTree + tree::PostorderTree # tree + ograph::OrderedGraph # ordered graph + representative::Vector{Int} # representative vertex + cardinality::Vector{Int} # supernode cardinality + ancestor::Vector{Int} # first ancestor + degrees::Vector{Int} # higher degrees +end + + +# Construct a supernodal elimination tree using an elimination algorithm. +function SupernodeTree( + graph::AbstractSymmetricGraph, + ealg::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, + stype::SupernodeType=DEFAULT_SUPERNODE_TYPE) + + SupernodeTree(EliminationTree(graph, ealg), stype) +end + + +# Construct a supernodal elimination tree. +function SupernodeTree(etree::EliminationTree, stype::SupernodeType=DEFAULT_SUPERNODE_TYPE) + degree = outdegrees(etree) + snode, parent, ancestor = stree(etree, degree, stype) + tree = Tree(parent) + sorder = postorder(tree) + tree = PostorderTree(tree, sorder) + + order = Order(vcat(snode[sorder]...)) + ograph = OrderedGraph(etree.ograph, order) + + n = length(tree) + representative = zeros(Int, n) + cardinality = zeros(Int, n) + _ancestor = zeros(Int, n) + _degree = zeros(Int, n) + + for i in 1:n - 1 + j = sorder[i] + representative[i] = inverse(order, snode[j][1]) + cardinality[i] = length(snode[j]) + _degree[i] = degree[j] + _ancestor[i] = inverse(order, ancestor[j]) + end + + representative[n] = inverse(order, snode[n][1]) + cardinality[n] = length(snode[n]) + _degree[n] = degree[n] + _ancestor[n] = ancestor[n] + + SupernodeTree(tree, ograph, representative, cardinality, _ancestor, _degree) +end + + +# Chordal Graphs and Semidefinite Optimization +# Vanderberghe and Andersen +# Algorithm 4.1: Maximal supernodes and supernodal elimination tree. +function stree(etree::EliminationTree, degree::AbstractVector, stype::SupernodeType) + n = length(etree) + index = zeros(Int, n) + snd = Vector{Int}[] + q = Int[] + a = Int[] + + for v in 1:n + ww = findchild(etree, degree, stype, v) + + if isnothing(ww) + i = length(snd) + 1 + index[v] = i + push!(snd, [v]) + push!(q, length(snd)) + push!(a, n + 1) + else + i = index[ww] + index[v] = i + push!(snd[i], v) + end + + for w in childindices(etree, v) + if w !== ww + j = index[w] + q[j] = i + a[j] = v + end + end + end + + snd, q, a +end + + +# Construct an elimination graph. +function eliminationgraph(stree::SupernodeTree) + ograph = deepcopy(stree.ograph) + n = length(stree) + + for i in 1:n - 1 + for u in supernode(stree, i)[1:end - 1] + v = u + 1 + + for w in outneighbors(ograph, u) + if v < w + add_edge!(ograph, v, w) + end + end + end + + u = last(supernode(stree, i)) + v = first(supernode(stree, parentindex(stree, i))) + + for w in outneighbors(ograph, u) + if v < w + add_edge!(ograph, v, w) + end + end + end + + ograph +end + + +# Compute the width of a supernodal elimination tree. +function width(stree::SupernodeTree) + maximum(stree.degree[stree.representative]) +end + + +# Get the (sorted) supernode at node i. +function supernode(stree::SupernodeTree, i::Integer) + v = stree.representative[i] + n = stree.cardinality[i] + v:v + n - 1 +end + + +# Compute the (unsorted) seperators of every node in T. +function seperators(stree::SupernodeTree) + n = length(stree) + seperator = Vector{Vector{Int}}(undef, n) + ograph = eliminationgraph(stree) + + for i in 1:n + clique = collect(outneighbors(ograph, stree.representative[i])) + filter!(j -> stree.ancestor[i] <= j, clique) + seperator[i] = clique + end + + seperator +end + + +# Get the number of nodes in T. +function Base.length(stree::SupernodeTree) + length(stree.tree) +end + + +# Get the level of node i. +function level(stree::SupernodeTree, i) + level(stree.tree, i) +end + + +# Get the vertex σ(i). +function order(stree::SupernodeTree, i) + order(stree.ograph, i) +end + + +# Get the index σ⁻¹(v), +function inverse(stree::SupernodeTree, i) + inverse(stree.ograph, i) +end + + +########################## +# Indexed Tree Interface # +########################## + + +function AbstractTrees.rootindex(stree::SupernodeTree) + rootindex(stree.tree) +end + + +function AbstractTrees.parentindex(stree::SupernodeTree, i) + parentindex(stree.tree, i) +end + + +function AbstractTrees.childindices(stree::SupernodeTree, i) + childindices(stree.tree, i) +end + + +function AbstractTrees.NodeType(::Type{IndexNode{SupernodeTree, Int}}) + HasNodeType() +end + + +function AbstractTrees.nodetype(::Type{IndexNode{SupernodeTree, Int}}) + IndexNode{SupernodeTree, Int} +end diff --git a/src/junction_trees/supernode_types.jl b/src/junction_trees/supernode_types.jl new file mode 100644 index 0000000..da92a40 --- /dev/null +++ b/src/junction_trees/supernode_types.jl @@ -0,0 +1,70 @@ +""" + SupernodeType + +A type of supernode. The options are +- [`Node`](@ref) +- [`MaximalSupernode`](@ref) +- [`FundamentalSupernode`](@ref) +""" +abstract type SupernodeType end + + +""" + Node <: Supernode + +A single-vertex supernode. +""" +struct Node <: SupernodeType end + + +""" + MaximalSupernode <: Supernode + +A maximal supernode. +""" +struct Maximal <: SupernodeType end + + +""" + FundamentalSupernode <: Supernode + +A fundamental supernode. +""" +struct Fundamental <: SupernodeType end + + +# Find a child w of v such that +# v ∈ snd(w). +# If no such child exists, return nothing. +function findchild(etree::EliminationTree, degree::AbstractVector, stype::Node, v::Integer) end + + +# Find a child w of v such that +# v ∈ snd(w). +# If no such child exists, return nothing. +function findchild(etree::EliminationTree, degree::AbstractVector, stype::Maximal, v::Integer) + for w in childindices(etree, v) + if degree[w] == degree[v] + 1 + return w + end + end +end + + +# Find a child w of v such that +# v ∈ snd(w). +# If no such child exists, return nothing. +function findchild(etree::EliminationTree, degree::AbstractVector, stype::Fundamental, v::Integer) + ws = childindices(etree, v) + + if length(ws) == 1 + w = only(ws) + + if degree[w] == degree[v] + 1 + return w + end + end +end + + +const DEFAULT_SUPERNODE_TYPE = Maximal() diff --git a/src/junction_trees/supernodes.jl b/src/junction_trees/supernodes.jl deleted file mode 100644 index af940bb..0000000 --- a/src/junction_trees/supernodes.jl +++ /dev/null @@ -1,36 +0,0 @@ -""" - Supernode - -A type of supernode. The options are -- [`Node`](@ref) -- [`MaximalSupernode`](@ref) -- [`FundamentalSupernode`](@ref) -""" -abstract type Supernode end - - -""" - Node <: Supernode - -A single-vertex supernode. -""" -struct Node <: Supernode end - - -""" - MaximalSupernode <: Supernode - -A maximal supernode. -""" -struct MaximalSupernode <: Supernode end - - -""" - FundamentalSupernode <: Supernode - -A fundamental supernode. -""" -struct FundamentalSupernode <: Supernode end - - -const DEFAULT_SUPERNODE = MaximalSupernode() diff --git a/src/junction_trees/trees.jl b/src/junction_trees/trees.jl index ab8655f..1b2aac8 100644 --- a/src/junction_trees/trees.jl +++ b/src/junction_trees/trees.jl @@ -1,202 +1,37 @@ # A rooted tree. -struct Tree +struct Tree <: AbstractTree root::Int - parentlist::Vector{Int} - childrenlist::Vector{Vector{Int}} - levellist::Vector{Int} - firstdescendantlist::Vector{Int} -end - - -# Orient a tree towards the given root. -function Tree(root::Integer, tree::Tree) - i = root - parent = parentindex(tree, i) - parentlist = copy(tree.parentlist) - childrenlist = deepcopy(tree.childrenlist) - - while !isnothing(parent) - parentlist[parent] = i - push!(childrenlist[i], parent) - deletesorted!(childrenlist[parent], i) - i = parent - parent = parentindex(tree, i) - end - - Tree(root, parentlist, childrenlist) -end - - -# Construct a tree from a list of parent and a list of children. -function Tree(root::Integer, parentlist::AbstractVector, childrenlist::AbstractVector) - n = length(parentlist) - levellist = Vector{Int}(undef, n) - firstdescendantlist = Vector{Int}(undef, n) - Tree(root, parentlist, childrenlist, levellist, firstdescendantlist) + parent::Vector{Int} + children::Vector{Vector{Int}} end # Construct a tree from a list of parents. -function Tree(root::Integer, parentlist::AbstractVector) - n = length(parentlist) - childrenlist = Vector{Vector{Int}}(undef, n) - childrenlist .= [[]] - +function Tree(parent::AbstractVector) + n = root = length(parent) + children = Vector{Vector{Int}}(undef, n) + for i in 1:n - if i != root - push!(childrenlist[parentlist[i]], i) - end + children[i] = [] end - Tree(root, parentlist, childrenlist) -end - - -# Construct an elimination tree. -function Tree(graph::AbstractSymmetricGraph, order::Order) - n = nv(graph) - parentlist = makeetree(graph, order) - @assert count(parentlist .== 0) == 1 - Tree(n, parentlist) -end - - -function Base.length(tree::Tree) - length(tree.parentlist) -end - - -# Compute the parent vector of the elimination tree of the elimination graph of a ordered -# graph. -# -# The complexity is -# 𝒪(m log(n)) -# where m = |E| and n = |V|. -# -# doi:10.1145/6497.6499 -# Algorithm 4.2: Elimination Tree by Path Compression -function makeetree(graph::AbstractSymmetricGraph, order::Order) - graph = Graph(graph, order) - - n = nv(graph) - parent = Vector{Int}(undef, n) - ancestor = Vector{Int}(undef, n) - for i in 1:n - parent[i] = 0 - ancestor[i] = 0 - - for k in inneighbors(graph, i) - r = k - - while ancestor[r] != 0 && ancestor[r] != i - t = ancestor[r] - ancestor[r] = i - r = t - end - - if ancestor[r] == 0 - ancestor[r] = i - parent[r] = i - end - end - end - - parent -end - - -# Given an ordered graph -# (G, σ), -# construct a directed graph by ordering the edges in G from lower to higher index. -# -# The complexity is -# 𝒪(m) -# where m = |E|. -function BasicGraphs.Graph(graph::AbstractSymmetricGraph, order::Order) - n = nv(graph) - digraph = Graph(n) - - for v in vertices(graph) - i = order.index[v] - - for w in neighbors(graph, v) - j = order.index[w] - - if i < j - add_edge!(digraph, i, j) - end + j = parent[i] + + if i == j + root = i + else + push!(children[j], i) end end - digraph -end - - -############## -# Postorders # -############## - - -# Get the level of node i. -# This function only works on postordered trees. -function getlevel(tree::Tree, i::Integer) - tree.levellist[i] -end - - -# Get the first descendant of node i. -# This function only works on postordered trees. -function getfirstdescendant(tree::Tree, i::Integer) - tree.firstdescendantlist[i] -end - - -# Evaluate whether node i₁ is a descendant of node i₂. -# This function only works on postordered trees. -function AbstractTrees.isdescendant(tree::Tree, i₁::Integer, i₂::Integer) - getfirstdescendant(tree, i₂) <= i₁ < i₂ + Tree(root, parent, children) end -# Compute a postordering of a tree. -# -# The complexity is -# 𝒪(n) -# where n = |V|. -function makepostorder(tree::Tree) - n = length(tree) - order = Order(n) - parentlist = Vector{Int}(undef, n) - childrenlist = Vector{Vector{Int}}(undef, n) - levellist = Vector{Int}(undef, n) - firstdescendantlist = Vector{Int}(undef, n) - - root, nodes... = PreOrderDFS(IndexNode(tree)) - - order[n] = root.index - parentlist[n] = 0 - childrenlist[n] = [] - levellist[n] = 0 - - for (i, node) in enumerate(nodes) - j = n - i - order[j] = node.index - - k = order.index[parentindex(tree, node.index)] - parentlist[j] = k - childrenlist[j] = [] - pushfirst!(childrenlist[k], j) - levellist[j] = 1 + levellist[k] - end - - for i in 1:n - init = i - firstdescendantlist[i] = minimum(firstdescendantlist[childrenlist[i]]; init) - end - - tree = Tree(n, parentlist, childrenlist, levellist, firstdescendantlist) - order, tree +# Get the number of nodes in a tree. +function Base.length(tree::Tree) + length(tree.parent) end @@ -212,13 +47,13 @@ end function AbstractTrees.parentindex(tree::Tree, i::Integer) if i != rootindex(tree) - tree.parentlist[i] + tree.parent[i] end end function AbstractTrees.childindices(tree::Tree, i::Integer) - tree.childrenlist[i] + tree.children[i] end diff --git a/src/nested_uwds/NestedUWDs.jl b/src/nested_uwds/NestedUWDs.jl deleted file mode 100644 index 5de7eca..0000000 --- a/src/nested_uwds/NestedUWDs.jl +++ /dev/null @@ -1,34 +0,0 @@ -module NestedUWDs - - -using AbstractTrees -using Catlab.ACSetInterface -using Catlab.BasicGraphs -using Catlab.DirectedWiringDiagrams -using Catlab.DirectedWiringDiagrams: WiringDiagramACSet -using Catlab.MonoidalUndirectedWiringDiagrams -using Catlab.MonoidalUndirectedWiringDiagrams: UntypedHypergraphDiagram -using Catlab.RelationalPrograms -using Catlab.RelationalPrograms: TypedUnnamedRelationDiagram -using Catlab.Theories -using Catlab.UndirectedWiringDiagrams -using Catlab.WiringDiagramAlgebras - -using ..JunctionTrees -using ..JunctionTrees: insertsorted!, DEFAULT_ELIMINATION_ALGORITHM, DEFAULT_SUPERNODE - -# Elimination Algorithms -export EliminationAlgorithm, AMDJL_AMD, CuthillMcKeeJL_RCM, MetisJL_ND, MCS - -# Supernodes -export Supernode, Node, MaximalSupernode, FundamentalSupernode - -# Nested UWDs -export NestedUWD -export evalschedule, makeschedule, makeoperations - - -include("nested_uwds.jl") - - -end diff --git a/src/nested_uwds/nested_uwds.jl b/src/nested_uwds/nested_uwds.jl deleted file mode 100644 index 54dd648..0000000 --- a/src/nested_uwds/nested_uwds.jl +++ /dev/null @@ -1,312 +0,0 @@ -""" - NestedUWD{T, B, V} - -An undirected wiring diagram, represented as a nested collected of undirected wiring -diagrams. -""" -struct NestedUWD{T, B, V} - diagram::TypedUnnamedRelationDiagram{T, B, V} - jtree::JunctionTree - assignmentlist::Vector{Int} - assignmentindex::Vector{Vector{Int}} -end - - -function NestedUWD( - diagram::D, - jtree::JunctionTree, - assignmentlist::AbstractVector, - assignmentindex::AbstractVector) where D <: UndirectedWiringDiagram - - T, B, V = getattributetypes(D) - relation = TypedUnnamedRelationDiagram{T, B, V}() - copy_parts!(relation, diagram) - NestedUWD{T, B, V}(relation, jtree, assignmentlist, assignmentindex) -end - - -function NestedUWD(diagram::UndirectedWiringDiagram, jtree::JunctionTree) - n = nparts(diagram, :Box) - m = length(jtree) - assignmentlist = Vector{Int}(undef, n) - assignmentindex = Vector{Vector{Int}}(undef, m) - assignmentindex .= [[]] - - for b in 1:n - i = getsubtree(jtree, diagram[incident(diagram, b, :box), :junction]) - assignmentlist[b] = i - push!(assignmentindex[i], b) - end - - NestedUWD(diagram, jtree, assignmentlist, assignmentindex) -end - - -""" - NestedUWD( - diagram::UndirectedWiringDiagram, - [, algorithm::Union{Order, EliminationAlgorithm}] - [, supernode::Supernode]) - -Construct a nested undirected wiring diagram. -""" -function NestedUWD( - diagram::UndirectedWiringDiagram, - algorithm::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, - supernode::Supernode=DEFAULT_SUPERNODE) - - jtree = JunctionTree(diagram, algorithm, supernode) - NestedUWD(diagram, jtree) -end - - -# Construct a tree decomposition of the line graph of an undirected wiring diagram. -function JunctionTree( - diagram::UndirectedWiringDiagram, - algorithm::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, - supernode::Supernode=DEFAULT_SUPERNODE) - - graph = makegraph(diagram) - jtree = JunctionTree(graph, algorithm, supernode) - - query = diagram[:outer_junction] - JunctionTree(getsubtree(jtree, query), jtree) -end - - -# Construct the line graph of an undirected wiring diagram. -function makegraph(diagram::UndirectedWiringDiagram) - n = nparts(diagram, :Junction) - m = nparts(diagram, :Box) - graph = SymmetricGraph(n) - - for i in 1:m - junctions = diagram[incident(diagram, i, :box), :junction] - l = length(junctions) - - for i₁ in 1:l - 1 - j₁ = junctions[i₁] - - for i₂ in i₁ + 1:l - j₂ = junctions[i₂] - - if !has_edge(graph, j₁, j₂) - add_edge!(graph, j₁, j₂) - end - end - end - end - - junctions = diagram[:, :outer_junction] - l = length(junctions) - - for i₁ in 1:l - 1 - j₁ = junctions[i₁] - - for i₂ in i₁ + 1:l - j₂ = junctions[i₂] - - if !has_edge(graph, j₁, j₂) - add_edge!(graph, j₁, j₂) - end - end - end - - graph -end - - -""" - makeschedule(nuwd::NestedUWD) - -Construct a directed wiring diagram that represents the nesting structure of a nested UWD. -""" -function makeschedule(nuwd::NestedUWD{<:Any, T}) where T - m = length(nuwd.assignmentlist) - n = length(nuwd.jtree) - - parents = map(1:n - 1) do i - parentindex(nuwd.jtree, i) - end - - costs = map(1:n) do i - length(getresidual(nuwd.jtree, i)) + length(getseperator(nuwd.jtree, i)) - end - - schedule = WiringDiagramACSet{T, Nothing, Union{Int, AbstractBox}, DataType}() - - add_parts!(schedule, :Box, n) - add_parts!(schedule, :Wire, n - 1) - add_parts!(schedule, :InPort, m + n - 1) - add_parts!(schedule, :InWire, m) - add_parts!(schedule, :OutPort, n) - add_parts!(schedule, :OutWire, 1) - add_parts!(schedule, :OuterInPort, m) - add_parts!(schedule, :OuterOutPort, 1) - - schedule[:, :src] = 1:n - 1 - schedule[:, :tgt] = m + 1:m + n - 1 - schedule[:, :in_src] = 1:m - schedule[:, :in_tgt] = 1:m - schedule[:, :out_src] = n:n - schedule[:, :out_tgt] = 1:1 - schedule[:, :in_port_box] = [nuwd.assignmentlist; parents] - schedule[:, :out_port_box] = 1:n - - schedule[:, :value] = costs - schedule[:, :box_type] = Box{Int} - schedule[:, :outer_in_port_type] = nuwd.diagram[:, :name] - - Theory = ThSymmetricMonoidalCategory.Meta.T - WiringDiagram{Theory, T, Nothing, Int}(schedule, nothing) -end - - -""" - function evalschedule( - f, - nuwd::NestedUWD, - generators::Union{AbstractVector, AbstractDict} - [, operations::AbstractVector]) - -Evaluate an undirected wiring diagrams given a set of generators for the boxes. The -optional first argument `f` should be callable with the signature -``` - f(diagram, generators) -``` -where `diagram` is an undirected wiring diagram, and `generators` is a vector. If `f` is not -specified, then it defaults to `oapply`. -""" -function evalschedule( - f, - nuwd::NestedUWD, - generators::AbstractVector{T}, - operations::AbstractVector=makeoperations(nuwd)) where T - - n = length(nuwd.jtree) - mailboxes = Vector{T}(undef, n) - - for i in 1:n - g₁ = generators[nuwd.assignmentindex[i]] - g₂ = mailboxes[childindices(nuwd.jtree, i)] - mailboxes[i] = f(operations[i], [g₁; g₂]) - end - - mailboxes[n] -end - - -function evalschedule( - f, - nuwd::NestedUWD, - generators::AbstractDict{<:Any, T}, - operations::AbstractVector=makeoperations(nuwd)) where T - - g = generators - n = nparts(nuwd.diagram, :Box) - generators = Vector{T}(undef, n) - - for i in 1:n - generators[i] = g[nuwd.diagram[i, :name]] - end - - evalschedule(f, nuwd, generators, operations) -end - - -function evalschedule( - nuwd::NestedUWD, - generators::Union{AbstractVector, AbstractDict}, - operations::AbstractVector=makeoperations(nuwd)) - - evalschedule(oapply, nuwd, generators, operations) -end - - -# For each node i of a nested UWD, construct the undirected wiring diagram corresponding to i. -function makeoperations(nuwd::NestedUWD) - m = length(nuwd.jtree) - - map(1:m) do i - makeoperation(nuwd, i) - end -end - - -# Construct the undirected wiring diagram corresponding to node i of a nested UWD. -function makeoperation(nuwd::NestedUWD{T, B, V}, i::Integer) where {T, B, V} - function findjunction(j::Integer) - v = nuwd.jtree.stree.order.index[j] - v₁ = nuwd.jtree.stree.firstsupernodelist[i] - v₂ = nuwd.jtree.stree.lastsupernodelist[i] - - if v <= v₂ - v - v₁ + 1 - else - v₂ - v₁ + 1 + searchsortedfirst(nuwd.jtree.seperatorlist[i], v) - end - end - - residual = getresidual(nuwd.jtree, i) - seperator = getseperator(nuwd.jtree, i) - m = length(residual) - n = length(seperator) - - operation = TypedUnnamedRelationDiagram{T, B, V}() - add_parts!(operation, :Junction, m + n) - - operation[1:m, :junction_type] = nuwd.diagram[residual, :junction_type] - operation[1:m, :variable] = nuwd.diagram[residual, :variable] - operation[m + 1:m + n, :junction_type] = nuwd.diagram[seperator, :junction_type] - operation[m + 1:m + n, :variable] = nuwd.diagram[seperator, :variable] - - if i < length(nuwd.jtree) - for j in seperator - p′ = add_part!(operation, :OuterPort) - operation[p′, :outer_junction] = m + p′ - operation[p′, :outer_port_type] = nuwd.diagram[j, :junction_type] - end - else - for j in nuwd.diagram[:outer_junction] - p′ = add_part!(operation, :OuterPort) - operation[p′, :outer_junction] = findjunction(j) - operation[p′, :outer_port_type] = nuwd.diagram[j, :junction_type] - end - end - - for b in nuwd.assignmentindex[i] - b′ = add_part!(operation, :Box) - operation[b′, :name] = nuwd.diagram[b, :name] - - for j in nuwd.diagram[incident(nuwd.diagram, b, :box), :junction] - p′ = add_part!(operation, :Port) - operation[p′, :box] = b′ - operation[p′, :junction] = findjunction(j) - operation[p′, :port_type] = nuwd.diagram[j, :junction_type] - end - end - - for b in childindices(nuwd.jtree, i) - b′ = add_part!(operation, :Box) - - for j in getseperator(nuwd.jtree, b) - p′ = add_part!(operation, :Port) - operation[p′, :box] = b′ - operation[p′, :junction] = findjunction(j) - operation[p′, :port_type] = nuwd.diagram[j, :junction_type] - end - end - - operation -end - - -# Get the attribute types of an undirected wiring diagram. -function getattributetypes(::Type{<:UntypedRelationDiagram{B, V}}) where {B, V} - Nothing, B, V -end - - -function getattributetypes(::Type{<:TypedRelationDiagram{T, B, V}}) where {T, B, V} - T, B, V -end diff --git a/test/JunctionTrees.jl b/test/JunctionTrees.jl index 9001215..763dab5 100644 --- a/test/JunctionTrees.jl +++ b/test/JunctionTrees.jl @@ -15,18 +15,32 @@ add_edges!(graph, [3, 4, 5, 15, 3, 4, 9, 16, 9, 16, 8, 9, 15, 11, 13, 14, 17, 13, 14, 16, 17, 17]) order = JunctionTrees.Order(graph, CuthillMcKeeJL_RCM()) -@test order == [2, 14, 13, 11, 4, 3, 12, 10, 16, 1, 17, 5, 6, 15, 9, 7, 8] +@test length(order) == 17 order = JunctionTrees.Order(graph, AMDJL_AMD()) -@test order == [8, 11, 7, 2, 4, 3, 1, 6, 13, 14, 10, 12, 17, 16, 5, 9, 15] +@test length(order) == 17 order = JunctionTrees.Order(graph, MetisJL_ND()) -@test order == [11, 17, 14, 13, 10, 12, 8, 6, 7, 5, 4, 3, 9, 2, 1, 16, 15] - -order = JunctionTrees.Order(graph, MCS()) -@test order == [2, 3, 4, 8, 1, 5, 6, 9, 7, 11, 13, 10, 14, 16, 12, 15, 17] +@test length(order) == 17 order = JunctionTrees.Order(1:17) +@test length(order) == 17 + + +# Figure 4.2 +etree = EliminationTree(graph, order) +@test parentindex.([etree], 1:17) == [3, 3, 4, 5, 9, 9, 8, 9, 15, 11, 13, 13, 14, 16, 16, 17, 17] + + +# Figure 4.3 +stree = SupernodeTree(etree, Node()) + + +#= +@test supernode.([stree], order(stree, 1:17)) + + + parent = JunctionTrees.makeetree(graph, order) # Figure 4.2 @@ -283,3 +297,4 @@ jtree = JunctionTree(graph, order, FundamentalSupernode()) @test !isdescendant(jtree, getsubtree(jtree, 15), getsubtree(jtree, 10)) @test !isdescendant(jtree, getsubtree(jtree, 1), getsubtree(jtree, 1)) @test getwidth(jtree) == 4 +=# diff --git a/test/NestedUWDs.jl b/test/NestedUWDs.jl deleted file mode 100644 index e156351..0000000 --- a/test/NestedUWDs.jl +++ /dev/null @@ -1,122 +0,0 @@ -using Catlab.RelationalPrograms -using Catlab.UndirectedWiringDiagrams -using LinearAlgebra -using StructuredDecompositions.NestedUWDs -using Test - - -# CategoricalTensorNetworks.jl -# https://github.com/AlgebraicJulia/CategoricalTensorNetworks.jl/ -function contract_tensor_network(d::UndirectedWiringDiagram, - tensors::AbstractVector{<:AbstractArray}) - @assert nboxes(d) == length(tensors) - juncs = [junction(d, ports(d, b)) for b in boxes(d)] - j_out = junction(d, ports(d, outer=true), outer=true) - contract_tensor_network(tensors, juncs, j_out) -end - - -function contract_tensor_network(tensors::AbstractVector{<:AbstractArray{T}}, - juncs::AbstractVector, j_out) where T - # Handle important binary case with specialized code. - if length(tensors) == 2 && length(juncs) == 2 - return contract_tensor_network(Tuple(tensors), Tuple(juncs), j_out) - end - - jsizes = Tuple(infer_junction_sizes(tensors, juncs, j_out)) - juncs, j_out = map(Tuple, juncs), Tuple(j_out) - C = zeros(T, Tuple(jsizes[j] for j in j_out)) - for index in CartesianIndices(jsizes) - x = one(T) - for (A, junc) in zip(tensors, juncs) - x *= A[(index[j] for j in junc)...] - end - C[(index[j] for j in j_out)...] += x - end - return C -end - - -function contract_tensor_network( # Binary case. - (A, B)::Tuple{<:AbstractArray{T},<:AbstractArray{T}}, - (jA, jB), j_out) where T - jsizes = Tuple(infer_junction_sizes((A, B), (jA, jB), j_out)) - jA, jB, j_out = Tuple(jA), Tuple(jB), Tuple(j_out) - C = zeros(T, Tuple(jsizes[j] for j in j_out)) - for index in CartesianIndices(jsizes) - C[(index[j] for j in j_out)...] += - A[(index[j] for j in jA)...] * B[(index[j] for j in jB)...] - end - return C -end - - -function infer_junction_sizes(tensors, juncs, j_out) - @assert length(tensors) == length(juncs) - njunc = maximum(Iterators.flatten((Iterators.flatten(juncs), j_out))) - jsizes = fill(-1, njunc) - for (A, junc) in zip(tensors, juncs) - for (i, j) in enumerate(junc) - if jsizes[j] == -1 - jsizes[j] = size(A, i) - else - @assert jsizes[j] == size(A, i) - end - end - end - @assert all(s >= 0 for s in jsizes) - jsizes -end - - -# out[v,z] = A[v,w] * B[w,x] * C[x,y] * D[y,z] -diagram = @relation (v, z) begin - A(v, w) - B(w, x) - C(x, y) - D(y, z) -end - -nuwd = NestedUWD(diagram) -A, B, C, D = map(randn, [(3, 4), (4, 5), (5, 6), (6, 7)]) -generators = Dict{Symbol, Array{Float64}}(:A => A, :B => B, :C => C, :D => D) -out = evalschedule(contract_tensor_network, nuwd, generators) -@test out ≈ A * B * C * D - -# out[] = A[w,x] * B[x,y] * C[y,z] * D[z,w] -diagram = @relation () begin - A(w, x) - B(x, y) - C(y, z) - D(z, w) -end - -nuwd = NestedUWD(diagram) -A, B, C, D = map(randn, [(10, 5), (5, 5), (5, 5), (5, 10)]) -generators = Dict{Symbol, Array{Float64}}(:A => A, :B => B, :C => C, :D => D) -out = evalschedule(contract_tensor_network, nuwd, generators) -@test out[] ≈ tr(A * B * C * D) - -# out[w,x,y,z] = A[w,x] * B[y,z] -diagram = @relation (w, x, y, z) begin - A(w, x) - B(y, z) -end - -nuwd = NestedUWD(diagram) -A, B = map(randn, [(3, 4), (5, 6)]) -generators = Dict{Symbol, Array{Float64}}(:A => A, :B => B) -out = evalschedule(contract_tensor_network, nuwd, generators) -@test out ≈ (reshape(A, (3, 4, 1, 1)) .* reshape(B, (1, 1, 5, 6))) - -# out[] = A[x,y] * B[x,y] -diagram = @relation () begin - A(x, y) - B(x, y) -end - -nuwd = NestedUWD(diagram) -A, B = map(randn, [(5, 5), (5, 5)]) -generators = Dict{Symbol, Array{Float64}}(:A => A, :B => B) -out = evalschedule(contract_tensor_network, nuwd, generators) -@test out[] ≈ dot(vec(A), vec(B)) diff --git a/test/runtests.jl b/test/runtests.jl index 6fdd6b8..fbb1669 100644 --- a/test/runtests.jl +++ b/test/runtests.jl @@ -15,7 +15,3 @@ end @testset "JunctionTrees" begin include("JunctionTrees.jl") end - -@testset "NestedUWDs" begin - include("NestedUWDs.jl") -end From 07fc95a407ae9e02925368c5540d17db64e93035 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sat, 16 Nov 2024 23:30:53 -0500 Subject: [PATCH 02/48] initial attempt at integration --- src/Decompositions.jl | 99 +++++++++++++++++++++++++ src/JunctionTrees.jl | 41 ++++++++++ src/StructuredDecompositions.jl | 2 +- src/junction_trees/JunctionTrees.jl | 41 ---------- src/junction_trees/elimination_trees.jl | 6 +- src/junction_trees/ordered_graphs.jl | 4 +- src/junction_trees/orders.jl | 4 +- src/junction_trees/supernode_trees.jl | 2 +- src/junction_trees/trees.jl | 6 +- 9 files changed, 152 insertions(+), 53 deletions(-) create mode 100644 src/JunctionTrees.jl delete mode 100644 src/junction_trees/JunctionTrees.jl diff --git a/src/Decompositions.jl b/src/Decompositions.jl index ef83687..2a9447b 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -5,15 +5,21 @@ export StructuredDecomposition, StrDecomp, 𝐃, bags, adhesions, adhesionSpans, ∫ +using ..JunctionTrees + using PartialFunctions using MLStyle +using AbstractTrees using Catlab using Catlab.CategoricalAlgebra using Catlab.Graphs using Catlab.ACSetInterface using Catlab.CategoricalAlgebra.Diagrams import Catlab.CategoricalAlgebra.Diagrams: ob_map, hom_map, colimit, limit +using Catlab.FinSets: FinSetInt, FinDomFunctionVector + + ##################### # DATA @@ -180,4 +186,97 @@ function 𝐃(f, d ::StructuredDecomposition, t::DecompType = d.decomp_type)::St StrDecomp(d.decomp_shape, Q, t) end + +################################## +# Integration with JunctionTrees # +################################## + + +const OTYPE = SymmetricGraph + + +const MTYPE = StructTightACSetTransformation{ + TypeLevelBasicSchema{ + Symbol, + Tuple{:V, :E}, + Tuple{(:src, :E, :V), (:tgt, :E, :V), (:inv, :E, :E)}, + Tuple{}, + Tuple{}, + Tuple{ + (nothing, :E, :E, ((:inv, :inv), ())), + (nothing, :E, :V, ((:inv, :src), (:tgt,))), + (nothing, :E, :V, ((:inv, :tgt), (:src,)))}}, + @NamedTuple{V::FinDomFunctionVector{Int64, Vector{Int64}, FinSetInt}, E::FinDomFunctionVector{Int64, Vector{Int64}, FinSetInt}}, + SymmetricGraph, + SymmetricGraph} + + +function Decompositions.StrDecomp(sgraph::AbstractSymmetricGraph, stree::SupernodeTree, seperator::AbstractVector) + n = length(stree) + graph = Graph(n) + objects = Vector{OTYPE}(undef, 2n - 1) + morphisms = Vector{MTYPE}(undef, 2n - 2) + + for i in 1:n + snd = supernode(stree, i) + sep = seperator[i] + objects[i] = induced_subgraph(sgraph, order(stree, [snd; sep])) + end + + for i in 1:n - 1 + sep = seperator[i] + objects[n + i] = induced_subgraph(sgraph, order(stree, sep)) + end + + for i in 1:n - 1 + j = parentindex(stree, i) + add_edge!(graph, i, j) + + sep_i = seperator[i] + sep_j = seperator[j] + snd_i = supernode(stree, i) + snd_j = supernode(stree, j) + + rep_j = first(snd_j) + len_j = length(snd_j) + + mapping = map(sep_i) do v + if v in snd_j + v - rep_j + 1 + else + len_j + searchsortedfirst(sep_j, v) + end + end + + morphisms[i] = induced_homomorphism(mapping, objects[n + i], objects[j]) + end + + for i in 1:n - 1 + sep = seperator[i] + snd = supernode(stree, i) + + mapping = length(snd) + 1:length(snd) + length(sep) + morphisms[n + i - 1] = induced_homomorphism(mapping, objects[n + i], objects[i]) + end + + StrDecomp(graph, FinDomFunctor(objects, morphisms, ∫(graph))) +end + + +function induced_homomorphism(vmapping, domain, codomain) + emapping = Vector{Int}(undef, ne(domain)) + index = Dict{Tuple{Int, Int}, Int}() + + for e in edges(codomain) + index[src(codomain, e), tgt(codomain, e)] = e + end + + for e in edges(domain) + emapping[e] = index[vmapping[src(domain, e)], vmapping[tgt(domain, e)]] + end + + ACSetTransformation(domain, codomain, V=vmapping, E=emapping) +end + + end diff --git a/src/JunctionTrees.jl b/src/JunctionTrees.jl new file mode 100644 index 0000000..1c756b4 --- /dev/null +++ b/src/JunctionTrees.jl @@ -0,0 +1,41 @@ +module JunctionTrees + + +import AMD +import CuthillMcKee +import Metis + +using AbstractTrees +using Catlab.BasicGraphs +using DataStructures +using SparseArrays + + +# Orders +export Order + + +# Elimination Algorithms +export AMDJL_AMD, CuthillMcKeeJL_RCM, MetisJL_ND + + +# Supernode Types +export Node, Maximal, Fundamental + + +# Supernode Trees +export SupernodeTree, width, height, seperators, supernode, order + + +include("junction_trees/orders.jl") +include("junction_trees/elimination_algorithms.jl") +include("junction_trees/ordered_graphs.jl") +include("junction_trees/abstract_trees.jl") +include("junction_trees/trees.jl") +include("junction_trees/postorder_trees.jl") +include("junction_trees/elimination_trees.jl") +include("junction_trees/supernode_types.jl") +include("junction_trees/supernode_trees.jl") + + +end diff --git a/src/StructuredDecompositions.jl b/src/StructuredDecompositions.jl index 2e82fef..2feb087 100644 --- a/src/StructuredDecompositions.jl +++ b/src/StructuredDecompositions.jl @@ -1,10 +1,10 @@ module StructuredDecompositions +include("JunctionTrees.jl") include("Decompositions.jl") include("FunctorUtils.jl") include("DecidingSheaves.jl") -include("junction_trees/JunctionTrees.jl") end diff --git a/src/junction_trees/JunctionTrees.jl b/src/junction_trees/JunctionTrees.jl deleted file mode 100644 index 8fbfe33..0000000 --- a/src/junction_trees/JunctionTrees.jl +++ /dev/null @@ -1,41 +0,0 @@ -module JunctionTrees - - -import AMD -import CuthillMcKee -import Metis - -using AbstractTrees -using Catlab.BasicGraphs -using DataStructures -using SparseArrays - - -# Orders -export Order - - -# Elimination Algorithms -export AMDJL_AMD, CuthillMcKeeJL_RCM, MetisJL_ND - - -# Supernode Types -export Node, Maximal, Fundamental - - -# Supernode Trees -export SupernodeTree, width, height - - -include("orders.jl") -include("elimination_algorithms.jl") -include("ordered_graphs.jl") -include("abstract_trees.jl") -include("trees.jl") -include("postorder_trees.jl") -include("elimination_trees.jl") -include("supernode_types.jl") -include("supernode_trees.jl") - - -end diff --git a/src/junction_trees/elimination_trees.jl b/src/junction_trees/elimination_trees.jl index 0102eff..97ddfb6 100644 --- a/src/junction_trees/elimination_trees.jl +++ b/src/junction_trees/elimination_trees.jl @@ -1,8 +1,8 @@ # An ordered graph (G, σ) equipped with the elimination tree T of its elimination graph. # Nodes i in T correspond to vertices σ(i) in G. struct EliminationTree{T <: AbstractTree} <: AbstractTree - tree::T - ograph::OrderedGraph + tree::T # elimination tree + ograph::OrderedGraph # ordered graph end @@ -30,7 +30,7 @@ end # Liu # Algorithm 4.2: Elimination Tree by Path Compression. function etree(ograph::OrderedGraph) - n = nv(graph) + n = nv(ograph) parent = collect(1:n) ancestor = collect(1:n) diff --git a/src/junction_trees/ordered_graphs.jl b/src/junction_trees/ordered_graphs.jl index 1915e69..2338126 100644 --- a/src/junction_trees/ordered_graphs.jl +++ b/src/junction_trees/ordered_graphs.jl @@ -1,7 +1,7 @@ # An ordered graph (G, σ). struct OrderedGraph - graph::Graph - order::Order + graph::Graph # graph + order::Order # permutation end diff --git a/src/junction_trees/orders.jl b/src/junction_trees/orders.jl index 121ed64..84d3f71 100644 --- a/src/junction_trees/orders.jl +++ b/src/junction_trees/orders.jl @@ -1,7 +1,7 @@ # A permutation σ of the set {1, ..., n}. struct Order <: AbstractVector{Int} - order::Vector{Int} - index::Vector{Int} + order::Vector{Int} # permutation + index::Vector{Int} # inverse permutation end diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl index 5bb5ba6..fee6476 100644 --- a/src/junction_trees/supernode_trees.jl +++ b/src/junction_trees/supernode_trees.jl @@ -1,6 +1,6 @@ # An ordered graph (G, σ) equipped with a supernodal elimination tree T. struct SupernodeTree <: AbstractTree - tree::PostorderTree # tree + tree::PostorderTree # supernodal elimination tree ograph::OrderedGraph # ordered graph representative::Vector{Int} # representative vertex cardinality::Vector{Int} # supernode cardinality diff --git a/src/junction_trees/trees.jl b/src/junction_trees/trees.jl index 1b2aac8..0b26371 100644 --- a/src/junction_trees/trees.jl +++ b/src/junction_trees/trees.jl @@ -1,8 +1,8 @@ # A rooted tree. struct Tree <: AbstractTree - root::Int - parent::Vector{Int} - children::Vector{Vector{Int}} + root::Int # root + parent::Vector{Int} # parent + children::Vector{Vector{Int}} # children end From f02b09c082ec9fe0da3c34a192e4c6221cefb0ff Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sat, 16 Nov 2024 23:49:31 -0500 Subject: [PATCH 03/48] Fixed a bug added another constructor. --- src/Decompositions.jl | 18 +++++++++++++++++- src/junction_trees/ordered_graphs.jl | 8 ++++++++ 2 files changed, 25 insertions(+), 1 deletion(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index 2a9447b..d4868a4 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -6,6 +6,7 @@ export StructuredDecomposition, StrDecomp, ∫ using ..JunctionTrees +using ..JunctionTrees: EliminationAlgorithm, SupernodeType, DEFAULT_ELIMINATION_ALGORITHM, DEFAULT_SUPERNODE_TYPE using PartialFunctions using MLStyle @@ -211,7 +212,22 @@ const MTYPE = StructTightACSetTransformation{ SymmetricGraph} -function Decompositions.StrDecomp(sgraph::AbstractSymmetricGraph, stree::SupernodeTree, seperator::AbstractVector) + +# Construct a tree decomposition of a graph. +function Decompositions.StrDecomp( + sgraph::AbstractSymmetricGraph, + ealg::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, + stype::SupernodeType=DEFAULT_SUPERNODE_TYPE) + + StrDecomp(sgraph, SupernodeTree(sgraph, ealg, stype)) +end + + +# Construct a tree decomposition of a graph. +function Decompositions.StrDecomp(sgraph::AbstractSymmetricGraph, stree::SupernodeTree) + seperator = seperators(stree) + foreach(sort!, seperator) + n = length(stree) graph = Graph(n) objects = Vector{OTYPE}(undef, 2n - 1) diff --git a/src/junction_trees/ordered_graphs.jl b/src/junction_trees/ordered_graphs.jl index 2338126..202e36a 100644 --- a/src/junction_trees/ordered_graphs.jl +++ b/src/junction_trees/ordered_graphs.jl @@ -5,6 +5,14 @@ struct OrderedGraph end +# Given a graph G, construct the ordered graph +# (G, σ), +# where the permutation σ is computed using an elimination algorithm. +function OrderedGraph(sgraph::AbstractSymmetricGraph, ealg::EliminationAlgorithm=DEFAULT_ELIMINATION_ALGORITHM) + OrderedGraph(sgraph, Order(sgraph, ealg)) +end + + # Given a graph G and permutation σ, construct the ordered graph # (G, σ). function OrderedGraph(sgraph::AbstractSymmetricGraph, order::Order) From a96febaa648e7243105dce20334fe517c3a75246 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sun, 17 Nov 2024 00:56:59 -0500 Subject: [PATCH 04/48] temp commit --- src/Decompositions.jl | 8 +- src/junction_trees/elimination_trees.jl | 103 +++--------------------- src/junction_trees/supernode_trees.jl | 66 ++------------- src/junction_trees/supernode_types.jl | 4 +- 4 files changed, 23 insertions(+), 158 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index d4868a4..e611ec2 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -228,7 +228,7 @@ function Decompositions.StrDecomp(sgraph::AbstractSymmetricGraph, stree::Superno seperator = seperators(stree) foreach(sort!, seperator) - n = length(stree) + n = length(stree.tree) graph = Graph(n) objects = Vector{OTYPE}(undef, 2n - 1) morphisms = Vector{MTYPE}(undef, 2n - 2) @@ -236,16 +236,16 @@ function Decompositions.StrDecomp(sgraph::AbstractSymmetricGraph, stree::Superno for i in 1:n snd = supernode(stree, i) sep = seperator[i] - objects[i] = induced_subgraph(sgraph, order(stree, [snd; sep])) + objects[i] = induced_subgraph(sgraph, order(stree.ograph, [snd; sep])) end for i in 1:n - 1 sep = seperator[i] - objects[n + i] = induced_subgraph(sgraph, order(stree, sep)) + objects[n + i] = induced_subgraph(sgraph, order(stree.ograph, sep)) end for i in 1:n - 1 - j = parentindex(stree, i) + j = parentindex(stree.tree, i) add_edge!(graph, i, j) sep_i = seperator[i] diff --git a/src/junction_trees/elimination_trees.jl b/src/junction_trees/elimination_trees.jl index 97ddfb6..ee76266 100644 --- a/src/junction_trees/elimination_trees.jl +++ b/src/junction_trees/elimination_trees.jl @@ -1,6 +1,6 @@ # An ordered graph (G, σ) equipped with the elimination tree T of its elimination graph. # Nodes i in T correspond to vertices σ(i) in G. -struct EliminationTree{T <: AbstractTree} <: AbstractTree +struct EliminationTree{T <: AbstractTree} tree::T # elimination tree ograph::OrderedGraph # ordered graph end @@ -59,7 +59,7 @@ end # Gilbert, Ng, and Peyton # Figure 3: Implementation of algorithm to compute row and column counts. function supcnt(etree::EliminationTree) - order = postorder(etree) + order = postorder(etree.tree) index = inverse(order) rc, cc = supcnt(EliminationTree{PostorderTree}(etree, order)) rc[index], cc[index] @@ -70,7 +70,7 @@ end # Gilbert, Ng, and Peyton # Figure 3: Implementation of algorithm to compute row and column counts. function supcnt(etree::EliminationTree{PostorderTree}) - n = length(etree) + n = length(etree.tree) #### Disjoint Set Union #### @@ -96,22 +96,22 @@ function supcnt(etree::EliminationTree{PostorderTree}) wt = ones(Int, n) for u in 1:n - 1 - wt[parentindex(etree, u)] = 0 + wt[parentindex(etree.tree, u)] = 0 end for p in 1:n - 1 - wt[parentindex(etree, p)] -= 1 + wt[parentindex(etree.tree, p)] -= 1 - for u in outneighbors(etree, p) - if firstdescendant(etree, p) > prev_nbr[u] + for u in outneighbors(etree.ograph, p) + if firstdescendant(etree.tree, p) > prev_nbr[u] wt[p] += 1 pp = prev_p[u] if iszero(pp) - rc[u] += level(etree, p) - level(etree, u) + rc[u] += level(etree.tree, p) - level(etree.tree, u) else q = find(pp) - rc[u] += level(etree, p) - level(etree, q) + rc[u] += level(etree.tree, p) - level(etree.tree, q) wt[q] -= 1 end @@ -121,13 +121,13 @@ function supcnt(etree::EliminationTree{PostorderTree}) prev_nbr[u] = p end - union(p, parentindex(etree, p)) + union(p, parentindex(etree.tree, p)) end cc = wt for v in 1:n - 1 - cc[parentindex(etree, v)] += cc[v] + cc[parentindex(etree.tree, v)] += cc[v] end rc, cc @@ -140,84 +140,3 @@ function outdegrees(etree::EliminationTree) rc, cc = supcnt(etree) cc .- 1 end - - -# Get the number of nodes in T. -function Base.length(etree::EliminationTree) - length(etree.tree) -end - - -# Get the level of a node i. -function level(etree::EliminationTree, i::Integer) - level(etree.tree, i) -end - - -# Get the first descendant of a node i. -function firstdescendant(etree::EliminationTree, i::Integer) - firstdescendant(etree.tree, i) -end - - -# Determine whether a vertex i is a descendant of a node j. -function isdescendant(etree::EliminationTree, i::Integer, j::Integer) - isdescendant(etree.tree, i, j) -end - - -# Get the vertex σ(i). -function order(etree::EliminationTree, i) - order(etree.ograph, i) -end - - -# Get the index σ⁻¹(v), -function inverse(etree::EliminationTree, i) - inverse(etree.ograph, i) -end - - -########################## -# Indexed Tree Interface # -########################## - - -function AbstractTrees.rootindex(etree::EliminationTree) - rootindex(etree.tree) -end - - -function AbstractTrees.parentindex(etree::EliminationTree, i::Integer) - parentindex(etree.tree, i) -end - - -function AbstractTrees.childindices(etree::EliminationTree, i::Integer) - childindices(etree.tree, i) -end - - -function AbstractTrees.NodeType(::Type{IndexNode{EliminationTree, Int}}) - HasNodeType() -end - - -function AbstractTrees.nodetype(::Type{IndexNode{EliminationTree{T}, Int}}) where T - IndexNode{EliminationTree{T}, Int} -end - - -############################ -# Abstract Graph Interface # -############################ - - -function BasicGraphs.inneighbors(etree::EliminationTree, i::Integer) - inneighbors(etree.ograph, i) -end - - -function BasicGraphs.outneighbors(etree::EliminationTree, i::Integer) - outneighbors(etree.ograph, i) -end diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl index fee6476..2402b54 100644 --- a/src/junction_trees/supernode_trees.jl +++ b/src/junction_trees/supernode_trees.jl @@ -1,5 +1,5 @@ # An ordered graph (G, σ) equipped with a supernodal elimination tree T. -struct SupernodeTree <: AbstractTree +struct SupernodeTree tree::PostorderTree # supernodal elimination tree ograph::OrderedGraph # ordered graph representative::Vector{Int} # representative vertex @@ -57,7 +57,7 @@ end # Vanderberghe and Andersen # Algorithm 4.1: Maximal supernodes and supernodal elimination tree. function stree(etree::EliminationTree, degree::AbstractVector, stype::SupernodeType) - n = length(etree) + n = length(etree.tree) index = zeros(Int, n) snd = Vector{Int}[] q = Int[] @@ -78,7 +78,7 @@ function stree(etree::EliminationTree, degree::AbstractVector, stype::SupernodeT push!(snd[i], v) end - for w in childindices(etree, v) + for w in childindices(etree.tree, v) if w !== ww j = index[w] q[j] = i @@ -94,7 +94,7 @@ end # Construct an elimination graph. function eliminationgraph(stree::SupernodeTree) ograph = deepcopy(stree.ograph) - n = length(stree) + n = length(stree.tree) for i in 1:n - 1 for u in supernode(stree, i)[1:end - 1] @@ -108,7 +108,7 @@ function eliminationgraph(stree::SupernodeTree) end u = last(supernode(stree, i)) - v = first(supernode(stree, parentindex(stree, i))) + v = first(supernode(stree, parentindex(stree.tree, i))) for w in outneighbors(ograph, u) if v < w @@ -137,7 +137,7 @@ end # Compute the (unsorted) seperators of every node in T. function seperators(stree::SupernodeTree) - n = length(stree) + n = length(stree.tree) seperator = Vector{Vector{Int}}(undef, n) ograph = eliminationgraph(stree) @@ -149,57 +149,3 @@ function seperators(stree::SupernodeTree) seperator end - - -# Get the number of nodes in T. -function Base.length(stree::SupernodeTree) - length(stree.tree) -end - - -# Get the level of node i. -function level(stree::SupernodeTree, i) - level(stree.tree, i) -end - - -# Get the vertex σ(i). -function order(stree::SupernodeTree, i) - order(stree.ograph, i) -end - - -# Get the index σ⁻¹(v), -function inverse(stree::SupernodeTree, i) - inverse(stree.ograph, i) -end - - -########################## -# Indexed Tree Interface # -########################## - - -function AbstractTrees.rootindex(stree::SupernodeTree) - rootindex(stree.tree) -end - - -function AbstractTrees.parentindex(stree::SupernodeTree, i) - parentindex(stree.tree, i) -end - - -function AbstractTrees.childindices(stree::SupernodeTree, i) - childindices(stree.tree, i) -end - - -function AbstractTrees.NodeType(::Type{IndexNode{SupernodeTree, Int}}) - HasNodeType() -end - - -function AbstractTrees.nodetype(::Type{IndexNode{SupernodeTree, Int}}) - IndexNode{SupernodeTree, Int} -end diff --git a/src/junction_trees/supernode_types.jl b/src/junction_trees/supernode_types.jl index da92a40..b0d3b94 100644 --- a/src/junction_trees/supernode_types.jl +++ b/src/junction_trees/supernode_types.jl @@ -43,7 +43,7 @@ function findchild(etree::EliminationTree, degree::AbstractVector, stype::Node, # v ∈ snd(w). # If no such child exists, return nothing. function findchild(etree::EliminationTree, degree::AbstractVector, stype::Maximal, v::Integer) - for w in childindices(etree, v) + for w in childindices(etree.tree, v) if degree[w] == degree[v] + 1 return w end @@ -55,7 +55,7 @@ end # v ∈ snd(w). # If no such child exists, return nothing. function findchild(etree::EliminationTree, degree::AbstractVector, stype::Fundamental, v::Integer) - ws = childindices(etree, v) + ws = childindices(etree.tree, v) if length(ws) == 1 w = only(ws) From a485644528d2ee904aaa6f6f371a00d44c99cc5a Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sun, 17 Nov 2024 01:02:14 -0500 Subject: [PATCH 05/48] temp commit --- src/JunctionTrees.jl | 1 - src/junction_trees/abstract_trees.jl | 25 ------------------------- src/junction_trees/elimination_trees.jl | 2 +- src/junction_trees/postorder_trees.jl | 10 ++++++++-- src/junction_trees/trees.jl | 18 +++++++++++++++++- 5 files changed, 26 insertions(+), 30 deletions(-) delete mode 100644 src/junction_trees/abstract_trees.jl diff --git a/src/JunctionTrees.jl b/src/JunctionTrees.jl index 1c756b4..e99a8a2 100644 --- a/src/JunctionTrees.jl +++ b/src/JunctionTrees.jl @@ -30,7 +30,6 @@ export SupernodeTree, width, height, seperators, supernode, order include("junction_trees/orders.jl") include("junction_trees/elimination_algorithms.jl") include("junction_trees/ordered_graphs.jl") -include("junction_trees/abstract_trees.jl") include("junction_trees/trees.jl") include("junction_trees/postorder_trees.jl") include("junction_trees/elimination_trees.jl") diff --git a/src/junction_trees/abstract_trees.jl b/src/junction_trees/abstract_trees.jl deleted file mode 100644 index 618e23d..0000000 --- a/src/junction_trees/abstract_trees.jl +++ /dev/null @@ -1,25 +0,0 @@ -# A rooted tree. -abstract type AbstractTree end - - -# Compute a postordering of tree's vertices. -function postorder(tree::AbstractTree) - n = length(tree) - order = Vector{Int}(undef, n) - index = Vector{Int}(undef, n) - - for node in PreOrderDFS(IndexNode(tree)) - order[n] = node.index - index[node.index] = n - n -= 1 - end - - Order(order, index) -end - - -# Get the height of a tree. -function height(tree::AbstractTree) - n = length(tree) - maximum(map(i -> level(tree, i), 1:n)) -end diff --git a/src/junction_trees/elimination_trees.jl b/src/junction_trees/elimination_trees.jl index ee76266..f6cbe18 100644 --- a/src/junction_trees/elimination_trees.jl +++ b/src/junction_trees/elimination_trees.jl @@ -1,6 +1,6 @@ # An ordered graph (G, σ) equipped with the elimination tree T of its elimination graph. # Nodes i in T correspond to vertices σ(i) in G. -struct EliminationTree{T <: AbstractTree} +struct EliminationTree{T <: Union{Tree, PostorderTree}} tree::T # elimination tree ograph::OrderedGraph # ordered graph end diff --git a/src/junction_trees/postorder_trees.jl b/src/junction_trees/postorder_trees.jl index 450282c..319454b 100644 --- a/src/junction_trees/postorder_trees.jl +++ b/src/junction_trees/postorder_trees.jl @@ -1,5 +1,5 @@ # A postordered rooted tree. -struct PostorderTree <: AbstractTree +struct PostorderTree parent::Vector{Int} # parent children::Vector{Vector{Int}} # children level::Vector{Int} # level @@ -36,7 +36,7 @@ end # Postorder a tree. -function PostorderTree(tree::AbstractTree, order::Order) +function PostorderTree(tree::Tree, order::Order) n = length(tree) parent = collect(1:n) @@ -72,6 +72,12 @@ function isdescendant(tree::PostorderTree, i::Integer, j::Integer) end +# Get the height of a tree. +function height(tree::PostorderTree) + maximum(tree.level) +end + + ########################## # Indexed Tree Interface # ########################## diff --git a/src/junction_trees/trees.jl b/src/junction_trees/trees.jl index 0b26371..fcc2069 100644 --- a/src/junction_trees/trees.jl +++ b/src/junction_trees/trees.jl @@ -1,5 +1,5 @@ # A rooted tree. -struct Tree <: AbstractTree +struct Tree root::Int # root parent::Vector{Int} # parent children::Vector{Vector{Int}} # children @@ -35,6 +35,22 @@ function Base.length(tree::Tree) end +# Compute a postordering of tree's vertices. +function postorder(tree::Tree) + n = length(tree) + order = Vector{Int}(undef, n) + index = Vector{Int}(undef, n) + + for node in PreOrderDFS(IndexNode(tree)) + order[n] = node.index + index[node.index] = n + n -= 1 + end + + Order(order, index) +end + + ########################## # Indexed Tree Interface # ########################## From 8ab81371f5f3cac69597756bf9c560bc73f0f1e8 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sun, 17 Nov 2024 01:04:30 -0500 Subject: [PATCH 06/48] Moved things around. --- src/Decompositions.jl | 4 ++-- src/junction_trees/elimination_trees.jl | 16 ++++++++-------- src/junction_trees/supernode_trees.jl | 22 +++++++++++----------- 3 files changed, 21 insertions(+), 21 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index e611ec2..a822164 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -236,12 +236,12 @@ function Decompositions.StrDecomp(sgraph::AbstractSymmetricGraph, stree::Superno for i in 1:n snd = supernode(stree, i) sep = seperator[i] - objects[i] = induced_subgraph(sgraph, order(stree.ograph, [snd; sep])) + objects[i] = induced_subgraph(sgraph, order(stree.graph, [snd; sep])) end for i in 1:n - 1 sep = seperator[i] - objects[n + i] = induced_subgraph(sgraph, order(stree.ograph, sep)) + objects[n + i] = induced_subgraph(sgraph, order(stree.graph, sep)) end for i in 1:n - 1 diff --git a/src/junction_trees/elimination_trees.jl b/src/junction_trees/elimination_trees.jl index f6cbe18..1a87c9a 100644 --- a/src/junction_trees/elimination_trees.jl +++ b/src/junction_trees/elimination_trees.jl @@ -2,7 +2,7 @@ # Nodes i in T correspond to vertices σ(i) in G. struct EliminationTree{T <: Union{Tree, PostorderTree}} tree::T # elimination tree - ograph::OrderedGraph # ordered graph + graph::OrderedGraph # ordered graph end @@ -15,27 +15,27 @@ end # Construct the elimination tree of an ordered graph. -function EliminationTree(ograph::OrderedGraph) - EliminationTree(Tree(etree(ograph)), ograph) +function EliminationTree(graph::OrderedGraph) + EliminationTree(Tree(etree(graph)), graph) end # Postorder an elimination tree. function EliminationTree{PostorderTree}(etree::EliminationTree, order::Order) - EliminationTree(PostorderTree(etree.tree, order), OrderedGraph(etree.ograph, order)) + EliminationTree(PostorderTree(etree.tree, order), OrderedGraph(etree.graph, order)) end # A Compact Row Storage Scheme for Cholesky Factors Using Elimination Trees # Liu # Algorithm 4.2: Elimination Tree by Path Compression. -function etree(ograph::OrderedGraph) - n = nv(ograph) +function etree(graph::OrderedGraph) + n = nv(graph) parent = collect(1:n) ancestor = collect(1:n) for i in 1:n - for k in inneighbors(ograph, i) + for k in inneighbors(graph, i) r = k while ancestor[r] != r && ancestor[r] != i @@ -102,7 +102,7 @@ function supcnt(etree::EliminationTree{PostorderTree}) for p in 1:n - 1 wt[parentindex(etree.tree, p)] -= 1 - for u in outneighbors(etree.ograph, p) + for u in outneighbors(etree.graph, p) if firstdescendant(etree.tree, p) > prev_nbr[u] wt[p] += 1 pp = prev_p[u] diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl index 2402b54..24dbc43 100644 --- a/src/junction_trees/supernode_trees.jl +++ b/src/junction_trees/supernode_trees.jl @@ -1,7 +1,7 @@ # An ordered graph (G, σ) equipped with a supernodal elimination tree T. struct SupernodeTree tree::PostorderTree # supernodal elimination tree - ograph::OrderedGraph # ordered graph + graph::OrderedGraph # ordered graph representative::Vector{Int} # representative vertex cardinality::Vector{Int} # supernode cardinality ancestor::Vector{Int} # first ancestor @@ -28,7 +28,7 @@ function SupernodeTree(etree::EliminationTree, stype::SupernodeType=DEFAULT_SUPE tree = PostorderTree(tree, sorder) order = Order(vcat(snode[sorder]...)) - ograph = OrderedGraph(etree.ograph, order) + graph = OrderedGraph(etree.graph, order) n = length(tree) representative = zeros(Int, n) @@ -49,7 +49,7 @@ function SupernodeTree(etree::EliminationTree, stype::SupernodeType=DEFAULT_SUPE _degree[n] = degree[n] _ancestor[n] = ancestor[n] - SupernodeTree(tree, ograph, representative, cardinality, _ancestor, _degree) + SupernodeTree(tree, graph, representative, cardinality, _ancestor, _degree) end @@ -93,16 +93,16 @@ end # Construct an elimination graph. function eliminationgraph(stree::SupernodeTree) - ograph = deepcopy(stree.ograph) + graph = deepcopy(stree.graph) n = length(stree.tree) for i in 1:n - 1 for u in supernode(stree, i)[1:end - 1] v = u + 1 - for w in outneighbors(ograph, u) + for w in outneighbors(graph, u) if v < w - add_edge!(ograph, v, w) + add_edge!(graph, v, w) end end end @@ -110,14 +110,14 @@ function eliminationgraph(stree::SupernodeTree) u = last(supernode(stree, i)) v = first(supernode(stree, parentindex(stree.tree, i))) - for w in outneighbors(ograph, u) + for w in outneighbors(graph, u) if v < w - add_edge!(ograph, v, w) + add_edge!(graph, v, w) end end end - ograph + graph end @@ -139,10 +139,10 @@ end function seperators(stree::SupernodeTree) n = length(stree.tree) seperator = Vector{Vector{Int}}(undef, n) - ograph = eliminationgraph(stree) + graph = eliminationgraph(stree) for i in 1:n - clique = collect(outneighbors(ograph, stree.representative[i])) + clique = collect(outneighbors(graph, stree.representative[i])) filter!(j -> stree.ancestor[i] <= j, clique) seperator[i] = clique end From 11aabcbef079cbcee809c2a2253fb5155589126b Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sun, 17 Nov 2024 01:13:41 -0500 Subject: [PATCH 07/48] alignment --- src/junction_trees/elimination_trees.jl | 2 +- src/junction_trees/supernode_trees.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/junction_trees/elimination_trees.jl b/src/junction_trees/elimination_trees.jl index 1a87c9a..0ef95dd 100644 --- a/src/junction_trees/elimination_trees.jl +++ b/src/junction_trees/elimination_trees.jl @@ -1,7 +1,7 @@ # An ordered graph (G, σ) equipped with the elimination tree T of its elimination graph. # Nodes i in T correspond to vertices σ(i) in G. struct EliminationTree{T <: Union{Tree, PostorderTree}} - tree::T # elimination tree + tree::T # elimination tree graph::OrderedGraph # ordered graph end diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl index 24dbc43..f47d867 100644 --- a/src/junction_trees/supernode_trees.jl +++ b/src/junction_trees/supernode_trees.jl @@ -1,7 +1,7 @@ # An ordered graph (G, σ) equipped with a supernodal elimination tree T. struct SupernodeTree tree::PostorderTree # supernodal elimination tree - graph::OrderedGraph # ordered graph + graph::OrderedGraph # ordered graph representative::Vector{Int} # representative vertex cardinality::Vector{Int} # supernode cardinality ancestor::Vector{Int} # first ancestor From f929e9e24b89174ed9775b2a9655edf4b607891a Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sun, 17 Nov 2024 10:54:05 -0500 Subject: [PATCH 08/48] Removed unused functions. --- src/JunctionTrees.jl | 2 +- src/junction_trees/ordered_graphs.jl | 26 -------------------------- 2 files changed, 1 insertion(+), 27 deletions(-) diff --git a/src/JunctionTrees.jl b/src/JunctionTrees.jl index e99a8a2..8a3aefc 100644 --- a/src/JunctionTrees.jl +++ b/src/JunctionTrees.jl @@ -24,7 +24,7 @@ export Node, Maximal, Fundamental # Supernode Trees -export SupernodeTree, width, height, seperators, supernode, order +export SupernodeTree, width, height, seperators, supernode, order, inverse include("junction_trees/orders.jl") diff --git a/src/junction_trees/ordered_graphs.jl b/src/junction_trees/ordered_graphs.jl index 202e36a..5400b1d 100644 --- a/src/junction_trees/ordered_graphs.jl +++ b/src/junction_trees/ordered_graphs.jl @@ -52,25 +52,6 @@ function OrderedGraph(ograph::OrderedGraph, order::Order) OrderedGraph(graph, compose(order, ograph.order)) end - -# Construct the elimination graph of an ordered graph. -function eliminationgraph(ograph::OrderedGraph) - ograph = deepcopy(ograph) - - for i in vertices(ograph) - ns = collect(outneighbors(ograph, i)) - n = length(ns) - - for j in 1:n - for k in j + 1:n - add_edge!(ograph, ns[j], ns[k]) - end - end - end - - ograph -end - function Base.deepcopy(ograph::OrderedGraph) order = deepcopy(ograph.order) @@ -116,13 +97,6 @@ function BasicGraphs.edges(ograph::OrderedGraph) end -function BasicGraphs.has_edge(ograph::OrderedGraph, i, j) - u = min(i, j) - v = max(i, j) - has_edge(ograph.graph, u, v) -end - - function BasicGraphs.add_edge!(ograph::OrderedGraph, i, j) u = min(i, j) v = max(i, j) From d3794e65ad8be7ef88faea954e423ef96de9704d Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sun, 17 Nov 2024 13:54:05 -0500 Subject: [PATCH 09/48] temp commit --- src/Decompositions.jl | 3 +- src/JunctionTrees.jl | 2 +- src/junction_trees/elimination_trees.jl | 29 ---- src/junction_trees/ordered_graphs.jl | 31 +++- src/junction_trees/supernode_trees.jl | 61 ++------ src/junction_trees/supernode_types.jl | 38 +++++ test/JunctionTrees.jl | 193 +++++++++++++++++++++++- 7 files changed, 267 insertions(+), 90 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index a822164..271aefc 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -207,12 +207,11 @@ const MTYPE = StructTightACSetTransformation{ (nothing, :E, :E, ((:inv, :inv), ())), (nothing, :E, :V, ((:inv, :src), (:tgt,))), (nothing, :E, :V, ((:inv, :tgt), (:src,)))}}, - @NamedTuple{V::FinDomFunctionVector{Int64, Vector{Int64}, FinSetInt}, E::FinDomFunctionVector{Int64, Vector{Int64}, FinSetInt}}, + @NamedTuple{V::FinDomFunctionVector{Int, Vector{Int}, FinSetInt}, E::FinDomFunctionVector{Int, Vector{Int}, FinSetInt}}, SymmetricGraph, SymmetricGraph} - # Construct a tree decomposition of a graph. function Decompositions.StrDecomp( sgraph::AbstractSymmetricGraph, diff --git a/src/JunctionTrees.jl b/src/JunctionTrees.jl index 8a3aefc..dac274c 100644 --- a/src/JunctionTrees.jl +++ b/src/JunctionTrees.jl @@ -24,7 +24,7 @@ export Node, Maximal, Fundamental # Supernode Trees -export SupernodeTree, width, height, seperators, supernode, order, inverse +export SupernodeTree, width, height, seperators, supernode, permutation, inverse include("junction_trees/orders.jl") diff --git a/src/junction_trees/elimination_trees.jl b/src/junction_trees/elimination_trees.jl index 0ef95dd..55c250e 100644 --- a/src/junction_trees/elimination_trees.jl +++ b/src/junction_trees/elimination_trees.jl @@ -26,35 +26,6 @@ function EliminationTree{PostorderTree}(etree::EliminationTree, order::Order) end -# A Compact Row Storage Scheme for Cholesky Factors Using Elimination Trees -# Liu -# Algorithm 4.2: Elimination Tree by Path Compression. -function etree(graph::OrderedGraph) - n = nv(graph) - parent = collect(1:n) - ancestor = collect(1:n) - - for i in 1:n - for k in inneighbors(graph, i) - r = k - - while ancestor[r] != r && ancestor[r] != i - t = ancestor[r] - ancestor[r] = i - r = t - end - - if ancestor[r] == r - ancestor[r] = i - parent[r] = i - end - end - end - - parent -end - - # An Efficient Algorithm to Compute Row and Column Counts for Sparse Cholesky Factorization # Gilbert, Ng, and Peyton # Figure 3: Implementation of algorithm to compute row and column counts. diff --git a/src/junction_trees/ordered_graphs.jl b/src/junction_trees/ordered_graphs.jl index 5400b1d..0b28da5 100644 --- a/src/junction_trees/ordered_graphs.jl +++ b/src/junction_trees/ordered_graphs.jl @@ -52,6 +52,35 @@ function OrderedGraph(ograph::OrderedGraph, order::Order) OrderedGraph(graph, compose(order, ograph.order)) end + +# A Compact Row Storage Scheme for Cholesky Factors Using Elimination Trees +# Liu +# Algorithm 4.2: Elimination Tree by Path Compression. +function etree(graph::OrderedGraph) + n = nv(graph) + parent = collect(1:n) + ancestor = collect(1:n) + + for i in 1:n + for k in inneighbors(graph, i) + r = k + + while ancestor[r] != r && ancestor[r] != i + t = ancestor[r] + ancestor[r] = i + r = t + end + + if ancestor[r] == r + ancestor[r] = i + parent[r] = i + end + end + end + + parent +end + function Base.deepcopy(ograph::OrderedGraph) order = deepcopy(ograph.order) @@ -61,7 +90,7 @@ end # Get the vertex σ(i). -function order(ograph::OrderedGraph, i) +function permutation(ograph::OrderedGraph, i) ograph.order[i] end diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl index f47d867..e1cfd15 100644 --- a/src/junction_trees/supernode_trees.jl +++ b/src/junction_trees/supernode_trees.jl @@ -5,7 +5,7 @@ struct SupernodeTree representative::Vector{Int} # representative vertex cardinality::Vector{Int} # supernode cardinality ancestor::Vector{Int} # first ancestor - degrees::Vector{Int} # higher degrees + degree::Vector{Int} # higher degrees end @@ -30,64 +30,25 @@ function SupernodeTree(etree::EliminationTree, stype::SupernodeType=DEFAULT_SUPE order = Order(vcat(snode[sorder]...)) graph = OrderedGraph(etree.graph, order) + permute!(snode, sorder) + permute!(ancestor, sorder) + permute!(degree, order) + n = length(tree) representative = zeros(Int, n) cardinality = zeros(Int, n) - _ancestor = zeros(Int, n) - _degree = zeros(Int, n) - + for i in 1:n - 1 - j = sorder[i] - representative[i] = inverse(order, snode[j][1]) - cardinality[i] = length(snode[j]) - _degree[i] = degree[j] - _ancestor[i] = inverse(order, ancestor[j]) + representative[i] = inverse(order, snode[i][1]) + cardinality[i] = length(snode[i]) + ancestor[i] = inverse(order, ancestor[i]) end representative[n] = inverse(order, snode[n][1]) cardinality[n] = length(snode[n]) - _degree[n] = degree[n] - _ancestor[n] = ancestor[n] - - SupernodeTree(tree, graph, representative, cardinality, _ancestor, _degree) -end - - -# Chordal Graphs and Semidefinite Optimization -# Vanderberghe and Andersen -# Algorithm 4.1: Maximal supernodes and supernodal elimination tree. -function stree(etree::EliminationTree, degree::AbstractVector, stype::SupernodeType) - n = length(etree.tree) - index = zeros(Int, n) - snd = Vector{Int}[] - q = Int[] - a = Int[] - - for v in 1:n - ww = findchild(etree, degree, stype, v) - - if isnothing(ww) - i = length(snd) + 1 - index[v] = i - push!(snd, [v]) - push!(q, length(snd)) - push!(a, n + 1) - else - i = index[ww] - index[v] = i - push!(snd[i], v) - end - - for w in childindices(etree.tree, v) - if w !== ww - j = index[w] - q[j] = i - a[j] = v - end - end - end + ancestor[n] = ancestor[n] - snd, q, a + SupernodeTree(tree, graph, representative, cardinality, ancestor, degree) end diff --git a/src/junction_trees/supernode_types.jl b/src/junction_trees/supernode_types.jl index b0d3b94..29158e6 100644 --- a/src/junction_trees/supernode_types.jl +++ b/src/junction_trees/supernode_types.jl @@ -33,6 +33,44 @@ A fundamental supernode. struct Fundamental <: SupernodeType end +# Chordal Graphs and Semidefinite Optimization +# Vanderberghe and Andersen +# Algorithm 4.1: Maximal supernodes and supernodal elimination tree. +function stree(etree::EliminationTree, degree::AbstractVector, stype::SupernodeType) + n = length(etree.tree) + index = zeros(Int, n) + snd = Vector{Int}[] + q = Int[] + a = Int[] + + for v in 1:n + ww = findchild(etree, degree, stype, v) + + if isnothing(ww) + i = length(snd) + 1 + index[v] = i + push!(snd, [v]) + push!(q, length(snd)) + push!(a, n + 1) + else + i = index[ww] + index[v] = i + push!(snd[i], v) + end + + for w in childindices(etree.tree, v) + if w !== ww + j = index[w] + q[j] = i + a[j] = v + end + end + end + + snd, q, a +end + + # Find a child w of v such that # v ∈ snd(w). # If no such child exists, return nothing. diff --git a/test/JunctionTrees.jl b/test/JunctionTrees.jl index 763dab5..715fdcd 100644 --- a/test/JunctionTrees.jl +++ b/test/JunctionTrees.jl @@ -1,13 +1,15 @@ +using StructuredDecompositions.JunctionTrees + + using AbstractTrees using Catlab.BasicGraphs using Catlab.RelationalPrograms using Catlab.UndirectedWiringDiagrams -using StructuredDecompositions.JunctionTrees using Test -# Vandenberghe and Andersen # Chordal Graphs and Semidefinite Optimization +# Vandenberghe and Andersen graph = SymmetricGraph(17) add_edges!(graph, @@ -26,14 +28,191 @@ order = JunctionTrees.Order(graph, MetisJL_ND()) order = JunctionTrees.Order(1:17) @test length(order) == 17 +# Figure 4.3 +stree = SupernodeTree(graph, order, Node()) +@test width(stree) == 4 +@test length(stree.tree) == 17 +@test height(stree.tree) == 7 + +@test map(i -> inverse(stree.graph, i), 1:17) == [ + 5, # a + 4, # b + 6, # c + 7, # d + 8, # e + 3, # f + 1, # g + 2, # h + 9, # i + 12, # j + 13, # k + 11, # l + 14, # m + 15, # n + 10, # o + 16, # p + 17, # q +] + +@test map(i -> parentindex(stree.tree, i), 1:17) == [ + 2, + 9, + 9, + 6, + 6, + 7, + 8, + 9, + 10, + 16, + 14, + 13, + 14, + 15, + 16, + 17, + nothing, +] + +@test map(i -> childindices(stree.tree, i), 1:17) == [ + [], + [1], + [], + [], + [], + [4, 5], + [6], + [7], + [2, 3, 8], + [9], + [], + [], + [12], + [11, 13], + [14], + [10, 15], + [16], +] + +@test map(i -> supernode(stree, i), 1:17) == [ + [1], # g + [2], # h + [3], # f + [4], # b + [5], # a + [6], # c + [7], # d + [8], # e + [9], # i + [10], # o + [11], # l + [12], # j + [13], # k + [14], # m + [15], # n + [16], # p + [17], # q +] + +@test map(sort, seperators(stree)) == [ + [2, 9, 10], # h i o + [9, 10], # i o + [9, 16], # i p + [6, 7], # c d + [6, 7, 8, 10], # c d e o + [7, 8, 10], # d e o + [8, 10], # e o + [9, 10, 16], # i o p + [10, 16], # o p + [16, 17], # p q + [14, 15, 16, 17], # m n p q + [13, 14, 15, 17], # k m n q + [14, 15, 17], # m n q + [15, 16, 17], # n p q + [16, 17], # p q + [17], # q + [], # +] -# Figure 4.2 -etree = EliminationTree(graph, order) -@test parentindex.([etree], 1:17) == [3, 3, 4, 5, 9, 9, 8, 9, 15, 11, 13, 13, 14, 16, 16, 17, 17] + +# Figure 4.7 (left) +stree = SupernodeTree(graph, order, Maximal()) +@test width(stree) == 4 +@test length(stree.tree) == 8 +@test height(stree.tree) == 4 + +@test map(i -> inverse(stree.graph, i), 1:17) == [ + 5, # a + 4, # b + 6, # c + 7, # d + 8, # e + 3, # f + 1, # g + 2, # h + 9, # i + 11, # j + 12, # k + 13, # l + 14, # m + 15, # n + 10, # o + 16, # p + 17, # q +] + +@test map(i -> parentindex(stree.tree, i), 1:8) == [ + 5, + 5, + 4, + 5, + 6, + 8, + 8, + nothing +] + +@test map(i -> childindices(stree.tree, i), 1:8) == [ + [], + [], + [], + [3], + [1, 2, 4], + [5], + [], + [6, 7], +] + + +@test map(i -> supernode(stree, i), 1:8) == [ + [1, 2], # g h + [3], # f + [4], # b + [5, 6, 7], # a c d + [8, 9], # e i + [10], # o + [11, 12], # j k + [13, 14, 15, 16, 17], # l m n p q +] + + +@test map(sort, seperators(stree)) == [ + [9, 10], # i o + [9, 16], # i p + [6, 7], # c d + [8, 10], # e o + [10, 16], # o p + [16, 17], # p q + [15, 15, 17], # m n q + [], # +] -# Figure 4.3 -stree = SupernodeTree(etree, Node()) +# Figure 4.9 +stree = SupernodeTree(graph, order, Fundamental()) +@test width(stree) == 4 +@test length(stree.tree) == 12 +@test height(stree.tree) == 5 #= From 0d008c8159bca357d2980d1e9e5776e5c5cc14bb Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sun, 17 Nov 2024 14:18:27 -0500 Subject: [PATCH 10/48] Tests passing. --- src/junction_trees/supernode_trees.jl | 24 +- test/JunctionTrees.jl | 318 ++++++-------------------- 2 files changed, 78 insertions(+), 264 deletions(-) diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl index e1cfd15..ca88be9 100644 --- a/src/junction_trees/supernode_trees.jl +++ b/src/junction_trees/supernode_trees.jl @@ -24,29 +24,29 @@ function SupernodeTree(etree::EliminationTree, stype::SupernodeType=DEFAULT_SUPE degree = outdegrees(etree) snode, parent, ancestor = stree(etree, degree, stype) tree = Tree(parent) - sorder = postorder(tree) - tree = PostorderTree(tree, sorder) - - order = Order(vcat(snode[sorder]...)) - graph = OrderedGraph(etree.graph, order) - permute!(snode, sorder) - permute!(ancestor, sorder) - permute!(degree, order) + treeorder = postorder(tree) + graphorder = Order(vcat(snode[treeorder]...)) + + tree = PostorderTree(tree, treeorder) + graph = OrderedGraph(etree.graph, graphorder) + + permute!(snode, treeorder) + permute!(ancestor, treeorder) + permute!(degree, graphorder) n = length(tree) representative = zeros(Int, n) cardinality = zeros(Int, n) for i in 1:n - 1 - representative[i] = inverse(order, snode[i][1]) + representative[i] = inverse(graphorder, snode[i][1]) cardinality[i] = length(snode[i]) - ancestor[i] = inverse(order, ancestor[i]) + ancestor[i] = inverse(graphorder, ancestor[i]) end - representative[n] = inverse(order, snode[n][1]) + representative[n] = inverse(graphorder, snode[n][1]) cardinality[n] = length(snode[n]) - ancestor[n] = ancestor[n] SupernodeTree(tree, graph, representative, cardinality, ancestor, degree) end diff --git a/test/JunctionTrees.jl b/test/JunctionTrees.jl index 715fdcd..f77e4e8 100644 --- a/test/JunctionTrees.jl +++ b/test/JunctionTrees.jl @@ -183,7 +183,6 @@ stree = SupernodeTree(graph, order, Maximal()) [6, 7], ] - @test map(i -> supernode(stree, i), 1:8) == [ [1, 2], # g h [3], # f @@ -195,7 +194,6 @@ stree = SupernodeTree(graph, order, Maximal()) [13, 14, 15, 16, 17], # l m n p q ] - @test map(sort, seperators(stree)) == [ [9, 10], # i o [9, 16], # i p @@ -203,11 +201,10 @@ stree = SupernodeTree(graph, order, Maximal()) [8, 10], # e o [10, 16], # o p [16, 17], # p q - [15, 15, 17], # m n q + [14, 15, 17], # m n q [], # ] - # Figure 4.9 stree = SupernodeTree(graph, order, Fundamental()) @test width(stree) == 4 @@ -215,265 +212,82 @@ stree = SupernodeTree(graph, order, Fundamental()) @test height(stree.tree) == 5 -#= -@test supernode.([stree], order(stree, 1:17)) - - - -parent = JunctionTrees.makeetree(graph, order) - -# Figure 4.2 -@test parent == [3, 3, 4, 5, 9, 9, 8, 9, 15, 11, 13, 13, 14, 16, 16, 17, 0] - -etree = JunctionTrees.Tree(17, parent) -indegrees, outdegrees = JunctionTrees.getdegrees(graph, order, etree) - -@test indegrees == [0, 0, 2, 3, 3, 0, 0, 1, 4, 0, 1, 0, 3, 4, 7, 7, 7] -@test outdegrees == [4, 2, 3, 2, 3, 2, 3, 2, 2, 4, 3, 4, 3, 2, 2, 1, 0] - -# Figure 4.3 -snd, sbt, q, a = JunctionTrees.makestree(etree, outdegrees, Node()) - -@test snd == [ - [1], - [2], - [3], - [4], - [5], - [6], - [7], - [8], - [9], - [10], - [11], - [12], - [13], - [14], - [15], - [16], - [17]] - -@test sbt == 1:17 -@test q == [3, 3, 4, 5, 9, 9, 8, 9, 15, 11, 13, 13, 14, 16, 16, 17, 0] -@test a == [3, 3, 4, 5, 9, 9, 8, 9, 15, 11, 13, 13, 14, 16, 16, 17, 0] - -# Figure 4.7 (left) -snd, sbt, q, a = JunctionTrees.makestree(etree, outdegrees, MaximalSupernode()) - -@test snd == [ - [1, 3, 4], - [2], - [5, 9], - [6], - [7, 8], - [10, 11], - [12, 13, 14, 16, 17], - [15] ] - -@test sbt == [1, 2, 1, 1, 3, 4, 5, 5, 3, 6, 6, 7, 7, 7, 8, 7, 7] -@test q == [3, 1, 8, 3, 3, 7, 0, 7] -@test a == [5, 3, 15, 9, 9, 13, 0, 16] - -# Figure 4.9 -snd, sbt, q, a = JunctionTrees.makestree(etree, outdegrees, FundamentalSupernode()) - -@test snd == [ - [1], - [2], - [3, 4], - [5], - [6], - [7, 8], - [9], - [10, 11], - [12], - [13, 14], - [15], - [16, 17] ] - -@test sbt == [1, 2, 3, 3, 4, 5, 6, 6, 7, 8, 8, 9, 10, 10, 11, 12, 12] -@test q == [3, 3, 4, 7, 7, 7, 11, 10, 10, 12, 12, 0] -@test a == [3, 3, 5, 9, 9, 9, 15, 13, 13, 16, 16, 0] - -# Figure 4.3 -jtree = JunctionTree(graph, order, Node()) - -@test getresidual.([jtree], getsubtree.([jtree], 1:17)) == [ - [1], - [2], - [3], - [4], - [5], - [6], - [7], - [8], - [9], - [10], - [11], - [12], - [13], - [14], - [15], - [16], - [17] ] - -@test getseperator.([jtree], getsubtree.([jtree], 1:17)) == [ - [3, 4, 5, 15], - [3, 4], - [4, 5, 15], - [5, 15], - [9, 15, 16], - [9, 16], - [8, 9, 15], - [9, 15], - [15, 16], - [11, 13, 14, 17], - [13, 14, 17], - [13, 14, 16, 17], - [14, 16, 17], - [16, 17], - [16, 17], - [17], - [] ] +@test map(i -> inverse(stree.graph, i), 1:17) == [ + 5, # a + 4, # b + 6, # c + 7, # d + 8, # e + 3, # f + 1, # g + 2, # h + 9, # i + 12, # j + 13, # k + 11, # l + 14, # m + 15, # n + 10, # o + 16, # p + 17, # q +] -@test getlevel.([jtree], getsubtree.([jtree], 1:17)) == [ +@test map(i -> parentindex(stree.tree, i), 1:12) == [ 7, 7, - 6, 5, - 4, - 4, - 5, - 4, - 3, 5, - 4, - 4, - 3, - 2, - 2, - 1, - 0 ] - -@test isdescendant(jtree, getsubtree(jtree, 5), getsubtree(jtree, 15)) -@test !isdescendant(jtree, getsubtree(jtree, 15), getsubtree(jtree, 5)) -@test !isdescendant(jtree, getsubtree(jtree, 10), getsubtree(jtree, 15)) -@test !isdescendant(jtree, getsubtree(jtree, 15), getsubtree(jtree, 10)) -@test !isdescendant(jtree, getsubtree(jtree, 1), getsubtree(jtree, 1)) -@test getwidth(jtree) == 4 - -# Figure 4.7 (left) -jtree = JunctionTree(graph, order, MaximalSupernode()) - -@test getresidual.([jtree], getsubtree.([jtree], 1:17)) == [ - [1, 3, 4], - [2], - [1, 3, 4], - [1, 3, 4], - [5, 9], - [6], - [7, 8], - [7, 8], - [5, 9], - [10, 11], - [10, 11], - [12, 13, 14, 16, 17], - [12, 13, 14, 16, 17], - [12, 13, 14, 16, 17], - [15], - [12, 13, 14, 16, 17], - [12, 13, 14, 16, 17]] + 6, + 7, + 8, + 12, + 11, + 11, + 12, + nothing, +] -@test getseperator.([jtree], getsubtree.([jtree], 1:17)) == [ - [5, 15], - [3, 4], - [5, 15], - [5, 15], - [15, 16], - [9, 16], - [9, 15], - [9, 15], - [15, 16], - [13, 14, 17], - [13, 14, 17], +@test map(i -> childindices(stree.tree, i), 1:12) == [ [], [], [], - [16, 17], [], - []] - -@test getlevel.([jtree], getsubtree.([jtree], 1:17)) == [ - 3, - 4, - 3, - 3, - 2, - 3, - 3, - 3, - 2, - 1, - 1, - 0, - 0, - 0, - 1, - 0, - 0 ] - -@test isdescendant(jtree, getsubtree(jtree, 5), getsubtree(jtree, 15)) -@test !isdescendant(jtree, getsubtree(jtree, 15), getsubtree(jtree, 5)) -@test !isdescendant(jtree, getsubtree(jtree, 10), getsubtree(jtree, 15)) -@test !isdescendant(jtree, getsubtree(jtree, 15), getsubtree(jtree, 10)) -@test !isdescendant(jtree, getsubtree(jtree, 1), getsubtree(jtree, 1)) -@test getwidth(jtree) == 4 - -# Figure 4.9 -jtree = JunctionTree(graph, order, FundamentalSupernode()) - -@test getresidual.([jtree], getsubtree.([jtree], 1:17)) == [ - [1], - [2], - [3, 4], [3, 4], [5], - [6], - [7, 8], - [7, 8], - [9], - [10, 11], - [10, 11], - [12], - [13, 14], - [13, 14], - [15], - [16, 17], - [16, 17]] - -@test getseperator.([jtree], getsubtree.([jtree], 1:17)) == [ - [3, 4, 5, 15], - [3, 4], - [5, 15], - [5, 15], - [9, 15, 16], - [9, 16], - [9, 15], - [9, 15], - [15, 16], - [13, 14, 17], - [13, 14, 17], - [13, 14, 16, 17], - [16, 17], - [16, 17], - [16, 17], + [1, 2, 6], + [7], [], - []] + [], + [9, 10], + [8, 11], +] + +@test map(i -> supernode(stree, i), 1:12) == [ + [1, 2], # g h + [3], # f + [4], # b + [5], # a + [6, 7], # c d + [8], # e + [9], # i + [10], # o + [11], # l + [12, 13], # j k + [14, 15], # m n + [16, 17], # p q +] -@test isdescendant(jtree, getsubtree(jtree, 5), getsubtree(jtree, 15)) -@test !isdescendant(jtree, getsubtree(jtree, 15), getsubtree(jtree, 5)) -@test !isdescendant(jtree, getsubtree(jtree, 10), getsubtree(jtree, 15)) -@test !isdescendant(jtree, getsubtree(jtree, 15), getsubtree(jtree, 10)) -@test !isdescendant(jtree, getsubtree(jtree, 1), getsubtree(jtree, 1)) -@test getwidth(jtree) == 4 -=# +@test map(sort, seperators(stree)) == [ + [9, 10], # i o + [9, 16], # i p + [6, 7], # c d + [6, 7, 8, 10], # c d e o + [8, 10], # e o + [9, 10, 16], # i o p + [10, 16], # o p + [16, 17], # p q + [14, 15, 16, 17], # m n p q + [14, 15, 17], # m n q + [16, 17], # p q + [], # +] From abf6ce051f09dd1105f59912a63ce22f824acd9f Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sun, 17 Nov 2024 14:20:25 -0500 Subject: [PATCH 11/48] Fixed bug in function name. --- src/Decompositions.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index 271aefc..1ec428a 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -214,11 +214,11 @@ const MTYPE = StructTightACSetTransformation{ # Construct a tree decomposition of a graph. function Decompositions.StrDecomp( - sgraph::AbstractSymmetricGraph, + graph::AbstractSymmetricGraph, ealg::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, stype::SupernodeType=DEFAULT_SUPERNODE_TYPE) - StrDecomp(sgraph, SupernodeTree(sgraph, ealg, stype)) + StrDecomp(graph, SupernodeTree(graph, ealg, stype)) end @@ -235,12 +235,12 @@ function Decompositions.StrDecomp(sgraph::AbstractSymmetricGraph, stree::Superno for i in 1:n snd = supernode(stree, i) sep = seperator[i] - objects[i] = induced_subgraph(sgraph, order(stree.graph, [snd; sep])) + objects[i] = induced_subgraph(sgraph, permutation(stree.graph, [snd; sep])) end for i in 1:n - 1 sep = seperator[i] - objects[n + i] = induced_subgraph(sgraph, order(stree.graph, sep)) + objects[n + i] = induced_subgraph(sgraph, permutation(stree.graph, sep)) end for i in 1:n - 1 From 41e6a4e0a1f88cad14b51ab924d654c4608c4daf Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sun, 17 Nov 2024 17:22:10 -0500 Subject: [PATCH 12/48] temp commit --- src/junction_trees/elimination_algorithms.jl | 12 ++++++-- src/junction_trees/supernode_trees.jl | 31 ++++++++++++++++++++ 2 files changed, 40 insertions(+), 3 deletions(-) diff --git a/src/junction_trees/elimination_algorithms.jl b/src/junction_trees/elimination_algorithms.jl index f87dc5b..9d273e9 100644 --- a/src/junction_trees/elimination_algorithms.jl +++ b/src/junction_trees/elimination_algorithms.jl @@ -34,9 +34,15 @@ The nested dissection heuristic. Uses Metis.jl. struct MetisJL_ND <: EliminationAlgorithm end +# Construct an order using the default graph elimination algorithm. +function Order(graph::AbstractSymmetricGraph) + Order(graph, DEFAULT_ELIMINATION_ALGORITHM) +end + + # Construct an order using the reverse Cuthill-McKee algorithm. Uses # CuthillMcKee.jl. -function Order(graph::AbstractSymmetricGraph, ::CuthillMcKeeJL_RCM) +function Order(graph::AbstractSymmetricGraph, ealg::CuthillMcKeeJL_RCM) order = CuthillMcKee.symrcm(adjacencymatrix(graph)) Order(order) end @@ -44,14 +50,14 @@ end # Construct an order using the approximate minimum degree algorithm. Uses # AMD.jl. -function Order(graph::AbstractSymmetricGraph, ::AMDJL_AMD) +function Order(graph::AbstractSymmetricGraph, ealg::AMDJL_AMD) order = AMD.symamd(adjacencymatrix(graph)) Order(order) end # Construct an order using the nested dissection heuristic. Uses Metis.jl. -function Order(graph::AbstractSymmetricGraph, ::MetisJL_ND) +function Order(graph::AbstractSymmetricGraph, ealg::MetisJL_ND) order, index = Metis.permutation(adjacencymatrix(graph)) Order(order, index) end diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl index ca88be9..4f88e87 100644 --- a/src/junction_trees/supernode_trees.jl +++ b/src/junction_trees/supernode_trees.jl @@ -110,3 +110,34 @@ function seperators(stree::SupernodeTree) seperator end + + +# Compute the (unsorted) seperators of every node in T. +function _seperators(stree::SupernodeTree) + n = length(stree.tree) + seperator = Vector{SparseIntSet}(undef, n) + + for i in 1:n - 1 + seperator[i] = SparseIntSet(stree.ancestor[i]) + end + + seperator[n] = SparseIntSet() + + for i in 1:n - 1 + for v in outneighbors(stree.graph, stree.representative[i]) + if stree.ancestor[i] < v + push!(seperator[i], v) + end + end + + j = parentindex(stree.tree, i) + + for v in seperator[i] + if stree.ancestor[j] < v + push!(seperator[j], v) + end + end + end + + map(collect, seperator) +end From 9fe79bc97352807b2dcf6c1f6333ed916e1ab783 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sun, 17 Nov 2024 17:36:10 -0500 Subject: [PATCH 13/48] Improved performance of `seperators` function. --- src/junction_trees/supernode_trees.jl | 60 ++++----------------------- test/JunctionTrees.jl | 6 +-- 2 files changed, 10 insertions(+), 56 deletions(-) diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl index 4f88e87..8fb3c98 100644 --- a/src/junction_trees/supernode_trees.jl +++ b/src/junction_trees/supernode_trees.jl @@ -52,36 +52,6 @@ function SupernodeTree(etree::EliminationTree, stype::SupernodeType=DEFAULT_SUPE end -# Construct an elimination graph. -function eliminationgraph(stree::SupernodeTree) - graph = deepcopy(stree.graph) - n = length(stree.tree) - - for i in 1:n - 1 - for u in supernode(stree, i)[1:end - 1] - v = u + 1 - - for w in outneighbors(graph, u) - if v < w - add_edge!(graph, v, w) - end - end - end - - u = last(supernode(stree, i)) - v = first(supernode(stree, parentindex(stree.tree, i))) - - for w in outneighbors(graph, u) - if v < w - add_edge!(graph, v, w) - end - end - end - - graph -end - - # Compute the width of a supernodal elimination tree. function width(stree::SupernodeTree) maximum(stree.degree[stree.representative]) @@ -99,37 +69,19 @@ end # Compute the (unsorted) seperators of every node in T. function seperators(stree::SupernodeTree) n = length(stree.tree) - seperator = Vector{Vector{Int}}(undef, n) - graph = eliminationgraph(stree) - - for i in 1:n - clique = collect(outneighbors(graph, stree.representative[i])) - filter!(j -> stree.ancestor[i] <= j, clique) - seperator[i] = clique - end - - seperator -end - - -# Compute the (unsorted) seperators of every node in T. -function _seperators(stree::SupernodeTree) - n = length(stree.tree) - seperator = Vector{SparseIntSet}(undef, n) + seperator = Vector{Set{Int}}(undef, n) for i in 1:n - 1 - seperator[i] = SparseIntSet(stree.ancestor[i]) - end + seperator[i] = Set(stree.ancestor[i]) - seperator[n] = SparseIntSet() - - for i in 1:n - 1 for v in outneighbors(stree.graph, stree.representative[i]) if stree.ancestor[i] < v push!(seperator[i], v) end end + end + for i in 1:n - 2 j = parentindex(stree.tree, i) for v in seperator[i] @@ -137,7 +89,9 @@ function _seperators(stree::SupernodeTree) push!(seperator[j], v) end end + end - map(collect, seperator) + seperator[n] = Set() + seperator end diff --git a/test/JunctionTrees.jl b/test/JunctionTrees.jl index f77e4e8..f3829a5 100644 --- a/test/JunctionTrees.jl +++ b/test/JunctionTrees.jl @@ -114,7 +114,7 @@ stree = SupernodeTree(graph, order, Node()) [17], # q ] -@test map(sort, seperators(stree)) == [ +@test map(sort ∘ collect, seperators(stree)) == [ [2, 9, 10], # h i o [9, 10], # i o [9, 16], # i p @@ -194,7 +194,7 @@ stree = SupernodeTree(graph, order, Maximal()) [13, 14, 15, 16, 17], # l m n p q ] -@test map(sort, seperators(stree)) == [ +@test map(sort ∘ collect, seperators(stree)) == [ [9, 10], # i o [9, 16], # i p [6, 7], # c d @@ -277,7 +277,7 @@ stree = SupernodeTree(graph, order, Fundamental()) [16, 17], # p q ] -@test map(sort, seperators(stree)) == [ +@test map(sort ∘ collect, seperators(stree)) == [ [9, 10], # i o [9, 16], # i p [6, 7], # c d From 8c48760e397479ce3f4eaf5bd5e02830f244271a Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sun, 17 Nov 2024 18:06:34 -0500 Subject: [PATCH 14/48] Added back MCS --- src/JunctionTrees.jl | 6 +- src/junction_trees/elimination_algorithms.jl | 65 +++++++++++++++++++- test/JunctionTrees.jl | 3 + 3 files changed, 70 insertions(+), 4 deletions(-) diff --git a/src/JunctionTrees.jl b/src/JunctionTrees.jl index dac274c..b200de6 100644 --- a/src/JunctionTrees.jl +++ b/src/JunctionTrees.jl @@ -12,11 +12,11 @@ using SparseArrays # Orders -export Order +export Order, inverse # Elimination Algorithms -export AMDJL_AMD, CuthillMcKeeJL_RCM, MetisJL_ND +export AMDJL_AMD, CuthillMcKeeJL_RCM, MetisJL_ND, MCS # Supernode Types @@ -24,7 +24,7 @@ export Node, Maximal, Fundamental # Supernode Trees -export SupernodeTree, width, height, seperators, supernode, permutation, inverse +export SupernodeTree, width, height, seperators, supernode, permutation include("junction_trees/orders.jl") diff --git a/src/junction_trees/elimination_algorithms.jl b/src/junction_trees/elimination_algorithms.jl index 9d273e9..79a2b77 100644 --- a/src/junction_trees/elimination_algorithms.jl +++ b/src/junction_trees/elimination_algorithms.jl @@ -5,7 +5,7 @@ A graph elimination algorithm. The options are - [`CuthillMcKeeJL_RCM`](@ref) - [`AMDJL_AMD`](@ref) - [`MetisJL_ND`](@ref) -- [`Order`](@ref) +- [`MCS`](@ref) """ abstract type EliminationAlgorithm end @@ -34,6 +34,14 @@ The nested dissection heuristic. Uses Metis.jl. struct MetisJL_ND <: EliminationAlgorithm end +""" + MCS <: EliminationAlgorithm + +The maximum cardinality search algorithm. +""" +struct MCS <: EliminationAlgorithm end + + # Construct an order using the default graph elimination algorithm. function Order(graph::AbstractSymmetricGraph) Order(graph, DEFAULT_ELIMINATION_ALGORITHM) @@ -63,6 +71,13 @@ function Order(graph::AbstractSymmetricGraph, ealg::MetisJL_ND) end +# Concstruct an order using the maximum cardinality search algorithm. +function Order(graph::AbstractSymmetricGraph, ealg::MCS) + order, index = mcs(graph) + Order(order, index) +end + + # Construct the adjacency matrix of a graph. function adjacencymatrix(graph::AbstractSymmetricGraph) m = ne(graph) @@ -83,4 +98,52 @@ function adjacencymatrix(graph::AbstractSymmetricGraph) 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{Set{Int}}(undef, n) + + for i in 1:n + size[i] = 1 + set[i] = Set() + end + + union!(set[1], 1:n) + + i = n + j = 1 + + while i >= 1 + v = pop!(set[j]) + α[v] = i + α⁻¹[i] = v + size[v] = 0 + + for w in neighbors(graph, v) + if size[w] >= 1 + delete!(set[size[w]], w) + size[w] += 1 + push!(set[size[w]], w) + end + end + + i -= 1 + j += 1 + + while j >= 1 && isempty(set[j]) + j -= 1 + end + end + + α⁻¹, α +end + const DEFAULT_ELIMINATION_ALGORITHM = AMDJL_AMD() diff --git a/test/JunctionTrees.jl b/test/JunctionTrees.jl index f3829a5..5216571 100644 --- a/test/JunctionTrees.jl +++ b/test/JunctionTrees.jl @@ -25,6 +25,9 @@ order = JunctionTrees.Order(graph, AMDJL_AMD()) order = JunctionTrees.Order(graph, MetisJL_ND()) @test length(order) == 17 +order = JunctionTrees.Order(graph, MCS()) +@test length(order) == 17 + order = JunctionTrees.Order(1:17) @test length(order) == 17 From 2da4596e40214c0594dc8160c9b9871150c8c519 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sun, 17 Nov 2024 18:07:11 -0500 Subject: [PATCH 15/48] Spacing --- src/junction_trees/elimination_algorithms.jl | 1 + 1 file changed, 1 insertion(+) diff --git a/src/junction_trees/elimination_algorithms.jl b/src/junction_trees/elimination_algorithms.jl index 79a2b77..0e1caeb 100644 --- a/src/junction_trees/elimination_algorithms.jl +++ b/src/junction_trees/elimination_algorithms.jl @@ -146,4 +146,5 @@ function mcs(graph::AbstractSymmetricGraph) α⁻¹, α end + const DEFAULT_ELIMINATION_ALGORITHM = AMDJL_AMD() From 7a7b30097dd7d820b099ff2df0c173bbdfa3644d Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Sun, 17 Nov 2024 18:10:09 -0500 Subject: [PATCH 16/48] Fixed bug --- src/Decompositions.jl | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index 1ec428a..8f2e4b3 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -224,8 +224,7 @@ end # Construct a tree decomposition of a graph. function Decompositions.StrDecomp(sgraph::AbstractSymmetricGraph, stree::SupernodeTree) - seperator = seperators(stree) - foreach(sort!, seperator) + seperator = map(sort ∘ collect, seperators(stree)) n = length(stree.tree) graph = Graph(n) From 156871bf8504f10eed2835cd62d5898352df8310 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 00:44:21 -0500 Subject: [PATCH 17/48] This may have done nothing at all. --- src/junction_trees/ordered_graphs.jl | 11 +++++++---- 1 file changed, 7 insertions(+), 4 deletions(-) diff --git a/src/junction_trees/ordered_graphs.jl b/src/junction_trees/ordered_graphs.jl index 0b28da5..5050bbf 100644 --- a/src/junction_trees/ordered_graphs.jl +++ b/src/junction_trees/ordered_graphs.jl @@ -58,20 +58,23 @@ end # Algorithm 4.2: Elimination Tree by Path Compression. function etree(graph::OrderedGraph) n = nv(graph) - parent = collect(1:n) - ancestor = collect(1:n) + parent = Vector{Int}(undef, n) + ancestor = Vector{Int}(undef, n) for i in 1:n + parent[i] = i + ancestor[i] = 0 + for k in inneighbors(graph, i) r = k - while ancestor[r] != r && ancestor[r] != i + while ancestor[r] != 0 && ancestor[r] != i t = ancestor[r] ancestor[r] = i r = t end - if ancestor[r] == r + if ancestor[r] == 0 ancestor[r] = i parent[r] = i end From ba2d6a23bd9563562c5889c0e930339911447f0d Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 00:51:51 -0500 Subject: [PATCH 18/48] Another meaningless change... --- src/junction_trees/ordered_graphs.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/junction_trees/ordered_graphs.jl b/src/junction_trees/ordered_graphs.jl index 5050bbf..bcbf02d 100644 --- a/src/junction_trees/ordered_graphs.jl +++ b/src/junction_trees/ordered_graphs.jl @@ -62,7 +62,7 @@ function etree(graph::OrderedGraph) ancestor = Vector{Int}(undef, n) for i in 1:n - parent[i] = i + parent[i] = 0 ancestor[i] = 0 for k in inneighbors(graph, i) @@ -81,6 +81,7 @@ function etree(graph::OrderedGraph) end end + parent[n] = n parent end From 6e1105b6b835ec85c8dc7a95e5e3c641fbbbffdc Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 01:49:44 -0500 Subject: [PATCH 19/48] Cleaned up the new code a bit. --- src/Decompositions.jl | 137 ++++++++++++++++++++++-------------------- 1 file changed, 73 insertions(+), 64 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index 8f2e4b3..980867a 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -193,27 +193,8 @@ end ################################## -const OTYPE = SymmetricGraph - - -const MTYPE = StructTightACSetTransformation{ - TypeLevelBasicSchema{ - Symbol, - Tuple{:V, :E}, - Tuple{(:src, :E, :V), (:tgt, :E, :V), (:inv, :E, :E)}, - Tuple{}, - Tuple{}, - Tuple{ - (nothing, :E, :E, ((:inv, :inv), ())), - (nothing, :E, :V, ((:inv, :src), (:tgt,))), - (nothing, :E, :V, ((:inv, :tgt), (:src,)))}}, - @NamedTuple{V::FinDomFunctionVector{Int, Vector{Int}, FinSetInt}, E::FinDomFunctionVector{Int, Vector{Int}, FinSetInt}}, - SymmetricGraph, - SymmetricGraph} - - # Construct a tree decomposition of a graph. -function Decompositions.StrDecomp( +function StrDecomp( graph::AbstractSymmetricGraph, ealg::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, stype::SupernodeType=DEFAULT_SUPERNODE_TYPE) @@ -223,73 +204,101 @@ end # Construct a tree decomposition of a graph. -function Decompositions.StrDecomp(sgraph::AbstractSymmetricGraph, stree::SupernodeTree) - seperator = map(sort ∘ collect, seperators(stree)) - +function StrDecomp(graph::AbstractSymmetricGraph, stree::SupernodeTree) n = length(stree.tree) - graph = Graph(n) - objects = Vector{OTYPE}(undef, 2n - 1) - morphisms = Vector{MTYPE}(undef, 2n - 2) - - for i in 1:n - snd = supernode(stree, i) - sep = seperator[i] - objects[i] = induced_subgraph(sgraph, permutation(stree.graph, [snd; sep])) - end + tree = Graph(n) for i in 1:n - 1 - sep = seperator[i] - objects[n + i] = induced_subgraph(sgraph, permutation(stree.graph, sep)) + add_edge!(tree, i, parentindex(stree.tree, i)) end - + + StrDecomp(tree, FinDomFunctor(homomorphisms(graph, stree)..., ∫(tree))) +end + + +function mappings(stree::SupernodeTree) + seperator = map(sort ∘ collect, seperators(stree)) + + n = length(stree.tree) + set = Vector{Vector{Int}}(undef, 2n - 1) + mapping = Vector{Vector{Int}}(undef, 2n - 2) + + for i in 1:n + set[i] = permutation(stree.graph, [supernode(stree, i); seperator[i]]) + end + + for i in 1:n - 1 + set[n + i] = permutation(stree.graph, seperator[i]) + end + for i in 1:n - 1 j = parentindex(stree.tree, i) - add_edge!(graph, i, j) - - sep_i = seperator[i] - sep_j = seperator[j] - snd_i = supernode(stree, i) - snd_j = supernode(stree, j) - - rep_j = first(snd_j) - len_j = length(snd_j) - - mapping = map(sep_i) do v - if v in snd_j - v - rep_j + 1 - else - len_j + searchsortedfirst(sep_j, v) - end + mapping[i] = induced_mapping(seperator[i], supernode(stree, j), seperator[j]) + end + + for i in 1:n - 1 + mapping[n + i - 1] = (1:length(seperator[i])) .+ length(supernode(stree, i)) + end + + set, mapping +end + + +function induced_mapping(sep::AbstractVector, psnd::AbstractUnitRange, psep::AbstractVector) + i = 0 + mapping = Vector{Int}(undef, length(sep)) + + for (j, v) in enumerate(sep) + if v in psnd + mapping[j] = v - first(psnd) + 1 + else + i = searchsortedfirst(view(psep, i + 1:length(psep)), v) + mapping[j] = length(psnd) + i end - - morphisms[i] = induced_homomorphism(mapping, objects[n + i], objects[j]) + end + + mapping +end + + +function homomorphisms(graph::AbstractSymmetricGraph, stree::SupernodeTree) + set, mapping = mappings(stree) + + n = length(stree.tree) + subgraph = Vector{Any}(undef, 2n - 1) + homomorphism = Vector{Any}(undef, 2n - 2) + + for i in 1:2n - 1 + subgraph[i] = induced_subgraph(graph, set[i]) end for i in 1:n - 1 - sep = seperator[i] - snd = supernode(stree, i) - - mapping = length(snd) + 1:length(snd) + length(sep) - morphisms[n + i - 1] = induced_homomorphism(mapping, objects[n + i], objects[i]) + homomorphism[i] = induced_homomorphism(subgraph[n + i], subgraph[parentindex(stree.tree, i)], mapping[i]) end - - StrDecomp(graph, FinDomFunctor(objects, morphisms, ∫(graph))) + + for i in 1:n - 1 + homomorphism[n + i - 1] = induced_homomorphism(subgraph[n + i], subgraph[i], mapping[n + i - 1]) + end + + subgraph, homomorphism end -function induced_homomorphism(vmapping, domain, codomain) - emapping = Vector{Int}(undef, ne(domain)) +function induced_homomorphism(domain::AbstractSymmetricGraph, codomain::AbstractSymmetricGraph, vmap::AbstractVector) index = Dict{Tuple{Int, Int}, Int}() + sizehint!(index, ne(codomain)) for e in edges(codomain) index[src(codomain, e), tgt(codomain, e)] = e end + emap = Vector{Int}(undef, ne(domain)) + for e in edges(domain) - emapping[e] = index[vmapping[src(domain, e)], vmapping[tgt(domain, e)]] + emap[e] = index[vmap[src(domain, e)], vmap[tgt(domain, e)]] end - ACSetTransformation(domain, codomain, V=vmapping, E=emapping) + ACSetTransformation(domain, codomain, V=vmap, E=emap) end From ac89572257a48b91fd2e8edc507ad47c3d38fc45 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 02:02:01 -0500 Subject: [PATCH 20/48] Removed unused import. --- src/Decompositions.jl | 2 -- 1 file changed, 2 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index 980867a..b17d6d4 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -18,8 +18,6 @@ using Catlab.Graphs using Catlab.ACSetInterface using Catlab.CategoricalAlgebra.Diagrams import Catlab.CategoricalAlgebra.Diagrams: ob_map, hom_map, colimit, limit -using Catlab.FinSets: FinSetInt, FinDomFunctionVector - ##################### From 4db55312fb41e656c9b87749e162e91229a98898 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 02:09:49 -0500 Subject: [PATCH 21/48] Avoided slow constructor. --- src/Decompositions.jl | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index b17d6d4..6443534 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -210,7 +210,8 @@ function StrDecomp(graph::AbstractSymmetricGraph, stree::SupernodeTree) add_edge!(tree, i, parentindex(stree.tree, i)) end - StrDecomp(tree, FinDomFunctor(homomorphisms(graph, stree)..., ∫(tree))) + diagram = FinDomFunctor(homomorphisms(graph, stree)..., ∫(tree)) + StrDecomp(tree, diagram, Decomposition, dom(diagram)) end From 26e4805e9bb6aab80538517260553c662d4e03ba Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 02:19:32 -0500 Subject: [PATCH 22/48] Removed unused function. --- src/junction_trees/ordered_graphs.jl | 13 ------------- 1 file changed, 13 deletions(-) diff --git a/src/junction_trees/ordered_graphs.jl b/src/junction_trees/ordered_graphs.jl index bcbf02d..f98385c 100644 --- a/src/junction_trees/ordered_graphs.jl +++ b/src/junction_trees/ordered_graphs.jl @@ -130,19 +130,6 @@ function BasicGraphs.edges(ograph::OrderedGraph) end -function BasicGraphs.add_edge!(ograph::OrderedGraph, i, j) - u = min(i, j) - v = max(i, j) - - if !has_edge(ograph.graph, u, v) - add_edge!(ograph.graph, u, v) - true - else - false - end -end - - function BasicGraphs.vertices(ograph::OrderedGraph) vertices(ograph.graph) end From fa6016cab8957495c1f814c9e1a6a56f2983c5c3 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 02:21:58 -0500 Subject: [PATCH 23/48] Implemented `ne`. --- src/junction_trees/ordered_graphs.jl | 5 +++++ 1 file changed, 5 insertions(+) diff --git a/src/junction_trees/ordered_graphs.jl b/src/junction_trees/ordered_graphs.jl index f98385c..b041ac2 100644 --- a/src/junction_trees/ordered_graphs.jl +++ b/src/junction_trees/ordered_graphs.jl @@ -110,6 +110,11 @@ end ############################ +function BasicGraphs.ne(ograph::OrderedGraph) + ne(ograph.graph) +end + + function BasicGraphs.nv(ograph::OrderedGraph) nv(ograph.graph) end From 39e6862a25ade3588f7c770b931d5e96011c9549 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 02:46:42 -0500 Subject: [PATCH 24/48] Fixed bug. --- src/Decompositions.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index 6443534..792dbfb 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -251,7 +251,7 @@ function induced_mapping(sep::AbstractVector, psnd::AbstractUnitRange, psep::Abs if v in psnd mapping[j] = v - first(psnd) + 1 else - i = searchsortedfirst(view(psep, i + 1:length(psep)), v) + i += searchsortedfirst(view(psep, i + 1:length(psep)), v) mapping[j] = length(psnd) + i end end From ff2030f08ff5eb55321b3375d5ee021b81574712 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 03:12:40 -0500 Subject: [PATCH 25/48] moved things around. --- src/Decompositions.jl | 58 ++++++++++++++++++++++++++----------------- 1 file changed, 35 insertions(+), 23 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index 792dbfb..b7c115f 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -223,19 +223,23 @@ function mappings(stree::SupernodeTree) mapping = Vector{Vector{Int}}(undef, 2n - 2) for i in 1:n + # clique(i) set[i] = permutation(stree.graph, [supernode(stree, i); seperator[i]]) end for i in 1:n - 1 + # seperator(i) set[n + i] = permutation(stree.graph, seperator[i]) end for i in 1:n - 1 + # seperator(i) → clique(parent(i)) j = parentindex(stree.tree, i) mapping[i] = induced_mapping(seperator[i], supernode(stree, j), seperator[j]) end for i in 1:n - 1 + # seperator(i) → clique(i) mapping[n + i - 1] = (1:length(seperator[i])) .+ length(supernode(stree, i)) end @@ -243,23 +247,6 @@ function mappings(stree::SupernodeTree) end -function induced_mapping(sep::AbstractVector, psnd::AbstractUnitRange, psep::AbstractVector) - i = 0 - mapping = Vector{Int}(undef, length(sep)) - - for (j, v) in enumerate(sep) - if v in psnd - mapping[j] = v - first(psnd) + 1 - else - i += searchsortedfirst(view(psep, i + 1:length(psep)), v) - mapping[j] = length(psnd) + i - end - end - - mapping -end - - function homomorphisms(graph::AbstractSymmetricGraph, stree::SupernodeTree) set, mapping = mappings(stree) @@ -267,15 +254,23 @@ function homomorphisms(graph::AbstractSymmetricGraph, stree::SupernodeTree) subgraph = Vector{Any}(undef, 2n - 1) homomorphism = Vector{Any}(undef, 2n - 2) - for i in 1:2n - 1 + for i in 1:n + # clique(i) subgraph[i] = induced_subgraph(graph, set[i]) end - + for i in 1:n - 1 + # seperator(i) + subgraph[n + i] = induced_subgraph(graph, set[n + i]) + end + + for i in 1:n - 1 + # seperator(i) → clique(parent(i)) homomorphism[i] = induced_homomorphism(subgraph[n + i], subgraph[parentindex(stree.tree, i)], mapping[i]) end for i in 1:n - 1 + # seperator(i) → clique(i) homomorphism[n + i - 1] = induced_homomorphism(subgraph[n + i], subgraph[i], mapping[n + i - 1]) end @@ -283,7 +278,24 @@ function homomorphisms(graph::AbstractSymmetricGraph, stree::SupernodeTree) end -function induced_homomorphism(domain::AbstractSymmetricGraph, codomain::AbstractSymmetricGraph, vmap::AbstractVector) +function induced_mapping(domain::AbstractVector, codomain_left::AbstractUnitRange, codomain_right::AbstractVector) + i = 0 + mapping = Vector{Int}(undef, length(domain)) + + for (j, v) in enumerate(domain) + if v in codomain_left + mapping[j] = v - first(codomain_left) + 1 + else + i += searchsortedfirst(view(codomain_right, i + 1:length(codomain_right)), v) + mapping[j] = length(codomain_left) + i + end + end + + mapping +end + + +function induced_homomorphism(domain::AbstractSymmetricGraph, codomain::AbstractSymmetricGraph, V::AbstractVector) index = Dict{Tuple{Int, Int}, Int}() sizehint!(index, ne(codomain)) @@ -291,13 +303,13 @@ function induced_homomorphism(domain::AbstractSymmetricGraph, codomain::Abstract index[src(codomain, e), tgt(codomain, e)] = e end - emap = Vector{Int}(undef, ne(domain)) + E = Vector{Int}(undef, ne(domain)) for e in edges(domain) - emap[e] = index[vmap[src(domain, e)], vmap[tgt(domain, e)]] + E[e] = index[V[src(domain, e)], V[tgt(domain, e)]] end - ACSetTransformation(domain, codomain, V=vmap, E=emap) + ACSetTransformation(domain, codomain; V, E) end From 68e148fe415b27d0a08ff35cbb43e344eb6d3504 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 10:07:33 -0500 Subject: [PATCH 26/48] Commented everything. --- src/Decompositions.jl | 9 +++++++++ src/junction_trees/elimination_trees.jl | 11 +++++++++++ src/junction_trees/ordered_graphs.jl | 12 ++++++++++++ src/junction_trees/orders.jl | 3 +++ src/junction_trees/postorder_trees.jl | 7 +++++++ src/junction_trees/supernode_trees.jl | 9 +++++++++ src/junction_trees/trees.jl | 3 +++ 7 files changed, 54 insertions(+) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index b7c115f..80881dd 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -192,6 +192,11 @@ end # Construct a tree decomposition of a graph. +# ---------------------------------------- +# graph simple connected graph +# ealg elimination algorithm +# stype supernode type +# ---------------------------------------- function StrDecomp( graph::AbstractSymmetricGraph, ealg::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, @@ -202,6 +207,10 @@ end # Construct a tree decomposition of a graph. +# ---------------------------------------- +# graph simple connected graph +# stree supernodal elimination tree +# ---------------------------------------- function StrDecomp(graph::AbstractSymmetricGraph, stree::SupernodeTree) n = length(stree.tree) tree = Graph(n) diff --git a/src/junction_trees/elimination_trees.jl b/src/junction_trees/elimination_trees.jl index 55c250e..7a1afc5 100644 --- a/src/junction_trees/elimination_trees.jl +++ b/src/junction_trees/elimination_trees.jl @@ -7,6 +7,10 @@ end # Construct an elimination tree using an elimination algorithm. +# ---------------------------------------- +# graph simple connected graph +# ealg elimination algorithm +# ---------------------------------------- function EliminationTree( graph::AbstractSymmetricGraph, ealg::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM) @@ -15,12 +19,19 @@ end # Construct the elimination tree of an ordered graph. +# ---------------------------------------- +# graph ordered graph +# ---------------------------------------- function EliminationTree(graph::OrderedGraph) EliminationTree(Tree(etree(graph)), graph) end # Postorder an elimination tree. +# ---------------------------------------- +# graph ordered graph +# order postorder +# ---------------------------------------- function EliminationTree{PostorderTree}(etree::EliminationTree, order::Order) EliminationTree(PostorderTree(etree.tree, order), OrderedGraph(etree.graph, order)) end diff --git a/src/junction_trees/ordered_graphs.jl b/src/junction_trees/ordered_graphs.jl index b041ac2..a983767 100644 --- a/src/junction_trees/ordered_graphs.jl +++ b/src/junction_trees/ordered_graphs.jl @@ -8,6 +8,10 @@ end # Given a graph G, construct the ordered graph # (G, σ), # where the permutation σ is computed using an elimination algorithm. +# ---------------------------------------- +# sgraph simple connected graph +# ealg elimination algorithm +# ---------------------------------------- function OrderedGraph(sgraph::AbstractSymmetricGraph, ealg::EliminationAlgorithm=DEFAULT_ELIMINATION_ALGORITHM) OrderedGraph(sgraph, Order(sgraph, ealg)) end @@ -15,6 +19,10 @@ end # Given a graph G and permutation σ, construct the ordered graph # (G, σ). +# ---------------------------------------- +# sgraph simple connected graph +# order vertex order +# ---------------------------------------- function OrderedGraph(sgraph::AbstractSymmetricGraph, order::Order) n = nv(sgraph) graph = Graph(n) @@ -34,6 +42,10 @@ end # Given an ordered graph (G, σ) and permutation μ, construct the ordered graph # (G, σ ∘ μ). +# ---------------------------------------- +# ograph ordered graph +# order permutation +# ---------------------------------------- function OrderedGraph(ograph::OrderedGraph, order::Order) n = nv(ograph) graph = Graph(n) diff --git a/src/junction_trees/orders.jl b/src/junction_trees/orders.jl index 84d3f71..ab9caab 100644 --- a/src/junction_trees/orders.jl +++ b/src/junction_trees/orders.jl @@ -7,6 +7,9 @@ end # Construct an permutation σ from a vector # (σ(1), ..., σ(n)). +# ---------------------------------------- +# order permutation +# ---------------------------------------- function Order(order::AbstractVector) n = length(order) index = Vector{Int}(undef, n) diff --git a/src/junction_trees/postorder_trees.jl b/src/junction_trees/postorder_trees.jl index 319454b..23ed2d8 100644 --- a/src/junction_trees/postorder_trees.jl +++ b/src/junction_trees/postorder_trees.jl @@ -8,6 +8,9 @@ end # Construct a tree from a postordered list of parents. +# ---------------------------------------- +# parent list of parents +# ---------------------------------------- function PostorderTree(parent::AbstractVector) n = length(parent) children = Vector{Vector{Int}}(undef, n) @@ -36,6 +39,10 @@ end # Postorder a tree. +# ---------------------------------------- +# tree tree +# order postorder +# ---------------------------------------- function PostorderTree(tree::Tree, order::Order) n = length(tree) parent = collect(1:n) diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl index 8fb3c98..830ef8c 100644 --- a/src/junction_trees/supernode_trees.jl +++ b/src/junction_trees/supernode_trees.jl @@ -10,6 +10,11 @@ end # Construct a supernodal elimination tree using an elimination algorithm. +# ---------------------------------------- +# graph simple connected graph +# ealg elimination algorithm +# stype supernode type +# ---------------------------------------- function SupernodeTree( graph::AbstractSymmetricGraph, ealg::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, @@ -20,6 +25,10 @@ end # Construct a supernodal elimination tree. +# ---------------------------------------- +# etree elimination tree +# stype supernode type +# ---------------------------------------- function SupernodeTree(etree::EliminationTree, stype::SupernodeType=DEFAULT_SUPERNODE_TYPE) degree = outdegrees(etree) snode, parent, ancestor = stree(etree, degree, stype) diff --git a/src/junction_trees/trees.jl b/src/junction_trees/trees.jl index fcc2069..759d19f 100644 --- a/src/junction_trees/trees.jl +++ b/src/junction_trees/trees.jl @@ -7,6 +7,9 @@ end # Construct a tree from a list of parents. +# ---------------------------------------- +# parent list of parents +# ---------------------------------------- function Tree(parent::AbstractVector) n = root = length(parent) children = Vector{Vector{Int}}(undef, n) From 7f13a30258891f9d1005a485c9acd11f30c50961 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 11:34:23 -0500 Subject: [PATCH 27/48] Spacing. --- src/junction_trees/elimination_algorithms.jl | 5 +---- 1 file changed, 1 insertion(+), 4 deletions(-) diff --git a/src/junction_trees/elimination_algorithms.jl b/src/junction_trees/elimination_algorithms.jl index 0e1caeb..816d5ca 100644 --- a/src/junction_trees/elimination_algorithms.jl +++ b/src/junction_trees/elimination_algorithms.jl @@ -98,11 +98,8 @@ function adjacencymatrix(graph::AbstractSymmetricGraph) end -# Simple Linear-Time Algorithms to Test Chordality of Graphs, Test Acyclicity -# of Hypergraphs, and Selectively Reduce Acyclic Hypergraphs -# +# 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) From 9e12333ac9b57baab2dbe77b2f5d335d24a35fd4 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 13:01:07 -0500 Subject: [PATCH 28/48] testing integration --- test/Decompositions.jl | 62 ++++++++++++++++++++++++++++++++++++++++-- 1 file changed, 60 insertions(+), 2 deletions(-) diff --git a/test/Decompositions.jl b/test/Decompositions.jl index 4761456..c9cb286 100644 --- a/test/Decompositions.jl +++ b/test/Decompositions.jl @@ -5,13 +5,14 @@ using PartialFunctions using StructuredDecompositions.Decompositions using StructuredDecompositions.FunctorUtils +using StructuredDecompositions.JunctionTrees: Order, Maximal using Catlab.Graphics using Catlab.Graphs using Catlab.ACSetInterface using Catlab.CategoricalAlgebra -#using Catlab.Graphi +#using Catlab.Graph #Define the instance####################### #bag 1 @@ -99,4 +100,61 @@ bigdecomp_skeleton = 𝐃ₛ(bigdecomp_to_sets) adhesionSpans(bigdecomp_skeleton) ) -end \ No newline at end of file + +################################## +# Integration with JunctionTrees # +################################## + + +graph = SymmetricGraph(17) + +add_edges!(graph, + [1, 1, 1, 1, 2, 2, 5, 5, 6, 6, 7, 7, 7, 10, 10, 10, 10, 12, 12, 12, 12, 15], + [3, 4, 5, 15, 3, 4, 9, 16, 9, 16, 8, 9, 15, 11, 13, 14, 17, 13, 14, 16, 17, 17]) + +decomposition = StrDecomp(graph, Order(1:17), Maximal()) + +@test decomposition.decomp_shape == @acset Graph begin + V = 8 + E = 7 + src = [1, 2, 3, 4, 5, 6, 7] + tgt = [5, 5, 4, 5, 6, 8, 8] +end + +@test map(i -> ob_map(decomposition.diagram, i), 1:15) == [ + induced_subgraph(graph, [7, 8, 9, 15]), # g h i o + induced_subgraph(graph, [6, 9, 16]), # f i p + induced_subgraph(graph, [2, 3, 4]), # b c d + induced_subgraph(graph, [1, 3, 4, 5, 15]), # a c d e o + induced_subgraph(graph, [5, 9, 15, 16]), # e i o p + induced_subgraph(graph, [15, 16, 17]), # o p q + induced_subgraph(graph, [10, 11, 13, 14, 17]), # j k m n q + induced_subgraph(graph, [12, 13, 14, 16, 17]), # l m n p q + induced_subgraph(graph, [9, 15]), # i o + induced_subgraph(graph, [9, 16]), # i p + induced_subgraph(graph, [3, 4]), # c d + induced_subgraph(graph, [5, 15]), # e o + induced_subgraph(graph, [15, 16]), # o p + induced_subgraph(graph, [16, 17]), # p q + induced_subgraph(graph, [13, 14, 17]), # m n q +] + +@test map(i -> hom_map(decomposition.diagram, i), 1:14) == [ + ACSetTransformation(induced_subgraph(graph, [9, 15]), induced_subgraph(graph, [5, 9, 15, 16]), V=[2, 3], E=Int[]), # i o → e i o p + ACSetTransformation(induced_subgraph(graph, [9, 16]), induced_subgraph(graph, [5, 9, 15, 16]), V=[2, 4], E=Int[]), # i p → e i o p + ACSetTransformation(induced_subgraph(graph, [3, 4]), induced_subgraph(graph, [1, 3, 4, 5, 15]), V=[2, 3], E=Int[]), # c d → a c d e o + ACSetTransformation(induced_subgraph(graph, [5, 15]), induced_subgraph(graph, [5, 9, 15, 16]), V=[1, 3], E=Int[]), # e o → e i o p + ACSetTransformation(induced_subgraph(graph, [15, 16]), induced_subgraph(graph, [15, 16, 17]), V=[1, 2], E=Int[]), # o p → o p q + ACSetTransformation(induced_subgraph(graph, [16, 17]), induced_subgraph(graph, [12, 13, 14, 16, 17]), V=[4, 5], E=Int[]), # p q → l m n p q + ACSetTransformation(induced_subgraph(graph, [13, 14, 17]), induced_subgraph(graph, [12, 13, 14, 16, 17]), V=[2, 3, 5], E=Int[]), # m n q → l m n p q + ACSetTransformation(induced_subgraph(graph, [9, 15]), induced_subgraph(graph, [7, 8, 9, 15]), V=[3, 4], E=Int[]), # i o → g h i o + ACSetTransformation(induced_subgraph(graph, [9, 16]), induced_subgraph(graph, [6, 9, 16]), V=[2, 3], E=Int[]), # i p → f i p + ACSetTransformation(induced_subgraph(graph, [3, 4]), induced_subgraph(graph, [2, 3, 4]), V=[2, 3], E=Int[]), # c d → b c d + ACSetTransformation(induced_subgraph(graph, [5, 15]), induced_subgraph(graph, [1, 3, 4, 5, 15]), V=[4, 5], E=Int[]), # e o → a c d e o + ACSetTransformation(induced_subgraph(graph, [15, 16]), induced_subgraph(graph, [5, 9, 15, 16]), V=[3, 4], E=Int[]), # o p → e i o p + ACSetTransformation(induced_subgraph(graph, [16, 17]), induced_subgraph(graph, [15, 16, 17]), V=[2, 3], E=Int[]), # p q → o p q + ACSetTransformation(induced_subgraph(graph, [13, 14, 17]), induced_subgraph(graph, [10, 11, 13, 14, 17]), V=[3, 4, 5], E=Int[]), # m n q → j k m n q +] + + +end From 5a0bf5cb3e76d11d212862608e5c5ba1d3f38609 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 21:39:27 -0500 Subject: [PATCH 29/48] Unconnected graphs. --- src/Decompositions.jl | 69 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 63 insertions(+), 6 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index 80881dd..c485769 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -191,22 +191,22 @@ end ################################## -# Construct a tree decomposition of a graph. +# Construct a tree decomposition. # ---------------------------------------- -# graph simple connected graph -# ealg elimination algorithm -# stype supernode type +# graph simple graph +# ealg elimination algorithm +# stype supernode type # ---------------------------------------- function StrDecomp( graph::AbstractSymmetricGraph, ealg::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, stype::SupernodeType=DEFAULT_SUPERNODE_TYPE) - StrDecomp(graph, SupernodeTree(graph, ealg, stype)) + merge_decompositions(decompositions(graph, ealg, stype)) end -# Construct a tree decomposition of a graph. +# Construct a tree decomposition. # ---------------------------------------- # graph simple connected graph # stree supernodal elimination tree @@ -224,6 +224,58 @@ function StrDecomp(graph::AbstractSymmetricGraph, stree::SupernodeTree) end +function merge_decompositions(decomposition::AbstractVector) + tree = map(decomposition) do d + d.decomp_shape + end + + domain = map(decomposition) do d + Catlab.graph(dom(d.diagram)) + end + + subgraph = map(decomposition) do d + ob_map(d.diagram) + end + + homomorphism = map(decomposition) do d + hom_map(d.diagram) + end + + diagram = FinDomFunctor(vcat(subgraph...), vcat(homomorphism...), FinCat(apex(coproduct(domain)))) + StrDecomp(apex(coproduct(tree)), diagram, Decomposition, dom(diagram)) +end + + +function decompositions(graph::AbstractSymmetricGraph, ealg::EliminationAlgorithm, stype::SupernodeType) + component = connected_components(graph) + + n = length(component) + decomposition = Vector(undef, n) + + for i in 1:n # TODO: construct decompositions in parallel + subgraph = induced_subgraph(graph, component[i]) + decomposition[i] = StrDecomp(subgraph, SupernodeTree(subgraph, ealg, stype)) + end + + decomposition +end + + +function decompositions(graph::AbstractSymmetricGraph, order::Order, stype::SupernodeType) + component = connected_components(graph) + + n = length(component) + decomposition = Vector(undef, n) + + for i in 1:n # TODO: construct decompositions in parallal + subgraph = induced_subgraph(graph, component[i]) + decomposition[i] = StrDecomp(subgraph, SupernodeTree(subgraph, induced_order(order, component[i]), stype)) + end + + decomposition +end + + function mappings(stree::SupernodeTree) seperator = map(sort ∘ collect, seperators(stree)) @@ -287,6 +339,11 @@ function homomorphisms(graph::AbstractSymmetricGraph, stree::SupernodeTree) end +function induced_order(order::Order, elements::AbstractVector) + Order(sortperm(inverse(order, elements))) +end + + function induced_mapping(domain::AbstractVector, codomain_left::AbstractUnitRange, codomain_right::AbstractVector) i = 0 mapping = Vector{Int}(undef, length(domain)) From f6e0d81b3d94efac1c935ffc6ea7a80bfde62d56 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 21:51:41 -0500 Subject: [PATCH 30/48] multi-threading --- src/Decompositions.jl | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index c485769..a29785a 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -252,7 +252,7 @@ function decompositions(graph::AbstractSymmetricGraph, ealg::EliminationAlgorith n = length(component) decomposition = Vector(undef, n) - for i in 1:n # TODO: construct decompositions in parallel + Threads.@threads for i in 1:n subgraph = induced_subgraph(graph, component[i]) decomposition[i] = StrDecomp(subgraph, SupernodeTree(subgraph, ealg, stype)) end @@ -267,7 +267,7 @@ function decompositions(graph::AbstractSymmetricGraph, order::Order, stype::Supe n = length(component) decomposition = Vector(undef, n) - for i in 1:n # TODO: construct decompositions in parallal + Threads.@threads for i in 1:n subgraph = induced_subgraph(graph, component[i]) decomposition[i] = StrDecomp(subgraph, SupernodeTree(subgraph, induced_order(order, component[i]), stype)) end From f520de56c2904ac5b54a78080116652c3666af5b Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Mon, 18 Nov 2024 21:57:29 -0500 Subject: [PATCH 31/48] moved things around --- src/Decompositions.jl | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index a29785a..dadc502 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -12,11 +12,13 @@ using PartialFunctions using MLStyle using AbstractTrees +using Base.Threads using Catlab using Catlab.CategoricalAlgebra using Catlab.Graphs using Catlab.ACSetInterface using Catlab.CategoricalAlgebra.Diagrams + import Catlab.CategoricalAlgebra.Diagrams: ob_map, hom_map, colimit, limit @@ -252,7 +254,7 @@ function decompositions(graph::AbstractSymmetricGraph, ealg::EliminationAlgorith n = length(component) decomposition = Vector(undef, n) - Threads.@threads for i in 1:n + @threads for i in 1:n subgraph = induced_subgraph(graph, component[i]) decomposition[i] = StrDecomp(subgraph, SupernodeTree(subgraph, ealg, stype)) end @@ -267,7 +269,7 @@ function decompositions(graph::AbstractSymmetricGraph, order::Order, stype::Supe n = length(component) decomposition = Vector(undef, n) - Threads.@threads for i in 1:n + @threads for i in 1:n subgraph = induced_subgraph(graph, component[i]) decomposition[i] = StrDecomp(subgraph, SupernodeTree(subgraph, induced_order(order, component[i]), stype)) end From c7d204c3f51038d9db7b6d00498e82df4b42cfee Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 11:11:26 -0500 Subject: [PATCH 32/48] printing --- src/junction_trees/supernode_trees.jl | 10 ++++++++++ 1 file changed, 10 insertions(+) diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl index 830ef8c..285b5f3 100644 --- a/src/junction_trees/supernode_trees.jl +++ b/src/junction_trees/supernode_trees.jl @@ -104,3 +104,13 @@ function seperators(stree::SupernodeTree) seperator[n] = Set() seperator end + + +function Base.show(io::IO, stree::SupernodeTree) + n = width(stree) + print(io, "width: $n\nsupernodal elimination tree:\n") + + print_tree(io, IndexNode(stree.tree)) do io, node + show(IOContext(io, :compact => true, :limit => true), permutation(stree.graph, supernode(stree, node.index))) + end +end From 562f2acee7e8da0b05b34c95b7b6ac07ce6b1719 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 12:46:52 -0500 Subject: [PATCH 33/48] Added another elimination algorithm: BT. --- Project.toml | 4 +++- src/JunctionTrees.jl | 4 +++- src/junction_trees/elimination_algorithms.jl | 22 ++++++++++++++++---- test/JunctionTrees.jl | 3 +++ 4 files changed, 27 insertions(+), 6 deletions(-) diff --git a/Project.toml b/Project.toml index b7cbce2..1797016 100644 --- a/Project.toml +++ b/Project.toml @@ -9,15 +9,17 @@ 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" MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078" Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" PartialFunctions = "570af359-4316-4cb7-8c74-252c00c2016b" SparseArrays = "2f01184e-e22b-5df5-ae63-d93ebab69eaf" +TreeWidthSolver = "7d267fc5-9ace-409f-a54c-cd2374872a55" [compat] -AbstractTrees = "0.4" AMD = "0.5" +AbstractTrees = "0.4" Catlab = "0.16" CuthillMcKee = "0.1" DataStructures = "0.18" diff --git a/src/JunctionTrees.jl b/src/JunctionTrees.jl index b200de6..c8a1a0f 100644 --- a/src/JunctionTrees.jl +++ b/src/JunctionTrees.jl @@ -3,7 +3,9 @@ module JunctionTrees import AMD import CuthillMcKee +import Graphs import Metis +import TreeWidthSolver using AbstractTrees using Catlab.BasicGraphs @@ -16,7 +18,7 @@ export Order, inverse # Elimination Algorithms -export AMDJL_AMD, CuthillMcKeeJL_RCM, MetisJL_ND, MCS +export AMDJL_AMD, CuthillMcKeeJL_RCM, MetisJL_ND, TreeWidthSolverJL_BT, MCS # Supernode Types diff --git a/src/junction_trees/elimination_algorithms.jl b/src/junction_trees/elimination_algorithms.jl index 816d5ca..f83c58d 100644 --- a/src/junction_trees/elimination_algorithms.jl +++ b/src/junction_trees/elimination_algorithms.jl @@ -5,6 +5,7 @@ A graph elimination algorithm. The options are - [`CuthillMcKeeJL_RCM`](@ref) - [`AMDJL_AMD`](@ref) - [`MetisJL_ND`](@ref) +- [`TreeWidthSolverJL_BT`](@ref) - [`MCS`](@ref) """ abstract type EliminationAlgorithm end @@ -34,6 +35,14 @@ The nested dissection heuristic. Uses Metis.jl. struct MetisJL_ND <: EliminationAlgorithm end +""" + TreeWidthSolverJL_BT <: EliminationAlgorithm + +The Bouchitte-Todinca algorithm. Uses TreeWidthSolver.jl. +""" +struct TreeWidthSolverJL_BT <: EliminationAlgorithm end + + """ MCS <: EliminationAlgorithm @@ -48,16 +57,14 @@ function Order(graph::AbstractSymmetricGraph) end -# Construct an order using the reverse Cuthill-McKee algorithm. Uses -# CuthillMcKee.jl. +# Construct an order using the reverse Cuthill-McKee algorithm. Uses CuthillMcKee.jl. function Order(graph::AbstractSymmetricGraph, ealg::CuthillMcKeeJL_RCM) order = CuthillMcKee.symrcm(adjacencymatrix(graph)) Order(order) end -# Construct an order using the approximate minimum degree algorithm. Uses -# AMD.jl. +# Construct an order using the approximate minimum degree algorithm. Uses AMD.jl. function Order(graph::AbstractSymmetricGraph, ealg::AMDJL_AMD) order = AMD.symamd(adjacencymatrix(graph)) Order(order) @@ -71,6 +78,13 @@ function Order(graph::AbstractSymmetricGraph, ealg::MetisJL_ND) end +# Construct an order using the Bouchitte-Todinca algorithm. Uses TreeWidthSolver.jl. +function Order(graph::AbstractSymmetricGraph, ealg::TreeWidthSolverJL_BT) + order = reverse(vcat(TreeWidthSolver.elimination_order(Graphs.Graph(adjacencymatrix(graph)))...)) + Order(order) +end + + # Concstruct an order using the maximum cardinality search algorithm. function Order(graph::AbstractSymmetricGraph, ealg::MCS) order, index = mcs(graph) diff --git a/test/JunctionTrees.jl b/test/JunctionTrees.jl index 5216571..b5fcf7d 100644 --- a/test/JunctionTrees.jl +++ b/test/JunctionTrees.jl @@ -25,6 +25,9 @@ order = JunctionTrees.Order(graph, AMDJL_AMD()) order = JunctionTrees.Order(graph, MetisJL_ND()) @test length(order) == 17 +order = JunctionTrees.Order(graph, TreeWidthSolverJL_BT()) +@test length(order) == 17 + order = JunctionTrees.Order(graph, MCS()) @test length(order) == 17 From 799365e676113ce4f5e388e4a848cc2c63cecd70 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 17:16:19 -0500 Subject: [PATCH 34/48] removed Graphs.jl dependency --- Project.toml | 1 - src/JunctionTrees.jl | 1 - src/junction_trees/elimination_algorithms.jl | 14 +++++++++++++- 3 files changed, 13 insertions(+), 3 deletions(-) diff --git a/Project.toml b/Project.toml index 1797016..08c9e46 100644 --- a/Project.toml +++ b/Project.toml @@ -9,7 +9,6 @@ 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" MLStyle = "d8e11817-5142-5d16-987a-aa16d5891078" Metis = "2679e427-3c69-5b7f-982b-ece356f1e94b" diff --git a/src/JunctionTrees.jl b/src/JunctionTrees.jl index c8a1a0f..c318e21 100644 --- a/src/JunctionTrees.jl +++ b/src/JunctionTrees.jl @@ -3,7 +3,6 @@ module JunctionTrees import AMD import CuthillMcKee -import Graphs import Metis import TreeWidthSolver diff --git a/src/junction_trees/elimination_algorithms.jl b/src/junction_trees/elimination_algorithms.jl index f83c58d..03d9ac4 100644 --- a/src/junction_trees/elimination_algorithms.jl +++ b/src/junction_trees/elimination_algorithms.jl @@ -80,7 +80,19 @@ end # Construct an order using the Bouchitte-Todinca algorithm. Uses TreeWidthSolver.jl. function Order(graph::AbstractSymmetricGraph, ealg::TreeWidthSolverJL_BT) - order = reverse(vcat(TreeWidthSolver.elimination_order(Graphs.Graph(adjacencymatrix(graph)))...)) + n = nv(graph) + 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))) + 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(order) end From 58153b2d59198d744f89a810c801bd61a5de6a3d Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 17:49:48 -0500 Subject: [PATCH 35/48] updated compat --- Project.toml | 1 + 1 file changed, 1 insertion(+) diff --git a/Project.toml b/Project.toml index 08c9e46..5dbb8a7 100644 --- a/Project.toml +++ b/Project.toml @@ -25,4 +25,5 @@ DataStructures = "0.18" MLStyle = "0.4" Metis = "1" PartialFunctions = "1.1" +TreeWidthSolver = "0.3" julia = "1.7" From e3c06db126c1fb345d4fcd0b9294e4d95916cccd Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 18:23:11 -0500 Subject: [PATCH 36/48] sped up mcs using linked lists --- Project.toml | 2 ++ src/JunctionTrees.jl | 1 + src/junction_trees/elimination_algorithms.jl | 21 ++++++++++---------- 3 files changed, 14 insertions(+), 10 deletions(-) diff --git a/Project.toml b/Project.toml index 5dbb8a7..4bdc28f 100644 --- a/Project.toml +++ b/Project.toml @@ -10,6 +10,7 @@ Catlab = "134e5e36-593f-5add-ad60-77f754baafbe" CuthillMcKee = "17f17636-5e38-52e3-a803-7ae3aaaf3da9" DataStructures = "864edb3b-99cc-5e75-8d2d-829cb0a9cfe8" 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" @@ -22,6 +23,7 @@ AbstractTrees = "0.4" Catlab = "0.16" CuthillMcKee = "0.1" DataStructures = "0.18" +LinkedLists = "0.1" MLStyle = "0.4" Metis = "1" PartialFunctions = "1.1" diff --git a/src/JunctionTrees.jl b/src/JunctionTrees.jl index c318e21..425b33f 100644 --- a/src/JunctionTrees.jl +++ b/src/JunctionTrees.jl @@ -3,6 +3,7 @@ module JunctionTrees import AMD import CuthillMcKee +import LinkedLists import Metis import TreeWidthSolver diff --git a/src/junction_trees/elimination_algorithms.jl b/src/junction_trees/elimination_algorithms.jl index 03d9ac4..c545791 100644 --- a/src/junction_trees/elimination_algorithms.jl +++ b/src/junction_trees/elimination_algorithms.jl @@ -130,31 +130,32 @@ end function mcs(graph::AbstractSymmetricGraph) n = nv(graph) α = Vector{Int}(undef, n) - α⁻¹ = Vector{Int}(undef, n) + β = Vector{Int}(undef, n) size = Vector{Int}(undef, n) - set = Vector{Set{Int}}(undef, n) + set = Vector{LinkedLists.LinkedList{Int}}(undef, n) + pointer = Vector{LinkedLists.ListNode{Int}}(undef, n) for i in 1:n size[i] = 1 - set[i] = Set() + set[i] = LinkedLists.LinkedList{Int}() + pointer[i] = push!(set[1], i) end - union!(set[1], 1:n) - i = n j = 1 while i >= 1 - v = pop!(set[j]) + v = first(set[j]) + deleteat!(set[j], pointer[v]) α[v] = i - α⁻¹[i] = v + β[i] = v size[v] = 0 for w in neighbors(graph, v) if size[w] >= 1 - delete!(set[size[w]], w) + deleteat!(set[size[w]], pointer[w]) size[w] += 1 - push!(set[size[w]], w) + pointer[w] = push!(set[size[w]], w) end end @@ -166,7 +167,7 @@ function mcs(graph::AbstractSymmetricGraph) end end - α⁻¹, α + β, α end From 1871eb6e5a7068517262ade2f914b2dd3517dbf2 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 18:37:45 -0500 Subject: [PATCH 37/48] moved things around --- src/junction_trees/supernode_trees.jl | 27 +++++++++++++-------------- 1 file changed, 13 insertions(+), 14 deletions(-) diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl index 285b5f3..11399fe 100644 --- a/src/junction_trees/supernode_trees.jl +++ b/src/junction_trees/supernode_trees.jl @@ -31,31 +31,30 @@ end # ---------------------------------------- function SupernodeTree(etree::EliminationTree, stype::SupernodeType=DEFAULT_SUPERNODE_TYPE) degree = outdegrees(etree) - snode, parent, ancestor = stree(etree, degree, stype) + supernode, parent, ancestor = stree(etree, degree, stype) tree = Tree(parent) - treeorder = postorder(tree) - graphorder = Order(vcat(snode[treeorder]...)) + order = postorder(tree) + tree = PostorderTree(tree, order) + permute!(supernode, order) + permute!(ancestor, order) - tree = PostorderTree(tree, treeorder) - graph = OrderedGraph(etree.graph, graphorder) - - permute!(snode, treeorder) - permute!(ancestor, treeorder) - permute!(degree, graphorder) + order = Order(vcat(supernode...)) + graph = OrderedGraph(etree.graph, order) + permute!(degree, order) n = length(tree) representative = zeros(Int, n) cardinality = zeros(Int, n) for i in 1:n - 1 - representative[i] = inverse(graphorder, snode[i][1]) - cardinality[i] = length(snode[i]) - ancestor[i] = inverse(graphorder, ancestor[i]) + representative[i] = inverse(order, supernode[i][1]) + cardinality[i] = length(supernode[i]) + ancestor[i] = inverse(order, ancestor[i]) end - representative[n] = inverse(graphorder, snode[n][1]) - cardinality[n] = length(snode[n]) + representative[n] = inverse(order, supernode[n][1]) + cardinality[n] = length(supernode[n]) SupernodeTree(tree, graph, representative, cardinality, ancestor, degree) end From 254a147914bf16116d71d89798e005dcd67172cf Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 19:12:53 -0500 Subject: [PATCH 38/48] moved things around --- src/junction_trees/supernode_trees.jl | 29 +++++++++------------------ src/junction_trees/supernode_types.jl | 11 +++++----- 2 files changed, 15 insertions(+), 25 deletions(-) diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl index 11399fe..0104dfb 100644 --- a/src/junction_trees/supernode_trees.jl +++ b/src/junction_trees/supernode_trees.jl @@ -43,18 +43,10 @@ function SupernodeTree(etree::EliminationTree, stype::SupernodeType=DEFAULT_SUPE graph = OrderedGraph(etree.graph, order) permute!(degree, order) - n = length(tree) - representative = zeros(Int, n) - cardinality = zeros(Int, n) - - for i in 1:n - 1 - representative[i] = inverse(order, supernode[i][1]) - cardinality[i] = length(supernode[i]) - ancestor[i] = inverse(order, ancestor[i]) - end - - representative[n] = inverse(order, supernode[n][1]) - cardinality[n] = length(supernode[n]) + representative = map(first, supernode) + cardinality = map(length, supernode) + map!(i -> inverse(order, i), ancestor, ancestor) + map!(i -> inverse(order, i), representative, representative) SupernodeTree(tree, graph, representative, cardinality, ancestor, degree) end @@ -89,15 +81,14 @@ function seperators(stree::SupernodeTree) end end - for i in 1:n - 2 - j = parentindex(stree.tree, i) - - for v in seperator[i] - if stree.ancestor[j] < v - push!(seperator[j], v) + for j in 1:n - 1 + for i in childindices(stree.tree, j) + for v in seperator[i] + if stree.ancestor[j] < v + push!(seperator[j], v) + end end end - end seperator[n] = Set() diff --git a/src/junction_trees/supernode_types.jl b/src/junction_trees/supernode_types.jl index 29158e6..0d58a81 100644 --- a/src/junction_trees/supernode_types.jl +++ b/src/junction_trees/supernode_types.jl @@ -37,6 +37,7 @@ struct Fundamental <: SupernodeType end # Vanderberghe and Andersen # Algorithm 4.1: Maximal supernodes and supernodal elimination tree. function stree(etree::EliminationTree, degree::AbstractVector, stype::SupernodeType) + m = 0 n = length(etree.tree) index = zeros(Int, n) snd = Vector{Int}[] @@ -47,14 +48,12 @@ function stree(etree::EliminationTree, degree::AbstractVector, stype::SupernodeT ww = findchild(etree, degree, stype, v) if isnothing(ww) - i = length(snd) + 1 - index[v] = i + index[v] = i = m += 1 push!(snd, [v]) - push!(q, length(snd)) - push!(a, n + 1) + push!(q, m) + push!(a, n) else - i = index[ww] - index[v] = i + index[v] = i = index[ww] push!(snd[i], v) end From f7f9af52647668340a0e7a5aabee68a9f32525f4 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 20:15:16 -0500 Subject: [PATCH 39/48] deleted pointless variable --- src/junction_trees/supernode_types.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/junction_trees/supernode_types.jl b/src/junction_trees/supernode_types.jl index 0d58a81..db46e37 100644 --- a/src/junction_trees/supernode_types.jl +++ b/src/junction_trees/supernode_types.jl @@ -48,19 +48,19 @@ function stree(etree::EliminationTree, degree::AbstractVector, stype::SupernodeT ww = findchild(etree, degree, stype, v) if isnothing(ww) - index[v] = i = m += 1 + index[v] = m += 1 push!(snd, [v]) push!(q, m) push!(a, n) else - index[v] = i = index[ww] - push!(snd[i], v) + index[v] = index[ww] + push!(snd[index[v]], v) end for w in childindices(etree.tree, v) if w !== ww j = index[w] - q[j] = i + q[j] = index[v] a[j] = v end end From c206be31fd6e361195704f8253b8cdb6c792b5af Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 20:39:27 -0500 Subject: [PATCH 40/48] renamed some variables --- src/junction_trees/supernode_types.jl | 46 +++++++++++++-------------- 1 file changed, 23 insertions(+), 23 deletions(-) diff --git a/src/junction_trees/supernode_types.jl b/src/junction_trees/supernode_types.jl index db46e37..97da65d 100644 --- a/src/junction_trees/supernode_types.jl +++ b/src/junction_trees/supernode_types.jl @@ -33,40 +33,40 @@ A fundamental supernode. struct Fundamental <: SupernodeType end -# Chordal Graphs and Semidefinite Optimization -# Vanderberghe and Andersen -# Algorithm 4.1: Maximal supernodes and supernodal elimination tree. +# Compact Clique Tree Data Structures in Sparse Matrix Factorizations +# Pothen and Sun +# Figure 4: The Clique Tree Algorithm 2 function stree(etree::EliminationTree, degree::AbstractVector, stype::SupernodeType) - m = 0 n = length(etree.tree) - index = zeros(Int, n) - snd = Vector{Int}[] - q = Int[] - a = Int[] + new_in_clique = Vector{Int}(undef, n) + new = Vector{Int}[] + parent = Int[] + first_anc = Int[] + + i = 0 for v in 1:n - ww = findchild(etree, degree, stype, v) - - if isnothing(ww) - index[v] = m += 1 - push!(snd, [v]) - push!(q, m) - push!(a, n) + u = findchild(etree, degree, stype, v) + + if !isnothing(u) + new_in_clique[v] = new_in_clique[u] + push!(new[new_in_clique[v]], v) else - index[v] = index[ww] - push!(snd[index[v]], v) + new_in_clique[v] = i += 1 + push!(new, [v]) + push!(parent, i) + push!(first_anc, n) end - for w in childindices(etree.tree, v) - if w !== ww - j = index[w] - q[j] = index[v] - a[j] = v + for s in childindices(etree.tree, v) + if s !== u + parent[new_in_clique[s]] = new_in_clique[v] + first_anc[new_in_clique[s]] = v end end end - snd, q, a + new, parent, first_anc end From e3f4dc854a230441edcf058127c895a4f1cebb15 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 20:49:38 -0500 Subject: [PATCH 41/48] renamed function --- src/junction_trees/supernode_types.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/junction_trees/supernode_types.jl b/src/junction_trees/supernode_types.jl index 97da65d..84ced33 100644 --- a/src/junction_trees/supernode_types.jl +++ b/src/junction_trees/supernode_types.jl @@ -46,7 +46,7 @@ function stree(etree::EliminationTree, degree::AbstractVector, stype::SupernodeT i = 0 for v in 1:n - u = findchild(etree, degree, stype, v) + u = child_in_supernode(etree, degree, stype, v) if !isnothing(u) new_in_clique[v] = new_in_clique[u] @@ -73,13 +73,13 @@ end # Find a child w of v such that # v ∈ snd(w). # If no such child exists, return nothing. -function findchild(etree::EliminationTree, degree::AbstractVector, stype::Node, v::Integer) end +function child_in_supernode(etree::EliminationTree, degree::AbstractVector, stype::Node, v::Integer) end # Find a child w of v such that # v ∈ snd(w). # If no such child exists, return nothing. -function findchild(etree::EliminationTree, degree::AbstractVector, stype::Maximal, v::Integer) +function child_in_supernode(etree::EliminationTree, degree::AbstractVector, stype::Maximal, v::Integer) for w in childindices(etree.tree, v) if degree[w] == degree[v] + 1 return w @@ -91,7 +91,7 @@ end # Find a child w of v such that # v ∈ snd(w). # If no such child exists, return nothing. -function findchild(etree::EliminationTree, degree::AbstractVector, stype::Fundamental, v::Integer) +function child_in_supernode(etree::EliminationTree, degree::AbstractVector, stype::Fundamental, v::Integer) ws = childindices(etree.tree, v) if length(ws) == 1 From 30eea14f0b5ee822a7a77672d65a8d28dcc1cb93 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 21:48:41 -0500 Subject: [PATCH 42/48] temp commit --- src/Decompositions.jl | 50 +++++++++++++++++++++++++++---------------- 1 file changed, 32 insertions(+), 18 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index dadc502..2d9fd62 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -227,24 +227,38 @@ end function merge_decompositions(decomposition::AbstractVector) - tree = map(decomposition) do d - d.decomp_shape - end - - domain = map(decomposition) do d - Catlab.graph(dom(d.diagram)) - end - - subgraph = map(decomposition) do d - ob_map(d.diagram) - end - - homomorphism = map(decomposition) do d - hom_map(d.diagram) - end - - diagram = FinDomFunctor(vcat(subgraph...), vcat(homomorphism...), FinCat(apex(coproduct(domain)))) - StrDecomp(apex(coproduct(tree)), diagram, Decomposition, dom(diagram)) + tree = apex(coproduct(map(d -> d.decomp_shape, decomposition))) + l = length(decomposition) + m = nv(tree) + subgraph = Vector(undef, 2m - l) + homomorphism = Vector(undef, 2m - 2l) + + i = 0 + + for j in 1:l + n = nv(decomposition[j].decomp_shape) + + for k in 1:n + subgraph[k + 2i - j + 1] = ob_map(decomposition[j].diagram, k) + end + + for k in 1:n - 1 + subgraph[k + 2i - j + 1 + m] = ob_map(decomposition[j].diagram, n + k) + end + + for k in 1:n - 1 + homomorphism[k + 2i - 2j + 2] = hom_map(decomposition[j].diagram, k) + end + + for k in 1:n - 1 + homomorphism[k + 2i - 2j + 2 + m - l] = hom_map(decomposition[j].diagram, n + k - 1) + end + + i += n + end + + diagram = FinDomFunctor(subgraph, homomorphism, ∫(tree)) + StrDecomp(tree, diagram, Decomposition, dom(diagram)) end From 70ef424d35b9455754fdced03e8b2dee994ef1b0 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 22:36:32 -0500 Subject: [PATCH 43/48] fixed bug in `merge_decompositions` --- src/Decompositions.jl | 8 ++++---- 1 file changed, 4 insertions(+), 4 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index 2d9fd62..60c111d 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -239,19 +239,19 @@ function merge_decompositions(decomposition::AbstractVector) n = nv(decomposition[j].decomp_shape) for k in 1:n - subgraph[k + 2i - j + 1] = ob_map(decomposition[j].diagram, k) + subgraph[i + k] = ob_map(decomposition[j].diagram, k) end for k in 1:n - 1 - subgraph[k + 2i - j + 1 + m] = ob_map(decomposition[j].diagram, n + k) + subgraph[i - j + k + m + 1] = ob_map(decomposition[j].diagram, k + n) end for k in 1:n - 1 - homomorphism[k + 2i - 2j + 2] = hom_map(decomposition[j].diagram, k) + homomorphism[i - j + k + 1] = hom_map(decomposition[j].diagram, k) end for k in 1:n - 1 - homomorphism[k + 2i - 2j + 2 + m - l] = hom_map(decomposition[j].diagram, n + k - 1) + homomorphism[i - j + k - l + m + 1] = hom_map(decomposition[j].diagram, k + n - 1) end i += n From 247cfaa1554962c7039172cb441d39b0b62daf3a Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 23:55:35 -0500 Subject: [PATCH 44/48] some docstrings --- src/Decompositions.jl | 12 ++++++------ src/junction_trees/elimination_algorithms.jl | 6 +++++- src/junction_trees/orders.jl | 16 ++++++++++------ 3 files changed, 21 insertions(+), 13 deletions(-) diff --git a/src/Decompositions.jl b/src/Decompositions.jl index 60c111d..44af30f 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -193,12 +193,12 @@ end ################################## -# Construct a tree decomposition. -# ---------------------------------------- -# graph simple graph -# ealg elimination algorithm -# stype supernode type -# ---------------------------------------- +""" + StrDecomp(graph::AbstractSymmetricGraph[, ealg::Union{Order, EliminationAlgorithm}[, stype::SupernodeType]]) + +Construct a structured decomposition of a simple graph, optionally specifying an elimination algorithm and +supernode type. +""" function StrDecomp( graph::AbstractSymmetricGraph, ealg::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, diff --git a/src/junction_trees/elimination_algorithms.jl b/src/junction_trees/elimination_algorithms.jl index c545791..d130c65 100644 --- a/src/junction_trees/elimination_algorithms.jl +++ b/src/junction_trees/elimination_algorithms.jl @@ -51,7 +51,11 @@ The maximum cardinality search algorithm. struct MCS <: EliminationAlgorithm end -# Construct an order using the default graph elimination algorithm. +""" + Order(graph::AbstractSymmetricGraph[, ealg::EliminationAlgorithm]) + +Construct an elimination order for a simple graph, optionally specifying an elimination algorithm. +""" function Order(graph::AbstractSymmetricGraph) Order(graph, DEFAULT_ELIMINATION_ALGORITHM) end diff --git a/src/junction_trees/orders.jl b/src/junction_trees/orders.jl index ab9caab..cfb4d71 100644 --- a/src/junction_trees/orders.jl +++ b/src/junction_trees/orders.jl @@ -1,15 +1,19 @@ -# A permutation σ of the set {1, ..., n}. +""" + Order <: AbstractVector{Int} + +A permutation of the set ``\\{1, \\dot, n\\}.`` +""" struct Order <: AbstractVector{Int} order::Vector{Int} # permutation index::Vector{Int} # inverse permutation end -# Construct an permutation σ from a vector -# (σ(1), ..., σ(n)). -# ---------------------------------------- -# order permutation -# ---------------------------------------- +""" + Order(order::AbstractVector) + +Construct a permutation ``\\sigma`` from a sequence ``(\\sigma(1), \\dots, \\sigma(n)).`` +""" function Order(order::AbstractVector) n = length(order) index = Vector{Int}(undef, n) From 9147bdcc37a58b7807424c9af71325049b177353 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Tue, 19 Nov 2024 23:58:37 -0500 Subject: [PATCH 45/48] mispellings --- src/junction_trees/elimination_algorithms.jl | 2 +- src/junction_trees/orders.jl | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/src/junction_trees/elimination_algorithms.jl b/src/junction_trees/elimination_algorithms.jl index d130c65..71d5d5c 100644 --- a/src/junction_trees/elimination_algorithms.jl +++ b/src/junction_trees/elimination_algorithms.jl @@ -101,7 +101,7 @@ function Order(graph::AbstractSymmetricGraph, ealg::TreeWidthSolverJL_BT) end -# Concstruct an order using the maximum cardinality search algorithm. +# Construct an order using the maximum cardinality search algorithm. function Order(graph::AbstractSymmetricGraph, ealg::MCS) order, index = mcs(graph) Order(order, index) diff --git a/src/junction_trees/orders.jl b/src/junction_trees/orders.jl index cfb4d71..849f06e 100644 --- a/src/junction_trees/orders.jl +++ b/src/junction_trees/orders.jl @@ -1,7 +1,7 @@ """ Order <: AbstractVector{Int} -A permutation of the set ``\\{1, \\dot, n\\}.`` +A permutation of the set ``\\{1, \\dots, n\\}.`` """ struct Order <: AbstractVector{Int} order::Vector{Int} # permutation From 27d80cef50a5bc018f1cf5eade5cd944773c6399 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Wed, 20 Nov 2024 11:08:06 -0500 Subject: [PATCH 46/48] punctuation --- src/junction_trees/ordered_graphs.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/junction_trees/ordered_graphs.jl b/src/junction_trees/ordered_graphs.jl index a983767..f06fa96 100644 --- a/src/junction_trees/ordered_graphs.jl +++ b/src/junction_trees/ordered_graphs.jl @@ -111,7 +111,7 @@ function permutation(ograph::OrderedGraph, i) end -# Get the index σ⁻¹(v), +# Get the index σ⁻¹(v). function inverse(ograph::OrderedGraph, v) inverse(ograph.order, v) end From 6c4dccd0c0f0341d28ece381656dc628a2bf9ad6 Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Wed, 20 Nov 2024 12:48:56 -0500 Subject: [PATCH 47/48] temp commit --- src/Decompositions.jl | 78 +++---------- src/JunctionTrees.jl | 5 + src/junction_trees/junction_trees.jl | 157 ++++++++++++++++++++++++++ src/junction_trees/supernode_trees.jl | 18 +-- 4 files changed, 180 insertions(+), 78 deletions(-) create mode 100644 src/junction_trees/junction_trees.jl diff --git a/src/Decompositions.jl b/src/Decompositions.jl index 44af30f..1ceb835 100644 --- a/src/Decompositions.jl +++ b/src/Decompositions.jl @@ -211,17 +211,17 @@ end # Construct a tree decomposition. # ---------------------------------------- # graph simple connected graph -# stree supernodal elimination tree +# jtree junction tree # ---------------------------------------- -function StrDecomp(graph::AbstractSymmetricGraph, stree::SupernodeTree) - n = length(stree.tree) +function StrDecomp(graph::AbstractSymmetricGraph, jtree::JunctionTree) + n = length(jtree) tree = Graph(n) for i in 1:n - 1 - add_edge!(tree, i, parentindex(stree.tree, i)) + add_edge!(tree, i, parentindex(jtree, i)) end - diagram = FinDomFunctor(homomorphisms(graph, stree)..., ∫(tree)) + diagram = FinDomFunctor(homomorphisms(graph, jtree)..., ∫(tree)) StrDecomp(tree, diagram, Decomposition, dom(diagram)) end @@ -270,7 +270,7 @@ function decompositions(graph::AbstractSymmetricGraph, ealg::EliminationAlgorith @threads for i in 1:n subgraph = induced_subgraph(graph, component[i]) - decomposition[i] = StrDecomp(subgraph, SupernodeTree(subgraph, ealg, stype)) + decomposition[i] = StrDecomp(subgraph, JunctionTree(subgraph, ealg, stype)) end decomposition @@ -285,70 +285,37 @@ function decompositions(graph::AbstractSymmetricGraph, order::Order, stype::Supe @threads for i in 1:n subgraph = induced_subgraph(graph, component[i]) - decomposition[i] = StrDecomp(subgraph, SupernodeTree(subgraph, induced_order(order, component[i]), stype)) + decomposition[i] = StrDecomp(subgraph, JunctionTree(subgraph, induced_order(order, component[i]), stype)) end decomposition end -function mappings(stree::SupernodeTree) - seperator = map(sort ∘ collect, seperators(stree)) - - n = length(stree.tree) - set = Vector{Vector{Int}}(undef, 2n - 1) - mapping = Vector{Vector{Int}}(undef, 2n - 2) - - for i in 1:n - # clique(i) - set[i] = permutation(stree.graph, [supernode(stree, i); seperator[i]]) - end - - for i in 1:n - 1 - # seperator(i) - set[n + i] = permutation(stree.graph, seperator[i]) - end - - for i in 1:n - 1 - # seperator(i) → clique(parent(i)) - j = parentindex(stree.tree, i) - mapping[i] = induced_mapping(seperator[i], supernode(stree, j), seperator[j]) - end - - for i in 1:n - 1 - # seperator(i) → clique(i) - mapping[n + i - 1] = (1:length(seperator[i])) .+ length(supernode(stree, i)) - end - - set, mapping -end - - -function homomorphisms(graph::AbstractSymmetricGraph, stree::SupernodeTree) - set, mapping = mappings(stree) - - n = length(stree.tree) +function homomorphisms(graph::AbstractSymmetricGraph, jtree::JunctionTree) + n = length(jtree) subgraph = Vector{Any}(undef, 2n - 1) homomorphism = Vector{Any}(undef, 2n - 2) for i in 1:n # clique(i) - subgraph[i] = induced_subgraph(graph, set[i]) + subgraph[i] = induced_subgraph(graph, clique(jtree, i)) end for i in 1:n - 1 # seperator(i) - subgraph[n + i] = induced_subgraph(graph, set[n + i]) + subgraph[n + i] = induced_subgraph(graph, seperator(jtree, i)) end for i in 1:n - 1 # seperator(i) → clique(parent(i)) - homomorphism[i] = induced_homomorphism(subgraph[n + i], subgraph[parentindex(stree.tree, i)], mapping[i]) + j = parentindex(jtree, i) + homomorphism[i] = induced_homomorphism(subgraph[n + i], subgraph[j], seperator_to_parent(jtree, i)) end for i in 1:n - 1 # seperator(i) → clique(i) - homomorphism[n + i - 1] = induced_homomorphism(subgraph[n + i], subgraph[i], mapping[n + i - 1]) + homomorphism[n + i - 1] = induced_homomorphism(subgraph[n + i], subgraph[i], seperator_to_self(jtree, i)) end subgraph, homomorphism @@ -360,23 +327,6 @@ function induced_order(order::Order, elements::AbstractVector) end -function induced_mapping(domain::AbstractVector, codomain_left::AbstractUnitRange, codomain_right::AbstractVector) - i = 0 - mapping = Vector{Int}(undef, length(domain)) - - for (j, v) in enumerate(domain) - if v in codomain_left - mapping[j] = v - first(codomain_left) + 1 - else - i += searchsortedfirst(view(codomain_right, i + 1:length(codomain_right)), v) - mapping[j] = length(codomain_left) + i - end - end - - mapping -end - - function induced_homomorphism(domain::AbstractSymmetricGraph, codomain::AbstractSymmetricGraph, V::AbstractVector) index = Dict{Tuple{Int, Int}, Int}() sizehint!(index, ne(codomain)) diff --git a/src/JunctionTrees.jl b/src/JunctionTrees.jl index 425b33f..4c126ca 100644 --- a/src/JunctionTrees.jl +++ b/src/JunctionTrees.jl @@ -29,6 +29,10 @@ export Node, Maximal, Fundamental export SupernodeTree, width, height, seperators, supernode, permutation +# Junction Trees +export JunctionTree, width, height, seperator, residual, clique, seperator_to_parent, seperator_to_self + + include("junction_trees/orders.jl") include("junction_trees/elimination_algorithms.jl") include("junction_trees/ordered_graphs.jl") @@ -37,6 +41,7 @@ include("junction_trees/postorder_trees.jl") include("junction_trees/elimination_trees.jl") include("junction_trees/supernode_types.jl") include("junction_trees/supernode_trees.jl") +include("junction_trees/junction_trees.jl") end diff --git a/src/junction_trees/junction_trees.jl b/src/junction_trees/junction_trees.jl new file mode 100644 index 0000000..403053b --- /dev/null +++ b/src/junction_trees/junction_trees.jl @@ -0,0 +1,157 @@ +# A junction tree. +struct JunctionTree + stree::SupernodeTree # supernodal elimination tree + seperator::Vector{Vector{Int}} # seperator +end + + +""" + JunctionTree(graph::AbstractSymmetricGraph[, ealg::Union{Order, EliminationAlgorithm}[, stype::SupernodeType]]) + +Construct a tree decomposition of a connected simple graph, optionally specifying an elimination algorithm and +a supernode type. +""" +function JunctionTree( + graph::AbstractSymmetricGraph, + ealg::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, + stype::SupernodeType=DEFAULT_SUPERNODE_TYPE) + + JunctionTree(SupernodeTree(graph, ealg, stype)) +end + + +# Construct a junction tree. +function JunctionTree(stree::SupernodeTree) + JunctionTree(stree, map(sort ∘ collect, seperators(stree))) +end + + +""" + clique(jtree::JunctionTree, i::Integer) + +Get the clique at node i. +""" +function clique(jtree::JunctionTree, i::Integer) + [residual(jtree, i); seperator(jtree, i)] +end + + +""" + seperator(jtree::JunctionTree, i::Integer) + +Get the seperator at node i. +""" +function seperator(jtree::JunctionTree, i::Integer) + permutation(jtree.stree.graph, jtree.seperator[i]) +end + + +""" + residual(jtree::JunctionTree, i::Integer) + +Get the residual at node i. +""" +function residual(jtree::JunctionTree, i::Integer) + permutation(jtree.stree.graph, supernode(jtree.stree, i)) +end + + +# Construct the inclusion seperator(i) → clique(parent(i)). +function seperator_to_parent(jtree::JunctionTree, i::Integer) + j = parentindex(jtree, i) + sep = jtree.seperator[i] + sep_parent = jtree.seperator[j] + res_parent = supernode(jtree.stree, j) + + i = 0 + index = Vector{Int}(undef, length(sep)) + + for (j, v) in enumerate(sep) + if v in res_parent + index[j] = v - first(res_parent) + 1 + else + i += searchsortedfirst(view(sep_parent, i + 1:length(sep_parent)), v) + index[j] = length(res_parent) + i + end + end + + index +end + + +# Construct the inclusion seperator(i) → clique(i). +function seperator_to_self(jtree::JunctionTree, i::Integer) + sep = jtree.seperator[i] + res = supernode(jtree.stree, i) + length(res) + 1:length(res) + length(sep) +end + + +""" + length(jtree::JunctionTree) + +Get the number of nodes in a junction tree. +""" +function Base.length(jtree::JunctionTree) + length(jtree.stree.tree) +end + + +""" + height(jtree::JunctionTree) + +Compute the height of a junction tree. +""" +function height(jtree::JunctionTree) + height(jtree.stree.tree) +end + + +""" + width(jtree::JunctionTree) + +Compute the width of a junction tree. +""" +function width(jtree::JunctionTree) + width(jtree.stree) +end + + +function Base.show(io::IO, jtree::JunctionTree) + n = width(jtree) + print(io, "width: $n\njunction tree:\n") + + print_tree(io, IndexNode(jtree)) do io, node + show(IOContext(io, :compact => true, :limit => true), clique(jtree, node.index)) + end +end + + +########################## +# Indexed Tree Interface # +########################## + + +function AbstractTrees.rootindex(jtree::JunctionTree) + rootindex(jtree.stree.tree) +end + + +function AbstractTrees.parentindex(jtree::JunctionTree, i::Integer) + parentindex(jtree.stree.tree, i) +end + + +function AbstractTrees.childindices(jtree::JunctionTree, i::Integer) + childindices(jtree.stree.tree, i) +end + + +function AbstractTrees.NodeType(::Type{IndexNode{JunctionTree, Int}}) + HasNodeType() +end + + +function AbstractTrees.nodetype(::Type{IndexNode{JunctionTree, Int}}) + IndexNode{JunctionTree, Int} +end diff --git a/src/junction_trees/supernode_trees.jl b/src/junction_trees/supernode_trees.jl index 0104dfb..cf7ff0e 100644 --- a/src/junction_trees/supernode_trees.jl +++ b/src/junction_trees/supernode_trees.jl @@ -16,11 +16,11 @@ end # stype supernode type # ---------------------------------------- function SupernodeTree( - graph::AbstractSymmetricGraph, - ealg::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, - stype::SupernodeType=DEFAULT_SUPERNODE_TYPE) + graph::AbstractSymmetricGraph, + ealg::Union{Order, EliminationAlgorithm}=DEFAULT_ELIMINATION_ALGORITHM, + stype::SupernodeType=DEFAULT_SUPERNODE_TYPE) - SupernodeTree(EliminationTree(graph, ealg), stype) + SupernodeTree(EliminationTree(graph, ealg), stype) end @@ -94,13 +94,3 @@ function seperators(stree::SupernodeTree) seperator[n] = Set() seperator end - - -function Base.show(io::IO, stree::SupernodeTree) - n = width(stree) - print(io, "width: $n\nsupernodal elimination tree:\n") - - print_tree(io, IndexNode(stree.tree)) do io, node - show(IOContext(io, :compact => true, :limit => true), permutation(stree.graph, supernode(stree, node.index))) - end -end From 078ace278c14ed5a391767b7e0d2e719d38281fb Mon Sep 17 00:00:00 2001 From: Richard Samuelson Date: Wed, 20 Nov 2024 13:05:07 -0500 Subject: [PATCH 48/48] Made a new type: JunctionTree. --- src/JunctionTrees.jl | 4 - src/junction_trees/junction_trees.jl | 9 +- test/JunctionTrees.jl | 223 ++++++++++----------------- 3 files changed, 89 insertions(+), 147 deletions(-) diff --git a/src/JunctionTrees.jl b/src/JunctionTrees.jl index 4c126ca..5437baa 100644 --- a/src/JunctionTrees.jl +++ b/src/JunctionTrees.jl @@ -25,10 +25,6 @@ export AMDJL_AMD, CuthillMcKeeJL_RCM, MetisJL_ND, TreeWidthSolverJL_BT, MCS export Node, Maximal, Fundamental -# Supernode Trees -export SupernodeTree, width, height, seperators, supernode, permutation - - # Junction Trees export JunctionTree, width, height, seperator, residual, clique, seperator_to_parent, seperator_to_self diff --git a/src/junction_trees/junction_trees.jl b/src/junction_trees/junction_trees.jl index 403053b..39f525a 100644 --- a/src/junction_trees/junction_trees.jl +++ b/src/junction_trees/junction_trees.jl @@ -1,4 +1,8 @@ -# A junction tree. +""" + JunctionTree + +A junction tree. +""" struct JunctionTree stree::SupernodeTree # supernodal elimination tree seperator::Vector{Vector{Int}} # seperator @@ -21,6 +25,9 @@ end # Construct a junction tree. +# ---------------------------------------- +# stree supernodal elimination tree +# ---------------------------------------- function JunctionTree(stree::SupernodeTree) JunctionTree(stree, map(sort ∘ collect, seperators(stree))) end diff --git a/test/JunctionTrees.jl b/test/JunctionTrees.jl index b5fcf7d..89ecf85 100644 --- a/test/JunctionTrees.jl +++ b/test/JunctionTrees.jl @@ -35,32 +35,12 @@ order = JunctionTrees.Order(1:17) @test length(order) == 17 # Figure 4.3 -stree = SupernodeTree(graph, order, Node()) -@test width(stree) == 4 -@test length(stree.tree) == 17 -@test height(stree.tree) == 7 +jtree = JunctionTree(graph, order, Node()) +@test width(jtree) == 4 +@test height(jtree) == 7 +@test length(jtree) == 17 -@test map(i -> inverse(stree.graph, i), 1:17) == [ - 5, # a - 4, # b - 6, # c - 7, # d - 8, # e - 3, # f - 1, # g - 2, # h - 9, # i - 12, # j - 13, # k - 11, # l - 14, # m - 15, # n - 10, # o - 16, # p - 17, # q -] - -@test map(i -> parentindex(stree.tree, i), 1:17) == [ +@test map(i -> parentindex(jtree, i), 1:17) == [ 2, 9, 9, @@ -80,7 +60,7 @@ stree = SupernodeTree(graph, order, Node()) nothing, ] -@test map(i -> childindices(stree.tree, i), 1:17) == [ +@test map(i -> childindices(jtree, i), 1:17) == [ [], [1], [], @@ -100,41 +80,41 @@ stree = SupernodeTree(graph, order, Node()) [16], ] -@test map(i -> supernode(stree, i), 1:17) == [ - [1], # g - [2], # h - [3], # f - [4], # b - [5], # a - [6], # c - [7], # d - [8], # e +@test map(i -> residual(jtree, i), 1:17) == [ + [7], # g + [8], # h + [6], # f + [2], # b + [1], # a + [3], # c + [4], # d + [5], # e [9], # i - [10], # o - [11], # l - [12], # j - [13], # k - [14], # m - [15], # n + [15], # o + [12], # l + [10], # j + [11], # k + [13], # m + [14], # n [16], # p [17], # q ] -@test map(sort ∘ collect, seperators(stree)) == [ - [2, 9, 10], # h i o - [9, 10], # i o +@test map(i -> seperator(jtree, i), 1:17) == [ + [8, 9, 15], # h i o + [9, 15], # i o [9, 16], # i p - [6, 7], # c d - [6, 7, 8, 10], # c d e o - [7, 8, 10], # d e o - [8, 10], # e o - [9, 10, 16], # i o p - [10, 16], # o p + [3, 4], # c d + [3, 4, 5, 15], # c d e o + [4, 5, 15], # d e o + [5, 15], # e o + [9, 15, 16], # i o p + [15, 16], # o p [16, 17], # p q - [14, 15, 16, 17], # m n p q - [13, 14, 15, 17], # k m n q - [14, 15, 17], # m n q - [15, 16, 17], # n p q + [13, 14, 16, 17], # m n p q + [11, 13, 14, 17], # k m n q + [13, 14, 17], # m n q + [14, 16, 17], # n p q [16, 17], # p q [17], # q [], # @@ -142,32 +122,12 @@ stree = SupernodeTree(graph, order, Node()) # Figure 4.7 (left) -stree = SupernodeTree(graph, order, Maximal()) -@test width(stree) == 4 -@test length(stree.tree) == 8 -@test height(stree.tree) == 4 - -@test map(i -> inverse(stree.graph, i), 1:17) == [ - 5, # a - 4, # b - 6, # c - 7, # d - 8, # e - 3, # f - 1, # g - 2, # h - 9, # i - 11, # j - 12, # k - 13, # l - 14, # m - 15, # n - 10, # o - 16, # p - 17, # q -] +jtree = JunctionTree(graph, order, Maximal()) +@test width(jtree) == 4 +@test height(jtree) == 4 +@test length(jtree) == 8 -@test map(i -> parentindex(stree.tree, i), 1:8) == [ +@test map(i -> parentindex(jtree, i), 1:8) == [ 5, 5, 4, @@ -178,7 +138,7 @@ stree = SupernodeTree(graph, order, Maximal()) nothing ] -@test map(i -> childindices(stree.tree, i), 1:8) == [ +@test map(i -> childindices(jtree, i), 1:8) == [ [], [], [], @@ -189,56 +149,35 @@ stree = SupernodeTree(graph, order, Maximal()) [6, 7], ] -@test map(i -> supernode(stree, i), 1:8) == [ - [1, 2], # g h - [3], # f - [4], # b - [5, 6, 7], # a c d - [8, 9], # e i - [10], # o - [11, 12], # j k - [13, 14, 15, 16, 17], # l m n p q +@test map(i -> residual(jtree, i), 1:8) == [ + [7, 8], # g h + [6], # f + [2], # b + [1, 3, 4], # a c d + [5, 9], # e i + [15], # o + [10, 11], # j k + [12, 13, 14, 16, 17], # l m n p q ] -@test map(sort ∘ collect, seperators(stree)) == [ - [9, 10], # i o +@test map(i -> seperator(jtree, i), 1:8) == [ + [9, 15], # i o [9, 16], # i p - [6, 7], # c d - [8, 10], # e o - [10, 16], # o p + [3, 4], # c d + [5, 15], # e o + [15, 16], # o p [16, 17], # p q - [14, 15, 17], # m n q + [13, 14, 17], # m n q [], # ] # Figure 4.9 -stree = SupernodeTree(graph, order, Fundamental()) -@test width(stree) == 4 -@test length(stree.tree) == 12 -@test height(stree.tree) == 5 - - -@test map(i -> inverse(stree.graph, i), 1:17) == [ - 5, # a - 4, # b - 6, # c - 7, # d - 8, # e - 3, # f - 1, # g - 2, # h - 9, # i - 12, # j - 13, # k - 11, # l - 14, # m - 15, # n - 10, # o - 16, # p - 17, # q -] +jtree = JunctionTree(graph, order, Fundamental()) +@test width(jtree) == 4 +@test height(jtree) == 5 +@test length(jtree) == 12 -@test map(i -> parentindex(stree.tree, i), 1:12) == [ +@test map(i -> parentindex(jtree, i), 1:12) == [ 7, 7, 5, @@ -253,7 +192,7 @@ stree = SupernodeTree(graph, order, Fundamental()) nothing, ] -@test map(i -> childindices(stree.tree, i), 1:12) == [ +@test map(i -> childindices(jtree, i), 1:12) == [ [], [], [], @@ -268,32 +207,32 @@ stree = SupernodeTree(graph, order, Fundamental()) [8, 11], ] -@test map(i -> supernode(stree, i), 1:12) == [ - [1, 2], # g h - [3], # f - [4], # b - [5], # a - [6, 7], # c d - [8], # e +@test map(i -> residual(jtree, i), 1:12) == [ + [7, 8], # g h + [6], # f + [2], # b + [1], # a + [3, 4], # c d + [5], # e [9], # i - [10], # o - [11], # l - [12, 13], # j k - [14, 15], # m n + [15], # o + [12], # l + [10, 11], # j k + [13, 14], # m n [16, 17], # p q ] -@test map(sort ∘ collect, seperators(stree)) == [ - [9, 10], # i o +@test map(i -> seperator(jtree, i), 1:12) == [ + [9, 15], # i o [9, 16], # i p - [6, 7], # c d - [6, 7, 8, 10], # c d e o - [8, 10], # e o - [9, 10, 16], # i o p - [10, 16], # o p + [3, 4], # c d + [3, 4, 5, 15], # c d e o + [5, 15], # e o + [9, 15, 16], # i o p + [15, 16], # o p [16, 17], # p q - [14, 15, 16, 17], # m n p q - [14, 15, 17], # m n q + [13, 14, 16, 17], # m n p q + [13, 14, 17], # m n q [16, 17], # p q [], # ]