Edmonds' blossom algorithm in Python
Finds maximum matchings in general graphs by contracting blossoms during augmenting-path search.
- File
edmonds_blossom.py, 141 lines- Level
- I0 not checked yet
- Toolchain
- Python 3.14, pytest 9
- Test framework
- pytest
- Run the tests
python3 -m pytest
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/edmonds-blossom/. The authors are listed in prov/6F205G.jsonl.
from collections import deque
from typing import Callable, NamedTuple, Optional
class Matching(NamedTuple):
size: int
mate: list[int] # mate[v] is the partner of v, or -1 if v is free
phases: int # number of forest searches, the last one included
def edmonds_blossom(
n: int,
edges: list[tuple[int, int]],
on_event: Optional[Callable[[str, tuple], None]] = None,
) -> Matching:
"""Maximum matching of a general graph by Edmonds' blossom algorithm.
Vertices are 0 .. n - 1 and each pair (u, v) in edges is an undirected edge;
the graph is unweighted and need not be bipartite. Self-loops are ignored and
parallel edges are allowed. The optional callback gets ("blossom", (base,
vertices)) for each contraction and ("augment", path) for each augmentation.
At most n/2 + 1 phases of O(n^2 + m) time each, so O(n^3) time for a simple
graph (O(n (n^2 + m)) in general); O(n + m) space.
"""
adj: list[list[int]] = [[] for _ in range(n)]
for u, v in edges:
if not (0 <= u < n and 0 <= v < n):
raise ValueError("edge end out of range")
if u != v:
adj[u].append(v)
adj[v].append(u)
mate = [-1] * n # @line L1
size = 0
phases = 0
while True: # @line L2
phases += 1
path = _phase(n, adj, mate, on_event)
if path is None: # @line L12
return Matching(size, mate, phases)
size += 1
def _phase(n, adj, mate, on_event):
"""Grow an alternating forest from all free vertices; augment once or prove maximality."""
base = list(range(n)) # @line L3
parent = [-1] * n # parent[y] of an inner vertex y, and links inside blossoms
outer = [False] * n
tree = [-1] * n
queue: deque[int] = deque()
for v in range(n):
if mate[v] < 0:
outer[v] = True
tree[v] = v
queue.append(v)
while queue: # @line L4
v = queue.popleft()
for w in adj[v]: # @line L5
if base[v] == base[w] or mate[v] == w: # @line L6
continue
if outer[w] and tree[w] != tree[v]: # @line L7
path = _root_path(v, mate, parent)[::-1] + _root_path(w, mate, parent)
_flip(v, mate, parent) # @line L8
_flip(w, mate, parent)
mate[v], mate[w] = w, v
if on_event is not None:
on_event("augment", tuple(path))
return path
if outer[w]: # @line L9: an edge inside one tree closes a blossom
b = _lca(v, w, base, mate, parent)
in_blossom = [False] * n
_mark_path(v, b, w, base, mate, parent, in_blossom) # @line L10
_mark_path(w, b, v, base, mate, parent, in_blossom)
members = [b]
for x in range(n):
if in_blossom[base[x]]:
if x != b:
members.append(x)
base[x] = b
if not outer[x]:
outer[x] = True
queue.append(x)
if on_event is not None:
on_event("blossom", (b, tuple(sorted(members))))
elif parent[w] < 0: # @line L11: w is not in the forest yet, so it is matched
parent[w] = v
tree[w] = tree[v]
u = mate[w]
outer[u] = True
tree[u] = tree[v]
queue.append(u)
return None
def _lca(a: int, b: int, base: list[int], mate: list[int], parent: list[int]) -> int:
"""Base of the blossom closed by an edge between two outer vertices of one tree."""
seen = set()
while True:
a = base[a]
seen.add(a)
if mate[a] < 0:
break
a = parent[mate[a]]
while True:
b = base[b]
if b in seen:
return b
b = parent[mate[b]]
def _mark_path(v, b, child, base, mate, parent, in_blossom):
"""Walk from v up to the base b, linking the outer vertices back towards child."""
while base[v] != b:
in_blossom[base[v]] = True
in_blossom[base[mate[v]]] = True
parent[v] = child
child = mate[v]
v = parent[mate[v]]
def _root_path(x: int, mate: list[int], parent: list[int]) -> list[int]:
"""The even alternating path from the outer vertex x to the root of its tree."""
path = [x]
y = mate[x]
while y >= 0:
z = parent[y]
path += [y, z]
y = mate[z]
return path
def _flip(x: int, mate: list[int], parent: list[int]) -> None:
"""Flip the path from the outer vertex x to its root; x is left free."""
y = mate[x]
while y >= 0:
z = parent[y]
nxt = mate[z]
mate[y], mate[z] = z, y
y = nxt
test_edmonds_blossom.py 148 lines
# SPDX-License-Identifier: MIT
# Part of https://algorithm.tamnd.com/edmonds-blossom/. The authors are listed in prov/6F205G.jsonl.
import random
import unittest
from functools import lru_cache
from edmonds_blossom import edmonds_blossom
def brute_nu(n: int, edges: list[tuple[int, int]]) -> int:
"""Largest matching size by dynamic programming over vertex subsets."""
adj = [0] * n
for u, v in edges:
if u != v:
adj[u] |= 1 << v
adj[v] |= 1 << u
@lru_cache(maxsize=None)
def best(mask: int) -> int:
if mask == 0:
return 0
v = (mask & -mask).bit_length() - 1
rest = mask & ~(1 << v)
res = best(rest)
cand = adj[v] & rest
while cand:
w = (cand & -cand).bit_length() - 1
cand &= cand - 1
res = max(res, 1 + best(rest & ~(1 << w)))
return res
return best((1 << n) - 1)
def run_checked(tc, n, edges):
"""Run with a callback that checks every reported path against the matching before it."""
es = {(u, v) for u, v in edges} | {(v, u) for u, v in edges}
mate = [-1] * n
events = []
def on_event(kind, data):
events.append((kind, data))
if kind != "augment":
return
path = list(data)
tc.assertEqual(len(set(path)), len(path))
tc.assertEqual(len(path) % 2, 0)
tc.assertEqual(mate[path[0]], -1)
tc.assertEqual(mate[path[-1]], -1)
for i in range(len(path) - 1):
a, b = path[i], path[i + 1]
tc.assertIn((a, b), es)
tc.assertEqual(mate[a] == b, i % 2 == 1)
for i in range(0, len(path), 2):
mate[path[i]], mate[path[i + 1]] = path[i + 1], path[i]
res = edmonds_blossom(n, edges, on_event)
tc.assertEqual(res.mate, mate)
count = 0
for v in range(n):
if res.mate[v] >= 0:
tc.assertEqual(res.mate[res.mate[v]], v)
tc.assertIn((v, res.mate[v]), es)
count += 1
tc.assertEqual(count, 2 * res.size)
return res, events
def cycle(k):
return [(i, (i + 1) % k) for i in range(k)]
PETERSEN = cycle(5) + [(i, i + 5) for i in range(5)] + [(5 + i, 5 + (i + 2) % 5) for i in range(5)]
class TestEdmondsBlossom(unittest.TestCase):
def test_worked_example(self):
# a, c, r, t, d, b = 0, 1, 2, 3, 4, 5; edges ab, cb, ar, cd, rd, at.
edges = [(0, 5), (1, 5), (0, 2), (1, 4), (2, 4), (0, 3)]
res, events = run_checked(self, 6, edges)
self.assertEqual(
events,
[
("augment", (0, 5)),
("augment", (1, 4)),
("blossom", (2, (0, 1, 2, 4, 5))),
("augment", (2, 4, 1, 5, 0, 3)),
],
)
self.assertEqual(res.size, 3)
self.assertEqual(res.phases, 4)
self.assertEqual(res.mate, [3, 5, 4, 0, 2, 1])
def test_empty(self):
self.assertEqual(edmonds_blossom(0, []), (0, [], 1))
self.assertEqual(edmonds_blossom(4, []).size, 0)
def test_odd_cycles(self):
for k in (3, 5, 7, 9):
res, _ = run_checked(self, k, cycle(k))
self.assertEqual(res.size, k // 2)
self.assertEqual(res.mate.count(-1), 1)
def test_no_perfect_matching(self):
# Star K_{1,3}, and two triangles joined at a cut vertex.
self.assertEqual(run_checked(self, 4, [(0, 1), (0, 2), (0, 3)])[0].size, 1)
bowtie = [(0, 1), (1, 2), (2, 0), (2, 3), (3, 4), (4, 2)]
self.assertEqual(run_checked(self, 5, bowtie)[0].size, 2)
# Three triangles hanging off one vertex: removing it leaves 3 odd components.
g = [(0, 1), (1, 2), (2, 0), (3, 4), (4, 5), (5, 3), (6, 7), (7, 8), (8, 6), (9, 0), (9, 3), (9, 6)]
self.assertEqual(run_checked(self, 10, g)[0].size, 4)
def test_petersen_and_complete(self):
self.assertEqual(run_checked(self, 10, PETERSEN)[0].size, 5)
k6 = [(u, v) for u in range(6) for v in range(u + 1, 6)]
self.assertEqual(run_checked(self, 6, k6)[0].size, 3)
def test_bipartite_and_loops(self):
self.assertEqual(run_checked(self, 6, cycle(6))[0].size, 3)
res = edmonds_blossom(3, [(0, 0), (0, 1), (0, 1), (2, 2)])
self.assertEqual(res.size, 1)
def test_nested_blossoms(self):
# A triangle inside a 5-cycle-like structure with a long stem.
edges = [(0, 1), (1, 2), (2, 3), (3, 4), (4, 2), (3, 5), (4, 6), (5, 6), (6, 7), (1, 8), (8, 9)]
res, _ = run_checked(self, 10, edges)
self.assertEqual(res.size, brute_nu(10, edges))
def test_bad_edge(self):
with self.assertRaises(ValueError):
edmonds_blossom(2, [(0, 2)])
def test_random_against_brute_force(self):
rng = random.Random(1965)
for _ in range(700):
n = rng.randint(0, 11)
p = rng.random() * 0.7
edges = [(u, v) for u in range(n) for v in range(u + 1, n) if rng.random() < p]
rng.shuffle(edges)
edges = [(v, u) if rng.random() < 0.5 else (u, v) for u, v in edges]
res, events = run_checked(self, n, edges)
self.assertEqual(res.size, brute_nu(n, edges))
self.assertEqual(res.phases, res.size + 1)
if __name__ == "__main__":
unittest.main()
What the next level needs
- the tests did not run in this build