Use Hopcroft-Karp matching algorithm

Project.toml

@@ -4,7 +4,6 @@ authors = ["Robert Parker and contributors"]
 version = "0.1.0"
 
 [deps]
-BipartiteMatching = "79040ab4-24c8-4c92-950c-d48b5991a0f6"
 Graphs = "86223c79-3864-5bf0-83f7-82e725a168b6"
 JuMP = "4076af6c-e467-56ae-b986-b466b2749572"
 MathOptInterface = "b8f27783-ece8-5eb3-8dc8-9495eed66fee"

src/dulmage_mendelsohn.jl

@@ -129,15 +129,15 @@ function dulmage_mendelsohn(graph::Graphs.Graph, set1::Set)
     matched_with_reachable1 = [matching[n] for n in reachable1]
     matched_with_reachable2 = [matching[n] for n in reachable2]
 
-    filter = cat(
+    filter = Set(cat(
         unmatched1,
         unmatched2,
         reachable1,
         reachable2,
         matched_with_reachable1,
         matched_with_reachable2;
         dims=1,
-    )
+    ))
     other1 = [n for n in nodes1 if !(n in filter)]
     other2 = [n for n in nodes2 if !(n in filter)]
 

src/maximum_matching.jl

@@ -18,8 +18,9 @@
 #  ___________________________________________________________________________
 
 import Graphs
-import BipartiteMatching as BM
 
+const UNMATCHED = nothing
+MatchedNodeType{T} = Union{T,typeof(UNMATCHED)}
 
 """
 TODO: This should probably be promoted to Graphs.jl
@@ -44,20 +45,147 @@ function _is_valid_bipartition(graph::Graphs.Graph, set1::Set)
     return true
 end
 
+# The following three functions are copied from the branch in PR #291
+# of Graphs.jl, https://github.com/JuliaGraphs/Graphs.jl/pull/291.
+# They will be removed when/if this PR is merged in favor of using the
+# Graphs.jl maximum_matching function.
+
+"""
+Determine whether an augmenting path exists and mark distances
+so we can compute shortest-length augmenting paths in the DFS.
+"""
+function _hk_augmenting_bfs!(
+    graph::Graphs.AbstractGraph{T},
+    set1::Vector{T},
+    matching::Dict{T,MatchedNodeType{T}},
+    distance::Dict{MatchedNodeType{T},Float64},
+)::Bool where {T<:Integer}
+    # Initialize queue with the unmatched nodes in set1
+    queue = Vector{MatchedNodeType{eltype(graph)}}([
+        n for n in set1 if matching[n] == UNMATCHED
+    ])
+
+    distance[UNMATCHED] = Inf
+    for n in set1
+        if matching[n] == UNMATCHED
+            distance[n] = 0.0
+        else
+            distance[n] = Inf
+        end
+    end
+
+    while !isempty(queue)
+        n1 = popfirst!(queue)
+
+        # If n1 is (a) matched or (b) in set1
+        if distance[n1] < Inf && n1 != UNMATCHED
+            for n2 in Graphs.neighbors(graph, n1)
+                # If n2 has not been encountered
+                if distance[matching[n2]] == Inf
+                    # Give it a distance
+                    distance[matching[n2]] = distance[n1] + 1
+
+                    # Note that n2 could be unmatched
+                    push!(queue, matching[n2])
+                end
+            end
+        end
+    end
+
+    found_augmenting_path = (distance[UNMATCHED] < Inf)
+    # The distance to UNMATCHED is the length of the shortest augmenting path
+    return found_augmenting_path
+end
+
+"""
+Compute augmenting paths and update the matching
+"""
+function _hk_augmenting_dfs!(
+    graph::Graphs.AbstractGraph{T},
+    root::MatchedNodeType{T},
+    matching::Dict{T,MatchedNodeType{T}},
+    distance::Dict{MatchedNodeType{T},Float64},
+)::Bool where {T<:Integer}
+    if root != UNMATCHED
+        for n in Graphs.neighbors(graph, root)
+            # Traverse edges of the minimum-length alternating path
+            if distance[matching[n]] == distance[root] + 1
+                if _hk_augmenting_dfs!(graph, matching[n], matching, distance)
+                    # If the edge is part of an augmenting path, update the
+                    # matching
+                    matching[root] = n
+                    matching[n] = root
+                    return true
+                end
+            end
+        end
+        # If we could not find a matched edge that was part of an augmenting
+        # path, we need to make sure we don't consider this vertex again
+        distance[root] = Inf
+        return false
+    else
+        # Return true to indicate that we are part of an augmenting path
+        return true
+    end
+end
+
+"""
+    hopcroft_karp_matching(graph::AbstractGraph)::Dict
+
+Compute a maximum-cardinality matching of a bipartite graph via the
+[Hopcroft-Karp algorithm](https://en.wikipedia.org/wiki/Hopcroft-Karp_algorithm).
+
+The return type is a dict mapping nodes to nodes. All matched nodes are included
+as keys. For example, if `i` is matched with `j`, `i => j` and `j => i` are both
+included in the returned dict.
+
+### Performance
+
+The algorithms runs in O((m + n)n^0.5), where n is the number of vertices and
+m is the number of edges. As it does not assume the number of edges is O(n^2),
+this algorithm is particularly effective for sparse bipartite graphs.
+
+### Arguments
+
+* `graph`: The bipartite `Graph` for which a maximum matching is computed
+
+### Exceptions
+
+* `ArgumentError`: The provided graph is not bipartite
+
+"""
+function hopcroft_karp_matching(graph::Graphs.AbstractGraph{T})::Dict{T,T} where {T<:Integer}
+    bmap = Graphs.bipartite_map(graph)
+    if length(bmap) != Graphs.nv(graph)
+        throw(ArgumentError("Provided graph is not bipartite"))
+    end
+    set1 = [n for n in Graphs.vertices(graph) if bmap[n] == 1]
+
+    # Initialize "state" that is modified during the algorithm
+    matching = Dict{eltype(graph),MatchedNodeType{eltype(graph)}}(
+        n => UNMATCHED for n in Graphs.vertices(graph)
+    )
+    distance = Dict{MatchedNodeType{eltype(graph)},Float64}()
+
+    # BFS to determine whether any augmenting paths exist
+    while _hk_augmenting_bfs!(graph, set1, matching, distance)
+        for n1 in set1
+            if matching[n1] == UNMATCHED
+                # DFS to update the matching along a minimum-length
+                # augmenting path
+                _hk_augmenting_dfs!(graph, n1, matching, distance)
+            end
+        end
+    end
+    matching = Dict(i => j for (i, j) in matching if j != UNMATCHED)
+    return matching
+end
 
 function maximum_matching(graph::Graphs.Graph, set1::Set)
     if !_is_valid_bipartition(graph, set1)
         throw(Exception)
     end
-    n_nodes = Graphs.nv(graph)
-    card1 = length(set1)
-    nodes1 = sort([node for node in set1])
-    set2 = setdiff(Set(1:n_nodes), set1)
-    nodes2 = sort([node for node in set2])
-    edge_set = Set((n1, n2) for n1 in nodes1 for n2 in Graphs.neighbors(graph, n1))
-    amat = BitArray{2}((r, c) in edge_set for r in nodes1, c in nodes2)
-    matching, _ = BM.findmaxcardinalitybipartitematching(amat)
-    # Translate row/column coordinates back into nodes of the graph
-    graph_matching = Dict(nodes1[r] => nodes2[c] for (r, c) in matching)
-    return graph_matching
+    matching = hopcroft_karp_matching(graph)
+    matching = Dict(i => j for (i, j) in matching if i in set1)
+    return matching
 end