ioftools / networkxMiCe / networkxmaster / networkx / algorithms / matching.py @ 5cef0f13
History  View  Annotate  Download (37 KB)
1 
# Copyright 2016 NetworkX developers.


2 
# Copyright (C) 20042019 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 
selfloops 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 maximumweighted 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 maximumcardinality matching

204 
with maximum weight among all maximumcardinality 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 "primaldual" 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 "primaldual" 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 nontrivial blossom or subblossom."""

260  
261 
__slots__ = ['childs', 'edges', 'mybestedges'] 
262  
263 
# b.childs is an ordered list of b's subblossoms, 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 toplevel Sblossom,

271 
# b.mybestedges is a list of leastslack edges to neighbouring

272 
# Sblossoms, 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 toplevel blossom,

305 
# label.get(b) is None if b is unlabeled (free),

306 
# 1 if b is an Sblossom,

307 
# 2 if b is a Tblossom.

308 
# The label of a vertex is found by looking at the label of its toplevel

309 
# containing blossom.

310 
# If v is a vertex inside a Tblossom, label[v] is 2 iff v is reachable

311 
# from an Svertex outside the blossom.

312 
# Labels are assigned during a stage and reset after each augmentation.

313 
label = {} 
314  
315 
# If b is a labeled toplevel 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 Tblossom 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 toplevel blossom to which v

324 
# belongs.

325 
# If v is a toplevel vertex, inblossom[v] == v since v is itself

326 
# a (trivial) toplevel blossom.

327 
# Initially all vertices are toplevel trivial blossoms.

328 
inblossom = dict(zip(gnodes, gnodes)) 
329  
330 
# If b is a subblossom,

331 
# blossomparent[b] is its immediate parent (sub)blossom.

332 
# If b is a toplevel 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 subblossom).

337 
blossombase = dict(zip(gnodes, gnodes)) 
338  
339 
# If w is a free vertex (or an unreached vertex inside a Tblossom),

340 
# bestedge[w] = (v, w) is the leastslack edge from an Svertex,

341 
# or None if there is no such edge.

342 
# If b is a (possibly trivial) toplevel Sblossom,

343 
# bestedge[b] = (v, w) is the leastslack edge to a different Sblossom

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 nontrivial 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 Svertices.

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 toplevel 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 Svertex/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 Tvertex/blossom; assign label S to its mate.

391 
# (If b is a nontrivial 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 Tblossom; 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 Svertices v and w.

434 
# Label the new blossom as S; set its dual variable to zero;

435 
# relabel its Tvertices 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 subblossoms 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 subblossom; 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 Tvertex now turns into an Svertex because it becomes

484 
# part of an Sblossom; 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 leastslack edges.

493 
nblist = bv.mybestedges 
494 
# The subblossom won't need this data again.

495 
bv.mybestedges = None

496 
else:

497 
# This subblossom does not have a list of leastslack

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 leastslack 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 toplevel blossom.

530 
def expandBlossom(b, endstage): 
531 
# Convert subblossoms into toplevel 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 subblossom.

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 Tblossom during a stage, its subblossoms must be

544 
# relabeled.

545 
if (not endstage) and label.get(b) == 2: 
546 
# Start at the subblossom through which the expanding

547 
# blossom obtained its label, and relabel subblossoms untili

548 
# we reach the base.

549 
# Figure out through which subblossom 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 Tsubblossom.

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 Ssubblossom 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 Tsubblossom.

580 
allowedge[(v, w)] = allowedge[(w, v)] = True

581 
j += jstep 
582 
# Relabel the base Tsubblossom 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 subblossom to see whether

592 
# it is reachable from a neighbouring Svertex outside the

593 
# expanding blossom.

594 
bv = b.childs[j] 
595 
if label.get(bv) == 1: 
596 
# This subblossom 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 subblossom contains a reachable vertex, assign

607 
# label T to the subblossom.

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 
# subblossom of b.

629 
t = v 
630 
while blossomparent[t] != b:

631 
t = blossomparent[t] 
632 
# Recursively deal with the first subblossom.

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 subblossom 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 subblossom 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 subblossoms.

661 
mate[w] = x 
662 
mate[x] = w 
663 
# Rotate the list of subblossoms 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 Svertices 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 Sblossom 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 Tblossom 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 nonnegative 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 nonnegative

712 
assert min(dualvar.values()) + vdualoffset >= 0 
713 
assert len(blossomdual) == 0 or min(blossomdual.values()) >= 0 
714 
# 0. all edges have nonnegative 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 selfloops 
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 toplevel blossoms/vertices.

756 
label.clear() 
757 
labeledge.clear() 
758  
759 
# Forget all about leastslack 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 
# primaldual 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 selfloops 
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 Svertex (not in the same blossom);

817 
# follow backlinks 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 Sblossom.

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 Tblossom, 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 Tblossom 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 leastslack nonallowable edge to

840 
# a different Sblossom.

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 Tblossom) but we can not reach it yet;

847 
# keep track of the leastslack 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 premultiplied 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 Svertex 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 Sblossoms.

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 Tblossom.

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; maxcardinality 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 
# Svertex: 2*u = 2*u  2*delta

914 
dualvar[v] = delta 
915 
elif label.get(inblossom[v]) == 2: 
916 
# Tvertex: 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 
# toplevel Sblossom: z = z + 2*delta

922 
blossomdual[b] += delta 
923 
elif label.get(b) == 2: 
924 
# toplevel Tblossom: 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 leastslack 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 leastslack 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 leastz 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 Sblossoms 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)
