Welcome, guest | Sign In | My Account | Store | Cart

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.

Python, 248 lines
 ``` 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.

 Created by David Eppstein on Fri, 12 Sep 2003 (PSF)