beam: add Hopcroft-Karp bipartite matching with brute-force-verified tests

wizard-work/beam/hopcroft_karp.py

@@ -0,0 +1,59 @@
+"""Hopcroft-Karp maximum bipartite matching (O(E sqrt(V))).
+
+Graph is given as adjacency from left vertices 0..nL-1 to right vertices
+0..nR-1. Returns (size, matchL, matchR) where matchL[u] = v or -1 and
+matchR[v] = u or -1.
+"""
+from collections import deque
+
+
+def hopcroft_karp(adj, n_left, n_right):
+    INF = 10 ** 9
+    pair_u = [-1] * n_left
+    pair_v = [-1] * n_right
+    dist = [0] * n_left
+
+    def bfs():
+        q = deque()
+        for u in range(n_left):
+            if pair_u[u] == -1:
+                dist[u] = 0
+                q.append(u)
+            else:
+                dist[u] = INF
+        found = False
+        while q:
+            u = q.popleft()
+            if dist[u] < INF:
+                for v in adj[u]:
+                    pu = pair_v[v]
+                    if pu != -1 and dist[pu] == INF:
+                        dist[pu] = dist[u] + 1
+                        q.append(pu)
+                    elif pu == -1:
+                        found = True
+        return found
+
+    def dfs(u):
+        if u != -1:
+            for v in adj[u]:
+                pu = pair_v[v]
+                if pu == -1 or (dist[pu] == dist[u] + 1 and dfs(pu)):
+                    pair_u[u] = v
+                    pair_v[v] = u
+                    return True
+            dist[u] = INF
+            return False
+        return True
+
+    matching = 0
+    while bfs():
+        for u in range(n_left):
+            if pair_u[u] == -1:
+                if dfs(u):
+                    matching += 1
+    return matching, pair_u, pair_v
+
+
+def max_bipartite_matching(adj, n_left, n_right):
+    return hopcroft_karp(adj, n_left, n_right)

wizard-work/beam/test_hopcroft_karp.py

@@ -0,0 +1,75 @@
+import random
+from hopcroft_karp import max_bipartite_matching
+
+
+def kuhn(adj, n_left, n_right):
+    match_r = [-1] * n_right
+
+    def try_k(u, seen):
+        for v in adj[u]:
+            if seen[v]:
+                continue
+            seen[v] = True
+            if match_r[v] == -1 or try_k(match_r[v], seen):
+                match_r[v] = u
+                return True
+        return False
+
+    size = 0
+    for u in range(n_left):
+        seen = [False] * n_right
+        if try_k(u, seen):
+            size += 1
+    return size
+
+
+def test_random():
+    rng = random.Random(42)
+    for _ in range(2000):
+        nL = rng.randint(0, 8)
+        nR = rng.randint(0, 8)
+        m = rng.randint(0, nL * nR if nL and nR else 0)
+        edges = []
+        seen = set()
+        for _ in range(m):
+            u = rng.randrange(nL) if nL else 0
+            v = rng.randrange(nR) if nR else 0
+            if (u, v) not in seen:
+                seen.add((u, v))
+                edges.append((u, v))
+        adj = [[] for _ in range(nL)]
+        for u, v in edges:
+            adj[u].append(v)
+        # shuffle adjacency order to stress
+        for u in range(nL):
+            rng.shuffle(adj[u])
+        size, ml, mr = max_bipartite_matching(adj, nL, nR)
+        ref = kuhn(adj, nL, nR)
+        assert size == ref, (size, ref, adj)
+        # verify matching consistency
+        assert all(ml[u] == -1 or mr[ml[u]] == u for u in range(nL))
+        assert sum(1 for u in range(nL) if ml[u] != -1) == size
+
+
+def test_empty():
+    size, ml, mr = max_bipartite_matching([], 0, 0)
+    assert size == 0 and ml == [] and mr == []
+
+
+def test_complete():
+    n = 5
+    adj = [[v for v in range(n)] for _ in range(n)]
+    size, _, _ = max_bipartite_matching(adj, n, n)
+    assert size == n
+
+
+def test_no_edges():
+    adj = [[] for _ in range(4)]
+    size, _, _ = max_bipartite_matching(adj, 4, 4)
+    assert size == 0
+
+
+def test_duplicates():
+    adj = [[0, 0, 1], [1]]
+    size, _, _ = max_bipartite_matching(adj, 2, 2)
+    assert size == 2