feat: add molzahn clique merging algorithm

src/OPFSDP.jl

@@ -65,6 +65,7 @@ export compute_powers, compute_power_from, compute_power_to, str_power_from, str
 export get_branches_in, get_branches_out
 export display_opf
 export chordal_extension
+export merge_cliques!, merge_molzahn!
 export solve!
 
 end

src/decompose/merge/merge.jl

@@ -3,3 +3,87 @@ mutable struct MolzahnMerge <: AbstractMerge
     L::Float64
     MolzahnMerge(L::Float64=0.1) = new(L)
 end
+
+# TODO: compute the merge cost only once and then update them in function the cliques merged
+"""
+    cost = compute_merge_cost(c1, c2)
+
+Return the cost of merging the clique `c1` and `c2`.
+"""
+function compute_merge_cost(c1::Vector{Int}, c2::Vector{Int})
+    di = length(c1)
+    dk = length(c2)
+    sik = length(intersect(c1, c2))
+    dik = di + dk - sik
+    delta_ik = dik * (2 * dik + 1) - di * (2 * di + 1) - dk * (2 * dk + 1) - sik * (2 * sik + 1)
+    return delta_ik
+end
+
+
+"""
+    costs = compute_merge_cost_all(maximal_cliques, clique_tree)
+
+Return the merge cost array of all the edges in `clique_tree`.
+"""
+function compute_merge_cost_all(maximal_cliques, clique_tree)
+    costs = []
+    for i in 1:clique_tree.n-1
+        for k in i+1:clique_tree.n
+            if clique_tree[i, k] != 0
+                push!(costs, [i, k, compute_merge_cost(maximal_cliques[i], maximal_cliques[k])])
+            end
+        end
+    end
+    return costs
+end
+
+
+"""
+    clique_tree = merge_clique!(i, k, maximal_cliques, clique_tree)
+
+Merge the cliques `maximal_cliques[i]` and `maximal_cliques[k]`.
+Modify `maximal_cliques` by deleting the cliques at index `i` and `k` and adding the new cliques formed by their merge.
+Return the modified clique_tree where the cliques at index `i` and `k` are merged.
+"""
+function merge_clique!(i, k, clique, maximal_cliques, clique_tree)
+	neighbors = union(_neighbors(clique_tree, i; exclude=[i, k]), neighbors(clique_tree, k; exclude=[i, k]))
+    neighbors = map(n -> n - (n > i) - (n > k), neighbors)
+
+    deleteat!(maximal_cliques, [i, k])
+    push!(maximal_cliques, clique)
+    # not a fan of using unicode char, but couldn't find the corresponding ascii char for ∉
+    # cannot modify a matrix in place in Julia ? So I have to return it instead
+    clique_tree = clique_tree[1:clique_tree.n .∉ [[i, k]], 1:clique_tree.n .∉ [[i, k]]]
+    clique_tree = [clique_tree; zeros(clique_tree.n)']
+    clique_tree = [clique_tree zeros(clique_tree.m)]
+    for n in neighbors
+		weight = length(intersect(maximal_cliques[n], maximal_cliques[end]))
+        clique_tree[n, end] = weight
+		clique_tree[end, n] = weight
+    end
+    return clique_tree
+end
+
+
+"""
+    merge_molzahn!(cadj, maximal_cliques, clique_tree, L=0.1)
+
+Apply the Molzahn et al. merging alogrithm to the chordal graph represented by `cadj`.
+`L` is a percentage of the number of cliques in `maximal_cliques`. It is used to stop the merging.
+"""
+function merge_molzahn!(cadj, maximal_cliques, clique_tree, L::Float64=0.1)
+	treshold = L * length(maximal_cliques)
+	if treshold < 2
+		treshold = 2
+	end
+    while length(maximal_cliques) > treshold
+        costs = compute_merge_cost_all(maximal_cliques, clique_tree)
+        i, k, cost = argmin(x -> x[3], costs)
+		clique = union(maximal_cliques[i], maximal_cliques[k])
+        clique_tree = merge_clique!(i, k, clique, maximal_cliques, clique_tree)
+    end
+	for c in maximal_cliques
+		_make_subgraph_complete!(cadj, c)
+	end
+end
+merge_cliques!(cadj, maximal_cliques, clique_tree, merge_alg::MolzahnMerge) = merge_molzahn!(cadj, maximal_cliques, clique_tree, merge_alg.L)

src/utils/graphs.jl

@@ -21,3 +21,191 @@ function edges(adj::SparseMatrixCSC)
 	end
 	return edges_list
 end
+
+
+"""
+    peo = mcs(A)
+Maximum cardinality search for graph adjacency matrix A.
+Returns a perfect elimination ordering for chordal graphs.
+"""
+function mcs(A)
+    n = size(A, 1)
+    w = zeros(Int, n)
+    peo = zeros(Int, n)
+    unnumbered = collect(1:n)
+
+    for i = n:-1:1
+        z = unnumbered[argmax(w[unnumbered])]
+        filter!(x -> x != z, unnumbered)
+        peo[i] = z
+
+        Nz = findall(x->x!=0, A[:, z])
+        for y in intersect(Nz, unnumbered)
+            w[y] += 1
+        end
+    end
+    return peo
+end
+
+
+"""
+    T = prim(A)
+Return minimum spanning tree adjacency matrix, given adjacency matrix.
+If minweight == false, return the *maximum* weight spanning tree.
+
+Convention: start with node 1.
+"""
+function prim(A)
+    n = size(A, 1)
+    candidate_edges = []
+    unvisited = collect(1:n)
+    next_node = 1 # convention
+    T = spzeros(Int, n, n)
+
+    while length(unvisited) > 1
+        current_node = next_node
+        filter!(node -> node != current_node, unvisited)
+        neighbors = intersect(findall(x->x!=0, A[:, current_node]), unvisited)
+        current_node_edges = [(current_node, i) for i in neighbors]
+        append!(candidate_edges, current_node_edges)
+        filter!(edge -> length(intersect(edge, unvisited)) == 1, candidate_edges)
+        weights = [A[edge...] for edge in candidate_edges]
+        next_edge = candidate_edges[argmax(weights)]
+        filter!(edge -> edge != next_edge, candidate_edges)
+        weight = maximum(weights)
+        T[next_edge[1], next_edge[2]] = weight
+        T[next_edge[2], next_edge[1]] = weight
+        next_node = intersect(next_edge, unvisited)[1]
+    end
+    return T
+end
+
+
+"""
+    A = overlap_graph(groups)
+Return adjacency matrix for overlap graph associated with `groups`.
+I.e. if `A[i, j] = k`, then `groups[i]` and `groups[j]` share `k` elements.
+"""
+function overlap_graph(groups)
+    n = length(groups)
+    I = Vector{Int}()
+    J = Vector{Int}()
+    V = Vector{Int}()
+    for (i, gi) in enumerate(groups)
+        for (j, gj) in enumerate(groups)
+            if gi != gj
+                overlap = length(intersect(gi, gj))
+                if overlap > 0
+                    push!(I, i)
+                    push!(J, j)
+                    push!(V, overlap)
+                end
+            end
+        end
+    end
+    return sparse(I, J, V, n, n)
+end
+
+
+"""
+    n = neighbors(adj, v; exclude)
+
+Return the neighbors of the vertex `v` in `adj` (undirected).
+"""
+function neighbors(adj::SparseMatrixCSC, v::Int; exclude=[v])
+	return [i for i in 1:adj.n if (adj[v, i] != 0 || adj[i, v] != 0) && i ∉ exclude]
+end
+
+
+"""
+    is_complete = check_neighboor_complete(adj, v, v_sub)
+
+Return `true` if the neighbors of `v` in the subgraph of `adj` formed by the vertices in `v_sub` is a complete graph.
+"""
+function check_neighboor_complete(adj::SparseMatrixCSC, v::Int, v_sub::Vector{Int}=1:adj.n)
+    neighbors = filter(n -> n in v_sub, neighbors(adj, v))
+    for i in 1:length(neighbors)-1
+        for j in i+1:length(neighbors)
+            if adj[neighbors[i], neighbors[j]] == 0
+                return false
+            end
+        end
+    end
+    return true
+end
+
+
+"""
+    is_chordal = check_chordal(adj)
+
+Return `true` if the adjacency matrix `adj` corresponds to a chordal graph.
+"""
+function check_chordal(adj::SparseMatrixCSC)
+    peo = mcs(adj)
+    for i in 1:length(peo)-1
+        if !_check_neighboor_complete(adj, peo[i], peo[i:end])
+            return false
+        end
+    end
+    return true
+end
+
+
+"""
+    make_neighborhood_complete!(adj, v)
+
+Make the neighborhood of `v` complete in `adj`.
+"""
+function make_neighborhood_complete!(adj::SparseMatrixCSC, v::Int)
+    neighbors = neighbors(adj, v)
+    for i in 1:length(neighbors) - 1
+        for j in i+1:length(neighbors)
+            adj[neighbors[i], neighbors[j]] = 1.0
+            adj[neighbors[j], neighbors[i]] = 1.0
+        end
+    end
+end
+
+
+"""
+    pmdo = perfect_minimum_degree_ordering(adj)
+
+Return a perfect minimum degree ordering of the graph represented by `adj`.
+"""
+function perfect_minimum_degree_ordering(adj::SparseMatrixCSC)
+    vertex_map = 1:adj.n
+    perm = Int[]
+    for i in 1:adj.n
+        vmindeg = argmin([sum(adj[1:end, j]) for j in 1:adj.n])
+        push!(perm, vertex_map[vmindeg])
+        make_neighborhood_complete!(adj, vmindeg)
+        vertex_map = [vertex_map[vmindeg] <= vertex_map[j] ? vertex_map[j + 1] : vertex_map[j] for j in 1:adj.n - 1]
+        adj = adj[1:end .!= vmindeg, 1:end .!= vmindeg]
+    end
+    return perm
+end
+
+
+"""
+    pmdo = perfect_minimum_degree_ordering(adj)
+
+Return a perfect minimum degree ordering of the graph build from the power model `pm`
+"""
+function perfect_minimum_degree_ordering(pm::AbstractSparseSDPWRMModel, nw::Int=nw_id_default)
+    adj, lookup_index = adjacency_matrix(pm, nw)
+    return perfect_minimum_degree_ordering(adj), lookup_index
+end
+
+"""
+    make_subgraph_complete!(adj, vertices)
+
+Make the subgraph of `adj` composed of the vertices in `vertices` complete.
+"""
+function make_subgraph_complete!(adj::SparseMatrixCSC, vertices::Vector{Int})
+    for i in 1:length(vertices)-1
+        for j in i+1:length(vertices)
+            adj[vertices[i], vertices[j]] = 1.0
+            adj[vertices[j], vertices[i]] = 1.0
+        end
+    end
+end