From e8f3a4948ceec561b13936ee164a4c89dfaf8a59 Mon Sep 17 00:00:00 2001 From: Etienne dg Date: Mon, 28 Aug 2023 14:50:26 +0200 Subject: [PATCH] Mincut (#105) * fixes mincut. fixes #64 * add check for non-negative edge weight * decrease recording of solutions * fix for directed graphs * disable for directed graphs * Revert "fix for directed graphs" This reverts commit 73104f86a82187d345390a55c36022ef01193ad9. * disable directed graphs * fix disable directed graphs * new heuristic; fixed many things * remove comment * fix test --- src/traversals/maxadjvisit.jl | 133 ++++++++++++++++++++++++--------- test/traversals/maxadjvisit.jl | 29 ++++++- 2 files changed, 123 insertions(+), 39 deletions(-) diff --git a/src/traversals/maxadjvisit.jl b/src/traversals/maxadjvisit.jl index 980f31734..0b3bbf5b4 100644 --- a/src/traversals/maxadjvisit.jl +++ b/src/traversals/maxadjvisit.jl @@ -12,60 +12,119 @@ Return a tuple `(parity, bestcut)`, where `parity` is a vector of integer values that determines the partition in `g` (1 or 2) and `bestcut` is the -weight of the cut that makes this partition. An optional `distmx` matrix may -be specified; if omitted, edge distances are assumed to be 1. +weight of the cut that makes this partition. An optional `distmx` matrix +of non-negative weights may be specified; if omitted, edge distances are +assumed to be 1. """ -function mincut(g::AbstractGraph, distmx::AbstractMatrix{T}=weights(g)) where {T<:Real} - U = eltype(g) - colormap = zeros(UInt8, nv(g)) ## 0 if unseen, 1 if processing and 2 if seen and closed - parities = falses(nv(g)) - bestweight = typemax(T) - cutweight = zero(T) - visited = zero(U) ## number of vertices visited - pq = PriorityQueue{U,T}(Base.Order.Reverse) +@traitfn function mincut(g::::(!IsDirected), distmx::AbstractMatrix{T}=weights(g)) where {T <: Real} - # Set number of visited neighbors for all vertices to 0 - for v in vertices(g) - pq[v] = zero(T) - end + nvg = nv(g) + U = eltype(g) # make sure we have at least two vertices, otherwise, there's nothing to cut, # in which case we'll return immediately. - (haskey(pq, one(U)) && nv(g) > one(U)) || return (Vector{Int8}([1]), cutweight) - - # Give the starting vertex high priority - pq[one(U)] = one(T) + (nvg > one(U)) || return (Vector{Int8}([1]), zero(T)) - while !isempty(pq) - u = dequeue!(pq) - colormap[u] = 1 + is_merged = falses(nvg) + merged_vertices = IntDisjointSets(U(nvg)) + graph_size = nvg + # We need to mutate the weight matrix, + # and we need it clean (0 for non edges) + w = zeros(T, nvg, nvg) + size(distmx) != (nvg, nvg) && throw(ArgumentError("Adjacency / distance matrix size should match the number of vertices")) + @inbounds for e in edges(g) + d = distmx[src(e), dst(e)] + (d < 0) && throw(DomainError(w, "weigths should be non-negative")) + w[src(e), dst(e)] = d + (d != distmx[dst(e), src(e)]) && throw(ArgumentError("Adjacency / distance matrix must be symmetric")) + w[dst(e), src(e)] = d + end + # we also need to mutate neighbors when merging vertices + fadjlist = [collect(outneighbors(g, v)) for v in vertices(g)] + parities = falses(nvg) + bestweight = typemax(T) + pq = PriorityQueue{U,T}(Base.Order.Reverse) + u = last_vertex = one(U) - for v in outneighbors(g, u) - # if the target of e is already marked then decrease cutweight - # otherwise, increase it - ew = distmx[u, v] - if colormap[v] != 0 - cutweight -= ew - else - cutweight += ew + is_processed = falses(nvg) + @inbounds while graph_size > 1 + cutweight = zero(T) + is_processed .= false + is_processed[u] = true + # initialize pq + for v in vertices(g) + is_merged[v] && continue + v == u && continue + pq[v] = zero(T) + end + for v in fadjlist[u] + (is_merged[v] || v == u ) && continue + pq[v] = w[u, v] + cutweight += w[u, v] + end + # Minimum cut phase + while true + last_vertex = u + u, adj_cost = first(pq) + dequeue!(pq) + isempty(pq) && break + for v in fadjlist[u] + (is_merged[v] || u == v) && continue + # if the target of e is already marked then decrease cutweight + # otherwise, increase it + ew = w[u, v] + if is_processed[v] + cutweight -= ew + else + cutweight += ew + pq[v] += ew + end end - if haskey(pq, v) - pq[v] += distmx[u, v] + is_processed[u] = true + # adj_cost is a lower bound on the cut separating the two last vertices + # encountered, so if adj_cost >= bestweight, we can already merge these + # vertices to save one phase. + if adj_cost >= bestweight + _merge_vertex!(merged_vertices, fadjlist, is_merged, w, u, last_vertex) + graph_size -= 1 end end - colormap[u] = 2 - visited += one(U) - if cutweight < bestweight && visited < nv(g) + # check if we improved the mincut + if cutweight < bestweight bestweight = cutweight - for u in vertices(g) - parities[u] = (colormap[u] == 2) + for v in vertices(g) + parities[v] = (find_root!(merged_vertices, v) == u) end end + + # merge u and last_vertex + root = _merge_vertex!(merged_vertices, fadjlist, is_merged, w, u, last_vertex) + graph_size -= 1 + u = root # we are sure this vertex was not merged, so the next phase start from it end return (convert(Vector{Int8}, parities) .+ one(Int8), bestweight) end +function _merge_vertex!(merged_vertices, fadjlist, is_merged, w, u, v) + root = union!(merged_vertices, u, v) + non_root = (root == u) ? v : u + is_merged[non_root] = true + # update weights + for v2 in fadjlist[non_root] + w[root, v2] += w[non_root, v2] + w[v2, root] = w[root, v2] + end + # update neighbors + fadjlist[root] = union(fadjlist[root], fadjlist[non_root]) + for v in fadjlist[non_root] + if root ∉ fadjlist[v] + push!(fadjlist[v], root) + end + end + return root +end + """ maximum_adjacency_visit(g[, distmx][, log][, io][, s]) maximum_adjacency_visit(g[, s]) @@ -106,7 +165,7 @@ function maximum_adjacency_visit( log && println(io, "discover vertex: $u") for v in outneighbors(g, u) log && println(io, " -- examine neighbor from $u to $v") - if has_key[v] + if has_key[v] && (u != v) ed = distmx[u, v] pq[v] += ed end diff --git a/test/traversals/maxadjvisit.jl b/test/traversals/maxadjvisit.jl index 286717ccd..e76359140 100644 --- a/test/traversals/maxadjvisit.jl +++ b/test/traversals/maxadjvisit.jl @@ -33,15 +33,17 @@ @test ne(g) == m parity, bestcut = @inferred(mincut(g, eweights)) + if parity[1] == 2 + parity .= 3 .- parity + end @test length(parity) == 8 - @test parity == [2, 2, 1, 1, 2, 2, 1, 1] + @test parity == [1, 1, 2, 2, 1, 1, 2, 2] @test bestcut == 4.0 parity, bestcut = @inferred(mincut(g)) @test length(parity) == 8 - @test parity == [2, 1, 1, 1, 1, 1, 1, 1] @test bestcut == 2.0 v = @inferred(maximum_adjacency_visit(g)) @@ -55,4 +57,27 @@ end @test maximum_adjacency_visit(gx, 1) == [1, 2, 5, 6, 3, 7, 4, 8] @test maximum_adjacency_visit(gx, 3) == [3, 2, 7, 4, 6, 8, 5, 1] + + # non regression test for #64 + g = clique_graph(4, 2) + w = zeros(Int, 8, 8) + for e in edges(g) + if src(e) in [1, 5] || dst(e) in [1, 5] + w[src(e), dst(e)] = 3 + else + w[src(e), dst(e)] = 2 + end + w[dst(e), src(e)] = w[src(e), dst(e)] + end + w[1, 5] = 6 + w[5, 1] = 6 + parity, bestcut = @inferred(mincut(g, w)) + if parity[1] == 2 + parity .= 3 .- parity + end + @test parity == [1, 1, 1, 1, 2, 2, 2, 2] + @test bestcut == 6 + + w[6,7] = -1 + @test_throws DomainError mincut(g, w) end