Coverage for /pythoncovmergedfiles/medio/medio/usr/local/lib/python3.9/dist-packages/networkx/algorithms/matching.py: 5%
511 statements
« prev ^ index » next coverage.py v7.3.2, created at 2023-10-20 07:00 +0000
1"""Functions for computing and verifying matchings in a graph."""
2from collections import Counter
3from itertools import combinations, repeat
5import networkx as nx
6from networkx.utils import not_implemented_for
8__all__ = [
9 "is_matching",
10 "is_maximal_matching",
11 "is_perfect_matching",
12 "max_weight_matching",
13 "min_weight_matching",
14 "maximal_matching",
15]
18@not_implemented_for("multigraph")
19@not_implemented_for("directed")
20@nx._dispatch
21def maximal_matching(G):
22 r"""Find a maximal matching in the graph.
24 A matching is a subset of edges in which no node occurs more than once.
25 A maximal matching cannot add more edges and still be a matching.
27 Parameters
28 ----------
29 G : NetworkX graph
30 Undirected graph
32 Returns
33 -------
34 matching : set
35 A maximal matching of the graph.
37 Examples
38 --------
39 >>> G = nx.Graph([(1, 2), (1, 3), (2, 3), (2, 4), (3, 5), (4, 5)])
40 >>> sorted(nx.maximal_matching(G))
41 [(1, 2), (3, 5)]
43 Notes
44 -----
45 The algorithm greedily selects a maximal matching M of the graph G
46 (i.e. no superset of M exists). It runs in $O(|E|)$ time.
47 """
48 matching = set()
49 nodes = set()
50 for edge in G.edges():
51 # If the edge isn't covered, add it to the matching
52 # then remove neighborhood of u and v from consideration.
53 u, v = edge
54 if u not in nodes and v not in nodes and u != v:
55 matching.add(edge)
56 nodes.update(edge)
57 return matching
60def matching_dict_to_set(matching):
61 """Converts matching dict format to matching set format
63 Converts a dictionary representing a matching (as returned by
64 :func:`max_weight_matching`) to a set representing a matching (as
65 returned by :func:`maximal_matching`).
67 In the definition of maximal matching adopted by NetworkX,
68 self-loops are not allowed, so the provided dictionary is expected
69 to never have any mapping from a key to itself. However, the
70 dictionary is expected to have mirrored key/value pairs, for
71 example, key ``u`` with value ``v`` and key ``v`` with value ``u``.
73 """
74 edges = set()
75 for edge in matching.items():
76 u, v = edge
77 if (v, u) in edges or edge in edges:
78 continue
79 if u == v:
80 raise nx.NetworkXError(f"Selfloops cannot appear in matchings {edge}")
81 edges.add(edge)
82 return edges
85@nx._dispatch
86def is_matching(G, matching):
87 """Return True if ``matching`` is a valid matching of ``G``
89 A *matching* in a graph is a set of edges in which no two distinct
90 edges share a common endpoint. Each node is incident to at most one
91 edge in the matching. The edges are said to be independent.
93 Parameters
94 ----------
95 G : NetworkX graph
97 matching : dict or set
98 A dictionary or set representing a matching. If a dictionary, it
99 must have ``matching[u] == v`` and ``matching[v] == u`` for each
100 edge ``(u, v)`` in the matching. If a set, it must have elements
101 of the form ``(u, v)``, where ``(u, v)`` is an edge in the
102 matching.
104 Returns
105 -------
106 bool
107 Whether the given set or dictionary represents a valid matching
108 in the graph.
110 Raises
111 ------
112 NetworkXError
113 If the proposed matching has an edge to a node not in G.
114 Or if the matching is not a collection of 2-tuple edges.
116 Examples
117 --------
118 >>> G = nx.Graph([(1, 2), (1, 3), (2, 3), (2, 4), (3, 5), (4, 5)])
119 >>> nx.is_maximal_matching(G, {1: 3, 2: 4}) # using dict to represent matching
120 True
122 >>> nx.is_matching(G, {(1, 3), (2, 4)}) # using set to represent matching
123 True
125 """
126 if isinstance(matching, dict):
127 matching = matching_dict_to_set(matching)
129 nodes = set()
130 for edge in matching:
131 if len(edge) != 2:
132 raise nx.NetworkXError(f"matching has non-2-tuple edge {edge}")
133 u, v = edge
134 if u not in G or v not in G:
135 raise nx.NetworkXError(f"matching contains edge {edge} with node not in G")
136 if u == v:
137 return False
138 if not G.has_edge(u, v):
139 return False
140 if u in nodes or v in nodes:
141 return False
142 nodes.update(edge)
143 return True
146@nx._dispatch
147def is_maximal_matching(G, matching):
148 """Return True if ``matching`` is a maximal matching of ``G``
150 A *maximal matching* in a graph is a matching in which adding any
151 edge would cause the set to no longer be a valid matching.
153 Parameters
154 ----------
155 G : NetworkX graph
157 matching : dict or set
158 A dictionary or set representing a matching. If a dictionary, it
159 must have ``matching[u] == v`` and ``matching[v] == u`` for each
160 edge ``(u, v)`` in the matching. If a set, it must have elements
161 of the form ``(u, v)``, where ``(u, v)`` is an edge in the
162 matching.
164 Returns
165 -------
166 bool
167 Whether the given set or dictionary represents a valid maximal
168 matching in the graph.
170 Examples
171 --------
172 >>> G = nx.Graph([(1, 2), (1, 3), (2, 3), (3, 4), (3, 5)])
173 >>> nx.is_maximal_matching(G, {(1, 2), (3, 4)})
174 True
176 """
177 if isinstance(matching, dict):
178 matching = matching_dict_to_set(matching)
179 # If the given set is not a matching, then it is not a maximal matching.
180 edges = set()
181 nodes = set()
182 for edge in matching:
183 if len(edge) != 2:
184 raise nx.NetworkXError(f"matching has non-2-tuple edge {edge}")
185 u, v = edge
186 if u not in G or v not in G:
187 raise nx.NetworkXError(f"matching contains edge {edge} with node not in G")
188 if u == v:
189 return False
190 if not G.has_edge(u, v):
191 return False
192 if u in nodes or v in nodes:
193 return False
194 nodes.update(edge)
195 edges.add(edge)
196 edges.add((v, u))
197 # A matching is maximal if adding any new edge from G to it
198 # causes the resulting set to match some node twice.
199 # Be careful to check for adding selfloops
200 for u, v in G.edges:
201 if (u, v) not in edges:
202 # could add edge (u, v) to edges and have a bigger matching
203 if u not in nodes and v not in nodes and u != v:
204 return False
205 return True
208@nx._dispatch
209def is_perfect_matching(G, matching):
210 """Return True if ``matching`` is a perfect matching for ``G``
212 A *perfect matching* in a graph is a matching in which exactly one edge
213 is incident upon each vertex.
215 Parameters
216 ----------
217 G : NetworkX graph
219 matching : dict or set
220 A dictionary or set representing a matching. If a dictionary, it
221 must have ``matching[u] == v`` and ``matching[v] == u`` for each
222 edge ``(u, v)`` in the matching. If a set, it must have elements
223 of the form ``(u, v)``, where ``(u, v)`` is an edge in the
224 matching.
226 Returns
227 -------
228 bool
229 Whether the given set or dictionary represents a valid perfect
230 matching in the graph.
232 Examples
233 --------
234 >>> G = nx.Graph([(1, 2), (1, 3), (2, 3), (2, 4), (3, 5), (4, 5), (4, 6)])
235 >>> my_match = {1: 2, 3: 5, 4: 6}
236 >>> nx.is_perfect_matching(G, my_match)
237 True
239 """
240 if isinstance(matching, dict):
241 matching = matching_dict_to_set(matching)
243 nodes = set()
244 for edge in matching:
245 if len(edge) != 2:
246 raise nx.NetworkXError(f"matching has non-2-tuple edge {edge}")
247 u, v = edge
248 if u not in G or v not in G:
249 raise nx.NetworkXError(f"matching contains edge {edge} with node not in G")
250 if u == v:
251 return False
252 if not G.has_edge(u, v):
253 return False
254 if u in nodes or v in nodes:
255 return False
256 nodes.update(edge)
257 return len(nodes) == len(G)
260@not_implemented_for("multigraph")
261@not_implemented_for("directed")
262@nx._dispatch(edge_attrs="weight")
263def min_weight_matching(G, weight="weight"):
264 """Computing a minimum-weight maximal matching of G.
266 Use the maximum-weight algorithm with edge weights subtracted
267 from the maximum weight of all edges.
269 A matching is a subset of edges in which no node occurs more than once.
270 The weight of a matching is the sum of the weights of its edges.
271 A maximal matching cannot add more edges and still be a matching.
272 The cardinality of a matching is the number of matched edges.
274 This method replaces the edge weights with 1 plus the maximum edge weight
275 minus the original edge weight.
277 new_weight = (max_weight + 1) - edge_weight
279 then runs :func:`max_weight_matching` with the new weights.
280 The max weight matching with these new weights corresponds
281 to the min weight matching using the original weights.
282 Adding 1 to the max edge weight keeps all edge weights positive
283 and as integers if they started as integers.
285 You might worry that adding 1 to each weight would make the algorithm
286 favor matchings with more edges. But we use the parameter
287 `maxcardinality=True` in `max_weight_matching` to ensure that the
288 number of edges in the competing matchings are the same and thus
289 the optimum does not change due to changes in the number of edges.
291 Read the documentation of `max_weight_matching` for more information.
293 Parameters
294 ----------
295 G : NetworkX graph
296 Undirected graph
298 weight: string, optional (default='weight')
299 Edge data key corresponding to the edge weight.
300 If key not found, uses 1 as weight.
302 Returns
303 -------
304 matching : set
305 A minimal weight matching of the graph.
307 See Also
308 --------
309 max_weight_matching
310 """
311 if len(G.edges) == 0:
312 return max_weight_matching(G, maxcardinality=True, weight=weight)
313 G_edges = G.edges(data=weight, default=1)
314 max_weight = 1 + max(w for _, _, w in G_edges)
315 InvG = nx.Graph()
316 edges = ((u, v, max_weight - w) for u, v, w in G_edges)
317 InvG.add_weighted_edges_from(edges, weight=weight)
318 return max_weight_matching(InvG, maxcardinality=True, weight=weight)
321@not_implemented_for("multigraph")
322@not_implemented_for("directed")
323@nx._dispatch(edge_attrs="weight")
324def max_weight_matching(G, maxcardinality=False, weight="weight"):
325 """Compute a maximum-weighted matching of G.
327 A matching is a subset of edges in which no node occurs more than once.
328 The weight of a matching is the sum of the weights of its edges.
329 A maximal matching cannot add more edges and still be a matching.
330 The cardinality of a matching is the number of matched edges.
332 Parameters
333 ----------
334 G : NetworkX graph
335 Undirected graph
337 maxcardinality: bool, optional (default=False)
338 If maxcardinality is True, compute the maximum-cardinality matching
339 with maximum weight among all maximum-cardinality matchings.
341 weight: string, optional (default='weight')
342 Edge data key corresponding to the edge weight.
343 If key not found, uses 1 as weight.
346 Returns
347 -------
348 matching : set
349 A maximal matching of the graph.
351 Examples
352 --------
353 >>> G = nx.Graph()
354 >>> edges = [(1, 2, 6), (1, 3, 2), (2, 3, 1), (2, 4, 7), (3, 5, 9), (4, 5, 3)]
355 >>> G.add_weighted_edges_from(edges)
356 >>> sorted(nx.max_weight_matching(G))
357 [(2, 4), (5, 3)]
359 Notes
360 -----
361 If G has edges with weight attributes the edge data are used as
362 weight values else the weights are assumed to be 1.
364 This function takes time O(number_of_nodes ** 3).
366 If all edge weights are integers, the algorithm uses only integer
367 computations. If floating point weights are used, the algorithm
368 could return a slightly suboptimal matching due to numeric
369 precision errors.
371 This method is based on the "blossom" method for finding augmenting
372 paths and the "primal-dual" method for finding a matching of maximum
373 weight, both methods invented by Jack Edmonds [1]_.
375 Bipartite graphs can also be matched using the functions present in
376 :mod:`networkx.algorithms.bipartite.matching`.
378 References
379 ----------
380 .. [1] "Efficient Algorithms for Finding Maximum Matching in Graphs",
381 Zvi Galil, ACM Computing Surveys, 1986.
382 """
383 #
384 # The algorithm is taken from "Efficient Algorithms for Finding Maximum
385 # Matching in Graphs" by Zvi Galil, ACM Computing Surveys, 1986.
386 # It is based on the "blossom" method for finding augmenting paths and
387 # the "primal-dual" method for finding a matching of maximum weight, both
388 # methods invented by Jack Edmonds.
389 #
390 # A C program for maximum weight matching by Ed Rothberg was used
391 # extensively to validate this new code.
392 #
393 # Many terms used in the code comments are explained in the paper
394 # by Galil. You will probably need the paper to make sense of this code.
395 #
397 class NoNode:
398 """Dummy value which is different from any node."""
400 class Blossom:
401 """Representation of a non-trivial blossom or sub-blossom."""
403 __slots__ = ["childs", "edges", "mybestedges"]
405 # b.childs is an ordered list of b's sub-blossoms, starting with
406 # the base and going round the blossom.
408 # b.edges is the list of b's connecting edges, such that
409 # b.edges[i] = (v, w) where v is a vertex in b.childs[i]
410 # and w is a vertex in b.childs[wrap(i+1)].
412 # If b is a top-level S-blossom,
413 # b.mybestedges is a list of least-slack edges to neighbouring
414 # S-blossoms, or None if no such list has been computed yet.
415 # This is used for efficient computation of delta3.
417 # Generate the blossom's leaf vertices.
418 def leaves(self):
419 stack = [*self.childs]
420 while stack:
421 t = stack.pop()
422 if isinstance(t, Blossom):
423 stack.extend(t.childs)
424 else:
425 yield t
427 # Get a list of vertices.
428 gnodes = list(G)
429 if not gnodes:
430 return set() # don't bother with empty graphs
432 # Find the maximum edge weight.
433 maxweight = 0
434 allinteger = True
435 for i, j, d in G.edges(data=True):
436 wt = d.get(weight, 1)
437 if i != j and wt > maxweight:
438 maxweight = wt
439 allinteger = allinteger and (str(type(wt)).split("'")[1] in ("int", "long"))
441 # If v is a matched vertex, mate[v] is its partner vertex.
442 # If v is a single vertex, v does not occur as a key in mate.
443 # Initially all vertices are single; updated during augmentation.
444 mate = {}
446 # If b is a top-level blossom,
447 # label.get(b) is None if b is unlabeled (free),
448 # 1 if b is an S-blossom,
449 # 2 if b is a T-blossom.
450 # The label of a vertex is found by looking at the label of its top-level
451 # containing blossom.
452 # If v is a vertex inside a T-blossom, label[v] is 2 iff v is reachable
453 # from an S-vertex outside the blossom.
454 # Labels are assigned during a stage and reset after each augmentation.
455 label = {}
457 # If b is a labeled top-level blossom,
458 # labeledge[b] = (v, w) is the edge through which b obtained its label
459 # such that w is a vertex in b, or None if b's base vertex is single.
460 # If w is a vertex inside a T-blossom and label[w] == 2,
461 # labeledge[w] = (v, w) is an edge through which w is reachable from
462 # outside the blossom.
463 labeledge = {}
465 # If v is a vertex, inblossom[v] is the top-level blossom to which v
466 # belongs.
467 # If v is a top-level vertex, inblossom[v] == v since v is itself
468 # a (trivial) top-level blossom.
469 # Initially all vertices are top-level trivial blossoms.
470 inblossom = dict(zip(gnodes, gnodes))
472 # If b is a sub-blossom,
473 # blossomparent[b] is its immediate parent (sub-)blossom.
474 # If b is a top-level blossom, blossomparent[b] is None.
475 blossomparent = dict(zip(gnodes, repeat(None)))
477 # If b is a (sub-)blossom,
478 # blossombase[b] is its base VERTEX (i.e. recursive sub-blossom).
479 blossombase = dict(zip(gnodes, gnodes))
481 # If w is a free vertex (or an unreached vertex inside a T-blossom),
482 # bestedge[w] = (v, w) is the least-slack edge from an S-vertex,
483 # or None if there is no such edge.
484 # If b is a (possibly trivial) top-level S-blossom,
485 # bestedge[b] = (v, w) is the least-slack edge to a different S-blossom
486 # (v inside b), or None if there is no such edge.
487 # This is used for efficient computation of delta2 and delta3.
488 bestedge = {}
490 # If v is a vertex,
491 # dualvar[v] = 2 * u(v) where u(v) is the v's variable in the dual
492 # optimization problem (if all edge weights are integers, multiplication
493 # by two ensures that all values remain integers throughout the algorithm).
494 # Initially, u(v) = maxweight / 2.
495 dualvar = dict(zip(gnodes, repeat(maxweight)))
497 # If b is a non-trivial blossom,
498 # blossomdual[b] = z(b) where z(b) is b's variable in the dual
499 # optimization problem.
500 blossomdual = {}
502 # If (v, w) in allowedge or (w, v) in allowedg, then the edge
503 # (v, w) is known to have zero slack in the optimization problem;
504 # otherwise the edge may or may not have zero slack.
505 allowedge = {}
507 # Queue of newly discovered S-vertices.
508 queue = []
510 # Return 2 * slack of edge (v, w) (does not work inside blossoms).
511 def slack(v, w):
512 return dualvar[v] + dualvar[w] - 2 * G[v][w].get(weight, 1)
514 # Assign label t to the top-level blossom containing vertex w,
515 # coming through an edge from vertex v.
516 def assignLabel(w, t, v):
517 b = inblossom[w]
518 assert label.get(w) is None and label.get(b) is None
519 label[w] = label[b] = t
520 if v is not None:
521 labeledge[w] = labeledge[b] = (v, w)
522 else:
523 labeledge[w] = labeledge[b] = None
524 bestedge[w] = bestedge[b] = None
525 if t == 1:
526 # b became an S-vertex/blossom; add it(s vertices) to the queue.
527 if isinstance(b, Blossom):
528 queue.extend(b.leaves())
529 else:
530 queue.append(b)
531 elif t == 2:
532 # b became a T-vertex/blossom; assign label S to its mate.
533 # (If b is a non-trivial blossom, its base is the only vertex
534 # with an external mate.)
535 base = blossombase[b]
536 assignLabel(mate[base], 1, base)
538 # Trace back from vertices v and w to discover either a new blossom
539 # or an augmenting path. Return the base vertex of the new blossom,
540 # or NoNode if an augmenting path was found.
541 def scanBlossom(v, w):
542 # Trace back from v and w, placing breadcrumbs as we go.
543 path = []
544 base = NoNode
545 while v is not NoNode:
546 # Look for a breadcrumb in v's blossom or put a new breadcrumb.
547 b = inblossom[v]
548 if label[b] & 4:
549 base = blossombase[b]
550 break
551 assert label[b] == 1
552 path.append(b)
553 label[b] = 5
554 # Trace one step back.
555 if labeledge[b] is None:
556 # The base of blossom b is single; stop tracing this path.
557 assert blossombase[b] not in mate
558 v = NoNode
559 else:
560 assert labeledge[b][0] == mate[blossombase[b]]
561 v = labeledge[b][0]
562 b = inblossom[v]
563 assert label[b] == 2
564 # b is a T-blossom; trace one more step back.
565 v = labeledge[b][0]
566 # Swap v and w so that we alternate between both paths.
567 if w is not NoNode:
568 v, w = w, v
569 # Remove breadcrumbs.
570 for b in path:
571 label[b] = 1
572 # Return base vertex, if we found one.
573 return base
575 # Construct a new blossom with given base, through S-vertices v and w.
576 # Label the new blossom as S; set its dual variable to zero;
577 # relabel its T-vertices to S and add them to the queue.
578 def addBlossom(base, v, w):
579 bb = inblossom[base]
580 bv = inblossom[v]
581 bw = inblossom[w]
582 # Create blossom.
583 b = Blossom()
584 blossombase[b] = base
585 blossomparent[b] = None
586 blossomparent[bb] = b
587 # Make list of sub-blossoms and their interconnecting edge endpoints.
588 b.childs = path = []
589 b.edges = edgs = [(v, w)]
590 # Trace back from v to base.
591 while bv != bb:
592 # Add bv to the new blossom.
593 blossomparent[bv] = b
594 path.append(bv)
595 edgs.append(labeledge[bv])
596 assert label[bv] == 2 or (
597 label[bv] == 1 and labeledge[bv][0] == mate[blossombase[bv]]
598 )
599 # Trace one step back.
600 v = labeledge[bv][0]
601 bv = inblossom[v]
602 # Add base sub-blossom; reverse lists.
603 path.append(bb)
604 path.reverse()
605 edgs.reverse()
606 # Trace back from w to base.
607 while bw != bb:
608 # Add bw to the new blossom.
609 blossomparent[bw] = b
610 path.append(bw)
611 edgs.append((labeledge[bw][1], labeledge[bw][0]))
612 assert label[bw] == 2 or (
613 label[bw] == 1 and labeledge[bw][0] == mate[blossombase[bw]]
614 )
615 # Trace one step back.
616 w = labeledge[bw][0]
617 bw = inblossom[w]
618 # Set label to S.
619 assert label[bb] == 1
620 label[b] = 1
621 labeledge[b] = labeledge[bb]
622 # Set dual variable to zero.
623 blossomdual[b] = 0
624 # Relabel vertices.
625 for v in b.leaves():
626 if label[inblossom[v]] == 2:
627 # This T-vertex now turns into an S-vertex because it becomes
628 # part of an S-blossom; add it to the queue.
629 queue.append(v)
630 inblossom[v] = b
631 # Compute b.mybestedges.
632 bestedgeto = {}
633 for bv in path:
634 if isinstance(bv, Blossom):
635 if bv.mybestedges is not None:
636 # Walk this subblossom's least-slack edges.
637 nblist = bv.mybestedges
638 # The sub-blossom won't need this data again.
639 bv.mybestedges = None
640 else:
641 # This subblossom does not have a list of least-slack
642 # edges; get the information from the vertices.
643 nblist = [
644 (v, w) for v in bv.leaves() for w in G.neighbors(v) if v != w
645 ]
646 else:
647 nblist = [(bv, w) for w in G.neighbors(bv) if bv != w]
648 for k in nblist:
649 (i, j) = k
650 if inblossom[j] == b:
651 i, j = j, i
652 bj = inblossom[j]
653 if (
654 bj != b
655 and label.get(bj) == 1
656 and ((bj not in bestedgeto) or slack(i, j) < slack(*bestedgeto[bj]))
657 ):
658 bestedgeto[bj] = k
659 # Forget about least-slack edge of the subblossom.
660 bestedge[bv] = None
661 b.mybestedges = list(bestedgeto.values())
662 # Select bestedge[b].
663 mybestedge = None
664 bestedge[b] = None
665 for k in b.mybestedges:
666 kslack = slack(*k)
667 if mybestedge is None or kslack < mybestslack:
668 mybestedge = k
669 mybestslack = kslack
670 bestedge[b] = mybestedge
672 # Expand the given top-level blossom.
673 def expandBlossom(b, endstage):
674 # This is an obnoxiously complicated recursive function for the sake of
675 # a stack-transformation. So, we hack around the complexity by using
676 # a trampoline pattern. By yielding the arguments to each recursive
677 # call, we keep the actual callstack flat.
679 def _recurse(b, endstage):
680 # Convert sub-blossoms into top-level blossoms.
681 for s in b.childs:
682 blossomparent[s] = None
683 if isinstance(s, Blossom):
684 if endstage and blossomdual[s] == 0:
685 # Recursively expand this sub-blossom.
686 yield s
687 else:
688 for v in s.leaves():
689 inblossom[v] = s
690 else:
691 inblossom[s] = s
692 # If we expand a T-blossom during a stage, its sub-blossoms must be
693 # relabeled.
694 if (not endstage) and label.get(b) == 2:
695 # Start at the sub-blossom through which the expanding
696 # blossom obtained its label, and relabel sub-blossoms untili
697 # we reach the base.
698 # Figure out through which sub-blossom the expanding blossom
699 # obtained its label initially.
700 entrychild = inblossom[labeledge[b][1]]
701 # Decide in which direction we will go round the blossom.
702 j = b.childs.index(entrychild)
703 if j & 1:
704 # Start index is odd; go forward and wrap.
705 j -= len(b.childs)
706 jstep = 1
707 else:
708 # Start index is even; go backward.
709 jstep = -1
710 # Move along the blossom until we get to the base.
711 v, w = labeledge[b]
712 while j != 0:
713 # Relabel the T-sub-blossom.
714 if jstep == 1:
715 p, q = b.edges[j]
716 else:
717 q, p = b.edges[j - 1]
718 label[w] = None
719 label[q] = None
720 assignLabel(w, 2, v)
721 # Step to the next S-sub-blossom and note its forward edge.
722 allowedge[(p, q)] = allowedge[(q, p)] = True
723 j += jstep
724 if jstep == 1:
725 v, w = b.edges[j]
726 else:
727 w, v = b.edges[j - 1]
728 # Step to the next T-sub-blossom.
729 allowedge[(v, w)] = allowedge[(w, v)] = True
730 j += jstep
731 # Relabel the base T-sub-blossom WITHOUT stepping through to
732 # its mate (so don't call assignLabel).
733 bw = b.childs[j]
734 label[w] = label[bw] = 2
735 labeledge[w] = labeledge[bw] = (v, w)
736 bestedge[bw] = None
737 # Continue along the blossom until we get back to entrychild.
738 j += jstep
739 while b.childs[j] != entrychild:
740 # Examine the vertices of the sub-blossom to see whether
741 # it is reachable from a neighbouring S-vertex outside the
742 # expanding blossom.
743 bv = b.childs[j]
744 if label.get(bv) == 1:
745 # This sub-blossom just got label S through one of its
746 # neighbours; leave it be.
747 j += jstep
748 continue
749 if isinstance(bv, Blossom):
750 for v in bv.leaves():
751 if label.get(v):
752 break
753 else:
754 v = bv
755 # If the sub-blossom contains a reachable vertex, assign
756 # label T to the sub-blossom.
757 if label.get(v):
758 assert label[v] == 2
759 assert inblossom[v] == bv
760 label[v] = None
761 label[mate[blossombase[bv]]] = None
762 assignLabel(v, 2, labeledge[v][0])
763 j += jstep
764 # Remove the expanded blossom entirely.
765 label.pop(b, None)
766 labeledge.pop(b, None)
767 bestedge.pop(b, None)
768 del blossomparent[b]
769 del blossombase[b]
770 del blossomdual[b]
772 # Now, we apply the trampoline pattern. We simulate a recursive
773 # callstack by maintaining a stack of generators, each yielding a
774 # sequence of function arguments. We grow the stack by appending a call
775 # to _recurse on each argument tuple, and shrink the stack whenever a
776 # generator is exhausted.
777 stack = [_recurse(b, endstage)]
778 while stack:
779 top = stack[-1]
780 for s in top:
781 stack.append(_recurse(s, endstage))
782 break
783 else:
784 stack.pop()
786 # Swap matched/unmatched edges over an alternating path through blossom b
787 # between vertex v and the base vertex. Keep blossom bookkeeping
788 # consistent.
789 def augmentBlossom(b, v):
790 # This is an obnoxiously complicated recursive function for the sake of
791 # a stack-transformation. So, we hack around the complexity by using
792 # a trampoline pattern. By yielding the arguments to each recursive
793 # call, we keep the actual callstack flat.
795 def _recurse(b, v):
796 # Bubble up through the blossom tree from vertex v to an immediate
797 # sub-blossom of b.
798 t = v
799 while blossomparent[t] != b:
800 t = blossomparent[t]
801 # Recursively deal with the first sub-blossom.
802 if isinstance(t, Blossom):
803 yield (t, v)
804 # Decide in which direction we will go round the blossom.
805 i = j = b.childs.index(t)
806 if i & 1:
807 # Start index is odd; go forward and wrap.
808 j -= len(b.childs)
809 jstep = 1
810 else:
811 # Start index is even; go backward.
812 jstep = -1
813 # Move along the blossom until we get to the base.
814 while j != 0:
815 # Step to the next sub-blossom and augment it recursively.
816 j += jstep
817 t = b.childs[j]
818 if jstep == 1:
819 w, x = b.edges[j]
820 else:
821 x, w = b.edges[j - 1]
822 if isinstance(t, Blossom):
823 yield (t, w)
824 # Step to the next sub-blossom and augment it recursively.
825 j += jstep
826 t = b.childs[j]
827 if isinstance(t, Blossom):
828 yield (t, x)
829 # Match the edge connecting those sub-blossoms.
830 mate[w] = x
831 mate[x] = w
832 # Rotate the list of sub-blossoms to put the new base at the front.
833 b.childs = b.childs[i:] + b.childs[:i]
834 b.edges = b.edges[i:] + b.edges[:i]
835 blossombase[b] = blossombase[b.childs[0]]
836 assert blossombase[b] == v
838 # Now, we apply the trampoline pattern. We simulate a recursive
839 # callstack by maintaining a stack of generators, each yielding a
840 # sequence of function arguments. We grow the stack by appending a call
841 # to _recurse on each argument tuple, and shrink the stack whenever a
842 # generator is exhausted.
843 stack = [_recurse(b, v)]
844 while stack:
845 top = stack[-1]
846 for args in top:
847 stack.append(_recurse(*args))
848 break
849 else:
850 stack.pop()
852 # Swap matched/unmatched edges over an alternating path between two
853 # single vertices. The augmenting path runs through S-vertices v and w.
854 def augmentMatching(v, w):
855 for s, j in ((v, w), (w, v)):
856 # Match vertex s to vertex j. Then trace back from s
857 # until we find a single vertex, swapping matched and unmatched
858 # edges as we go.
859 while 1:
860 bs = inblossom[s]
861 assert label[bs] == 1
862 assert (labeledge[bs] is None and blossombase[bs] not in mate) or (
863 labeledge[bs][0] == mate[blossombase[bs]]
864 )
865 # Augment through the S-blossom from s to base.
866 if isinstance(bs, Blossom):
867 augmentBlossom(bs, s)
868 # Update mate[s]
869 mate[s] = j
870 # Trace one step back.
871 if labeledge[bs] is None:
872 # Reached single vertex; stop.
873 break
874 t = labeledge[bs][0]
875 bt = inblossom[t]
876 assert label[bt] == 2
877 # Trace one more step back.
878 s, j = labeledge[bt]
879 # Augment through the T-blossom from j to base.
880 assert blossombase[bt] == t
881 if isinstance(bt, Blossom):
882 augmentBlossom(bt, j)
883 # Update mate[j]
884 mate[j] = s
886 # Verify that the optimum solution has been reached.
887 def verifyOptimum():
888 if maxcardinality:
889 # Vertices may have negative dual;
890 # find a constant non-negative number to add to all vertex duals.
891 vdualoffset = max(0, -min(dualvar.values()))
892 else:
893 vdualoffset = 0
894 # 0. all dual variables are non-negative
895 assert min(dualvar.values()) + vdualoffset >= 0
896 assert len(blossomdual) == 0 or min(blossomdual.values()) >= 0
897 # 0. all edges have non-negative slack and
898 # 1. all matched edges have zero slack;
899 for i, j, d in G.edges(data=True):
900 wt = d.get(weight, 1)
901 if i == j:
902 continue # ignore self-loops
903 s = dualvar[i] + dualvar[j] - 2 * wt
904 iblossoms = [i]
905 jblossoms = [j]
906 while blossomparent[iblossoms[-1]] is not None:
907 iblossoms.append(blossomparent[iblossoms[-1]])
908 while blossomparent[jblossoms[-1]] is not None:
909 jblossoms.append(blossomparent[jblossoms[-1]])
910 iblossoms.reverse()
911 jblossoms.reverse()
912 for bi, bj in zip(iblossoms, jblossoms):
913 if bi != bj:
914 break
915 s += 2 * blossomdual[bi]
916 assert s >= 0
917 if mate.get(i) == j or mate.get(j) == i:
918 assert mate[i] == j and mate[j] == i
919 assert s == 0
920 # 2. all single vertices have zero dual value;
921 for v in gnodes:
922 assert (v in mate) or dualvar[v] + vdualoffset == 0
923 # 3. all blossoms with positive dual value are full.
924 for b in blossomdual:
925 if blossomdual[b] > 0:
926 assert len(b.edges) % 2 == 1
927 for i, j in b.edges[1::2]:
928 assert mate[i] == j and mate[j] == i
929 # Ok.
931 # Main loop: continue until no further improvement is possible.
932 while 1:
933 # Each iteration of this loop is a "stage".
934 # A stage finds an augmenting path and uses that to improve
935 # the matching.
937 # Remove labels from top-level blossoms/vertices.
938 label.clear()
939 labeledge.clear()
941 # Forget all about least-slack edges.
942 bestedge.clear()
943 for b in blossomdual:
944 b.mybestedges = None
946 # Loss of labeling means that we can not be sure that currently
947 # allowable edges remain allowable throughout this stage.
948 allowedge.clear()
950 # Make queue empty.
951 queue[:] = []
953 # Label single blossoms/vertices with S and put them in the queue.
954 for v in gnodes:
955 if (v not in mate) and label.get(inblossom[v]) is None:
956 assignLabel(v, 1, None)
958 # Loop until we succeed in augmenting the matching.
959 augmented = 0
960 while 1:
961 # Each iteration of this loop is a "substage".
962 # A substage tries to find an augmenting path;
963 # if found, the path is used to improve the matching and
964 # the stage ends. If there is no augmenting path, the
965 # primal-dual method is used to pump some slack out of
966 # the dual variables.
968 # Continue labeling until all vertices which are reachable
969 # through an alternating path have got a label.
970 while queue and not augmented:
971 # Take an S vertex from the queue.
972 v = queue.pop()
973 assert label[inblossom[v]] == 1
975 # Scan its neighbours:
976 for w in G.neighbors(v):
977 if w == v:
978 continue # ignore self-loops
979 # w is a neighbour to v
980 bv = inblossom[v]
981 bw = inblossom[w]
982 if bv == bw:
983 # this edge is internal to a blossom; ignore it
984 continue
985 if (v, w) not in allowedge:
986 kslack = slack(v, w)
987 if kslack <= 0:
988 # edge k has zero slack => it is allowable
989 allowedge[(v, w)] = allowedge[(w, v)] = True
990 if (v, w) in allowedge:
991 if label.get(bw) is None:
992 # (C1) w is a free vertex;
993 # label w with T and label its mate with S (R12).
994 assignLabel(w, 2, v)
995 elif label.get(bw) == 1:
996 # (C2) w is an S-vertex (not in the same blossom);
997 # follow back-links to discover either an
998 # augmenting path or a new blossom.
999 base = scanBlossom(v, w)
1000 if base is not NoNode:
1001 # Found a new blossom; add it to the blossom
1002 # bookkeeping and turn it into an S-blossom.
1003 addBlossom(base, v, w)
1004 else:
1005 # Found an augmenting path; augment the
1006 # matching and end this stage.
1007 augmentMatching(v, w)
1008 augmented = 1
1009 break
1010 elif label.get(w) is None:
1011 # w is inside a T-blossom, but w itself has not
1012 # yet been reached from outside the blossom;
1013 # mark it as reached (we need this to relabel
1014 # during T-blossom expansion).
1015 assert label[bw] == 2
1016 label[w] = 2
1017 labeledge[w] = (v, w)
1018 elif label.get(bw) == 1:
1019 # keep track of the least-slack non-allowable edge to
1020 # a different S-blossom.
1021 if bestedge.get(bv) is None or kslack < slack(*bestedge[bv]):
1022 bestedge[bv] = (v, w)
1023 elif label.get(w) is None:
1024 # w is a free vertex (or an unreached vertex inside
1025 # a T-blossom) but we can not reach it yet;
1026 # keep track of the least-slack edge that reaches w.
1027 if bestedge.get(w) is None or kslack < slack(*bestedge[w]):
1028 bestedge[w] = (v, w)
1030 if augmented:
1031 break
1033 # There is no augmenting path under these constraints;
1034 # compute delta and reduce slack in the optimization problem.
1035 # (Note that our vertex dual variables, edge slacks and delta's
1036 # are pre-multiplied by two.)
1037 deltatype = -1
1038 delta = deltaedge = deltablossom = None
1040 # Compute delta1: the minimum value of any vertex dual.
1041 if not maxcardinality:
1042 deltatype = 1
1043 delta = min(dualvar.values())
1045 # Compute delta2: the minimum slack on any edge between
1046 # an S-vertex and a free vertex.
1047 for v in G.nodes():
1048 if label.get(inblossom[v]) is None and bestedge.get(v) is not None:
1049 d = slack(*bestedge[v])
1050 if deltatype == -1 or d < delta:
1051 delta = d
1052 deltatype = 2
1053 deltaedge = bestedge[v]
1055 # Compute delta3: half the minimum slack on any edge between
1056 # a pair of S-blossoms.
1057 for b in blossomparent:
1058 if (
1059 blossomparent[b] is None
1060 and label.get(b) == 1
1061 and bestedge.get(b) is not None
1062 ):
1063 kslack = slack(*bestedge[b])
1064 if allinteger:
1065 assert (kslack % 2) == 0
1066 d = kslack // 2
1067 else:
1068 d = kslack / 2.0
1069 if deltatype == -1 or d < delta:
1070 delta = d
1071 deltatype = 3
1072 deltaedge = bestedge[b]
1074 # Compute delta4: minimum z variable of any T-blossom.
1075 for b in blossomdual:
1076 if (
1077 blossomparent[b] is None
1078 and label.get(b) == 2
1079 and (deltatype == -1 or blossomdual[b] < delta)
1080 ):
1081 delta = blossomdual[b]
1082 deltatype = 4
1083 deltablossom = b
1085 if deltatype == -1:
1086 # No further improvement possible; max-cardinality optimum
1087 # reached. Do a final delta update to make the optimum
1088 # verifiable.
1089 assert maxcardinality
1090 deltatype = 1
1091 delta = max(0, min(dualvar.values()))
1093 # Update dual variables according to delta.
1094 for v in gnodes:
1095 if label.get(inblossom[v]) == 1:
1096 # S-vertex: 2*u = 2*u - 2*delta
1097 dualvar[v] -= delta
1098 elif label.get(inblossom[v]) == 2:
1099 # T-vertex: 2*u = 2*u + 2*delta
1100 dualvar[v] += delta
1101 for b in blossomdual:
1102 if blossomparent[b] is None:
1103 if label.get(b) == 1:
1104 # top-level S-blossom: z = z + 2*delta
1105 blossomdual[b] += delta
1106 elif label.get(b) == 2:
1107 # top-level T-blossom: z = z - 2*delta
1108 blossomdual[b] -= delta
1110 # Take action at the point where minimum delta occurred.
1111 if deltatype == 1:
1112 # No further improvement possible; optimum reached.
1113 break
1114 elif deltatype == 2:
1115 # Use the least-slack edge to continue the search.
1116 (v, w) = deltaedge
1117 assert label[inblossom[v]] == 1
1118 allowedge[(v, w)] = allowedge[(w, v)] = True
1119 queue.append(v)
1120 elif deltatype == 3:
1121 # Use the least-slack edge to continue the search.
1122 (v, w) = deltaedge
1123 allowedge[(v, w)] = allowedge[(w, v)] = True
1124 assert label[inblossom[v]] == 1
1125 queue.append(v)
1126 elif deltatype == 4:
1127 # Expand the least-z blossom.
1128 expandBlossom(deltablossom, False)
1130 # End of a this substage.
1132 # Paranoia check that the matching is symmetric.
1133 for v in mate:
1134 assert mate[mate[v]] == v
1136 # Stop when no more augmenting path can be found.
1137 if not augmented:
1138 break
1140 # End of a stage; expand all S-blossoms which have zero dual.
1141 for b in list(blossomdual.keys()):
1142 if b not in blossomdual:
1143 continue # already expanded
1144 if blossomparent[b] is None and label.get(b) == 1 and blossomdual[b] == 0:
1145 expandBlossom(b, True)
1147 # Verify that we reached the optimum solution (only for integer weights).
1148 if allinteger:
1149 verifyOptimum()
1151 return matching_dict_to_set(mate)