matching.py 43 KB

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