matching.py 41 KB

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