Statistics
| Branch: | Revision:

iof-tools / networkxMiCe / networkx-master / networkx / algorithms / matching.py @ 5cef0f13

History | View | Annotate | Download (37 KB)

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)