Shortcuts on this page
r m x toggle line displays
j k next/prev highlighted chunk
0 (zero) top of page
1 (one) first highlighted chunk
1"""Functions for computing and verifying matchings in a graph."""
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._dispatchable
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._dispatchable
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._dispatchable
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._dispatchable
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._dispatchable(edge_attrs="weight")
263def min_weight_matching(G, weight="weight"):
264 """Compute a minimum-weight maximum-cardinality matching of `G`.
266 The minimum-weight maximum-cardinality matching is the matching
267 that has the minimum weight among all maximum-cardinality matchings.
269 Use the maximum-weight algorithm with edge weights subtracted
270 from the maximum weight of all edges.
272 A matching is a subset of edges in which no node occurs more than once.
273 The weight of a matching is the sum of the weights of its edges.
274 A maximal matching cannot add more edges and still be a matching.
275 The cardinality of a matching is the number of matched edges.
277 This method replaces the edge weights with 1 plus the maximum edge weight
278 minus the original edge weight.
280 new_weight = (max_weight + 1) - edge_weight
282 then runs :func:`max_weight_matching` with the new weights.
283 The max weight matching with these new weights corresponds
284 to the min weight matching using the original weights.
285 Adding 1 to the max edge weight keeps all edge weights positive
286 and as integers if they started as integers.
288 Read the documentation of `max_weight_matching` for more information.
290 Parameters
291 ----------
292 G : NetworkX graph
293 Undirected graph
295 weight: string, optional (default='weight')
296 Edge data key corresponding to the edge weight.
297 If key not found, uses 1 as weight.
299 Returns
300 -------
301 matching : set
302 A minimal weight matching of the graph.
304 See Also
305 --------
306 max_weight_matching
307 """
308 if len(G.edges) == 0:
309 return max_weight_matching(G, maxcardinality=True, weight=weight)
310 G_edges = G.edges(data=weight, default=1)
311 max_weight = 1 + max(w for _, _, w in G_edges)
312 InvG = nx.Graph()
313 edges = ((u, v, max_weight - w) for u, v, w in G_edges)
314 InvG.add_weighted_edges_from(edges, weight=weight)
315 return max_weight_matching(InvG, maxcardinality=True, weight=weight)
318@not_implemented_for("multigraph")
319@not_implemented_for("directed")
320@nx._dispatchable(edge_attrs="weight")
321def max_weight_matching(G, maxcardinality=False, weight="weight"):
322 """Compute a maximum-weighted matching of G.
324 A matching is a subset of edges in which no node occurs more than once.
325 The weight of a matching is the sum of the weights of its edges.
326 A maximal matching cannot add more edges and still be a matching.
327 The cardinality of a matching is the number of matched edges.
329 Parameters
330 ----------
331 G : NetworkX graph
332 Undirected graph
334 maxcardinality: bool, optional (default=False)
335 If maxcardinality is True, compute the maximum-cardinality matching
336 with maximum weight among all maximum-cardinality matchings.
338 weight: string, optional (default='weight')
339 Edge data key corresponding to the edge weight.
340 If key not found, uses 1 as weight.
343 Returns
344 -------
345 matching : set
346 A maximal matching of the graph.
348 Examples
349 --------
350 >>> G = nx.Graph()
351 >>> edges = [(1, 2, 6), (1, 3, 2), (2, 3, 1), (2, 4, 7), (3, 5, 9), (4, 5, 3)]
352 >>> G.add_weighted_edges_from(edges)
353 >>> sorted(nx.max_weight_matching(G))
354 [(2, 4), (5, 3)]
356 Notes
357 -----
358 If G has edges with weight attributes the edge data are used as
359 weight values else the weights are assumed to be 1.
361 This function takes time O(number_of_nodes ** 3).
363 If all edge weights are integers, the algorithm uses only integer
364 computations. If floating point weights are used, the algorithm
365 could return a slightly suboptimal matching due to numeric
366 precision errors.
368 This method is based on the "blossom" method for finding augmenting
369 paths and the "primal-dual" method for finding a matching of maximum
370 weight, both methods invented by Jack Edmonds [1]_.
372 Bipartite graphs can also be matched using the functions present in
373 :mod:`networkx.algorithms.bipartite.matching`.
375 References
376 ----------
377 .. [1] "Efficient Algorithms for Finding Maximum Matching in Graphs",
378 Zvi Galil, ACM Computing Surveys, 1986.
379 """
380 #
381 # The algorithm is taken from "Efficient Algorithms for Finding Maximum
382 # Matching in Graphs" by Zvi Galil, ACM Computing Surveys, 1986.
383 # It is based on the "blossom" method for finding augmenting paths and
384 # the "primal-dual" method for finding a matching of maximum weight, both
385 # methods invented by Jack Edmonds.
386 #
387 # A C program for maximum weight matching by Ed Rothberg was used
388 # extensively to validate this new code.
389 #
390 # Many terms used in the code comments are explained in the paper
391 # by Galil. You will probably need the paper to make sense of this code.
392 #
394 class NoNode:
395 """Dummy value which is different from any node."""
397 class Blossom:
398 """Representation of a non-trivial blossom or sub-blossom."""
400 __slots__ = ["childs", "edges", "mybestedges"]
402 # b.childs is an ordered list of b's sub-blossoms, starting with
403 # the base and going round the blossom.
405 # b.edges is the list of b's connecting edges, such that
406 # b.edges[i] = (v, w) where v is a vertex in b.childs[i]
407 # and w is a vertex in b.childs[wrap(i+1)].
409 # If b is a top-level S-blossom,
410 # b.mybestedges is a list of least-slack edges to neighboring
411 # S-blossoms, or None if no such list has been computed yet.
412 # This is used for efficient computation of delta3.
414 # Generate the blossom's leaf vertices.
415 def leaves(self):
416 stack = [*self.childs]
417 while stack:
418 t = stack.pop()
419 if isinstance(t, Blossom):
420 stack.extend(t.childs)
421 else:
422 yield t
424 # Get a list of vertices.
425 gnodes = list(G)
426 if not gnodes:
427 return set() # don't bother with empty graphs
429 # Find the maximum edge weight.
430 maxweight = 0
431 allinteger = True
432 for i, j, d in G.edges(data=True):
433 wt = d.get(weight, 1)
434 if i != j and wt > maxweight:
435 maxweight = wt
436 allinteger = allinteger and (str(type(wt)).split("'")[1] in ("int", "long"))
438 # If v is a matched vertex, mate[v] is its partner vertex.
439 # If v is a single vertex, v does not occur as a key in mate.
440 # Initially all vertices are single; updated during augmentation.
441 mate = {}
443 # If b is a top-level blossom,
444 # label.get(b) is None if b is unlabeled (free),
445 # 1 if b is an S-blossom,
446 # 2 if b is a T-blossom.
447 # The label of a vertex is found by looking at the label of its top-level
448 # containing blossom.
449 # If v is a vertex inside a T-blossom, label[v] is 2 iff v is reachable
450 # from an S-vertex outside the blossom.
451 # Labels are assigned during a stage and reset after each augmentation.
452 label = {}
454 # If b is a labeled top-level blossom,
455 # labeledge[b] = (v, w) is the edge through which b obtained its label
456 # such that w is a vertex in b, or None if b's base vertex is single.
457 # If w is a vertex inside a T-blossom and label[w] == 2,
458 # labeledge[w] = (v, w) is an edge through which w is reachable from
459 # outside the blossom.
460 labeledge = {}
462 # If v is a vertex, inblossom[v] is the top-level blossom to which v
463 # belongs.
464 # If v is a top-level vertex, inblossom[v] == v since v is itself
465 # a (trivial) top-level blossom.
466 # Initially all vertices are top-level trivial blossoms.
467 inblossom = dict(zip(gnodes, gnodes))
469 # If b is a sub-blossom,
470 # blossomparent[b] is its immediate parent (sub-)blossom.
471 # If b is a top-level blossom, blossomparent[b] is None.
472 blossomparent = dict(zip(gnodes, repeat(None)))
474 # If b is a (sub-)blossom,
475 # blossombase[b] is its base VERTEX (i.e. recursive sub-blossom).
476 blossombase = dict(zip(gnodes, gnodes))
478 # If w is a free vertex (or an unreached vertex inside a T-blossom),
479 # bestedge[w] = (v, w) is the least-slack edge from an S-vertex,
480 # or None if there is no such edge.
481 # If b is a (possibly trivial) top-level S-blossom,
482 # bestedge[b] = (v, w) is the least-slack edge to a different S-blossom
483 # (v inside b), or None if there is no such edge.
484 # This is used for efficient computation of delta2 and delta3.
485 bestedge = {}
487 # If v is a vertex,
488 # dualvar[v] = 2 * u(v) where u(v) is the v's variable in the dual
489 # optimization problem (if all edge weights are integers, multiplication
490 # by two ensures that all values remain integers throughout the algorithm).
491 # Initially, u(v) = maxweight / 2.
492 dualvar = dict(zip(gnodes, repeat(maxweight)))
494 # If b is a non-trivial blossom,
495 # blossomdual[b] = z(b) where z(b) is b's variable in the dual
496 # optimization problem.
497 blossomdual = {}
499 # If (v, w) in allowedge or (w, v) in allowedg, then the edge
500 # (v, w) is known to have zero slack in the optimization problem;
501 # otherwise the edge may or may not have zero slack.
502 allowedge = {}
504 # Queue of newly discovered S-vertices.
505 queue = []
507 # Return 2 * slack of edge (v, w) (does not work inside blossoms).
508 def slack(v, w):
509 return dualvar[v] + dualvar[w] - 2 * G[v][w].get(weight, 1)
511 # Assign label t to the top-level blossom containing vertex w,
512 # coming through an edge from vertex v.
513 def assignLabel(w, t, v):
514 b = inblossom[w]
515 assert label.get(w) is None and label.get(b) is None
516 label[w] = label[b] = t
517 if v is not None:
518 labeledge[w] = labeledge[b] = (v, w)
519 else:
520 labeledge[w] = labeledge[b] = None
521 bestedge[w] = bestedge[b] = None
522 if t == 1:
523 # b became an S-vertex/blossom; add it(s vertices) to the queue.
524 if isinstance(b, Blossom):
525 queue.extend(b.leaves())
526 else:
527 queue.append(b)
528 elif t == 2:
529 # b became a T-vertex/blossom; assign label S to its mate.
530 # (If b is a non-trivial blossom, its base is the only vertex
531 # with an external mate.)
532 base = blossombase[b]
533 assignLabel(mate[base], 1, base)
535 # Trace back from vertices v and w to discover either a new blossom
536 # or an augmenting path. Return the base vertex of the new blossom,
537 # or NoNode if an augmenting path was found.
538 def scanBlossom(v, w):
539 # Trace back from v and w, placing breadcrumbs as we go.
540 path = []
541 base = NoNode
542 while v is not NoNode:
543 # Look for a breadcrumb in v's blossom or put a new breadcrumb.
544 b = inblossom[v]
545 if label[b] & 4:
546 base = blossombase[b]
547 break
548 assert label[b] == 1
549 path.append(b)
550 label[b] = 5
551 # Trace one step back.
552 if labeledge[b] is None:
553 # The base of blossom b is single; stop tracing this path.
554 assert blossombase[b] not in mate
555 v = NoNode
556 else:
557 assert labeledge[b][0] == mate[blossombase[b]]
558 v = labeledge[b][0]
559 b = inblossom[v]
560 assert label[b] == 2
561 # b is a T-blossom; trace one more step back.
562 v = labeledge[b][0]
563 # Swap v and w so that we alternate between both paths.
564 if w is not NoNode:
565 v, w = w, v
566 # Remove breadcrumbs.
567 for b in path:
568 label[b] = 1
569 # Return base vertex, if we found one.
570 return base
572 # Construct a new blossom with given base, through S-vertices v and w.
573 # Label the new blossom as S; set its dual variable to zero;
574 # relabel its T-vertices to S and add them to the queue.
575 def addBlossom(base, v, w):
576 bb = inblossom[base]
577 bv = inblossom[v]
578 bw = inblossom[w]
579 # Create blossom.
580 b = Blossom()
581 blossombase[b] = base
582 blossomparent[b] = None
583 blossomparent[bb] = b
584 # Make list of sub-blossoms and their interconnecting edge endpoints.
585 b.childs = path = []
586 b.edges = edgs = [(v, w)]
587 # Trace back from v to base.
588 while bv != bb:
589 # Add bv to the new blossom.
590 blossomparent[bv] = b
591 path.append(bv)
592 edgs.append(labeledge[bv])
593 assert label[bv] == 2 or (
594 label[bv] == 1 and labeledge[bv][0] == mate[blossombase[bv]]
595 )
596 # Trace one step back.
597 v = labeledge[bv][0]
598 bv = inblossom[v]
599 # Add base sub-blossom; reverse lists.
600 path.append(bb)
601 path.reverse()
602 edgs.reverse()
603 # Trace back from w to base.
604 while bw != bb:
605 # Add bw to the new blossom.
606 blossomparent[bw] = b
607 path.append(bw)
608 edgs.append((labeledge[bw][1], labeledge[bw][0]))
609 assert label[bw] == 2 or (
610 label[bw] == 1 and labeledge[bw][0] == mate[blossombase[bw]]
611 )
612 # Trace one step back.
613 w = labeledge[bw][0]
614 bw = inblossom[w]
615 # Set label to S.
616 assert label[bb] == 1
617 label[b] = 1
618 labeledge[b] = labeledge[bb]
619 # Set dual variable to zero.
620 blossomdual[b] = 0
621 # Relabel vertices.
622 for v in b.leaves():
623 if label[inblossom[v]] == 2:
624 # This T-vertex now turns into an S-vertex because it becomes
625 # part of an S-blossom; add it to the queue.
626 queue.append(v)
627 inblossom[v] = b
628 # Compute b.mybestedges.
629 bestedgeto = {}
630 for bv in path:
631 if isinstance(bv, Blossom):
632 if bv.mybestedges is not None:
633 # Walk this subblossom's least-slack edges.
634 nblist = bv.mybestedges
635 # The sub-blossom won't need this data again.
636 bv.mybestedges = None
637 else:
638 # This subblossom does not have a list of least-slack
639 # edges; get the information from the vertices.
640 nblist = [
641 (v, w) for v in bv.leaves() for w in G.neighbors(v) if v != w
642 ]
643 else:
644 nblist = [(bv, w) for w in G.neighbors(bv) if bv != w]
645 for k in nblist:
646 (i, j) = k
647 if inblossom[j] == b:
648 i, j = j, i
649 bj = inblossom[j]
650 if (
651 bj != b
652 and label.get(bj) == 1
653 and ((bj not in bestedgeto) or slack(i, j) < slack(*bestedgeto[bj]))
654 ):
655 bestedgeto[bj] = k
656 # Forget about least-slack edge of the subblossom.
657 bestedge[bv] = None
658 b.mybestedges = list(bestedgeto.values())
659 # Select bestedge[b].
660 mybestedge = None
661 bestedge[b] = None
662 for k in b.mybestedges:
663 kslack = slack(*k)
664 if mybestedge is None or kslack < mybestslack:
665 mybestedge = k
666 mybestslack = kslack
667 bestedge[b] = mybestedge
669 # Expand the given top-level blossom.
670 def expandBlossom(b, endstage):
671 # This is an obnoxiously complicated recursive function for the sake of
672 # a stack-transformation. So, we hack around the complexity by using
673 # a trampoline pattern. By yielding the arguments to each recursive
674 # call, we keep the actual callstack flat.
676 def _recurse(b, endstage):
677 # Convert sub-blossoms into top-level blossoms.
678 for s in b.childs:
679 blossomparent[s] = None
680 if isinstance(s, Blossom):
681 if endstage and blossomdual[s] == 0:
682 # Recursively expand this sub-blossom.
683 yield s
684 else:
685 for v in s.leaves():
686 inblossom[v] = s
687 else:
688 inblossom[s] = s
689 # If we expand a T-blossom during a stage, its sub-blossoms must be
690 # relabeled.
691 if (not endstage) and label.get(b) == 2:
692 # Start at the sub-blossom through which the expanding
693 # blossom obtained its label, and relabel sub-blossoms untili
694 # we reach the base.
695 # Figure out through which sub-blossom the expanding blossom
696 # obtained its label initially.
697 entrychild = inblossom[labeledge[b][1]]
698 # Decide in which direction we will go round the blossom.
699 j = b.childs.index(entrychild)
700 if j & 1:
701 # Start index is odd; go forward and wrap.
702 j -= len(b.childs)
703 jstep = 1
704 else:
705 # Start index is even; go backward.
706 jstep = -1
707 # Move along the blossom until we get to the base.
708 v, w = labeledge[b]
709 while j != 0:
710 # Relabel the T-sub-blossom.
711 if jstep == 1:
712 p, q = b.edges[j]
713 else:
714 q, p = b.edges[j - 1]
715 label[w] = None
716 label[q] = None
717 assignLabel(w, 2, v)
718 # Step to the next S-sub-blossom and note its forward edge.
719 allowedge[(p, q)] = allowedge[(q, p)] = True
720 j += jstep
721 if jstep == 1:
722 v, w = b.edges[j]
723 else:
724 w, v = b.edges[j - 1]
725 # Step to the next T-sub-blossom.
726 allowedge[(v, w)] = allowedge[(w, v)] = True
727 j += jstep
728 # Relabel the base T-sub-blossom WITHOUT stepping through to
729 # its mate (so don't call assignLabel).
730 bw = b.childs[j]
731 label[w] = label[bw] = 2
732 labeledge[w] = labeledge[bw] = (v, w)
733 bestedge[bw] = None
734 # Continue along the blossom until we get back to entrychild.
735 j += jstep
736 while b.childs[j] != entrychild:
737 # Examine the vertices of the sub-blossom to see whether
738 # it is reachable from a neighboring S-vertex outside the
739 # expanding blossom.
740 bv = b.childs[j]
741 if label.get(bv) == 1:
742 # This sub-blossom just got label S through one of its
743 # neighbors; leave it be.
744 j += jstep
745 continue
746 if isinstance(bv, Blossom):
747 for v in bv.leaves():
748 if label.get(v):
749 break
750 else:
751 v = bv
752 # If the sub-blossom contains a reachable vertex, assign
753 # label T to the sub-blossom.
754 if label.get(v):
755 assert label[v] == 2
756 assert inblossom[v] == bv
757 label[v] = None
758 label[mate[blossombase[bv]]] = None
759 assignLabel(v, 2, labeledge[v][0])
760 j += jstep
761 # Remove the expanded blossom entirely.
762 label.pop(b, None)
763 labeledge.pop(b, None)
764 bestedge.pop(b, None)
765 del blossomparent[b]
766 del blossombase[b]
767 del blossomdual[b]
769 # Now, we apply the trampoline pattern. We simulate a recursive
770 # callstack by maintaining a stack of generators, each yielding a
771 # sequence of function arguments. We grow the stack by appending a call
772 # to _recurse on each argument tuple, and shrink the stack whenever a
773 # generator is exhausted.
774 stack = [_recurse(b, endstage)]
775 while stack:
776 top = stack[-1]
777 for s in top:
778 stack.append(_recurse(s, endstage))
779 break
780 else:
781 stack.pop()
783 # Swap matched/unmatched edges over an alternating path through blossom b
784 # between vertex v and the base vertex. Keep blossom bookkeeping
785 # consistent.
786 def augmentBlossom(b, v):
787 # This is an obnoxiously complicated recursive function for the sake of
788 # a stack-transformation. So, we hack around the complexity by using
789 # a trampoline pattern. By yielding the arguments to each recursive
790 # call, we keep the actual callstack flat.
792 def _recurse(b, v):
793 # Bubble up through the blossom tree from vertex v to an immediate
794 # sub-blossom of b.
795 t = v
796 while blossomparent[t] != b:
797 t = blossomparent[t]
798 # Recursively deal with the first sub-blossom.
799 if isinstance(t, Blossom):
800 yield (t, v)
801 # Decide in which direction we will go round the blossom.
802 i = j = b.childs.index(t)
803 if i & 1:
804 # Start index is odd; go forward and wrap.
805 j -= len(b.childs)
806 jstep = 1
807 else:
808 # Start index is even; go backward.
809 jstep = -1
810 # Move along the blossom until we get to the base.
811 while j != 0:
812 # Step to the next sub-blossom and augment it recursively.
813 j += jstep
814 t = b.childs[j]
815 if jstep == 1:
816 w, x = b.edges[j]
817 else:
818 x, w = b.edges[j - 1]
819 if isinstance(t, Blossom):
820 yield (t, w)
821 # Step to the next sub-blossom and augment it recursively.
822 j += jstep
823 t = b.childs[j]
824 if isinstance(t, Blossom):
825 yield (t, x)
826 # Match the edge connecting those sub-blossoms.
827 mate[w] = x
828 mate[x] = w
829 # Rotate the list of sub-blossoms to put the new base at the front.
830 b.childs = b.childs[i:] + b.childs[:i]
831 b.edges = b.edges[i:] + b.edges[:i]
832 blossombase[b] = blossombase[b.childs[0]]
833 assert blossombase[b] == v
835 # Now, we apply the trampoline pattern. We simulate a recursive
836 # callstack by maintaining a stack of generators, each yielding a
837 # sequence of function arguments. We grow the stack by appending a call
838 # to _recurse on each argument tuple, and shrink the stack whenever a
839 # generator is exhausted.
840 stack = [_recurse(b, v)]
841 while stack:
842 top = stack[-1]
843 for args in top:
844 stack.append(_recurse(*args))
845 break
846 else:
847 stack.pop()
849 # Swap matched/unmatched edges over an alternating path between two
850 # single vertices. The augmenting path runs through S-vertices v and w.
851 def augmentMatching(v, w):
852 for s, j in ((v, w), (w, v)):
853 # Match vertex s to vertex j. Then trace back from s
854 # until we find a single vertex, swapping matched and unmatched
855 # edges as we go.
856 while 1:
857 bs = inblossom[s]
858 assert label[bs] == 1
859 assert (labeledge[bs] is None and blossombase[bs] not in mate) or (
860 labeledge[bs][0] == mate[blossombase[bs]]
861 )
862 # Augment through the S-blossom from s to base.
863 if isinstance(bs, Blossom):
864 augmentBlossom(bs, s)
865 # Update mate[s]
866 mate[s] = j
867 # Trace one step back.
868 if labeledge[bs] is None:
869 # Reached single vertex; stop.
870 break
871 t = labeledge[bs][0]
872 bt = inblossom[t]
873 assert label[bt] == 2
874 # Trace one more step back.
875 s, j = labeledge[bt]
876 # Augment through the T-blossom from j to base.
877 assert blossombase[bt] == t
878 if isinstance(bt, Blossom):
879 augmentBlossom(bt, j)
880 # Update mate[j]
881 mate[j] = s
883 # Verify that the optimum solution has been reached.
884 def verifyOptimum():
885 if maxcardinality:
886 # Vertices may have negative dual;
887 # find a constant non-negative number to add to all vertex duals.
888 vdualoffset = max(0, -min(dualvar.values()))
889 else:
890 vdualoffset = 0
891 # 0. all dual variables are non-negative
892 assert min(dualvar.values()) + vdualoffset >= 0
893 assert len(blossomdual) == 0 or min(blossomdual.values()) >= 0
894 # 0. all edges have non-negative slack and
895 # 1. all matched edges have zero slack;
896 for i, j, d in G.edges(data=True):
897 wt = d.get(weight, 1)
898 if i == j:
899 continue # ignore self-loops
900 s = dualvar[i] + dualvar[j] - 2 * wt
901 iblossoms = [i]
902 jblossoms = [j]
903 while blossomparent[iblossoms[-1]] is not None:
904 iblossoms.append(blossomparent[iblossoms[-1]])
905 while blossomparent[jblossoms[-1]] is not None:
906 jblossoms.append(blossomparent[jblossoms[-1]])
907 iblossoms.reverse()
908 jblossoms.reverse()
909 for bi, bj in zip(iblossoms, jblossoms):
910 if bi != bj:
911 break
912 s += 2 * blossomdual[bi]
913 assert s >= 0
914 if mate.get(i) == j or mate.get(j) == i:
915 assert mate[i] == j and mate[j] == i
916 assert s == 0
917 # 2. all single vertices have zero dual value;
918 for v in gnodes:
919 assert (v in mate) or dualvar[v] + vdualoffset == 0
920 # 3. all blossoms with positive dual value are full.
921 for b in blossomdual:
922 if blossomdual[b] > 0:
923 assert len(b.edges) % 2 == 1
924 for i, j in b.edges[1::2]:
925 assert mate[i] == j and mate[j] == i
926 # Ok.
928 # Main loop: continue until no further improvement is possible.
929 while 1:
930 # Each iteration of this loop is a "stage".
931 # A stage finds an augmenting path and uses that to improve
932 # the matching.
934 # Remove labels from top-level blossoms/vertices.
935 label.clear()
936 labeledge.clear()
938 # Forget all about least-slack edges.
939 bestedge.clear()
940 for b in blossomdual:
941 b.mybestedges = None
943 # Loss of labeling means that we can not be sure that currently
944 # allowable edges remain allowable throughout this stage.
945 allowedge.clear()
947 # Make queue empty.
948 queue[:] = []
950 # Label single blossoms/vertices with S and put them in the queue.
951 for v in gnodes:
952 if (v not in mate) and label.get(inblossom[v]) is None:
953 assignLabel(v, 1, None)
955 # Loop until we succeed in augmenting the matching.
956 augmented = 0
957 while 1:
958 # Each iteration of this loop is a "substage".
959 # A substage tries to find an augmenting path;
960 # if found, the path is used to improve the matching and
961 # the stage ends. If there is no augmenting path, the
962 # primal-dual method is used to pump some slack out of
963 # the dual variables.
965 # Continue labeling until all vertices which are reachable
966 # through an alternating path have got a label.
967 while queue and not augmented:
968 # Take an S vertex from the queue.
969 v = queue.pop()
970 assert label[inblossom[v]] == 1
972 # Scan its neighbors:
973 for w in G.neighbors(v):
974 if w == v:
975 continue # ignore self-loops
976 # w is a neighbor to v
977 bv = inblossom[v]
978 bw = inblossom[w]
979 if bv == bw:
980 # this edge is internal to a blossom; ignore it
981 continue
982 if (v, w) not in allowedge:
983 kslack = slack(v, w)
984 if kslack <= 0:
985 # edge k has zero slack => it is allowable
986 allowedge[(v, w)] = allowedge[(w, v)] = True
987 if (v, w) in allowedge:
988 if label.get(bw) is None:
989 # (C1) w is a free vertex;
990 # label w with T and label its mate with S (R12).
991 assignLabel(w, 2, v)
992 elif label.get(bw) == 1:
993 # (C2) w is an S-vertex (not in the same blossom);
994 # follow back-links to discover either an
995 # augmenting path or a new blossom.
996 base = scanBlossom(v, w)
997 if base is not NoNode:
998 # Found a new blossom; add it to the blossom
999 # bookkeeping and turn it into an S-blossom.
1000 addBlossom(base, v, w)
1001 else:
1002 # Found an augmenting path; augment the
1003 # matching and end this stage.
1004 augmentMatching(v, w)
1005 augmented = 1
1006 break
1007 elif label.get(w) is None:
1008 # w is inside a T-blossom, but w itself has not
1009 # yet been reached from outside the blossom;
1010 # mark it as reached (we need this to relabel
1011 # during T-blossom expansion).
1012 assert label[bw] == 2
1013 label[w] = 2
1014 labeledge[w] = (v, w)
1015 elif label.get(bw) == 1:
1016 # keep track of the least-slack non-allowable edge to
1017 # a different S-blossom.
1018 if bestedge.get(bv) is None or kslack < slack(*bestedge[bv]):
1019 bestedge[bv] = (v, w)
1020 elif label.get(w) is None:
1021 # w is a free vertex (or an unreached vertex inside
1022 # a T-blossom) but we can not reach it yet;
1023 # keep track of the least-slack edge that reaches w.
1024 if bestedge.get(w) is None or kslack < slack(*bestedge[w]):
1025 bestedge[w] = (v, w)
1027 if augmented:
1028 break
1030 # There is no augmenting path under these constraints;
1031 # compute delta and reduce slack in the optimization problem.
1032 # (Note that our vertex dual variables, edge slacks and delta's
1033 # are pre-multiplied by two.)
1034 deltatype = -1
1035 delta = deltaedge = deltablossom = None
1037 # Compute delta1: the minimum value of any vertex dual.
1038 if not maxcardinality:
1039 deltatype = 1
1040 delta = min(dualvar.values())
1042 # Compute delta2: the minimum slack on any edge between
1043 # an S-vertex and a free vertex.
1044 for v in G.nodes():
1045 if label.get(inblossom[v]) is None and bestedge.get(v) is not None:
1046 d = slack(*bestedge[v])
1047 if deltatype == -1 or d < delta:
1048 delta = d
1049 deltatype = 2
1050 deltaedge = bestedge[v]
1052 # Compute delta3: half the minimum slack on any edge between
1053 # a pair of S-blossoms.
1054 for b in blossomparent:
1055 if (
1056 blossomparent[b] is None
1057 and label.get(b) == 1
1058 and bestedge.get(b) is not None
1059 ):
1060 kslack = slack(*bestedge[b])
1061 if allinteger:
1062 assert (kslack % 2) == 0
1063 d = kslack // 2
1064 else:
1065 d = kslack / 2.0
1066 if deltatype == -1 or d < delta:
1067 delta = d
1068 deltatype = 3
1069 deltaedge = bestedge[b]
1071 # Compute delta4: minimum z variable of any T-blossom.
1072 for b in blossomdual:
1073 if (
1074 blossomparent[b] is None
1075 and label.get(b) == 2
1076 and (deltatype == -1 or blossomdual[b] < delta)
1077 ):
1078 delta = blossomdual[b]
1079 deltatype = 4
1080 deltablossom = b
1082 if deltatype == -1:
1083 # No further improvement possible; max-cardinality optimum
1084 # reached. Do a final delta update to make the optimum
1085 # verifiable.
1086 assert maxcardinality
1087 deltatype = 1
1088 delta = max(0, min(dualvar.values()))
1090 # Update dual variables according to delta.
1091 for v in gnodes:
1092 if label.get(inblossom[v]) == 1:
1093 # S-vertex: 2*u = 2*u - 2*delta
1094 dualvar[v] -= delta
1095 elif label.get(inblossom[v]) == 2:
1096 # T-vertex: 2*u = 2*u + 2*delta
1097 dualvar[v] += delta
1098 for b in blossomdual:
1099 if blossomparent[b] is None:
1100 if label.get(b) == 1:
1101 # top-level S-blossom: z = z + 2*delta
1102 blossomdual[b] += delta
1103 elif label.get(b) == 2:
1104 # top-level T-blossom: z = z - 2*delta
1105 blossomdual[b] -= delta
1107 # Take action at the point where minimum delta occurred.
1108 if deltatype == 1:
1109 # No further improvement possible; optimum reached.
1110 break
1111 elif deltatype == 2:
1112 # Use the least-slack edge to continue the search.
1113 (v, w) = deltaedge
1114 assert label[inblossom[v]] == 1
1115 allowedge[(v, w)] = allowedge[(w, v)] = True
1116 queue.append(v)
1117 elif deltatype == 3:
1118 # Use the least-slack edge to continue the search.
1119 (v, w) = deltaedge
1120 allowedge[(v, w)] = allowedge[(w, v)] = True
1121 assert label[inblossom[v]] == 1
1122 queue.append(v)
1123 elif deltatype == 4:
1124 # Expand the least-z blossom.
1125 expandBlossom(deltablossom, False)
1127 # End of a this substage.
1129 # Paranoia check that the matching is symmetric.
1130 for v in mate:
1131 assert mate[mate[v]] == v
1133 # Stop when no more augmenting path can be found.
1134 if not augmented:
1135 break
1137 # End of a stage; expand all S-blossoms which have zero dual.
1138 for b in list(blossomdual.keys()):
1139 if b not in blossomdual:
1140 continue # already expanded
1141 if blossomparent[b] is None and label.get(b) == 1 and blossomdual[b] == 0:
1142 expandBlossom(b, True)
1144 # Verify that we reached the optimum solution (only for integer weights).
1145 if allinteger:
1146 verifyOptimum()
1148 return matching_dict_to_set(mate)