This is an implementation of Edmonds' blossom-contraction algorithm for maximum cardinality matching in general graphs. It's maybe a little long and complex for the recipe book, but I hope it will spare someone else the agony of implementing it themselves.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 | # Find maximum cardinality matching in general undirected graph
# D. Eppstein, UC Irvine, 6 Sep 2003
from __future__ import generators
if 'True' not in globals():
globals()['True'] = not None
globals()['False'] = not True
class unionFind:
'''Union Find data structure. Modified from Josiah Carlson's code,
http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/215912
to allow arbitrarily many arguments in unions, use [] syntax for finds,
and eliminate unnecessary code.'''
def __init__(self):
self.weights = {}
self.parents = {}
def __getitem__(self, object):
'''Find the root of the set that an object is in.
Object must be hashable; previously unknown objects become new singleton sets.'''
# check for previously unknown object
if object not in self.parents:
self.parents[object] = object
self.weights[object] = 1
return object
# find path of objects leading to the root
path = [object]
root = self.parents[object]
while root != path[-1]:
path.append(root)
root = self.parents[root]
# compress the path and return
for ancestor in path:
self.parents[ancestor] = root
return root
def union(self, *objects):
'''Find the sets containing the given objects and merge them all.'''
roots = [self[x] for x in objects]
heaviest = max([(self.weights[r],r) for r in roots])[1]
for r in roots:
if r != heaviest:
self.weights[heaviest] += self.weights[r]
self.parents[r] = heaviest
def matching(G, initialMatching = {}):
'''Find a maximum cardinality matching in a graph G.
G is represented in modified GvR form: iter(G) lists its vertices;
iter(G[v]) lists the neighbors of v; w in G[v] tests adjacency.
The output is a dictionary mapping vertices to their matches;
unmatched vertices are omitted from the dictionary.
We use Edmonds' blossom-contraction algorithm, as described e.g.
in Galil's 1986 Computing Surveys paper.'''
# Copy initial matching so we can use it nondestructively
matching = {}
for x in initialMatching:
matching[x] = initialMatching[x]
# Form greedy matching to avoid some iterations of augmentation
for v in G:
if v not in matching:
for w in G[v]:
if w not in matching:
matching[v] = w
matching[w] = v
break
def augment():
'''Search for a single augmenting path.
Return value is true if the matching size was increased, false otherwise.'''
# Data structures for augmenting path search:
#
# leader: union-find structure; the leader of a blossom is one
# of its vertices (not necessarily topmost), and leader[v] always
# points to the leader of the largest blossom containing v
#
# S: dictionary of leader at even levels of the structure tree.
# Dictionary keys are names of leader (as returned by the union-find
# data structure) and values are the structure tree parent of the blossom
# (a T-node, or the top vertex if the blossom is a root of a structure tree).
#
# T: dictionary of vertices at odd levels of the structure tree.
# Dictionary keys are the vertices; T[x] is a vertex with an unmatched
# edge to x. To find the parent in the structure tree, use leader[T[x]].
#
# unexplored: collection of unexplored vertices within leader of S
#
# base: if x was originally a T-vertex, but becomes part of a blossom,
# base[t] will be the pair (v,w) at the base of the blossom, where v and t
# are on the same side of the blossom and w is on the other side.
leader = unionFind()
S = {}
T = {}
unexplored = []
base = {}
# Subroutines for augmenting path search.
# Many of these are called only from one place, but are split out
# as subroutines to improve modularization and readability.
def blossom(v,w,a):
'''Create a new blossom from edge v-w with common ancestor a.'''
def findSide(v,w):
path = [leader[v]]
b = (v,w) # new base for all T nodes found on the path
while path[-1] != a:
tnode = S[path[-1]]
path.append(tnode)
base[tnode] = b
unexplored.append(tnode)
path.append(leader[T[tnode]])
return path
a = leader[a] # sanity check
path1,path2 = findSide(v,w), findSide(w,v)
leader.union(*path1)
leader.union(*path2)
S[leader[a]] = S[a] # update structure tree
topless = object() # should be unequal to any graph vertex
def alternatingPath(start, goal = topless):
'''Return sequence of vertices on alternating path from start to goal.
Goal must be a T node along the path from the start to the root of the structure tree.
If goal is omitted, we find an alternating path to the structure tree root.'''
path = []
while 1:
while start in T:
v, w = base[start]
vs = alternatingPath(v, start)
vs.reverse()
path += vs
start = w
path.append(start)
if start not in matching:
return path # reached top of structure tree, done!
tnode = matching[start]
path.append(tnode)
if tnode == goal:
return path # finished recursive subpath
start = T[tnode]
def pairs(L):
'''Utility to partition list into pairs of items.
If list has odd length, the final pair is omitted silently.'''
i = 0
while i < len(L) - 1:
yield L[i],L[i+1]
i += 2
def alternate(v):
'''Make v unmatched by alternating the path to the root of its structure tree.'''
path = alternatingPath(v)
path.reverse()
for x,y in pairs(path):
matching[x] = y
matching[y] = x
def addMatch(v, w):
'''Here with an S-S edge vw connecting vertices in different structure trees.
Find the corresponding augmenting path and use it to augment the matching.'''
alternate(v)
alternate(w)
matching[v] = w
matching[w] = v
def ss(v,w):
'''Handle detection of an S-S edge in augmenting path search.
Like augment(), returns true iff the matching size was increased.'''
if leader[v] == leader[w]:
return False # self-loop within blossom, ignore
# parallel search up two branches of structure tree
# until we find a common ancestor of v and w
path1, head1 = {}, v
path2, head2 = {}, w
def step(path, head):
head = leader[head]
parent = leader[S[head]]
if parent == head:
return head # found root of structure tree
path[head] = parent
path[parent] = leader[T[parent]]
return path[parent]
while 1:
head1 = step(path1, head1)
head2 = step(path2, head2)
if head1 == head2:
blossom(v, w, head1)
return False
if leader[S[head1]] == head1 and leader[S[head2]] == head2:
addMatch(v, w)
return True
if head1 in path2:
blossom(v, w, head1)
return False
if head2 in path1:
blossom(v, w, head2)
return False
# Start of main augmenting path search code.
for v in G:
if v not in matching:
S[v] = v
unexplored.append(v)
current = 0 # index into unexplored, in FIFO order so we get short paths
while current < len(unexplored):
v = unexplored[current]
current += 1
for w in G[v]:
if leader[w] in S: # S-S edge: blossom or augmenting path
if ss(v,w):
return True
elif w not in T: # previously unexplored node, add as T-node
T[w] = v
u = matching[w]
if leader[u] not in S:
S[u] = w # and add its match as an S-node
unexplored.append(u)
return False # ran out of graph without finding an augmenting path
# augment the matching until it is maximum
while augment():
pass
return matching
|
My earlier recipe at http://aspn.activestate.com/ASPN/Cookbook/Python/Recipe/123641 solves the same type of matching problem, but is limited to bipartite graphs. I needed a version without that limitation for some graph drawing code. The Edmonds algorithm takes time O(mn alpha(m,n)) where alpha is the inverse Ackermann function; a theoretically faster O(m sqrt(n) algorithm is known but I'm not sure it has ever been implemented.