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 73104f8.

* disable directed graphs

* fix disable directed graphs

* new heuristic; fixed many things

* remove comment

* fix test

Author: etiennedeg

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

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