Mercurial > repos > guerler > springsuite
comparison planemo/lib/python3.7/site-packages/networkx/algorithms/matching.py @ 1:56ad4e20f292 draft
"planemo upload commit 6eee67778febed82ddd413c3ca40b3183a3898f1"
author | guerler |
---|---|
date | Fri, 31 Jul 2020 00:32:28 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
0:d30785e31577 | 1:56ad4e20f292 |
---|---|
1 # Copyright 2016 NetworkX developers. | |
2 # Copyright (C) 2004-2019 by | |
3 # Aric Hagberg <hagberg@lanl.gov> | |
4 # Dan Schult <dschult@colgate.edu> | |
5 # Pieter Swart <swart@lanl.gov> | |
6 # | |
7 # Copyright (C) 2008 by | |
8 # Joris van Rantwijk. | |
9 # | |
10 # Copyright (C) 2011 by | |
11 # Nicholas Mancuso <nick.mancuso@gmail.com> | |
12 # | |
13 # All rights reserved. | |
14 # BSD license. | |
15 """Functions for computing and verifying matchings in a graph.""" | |
16 from collections import Counter | |
17 from itertools import combinations | |
18 from itertools import repeat | |
19 | |
20 __all__ = ['is_matching', 'is_maximal_matching', 'is_perfect_matching', | |
21 'max_weight_matching', 'maximal_matching'] | |
22 | |
23 | |
24 def maximal_matching(G): | |
25 r"""Find a maximal matching in the graph. | |
26 | |
27 A matching is a subset of edges in which no node occurs more than once. | |
28 A maximal matching cannot add more edges and still be a matching. | |
29 | |
30 Parameters | |
31 ---------- | |
32 G : NetworkX graph | |
33 Undirected graph | |
34 | |
35 Returns | |
36 ------- | |
37 matching : set | |
38 A maximal matching of the graph. | |
39 | |
40 Notes | |
41 ----- | |
42 The algorithm greedily selects a maximal matching M of the graph G | |
43 (i.e. no superset of M exists). It runs in $O(|E|)$ time. | |
44 """ | |
45 matching = set() | |
46 nodes = set() | |
47 for u, v in G.edges(): | |
48 # If the edge isn't covered, add it to the matching | |
49 # then remove neighborhood of u and v from consideration. | |
50 if u not in nodes and v not in nodes and u != v: | |
51 matching.add((u, v)) | |
52 nodes.add(u) | |
53 nodes.add(v) | |
54 return matching | |
55 | |
56 | |
57 def matching_dict_to_set(matching): | |
58 """Converts a dictionary representing a matching (as returned by | |
59 :func:`max_weight_matching`) to a set representing a matching (as | |
60 returned by :func:`maximal_matching`). | |
61 | |
62 In the definition of maximal matching adopted by NetworkX, | |
63 self-loops are not allowed, so the provided dictionary is expected | |
64 to never have any mapping from a key to itself. However, the | |
65 dictionary is expected to have mirrored key/value pairs, for | |
66 example, key ``u`` with value ``v`` and key ``v`` with value ``u``. | |
67 | |
68 """ | |
69 # Need to compensate for the fact that both pairs (u, v) and (v, u) | |
70 # appear in `matching.items()`, so we use a set of sets. This way, | |
71 # only the (frozen)set `{u, v}` appears as an element in the | |
72 # returned set. | |
73 | |
74 return set((u, v) for (u, v) in set(map(frozenset, matching.items()))) | |
75 | |
76 | |
77 def is_matching(G, matching): | |
78 """Decides whether the given set or dictionary represents a valid | |
79 matching in ``G``. | |
80 | |
81 A *matching* in a graph is a set of edges in which no two distinct | |
82 edges share a common endpoint. | |
83 | |
84 Parameters | |
85 ---------- | |
86 G : NetworkX graph | |
87 | |
88 matching : dict or set | |
89 A dictionary or set representing a matching. If a dictionary, it | |
90 must have ``matching[u] == v`` and ``matching[v] == u`` for each | |
91 edge ``(u, v)`` in the matching. If a set, it must have elements | |
92 of the form ``(u, v)``, where ``(u, v)`` is an edge in the | |
93 matching. | |
94 | |
95 Returns | |
96 ------- | |
97 bool | |
98 Whether the given set or dictionary represents a valid matching | |
99 in the graph. | |
100 | |
101 """ | |
102 if isinstance(matching, dict): | |
103 matching = matching_dict_to_set(matching) | |
104 # TODO This is parallelizable. | |
105 return all(len(set(e1) & set(e2)) == 0 | |
106 for e1, e2 in combinations(matching, 2)) | |
107 | |
108 | |
109 def is_maximal_matching(G, matching): | |
110 """Decides whether the given set or dictionary represents a valid | |
111 maximal matching in ``G``. | |
112 | |
113 A *maximal matching* in a graph is a matching in which adding any | |
114 edge would cause the set to no longer be a valid matching. | |
115 | |
116 Parameters | |
117 ---------- | |
118 G : NetworkX graph | |
119 | |
120 matching : dict or set | |
121 A dictionary or set representing a matching. If a dictionary, it | |
122 must have ``matching[u] == v`` and ``matching[v] == u`` for each | |
123 edge ``(u, v)`` in the matching. If a set, it must have elements | |
124 of the form ``(u, v)``, where ``(u, v)`` is an edge in the | |
125 matching. | |
126 | |
127 Returns | |
128 ------- | |
129 bool | |
130 Whether the given set or dictionary represents a valid maximal | |
131 matching in the graph. | |
132 | |
133 """ | |
134 if isinstance(matching, dict): | |
135 matching = matching_dict_to_set(matching) | |
136 # If the given set is not a matching, then it is not a maximal matching. | |
137 if not is_matching(G, matching): | |
138 return False | |
139 # A matching is maximal if adding any unmatched edge to it causes | |
140 # the resulting set to *not* be a valid matching. | |
141 # | |
142 # HACK Since the `matching_dict_to_set` function returns a set of | |
143 # sets, we need to convert the list of edges to a set of sets in | |
144 # order for the set difference function to work. Ideally, we would | |
145 # just be able to do `set(G.edges()) - matching`. | |
146 all_edges = set(map(frozenset, G.edges())) | |
147 matched_edges = set(map(frozenset, matching)) | |
148 unmatched_edges = all_edges - matched_edges | |
149 # TODO This is parallelizable. | |
150 return all(not is_matching(G, matching | {e}) for e in unmatched_edges) | |
151 | |
152 | |
153 def is_perfect_matching(G, matching): | |
154 """Decides whether the given set represents a valid perfect matching in | |
155 ``G``. | |
156 | |
157 A *perfect matching* in a graph is a matching in which exactly one edge | |
158 is incident upon each vertex. | |
159 | |
160 Parameters | |
161 ---------- | |
162 G : NetworkX graph | |
163 | |
164 matching : dict or set | |
165 A dictionary or set representing a matching. If a dictionary, it | |
166 must have ``matching[u] == v`` and ``matching[v] == u`` for each | |
167 edge ``(u, v)`` in the matching. If a set, it must have elements | |
168 of the form ``(u, v)``, where ``(u, v)`` is an edge in the | |
169 matching. | |
170 | |
171 Returns | |
172 ------- | |
173 bool | |
174 Whether the given set or dictionary represents a valid perfect | |
175 matching in the graph. | |
176 | |
177 """ | |
178 if isinstance(matching, dict): | |
179 matching = matching_dict_to_set(matching) | |
180 | |
181 if not is_matching(G, matching): | |
182 return False | |
183 | |
184 counts = Counter(sum(matching, ())) | |
185 | |
186 return all(counts[v] == 1 for v in G) | |
187 | |
188 | |
189 def max_weight_matching(G, maxcardinality=False, weight='weight'): | |
190 """Compute a maximum-weighted matching of G. | |
191 | |
192 A matching is a subset of edges in which no node occurs more than once. | |
193 The weight of a matching is the sum of the weights of its edges. | |
194 A maximal matching cannot add more edges and still be a matching. | |
195 The cardinality of a matching is the number of matched edges. | |
196 | |
197 Parameters | |
198 ---------- | |
199 G : NetworkX graph | |
200 Undirected graph | |
201 | |
202 maxcardinality: bool, optional (default=False) | |
203 If maxcardinality is True, compute the maximum-cardinality matching | |
204 with maximum weight among all maximum-cardinality matchings. | |
205 | |
206 weight: string, optional (default='weight') | |
207 Edge data key corresponding to the edge weight. | |
208 If key not found, uses 1 as weight. | |
209 | |
210 | |
211 Returns | |
212 ------- | |
213 matching : set | |
214 A maximal matching of the graph. | |
215 | |
216 Notes | |
217 ----- | |
218 If G has edges with weight attributes the edge data are used as | |
219 weight values else the weights are assumed to be 1. | |
220 | |
221 This function takes time O(number_of_nodes ** 3). | |
222 | |
223 If all edge weights are integers, the algorithm uses only integer | |
224 computations. If floating point weights are used, the algorithm | |
225 could return a slightly suboptimal matching due to numeric | |
226 precision errors. | |
227 | |
228 This method is based on the "blossom" method for finding augmenting | |
229 paths and the "primal-dual" method for finding a matching of maximum | |
230 weight, both methods invented by Jack Edmonds [1]_. | |
231 | |
232 Bipartite graphs can also be matched using the functions present in | |
233 :mod:`networkx.algorithms.bipartite.matching`. | |
234 | |
235 References | |
236 ---------- | |
237 .. [1] "Efficient Algorithms for Finding Maximum Matching in Graphs", | |
238 Zvi Galil, ACM Computing Surveys, 1986. | |
239 """ | |
240 # | |
241 # The algorithm is taken from "Efficient Algorithms for Finding Maximum | |
242 # Matching in Graphs" by Zvi Galil, ACM Computing Surveys, 1986. | |
243 # It is based on the "blossom" method for finding augmenting paths and | |
244 # the "primal-dual" method for finding a matching of maximum weight, both | |
245 # methods invented by Jack Edmonds. | |
246 # | |
247 # A C program for maximum weight matching by Ed Rothberg was used | |
248 # extensively to validate this new code. | |
249 # | |
250 # Many terms used in the code comments are explained in the paper | |
251 # by Galil. You will probably need the paper to make sense of this code. | |
252 # | |
253 | |
254 class NoNode: | |
255 """Dummy value which is different from any node.""" | |
256 pass | |
257 | |
258 class Blossom: | |
259 """Representation of a non-trivial blossom or sub-blossom.""" | |
260 | |
261 __slots__ = ['childs', 'edges', 'mybestedges'] | |
262 | |
263 # b.childs is an ordered list of b's sub-blossoms, starting with | |
264 # the base and going round the blossom. | |
265 | |
266 # b.edges is the list of b's connecting edges, such that | |
267 # b.edges[i] = (v, w) where v is a vertex in b.childs[i] | |
268 # and w is a vertex in b.childs[wrap(i+1)]. | |
269 | |
270 # If b is a top-level S-blossom, | |
271 # b.mybestedges is a list of least-slack edges to neighbouring | |
272 # S-blossoms, or None if no such list has been computed yet. | |
273 # This is used for efficient computation of delta3. | |
274 | |
275 # Generate the blossom's leaf vertices. | |
276 def leaves(self): | |
277 for t in self.childs: | |
278 if isinstance(t, Blossom): | |
279 for v in t.leaves(): | |
280 yield v | |
281 else: | |
282 yield t | |
283 | |
284 # Get a list of vertices. | |
285 gnodes = list(G) | |
286 if not gnodes: | |
287 return set() # don't bother with empty graphs | |
288 | |
289 # Find the maximum edge weight. | |
290 maxweight = 0 | |
291 allinteger = True | |
292 for i, j, d in G.edges(data=True): | |
293 wt = d.get(weight, 1) | |
294 if i != j and wt > maxweight: | |
295 maxweight = wt | |
296 allinteger = allinteger and (str(type(wt)).split("'")[1] | |
297 in ('int', 'long')) | |
298 | |
299 # If v is a matched vertex, mate[v] is its partner vertex. | |
300 # If v is a single vertex, v does not occur as a key in mate. | |
301 # Initially all vertices are single; updated during augmentation. | |
302 mate = {} | |
303 | |
304 # If b is a top-level blossom, | |
305 # label.get(b) is None if b is unlabeled (free), | |
306 # 1 if b is an S-blossom, | |
307 # 2 if b is a T-blossom. | |
308 # The label of a vertex is found by looking at the label of its top-level | |
309 # containing blossom. | |
310 # If v is a vertex inside a T-blossom, label[v] is 2 iff v is reachable | |
311 # from an S-vertex outside the blossom. | |
312 # Labels are assigned during a stage and reset after each augmentation. | |
313 label = {} | |
314 | |
315 # If b is a labeled top-level blossom, | |
316 # labeledge[b] = (v, w) is the edge through which b obtained its label | |
317 # such that w is a vertex in b, or None if b's base vertex is single. | |
318 # If w is a vertex inside a T-blossom and label[w] == 2, | |
319 # labeledge[w] = (v, w) is an edge through which w is reachable from | |
320 # outside the blossom. | |
321 labeledge = {} | |
322 | |
323 # If v is a vertex, inblossom[v] is the top-level blossom to which v | |
324 # belongs. | |
325 # If v is a top-level vertex, inblossom[v] == v since v is itself | |
326 # a (trivial) top-level blossom. | |
327 # Initially all vertices are top-level trivial blossoms. | |
328 inblossom = dict(zip(gnodes, gnodes)) | |
329 | |
330 # If b is a sub-blossom, | |
331 # blossomparent[b] is its immediate parent (sub-)blossom. | |
332 # If b is a top-level blossom, blossomparent[b] is None. | |
333 blossomparent = dict(zip(gnodes, repeat(None))) | |
334 | |
335 # If b is a (sub-)blossom, | |
336 # blossombase[b] is its base VERTEX (i.e. recursive sub-blossom). | |
337 blossombase = dict(zip(gnodes, gnodes)) | |
338 | |
339 # If w is a free vertex (or an unreached vertex inside a T-blossom), | |
340 # bestedge[w] = (v, w) is the least-slack edge from an S-vertex, | |
341 # or None if there is no such edge. | |
342 # If b is a (possibly trivial) top-level S-blossom, | |
343 # bestedge[b] = (v, w) is the least-slack edge to a different S-blossom | |
344 # (v inside b), or None if there is no such edge. | |
345 # This is used for efficient computation of delta2 and delta3. | |
346 bestedge = {} | |
347 | |
348 # If v is a vertex, | |
349 # dualvar[v] = 2 * u(v) where u(v) is the v's variable in the dual | |
350 # optimization problem (if all edge weights are integers, multiplication | |
351 # by two ensures that all values remain integers throughout the algorithm). | |
352 # Initially, u(v) = maxweight / 2. | |
353 dualvar = dict(zip(gnodes, repeat(maxweight))) | |
354 | |
355 # If b is a non-trivial blossom, | |
356 # blossomdual[b] = z(b) where z(b) is b's variable in the dual | |
357 # optimization problem. | |
358 blossomdual = {} | |
359 | |
360 # If (v, w) in allowedge or (w, v) in allowedg, then the edge | |
361 # (v, w) is known to have zero slack in the optimization problem; | |
362 # otherwise the edge may or may not have zero slack. | |
363 allowedge = {} | |
364 | |
365 # Queue of newly discovered S-vertices. | |
366 queue = [] | |
367 | |
368 # Return 2 * slack of edge (v, w) (does not work inside blossoms). | |
369 def slack(v, w): | |
370 return dualvar[v] + dualvar[w] - 2 * G[v][w].get(weight, 1) | |
371 | |
372 # Assign label t to the top-level blossom containing vertex w, | |
373 # coming through an edge from vertex v. | |
374 def assignLabel(w, t, v): | |
375 b = inblossom[w] | |
376 assert label.get(w) is None and label.get(b) is None | |
377 label[w] = label[b] = t | |
378 if v is not None: | |
379 labeledge[w] = labeledge[b] = (v, w) | |
380 else: | |
381 labeledge[w] = labeledge[b] = None | |
382 bestedge[w] = bestedge[b] = None | |
383 if t == 1: | |
384 # b became an S-vertex/blossom; add it(s vertices) to the queue. | |
385 if isinstance(b, Blossom): | |
386 queue.extend(b.leaves()) | |
387 else: | |
388 queue.append(b) | |
389 elif t == 2: | |
390 # b became a T-vertex/blossom; assign label S to its mate. | |
391 # (If b is a non-trivial blossom, its base is the only vertex | |
392 # with an external mate.) | |
393 base = blossombase[b] | |
394 assignLabel(mate[base], 1, base) | |
395 | |
396 # Trace back from vertices v and w to discover either a new blossom | |
397 # or an augmenting path. Return the base vertex of the new blossom, | |
398 # or NoNode if an augmenting path was found. | |
399 def scanBlossom(v, w): | |
400 # Trace back from v and w, placing breadcrumbs as we go. | |
401 path = [] | |
402 base = NoNode | |
403 while v is not NoNode: | |
404 # Look for a breadcrumb in v's blossom or put a new breadcrumb. | |
405 b = inblossom[v] | |
406 if label[b] & 4: | |
407 base = blossombase[b] | |
408 break | |
409 assert label[b] == 1 | |
410 path.append(b) | |
411 label[b] = 5 | |
412 # Trace one step back. | |
413 if labeledge[b] is None: | |
414 # The base of blossom b is single; stop tracing this path. | |
415 assert blossombase[b] not in mate | |
416 v = NoNode | |
417 else: | |
418 assert labeledge[b][0] == mate[blossombase[b]] | |
419 v = labeledge[b][0] | |
420 b = inblossom[v] | |
421 assert label[b] == 2 | |
422 # b is a T-blossom; trace one more step back. | |
423 v = labeledge[b][0] | |
424 # Swap v and w so that we alternate between both paths. | |
425 if w is not NoNode: | |
426 v, w = w, v | |
427 # Remove breadcrumbs. | |
428 for b in path: | |
429 label[b] = 1 | |
430 # Return base vertex, if we found one. | |
431 return base | |
432 | |
433 # Construct a new blossom with given base, through S-vertices v and w. | |
434 # Label the new blossom as S; set its dual variable to zero; | |
435 # relabel its T-vertices to S and add them to the queue. | |
436 def addBlossom(base, v, w): | |
437 bb = inblossom[base] | |
438 bv = inblossom[v] | |
439 bw = inblossom[w] | |
440 # Create blossom. | |
441 b = Blossom() | |
442 blossombase[b] = base | |
443 blossomparent[b] = None | |
444 blossomparent[bb] = b | |
445 # Make list of sub-blossoms and their interconnecting edge endpoints. | |
446 b.childs = path = [] | |
447 b.edges = edgs = [(v, w)] | |
448 # Trace back from v to base. | |
449 while bv != bb: | |
450 # Add bv to the new blossom. | |
451 blossomparent[bv] = b | |
452 path.append(bv) | |
453 edgs.append(labeledge[bv]) | |
454 assert label[bv] == 2 or (label[bv] == 1 and labeledge[ | |
455 bv][0] == mate[blossombase[bv]]) | |
456 # Trace one step back. | |
457 v = labeledge[bv][0] | |
458 bv = inblossom[v] | |
459 # Add base sub-blossom; reverse lists. | |
460 path.append(bb) | |
461 path.reverse() | |
462 edgs.reverse() | |
463 # Trace back from w to base. | |
464 while bw != bb: | |
465 # Add bw to the new blossom. | |
466 blossomparent[bw] = b | |
467 path.append(bw) | |
468 edgs.append((labeledge[bw][1], labeledge[bw][0])) | |
469 assert label[bw] == 2 or (label[bw] == 1 and labeledge[ | |
470 bw][0] == mate[blossombase[bw]]) | |
471 # Trace one step back. | |
472 w = labeledge[bw][0] | |
473 bw = inblossom[w] | |
474 # Set label to S. | |
475 assert label[bb] == 1 | |
476 label[b] = 1 | |
477 labeledge[b] = labeledge[bb] | |
478 # Set dual variable to zero. | |
479 blossomdual[b] = 0 | |
480 # Relabel vertices. | |
481 for v in b.leaves(): | |
482 if label[inblossom[v]] == 2: | |
483 # This T-vertex now turns into an S-vertex because it becomes | |
484 # part of an S-blossom; add it to the queue. | |
485 queue.append(v) | |
486 inblossom[v] = b | |
487 # Compute b.mybestedges. | |
488 bestedgeto = {} | |
489 for bv in path: | |
490 if isinstance(bv, Blossom): | |
491 if bv.mybestedges is not None: | |
492 # Walk this subblossom's least-slack edges. | |
493 nblist = bv.mybestedges | |
494 # The sub-blossom won't need this data again. | |
495 bv.mybestedges = None | |
496 else: | |
497 # This subblossom does not have a list of least-slack | |
498 # edges; get the information from the vertices. | |
499 nblist = [(v, w) | |
500 for v in bv.leaves() | |
501 for w in G.neighbors(v) | |
502 if v != w] | |
503 else: | |
504 nblist = [(bv, w) | |
505 for w in G.neighbors(bv) | |
506 if bv != w] | |
507 for k in nblist: | |
508 (i, j) = k | |
509 if inblossom[j] == b: | |
510 i, j = j, i | |
511 bj = inblossom[j] | |
512 if (bj != b and label.get(bj) == 1 and | |
513 ((bj not in bestedgeto) or | |
514 slack(i, j) < slack(*bestedgeto[bj]))): | |
515 bestedgeto[bj] = k | |
516 # Forget about least-slack edge of the subblossom. | |
517 bestedge[bv] = None | |
518 b.mybestedges = list(bestedgeto.values()) | |
519 # Select bestedge[b]. | |
520 mybestedge = None | |
521 bestedge[b] = None | |
522 for k in b.mybestedges: | |
523 kslack = slack(*k) | |
524 if mybestedge is None or kslack < mybestslack: | |
525 mybestedge = k | |
526 mybestslack = kslack | |
527 bestedge[b] = mybestedge | |
528 | |
529 # Expand the given top-level blossom. | |
530 def expandBlossom(b, endstage): | |
531 # Convert sub-blossoms into top-level blossoms. | |
532 for s in b.childs: | |
533 blossomparent[s] = None | |
534 if isinstance(s, Blossom): | |
535 if endstage and blossomdual[s] == 0: | |
536 # Recursively expand this sub-blossom. | |
537 expandBlossom(s, endstage) | |
538 else: | |
539 for v in s.leaves(): | |
540 inblossom[v] = s | |
541 else: | |
542 inblossom[s] = s | |
543 # If we expand a T-blossom during a stage, its sub-blossoms must be | |
544 # relabeled. | |
545 if (not endstage) and label.get(b) == 2: | |
546 # Start at the sub-blossom through which the expanding | |
547 # blossom obtained its label, and relabel sub-blossoms untili | |
548 # we reach the base. | |
549 # Figure out through which sub-blossom the expanding blossom | |
550 # obtained its label initially. | |
551 entrychild = inblossom[labeledge[b][1]] | |
552 # Decide in which direction we will go round the blossom. | |
553 j = b.childs.index(entrychild) | |
554 if j & 1: | |
555 # Start index is odd; go forward and wrap. | |
556 j -= len(b.childs) | |
557 jstep = 1 | |
558 else: | |
559 # Start index is even; go backward. | |
560 jstep = -1 | |
561 # Move along the blossom until we get to the base. | |
562 v, w = labeledge[b] | |
563 while j != 0: | |
564 # Relabel the T-sub-blossom. | |
565 if jstep == 1: | |
566 p, q = b.edges[j] | |
567 else: | |
568 q, p = b.edges[j - 1] | |
569 label[w] = None | |
570 label[q] = None | |
571 assignLabel(w, 2, v) | |
572 # Step to the next S-sub-blossom and note its forward edge. | |
573 allowedge[(p, q)] = allowedge[(q, p)] = True | |
574 j += jstep | |
575 if jstep == 1: | |
576 v, w = b.edges[j] | |
577 else: | |
578 w, v = b.edges[j - 1] | |
579 # Step to the next T-sub-blossom. | |
580 allowedge[(v, w)] = allowedge[(w, v)] = True | |
581 j += jstep | |
582 # Relabel the base T-sub-blossom WITHOUT stepping through to | |
583 # its mate (so don't call assignLabel). | |
584 bw = b.childs[j] | |
585 label[w] = label[bw] = 2 | |
586 labeledge[w] = labeledge[bw] = (v, w) | |
587 bestedge[bw] = None | |
588 # Continue along the blossom until we get back to entrychild. | |
589 j += jstep | |
590 while b.childs[j] != entrychild: | |
591 # Examine the vertices of the sub-blossom to see whether | |
592 # it is reachable from a neighbouring S-vertex outside the | |
593 # expanding blossom. | |
594 bv = b.childs[j] | |
595 if label.get(bv) == 1: | |
596 # This sub-blossom just got label S through one of its | |
597 # neighbours; leave it be. | |
598 j += jstep | |
599 continue | |
600 if isinstance(bv, Blossom): | |
601 for v in bv.leaves(): | |
602 if label.get(v): | |
603 break | |
604 else: | |
605 v = bv | |
606 # If the sub-blossom contains a reachable vertex, assign | |
607 # label T to the sub-blossom. | |
608 if label.get(v): | |
609 assert label[v] == 2 | |
610 assert inblossom[v] == bv | |
611 label[v] = None | |
612 label[mate[blossombase[bv]]] = None | |
613 assignLabel(v, 2, labeledge[v][0]) | |
614 j += jstep | |
615 # Remove the expanded blossom entirely. | |
616 label.pop(b, None) | |
617 labeledge.pop(b, None) | |
618 bestedge.pop(b, None) | |
619 del blossomparent[b] | |
620 del blossombase[b] | |
621 del blossomdual[b] | |
622 | |
623 # Swap matched/unmatched edges over an alternating path through blossom b | |
624 # between vertex v and the base vertex. Keep blossom bookkeeping | |
625 # consistent. | |
626 def augmentBlossom(b, v): | |
627 # Bubble up through the blossom tree from vertex v to an immediate | |
628 # sub-blossom of b. | |
629 t = v | |
630 while blossomparent[t] != b: | |
631 t = blossomparent[t] | |
632 # Recursively deal with the first sub-blossom. | |
633 if isinstance(t, Blossom): | |
634 augmentBlossom(t, v) | |
635 # Decide in which direction we will go round the blossom. | |
636 i = j = b.childs.index(t) | |
637 if i & 1: | |
638 # Start index is odd; go forward and wrap. | |
639 j -= len(b.childs) | |
640 jstep = 1 | |
641 else: | |
642 # Start index is even; go backward. | |
643 jstep = -1 | |
644 # Move along the blossom until we get to the base. | |
645 while j != 0: | |
646 # Step to the next sub-blossom and augment it recursively. | |
647 j += jstep | |
648 t = b.childs[j] | |
649 if jstep == 1: | |
650 w, x = b.edges[j] | |
651 else: | |
652 x, w = b.edges[j - 1] | |
653 if isinstance(t, Blossom): | |
654 augmentBlossom(t, w) | |
655 # Step to the next sub-blossom and augment it recursively. | |
656 j += jstep | |
657 t = b.childs[j] | |
658 if isinstance(t, Blossom): | |
659 augmentBlossom(t, x) | |
660 # Match the edge connecting those sub-blossoms. | |
661 mate[w] = x | |
662 mate[x] = w | |
663 # Rotate the list of sub-blossoms to put the new base at the front. | |
664 b.childs = b.childs[i:] + b.childs[:i] | |
665 b.edges = b.edges[i:] + b.edges[:i] | |
666 blossombase[b] = blossombase[b.childs[0]] | |
667 assert blossombase[b] == v | |
668 | |
669 # Swap matched/unmatched edges over an alternating path between two | |
670 # single vertices. The augmenting path runs through S-vertices v and w. | |
671 def augmentMatching(v, w): | |
672 for (s, j) in ((v, w), (w, v)): | |
673 # Match vertex s to vertex j. Then trace back from s | |
674 # until we find a single vertex, swapping matched and unmatched | |
675 # edges as we go. | |
676 while 1: | |
677 bs = inblossom[s] | |
678 assert label[bs] == 1 | |
679 assert ( | |
680 labeledge[bs] is None and blossombase[bs] not in mate)\ | |
681 or (labeledge[bs][0] == mate[blossombase[bs]]) | |
682 # Augment through the S-blossom from s to base. | |
683 if isinstance(bs, Blossom): | |
684 augmentBlossom(bs, s) | |
685 # Update mate[s] | |
686 mate[s] = j | |
687 # Trace one step back. | |
688 if labeledge[bs] is None: | |
689 # Reached single vertex; stop. | |
690 break | |
691 t = labeledge[bs][0] | |
692 bt = inblossom[t] | |
693 assert label[bt] == 2 | |
694 # Trace one more step back. | |
695 s, j = labeledge[bt] | |
696 # Augment through the T-blossom from j to base. | |
697 assert blossombase[bt] == t | |
698 if isinstance(bt, Blossom): | |
699 augmentBlossom(bt, j) | |
700 # Update mate[j] | |
701 mate[j] = s | |
702 | |
703 # Verify that the optimum solution has been reached. | |
704 def verifyOptimum(): | |
705 if maxcardinality: | |
706 # Vertices may have negative dual; | |
707 # find a constant non-negative number to add to all vertex duals. | |
708 vdualoffset = max(0, -min(dualvar.values())) | |
709 else: | |
710 vdualoffset = 0 | |
711 # 0. all dual variables are non-negative | |
712 assert min(dualvar.values()) + vdualoffset >= 0 | |
713 assert len(blossomdual) == 0 or min(blossomdual.values()) >= 0 | |
714 # 0. all edges have non-negative slack and | |
715 # 1. all matched edges have zero slack; | |
716 for i, j, d in G.edges(data=True): | |
717 wt = d.get(weight, 1) | |
718 if i == j: | |
719 continue # ignore self-loops | |
720 s = dualvar[i] + dualvar[j] - 2 * wt | |
721 iblossoms = [i] | |
722 jblossoms = [j] | |
723 while blossomparent[iblossoms[-1]] is not None: | |
724 iblossoms.append(blossomparent[iblossoms[-1]]) | |
725 while blossomparent[jblossoms[-1]] is not None: | |
726 jblossoms.append(blossomparent[jblossoms[-1]]) | |
727 iblossoms.reverse() | |
728 jblossoms.reverse() | |
729 for (bi, bj) in zip(iblossoms, jblossoms): | |
730 if bi != bj: | |
731 break | |
732 s += 2 * blossomdual[bi] | |
733 assert s >= 0 | |
734 if mate.get(i) == j or mate.get(j) == i: | |
735 assert mate[i] == j and mate[j] == i | |
736 assert s == 0 | |
737 # 2. all single vertices have zero dual value; | |
738 for v in gnodes: | |
739 assert (v in mate) or dualvar[v] + vdualoffset == 0 | |
740 # 3. all blossoms with positive dual value are full. | |
741 for b in blossomdual: | |
742 if blossomdual[b] > 0: | |
743 assert len(b.edges) % 2 == 1 | |
744 for (i, j) in b.edges[1::2]: | |
745 assert mate[i] == j and mate[j] == i | |
746 # Ok. | |
747 | |
748 # Main loop: continue until no further improvement is possible. | |
749 while 1: | |
750 | |
751 # Each iteration of this loop is a "stage". | |
752 # A stage finds an augmenting path and uses that to improve | |
753 # the matching. | |
754 | |
755 # Remove labels from top-level blossoms/vertices. | |
756 label.clear() | |
757 labeledge.clear() | |
758 | |
759 # Forget all about least-slack edges. | |
760 bestedge.clear() | |
761 for b in blossomdual: | |
762 b.mybestedges = None | |
763 | |
764 # Loss of labeling means that we can not be sure that currently | |
765 # allowable edges remain allowable throughout this stage. | |
766 allowedge.clear() | |
767 | |
768 # Make queue empty. | |
769 queue[:] = [] | |
770 | |
771 # Label single blossoms/vertices with S and put them in the queue. | |
772 for v in gnodes: | |
773 if (v not in mate) and label.get(inblossom[v]) is None: | |
774 assignLabel(v, 1, None) | |
775 | |
776 # Loop until we succeed in augmenting the matching. | |
777 augmented = 0 | |
778 while 1: | |
779 | |
780 # Each iteration of this loop is a "substage". | |
781 # A substage tries to find an augmenting path; | |
782 # if found, the path is used to improve the matching and | |
783 # the stage ends. If there is no augmenting path, the | |
784 # primal-dual method is used to pump some slack out of | |
785 # the dual variables. | |
786 | |
787 # Continue labeling until all vertices which are reachable | |
788 # through an alternating path have got a label. | |
789 while queue and not augmented: | |
790 | |
791 # Take an S vertex from the queue. | |
792 v = queue.pop() | |
793 assert label[inblossom[v]] == 1 | |
794 | |
795 # Scan its neighbours: | |
796 for w in G.neighbors(v): | |
797 if w == v: | |
798 continue # ignore self-loops | |
799 # w is a neighbour to v | |
800 bv = inblossom[v] | |
801 bw = inblossom[w] | |
802 if bv == bw: | |
803 # this edge is internal to a blossom; ignore it | |
804 continue | |
805 if (v, w) not in allowedge: | |
806 kslack = slack(v, w) | |
807 if kslack <= 0: | |
808 # edge k has zero slack => it is allowable | |
809 allowedge[(v, w)] = allowedge[(w, v)] = True | |
810 if (v, w) in allowedge: | |
811 if label.get(bw) is None: | |
812 # (C1) w is a free vertex; | |
813 # label w with T and label its mate with S (R12). | |
814 assignLabel(w, 2, v) | |
815 elif label.get(bw) == 1: | |
816 # (C2) w is an S-vertex (not in the same blossom); | |
817 # follow back-links to discover either an | |
818 # augmenting path or a new blossom. | |
819 base = scanBlossom(v, w) | |
820 if base is not NoNode: | |
821 # Found a new blossom; add it to the blossom | |
822 # bookkeeping and turn it into an S-blossom. | |
823 addBlossom(base, v, w) | |
824 else: | |
825 # Found an augmenting path; augment the | |
826 # matching and end this stage. | |
827 augmentMatching(v, w) | |
828 augmented = 1 | |
829 break | |
830 elif label.get(w) is None: | |
831 # w is inside a T-blossom, but w itself has not | |
832 # yet been reached from outside the blossom; | |
833 # mark it as reached (we need this to relabel | |
834 # during T-blossom expansion). | |
835 assert label[bw] == 2 | |
836 label[w] = 2 | |
837 labeledge[w] = (v, w) | |
838 elif label.get(bw) == 1: | |
839 # keep track of the least-slack non-allowable edge to | |
840 # a different S-blossom. | |
841 if bestedge.get(bv) is None or \ | |
842 kslack < slack(*bestedge[bv]): | |
843 bestedge[bv] = (v, w) | |
844 elif label.get(w) is None: | |
845 # w is a free vertex (or an unreached vertex inside | |
846 # a T-blossom) but we can not reach it yet; | |
847 # keep track of the least-slack edge that reaches w. | |
848 if bestedge.get(w) is None or \ | |
849 kslack < slack(*bestedge[w]): | |
850 bestedge[w] = (v, w) | |
851 | |
852 if augmented: | |
853 break | |
854 | |
855 # There is no augmenting path under these constraints; | |
856 # compute delta and reduce slack in the optimization problem. | |
857 # (Note that our vertex dual variables, edge slacks and delta's | |
858 # are pre-multiplied by two.) | |
859 deltatype = -1 | |
860 delta = deltaedge = deltablossom = None | |
861 | |
862 # Compute delta1: the minimum value of any vertex dual. | |
863 if not maxcardinality: | |
864 deltatype = 1 | |
865 delta = min(dualvar.values()) | |
866 | |
867 # Compute delta2: the minimum slack on any edge between | |
868 # an S-vertex and a free vertex. | |
869 for v in G.nodes(): | |
870 if label.get(inblossom[v]) is None and \ | |
871 bestedge.get(v) is not None: | |
872 d = slack(*bestedge[v]) | |
873 if deltatype == -1 or d < delta: | |
874 delta = d | |
875 deltatype = 2 | |
876 deltaedge = bestedge[v] | |
877 | |
878 # Compute delta3: half the minimum slack on any edge between | |
879 # a pair of S-blossoms. | |
880 for b in blossomparent: | |
881 if (blossomparent[b] is None and label.get(b) == 1 and | |
882 bestedge.get(b) is not None): | |
883 kslack = slack(*bestedge[b]) | |
884 if allinteger: | |
885 assert (kslack % 2) == 0 | |
886 d = kslack // 2 | |
887 else: | |
888 d = kslack / 2.0 | |
889 if deltatype == -1 or d < delta: | |
890 delta = d | |
891 deltatype = 3 | |
892 deltaedge = bestedge[b] | |
893 | |
894 # Compute delta4: minimum z variable of any T-blossom. | |
895 for b in blossomdual: | |
896 if (blossomparent[b] is None and label.get(b) == 2 and | |
897 (deltatype == -1 or blossomdual[b] < delta)): | |
898 delta = blossomdual[b] | |
899 deltatype = 4 | |
900 deltablossom = b | |
901 | |
902 if deltatype == -1: | |
903 # No further improvement possible; max-cardinality optimum | |
904 # reached. Do a final delta update to make the optimum | |
905 # verifyable. | |
906 assert maxcardinality | |
907 deltatype = 1 | |
908 delta = max(0, min(dualvar.values())) | |
909 | |
910 # Update dual variables according to delta. | |
911 for v in gnodes: | |
912 if label.get(inblossom[v]) == 1: | |
913 # S-vertex: 2*u = 2*u - 2*delta | |
914 dualvar[v] -= delta | |
915 elif label.get(inblossom[v]) == 2: | |
916 # T-vertex: 2*u = 2*u + 2*delta | |
917 dualvar[v] += delta | |
918 for b in blossomdual: | |
919 if blossomparent[b] is None: | |
920 if label.get(b) == 1: | |
921 # top-level S-blossom: z = z + 2*delta | |
922 blossomdual[b] += delta | |
923 elif label.get(b) == 2: | |
924 # top-level T-blossom: z = z - 2*delta | |
925 blossomdual[b] -= delta | |
926 | |
927 # Take action at the point where minimum delta occurred. | |
928 if deltatype == 1: | |
929 # No further improvement possible; optimum reached. | |
930 break | |
931 elif deltatype == 2: | |
932 # Use the least-slack edge to continue the search. | |
933 (v, w) = deltaedge | |
934 assert label[inblossom[v]] == 1 | |
935 allowedge[(v, w)] = allowedge[(w, v)] = True | |
936 queue.append(v) | |
937 elif deltatype == 3: | |
938 # Use the least-slack edge to continue the search. | |
939 (v, w) = deltaedge | |
940 allowedge[(v, w)] = allowedge[(w, v)] = True | |
941 assert label[inblossom[v]] == 1 | |
942 queue.append(v) | |
943 elif deltatype == 4: | |
944 # Expand the least-z blossom. | |
945 expandBlossom(deltablossom, False) | |
946 | |
947 # End of a this substage. | |
948 | |
949 # Paranoia check that the matching is symmetric. | |
950 for v in mate: | |
951 assert mate[mate[v]] == v | |
952 | |
953 # Stop when no more augmenting path can be found. | |
954 if not augmented: | |
955 break | |
956 | |
957 # End of a stage; expand all S-blossoms which have zero dual. | |
958 for b in list(blossomdual.keys()): | |
959 if b not in blossomdual: | |
960 continue # already expanded | |
961 if (blossomparent[b] is None and label.get(b) == 1 and | |
962 blossomdual[b] == 0): | |
963 expandBlossom(b, True) | |
964 | |
965 # Verify that we reached the optimum solution (only for integer weights). | |
966 if allinteger: | |
967 verifyOptimum() | |
968 | |
969 return matching_dict_to_set(mate) |