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)