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