Statistics
| Branch: | Revision:

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

 1 ```# Copyright 2016 NetworkX developers. ``` ```# Copyright (C) 2004-2019 by ``` ```# Aric Hagberg ``` ```# Dan Schult ``` ```# Pieter Swart ``` ```# ``` ```# Copyright (C) 2008 by ``` ```# Joris van Rantwijk. ``` ```# ``` ```# Copyright (C) 2011 by ``` ```# Nicholas Mancuso ``` ```# ``` ```# All rights reserved. ``` ```# BSD license. ``` ```"""Functions for computing and verifying matchings in a graph.""" ``` ```from collections import Counter ``` ```from itertools import combinations ``` ```from itertools import repeat ``` ```__all__ = ['is_matching', 'is_maximal_matching', 'is_perfect_matching', ``` ``` 'max_weight_matching', 'maximal_matching'] ``` ```def maximal_matching(G): ``` ``` r"""Find a maximal matching in the graph. ``` ``` ``` ``` A matching is a subset of edges in which no node occurs more than once. ``` ``` A maximal matching cannot add more edges and still be a matching. ``` ``` ``` ``` Parameters ``` ``` ---------- ``` ``` G : NetworkX graph ``` ``` Undirected graph ``` ``` ``` ``` Returns ``` ``` ------- ``` ``` matching : set ``` ``` A maximal matching of the graph. ``` ``` ``` ``` Notes ``` ``` ----- ``` ``` The algorithm greedily selects a maximal matching M of the graph G ``` ``` (i.e. no superset of M exists). It runs in \$O(|E|)\$ time. ``` ``` """ ``` ``` matching = set() ``` ``` nodes = set() ``` ``` for u, v in G.edges(): ``` ``` # If the edge isn't covered, add it to the matching ``` ``` # then remove neighborhood of u and v from consideration. ``` ``` if u not in nodes and v not in nodes and u != v: ``` ``` matching.add((u, v)) ``` ``` nodes.add(u) ``` ``` nodes.add(v) ``` ``` return matching ``` ```def matching_dict_to_set(matching): ``` ``` """Converts a dictionary representing a matching (as returned by ``` ``` :func:`max_weight_matching`) to a set representing a matching (as ``` ``` returned by :func:`maximal_matching`). ``` ``` ``` ``` In the definition of maximal matching adopted by NetworkX, ``` ``` self-loops are not allowed, so the provided dictionary is expected ``` ``` to never have any mapping from a key to itself. However, the ``` ``` dictionary is expected to have mirrored key/value pairs, for ``` ``` example, key ``u`` with value ``v`` and key ``v`` with value ``u``. ``` ``` ``` ``` """ ``` ``` # Need to compensate for the fact that both pairs (u, v) and (v, u) ``` ``` # appear in `matching.items()`, so we use a set of sets. This way, ``` ``` # only the (frozen)set `{u, v}` appears as an element in the ``` ``` # returned set. ``` ``` return set((u, v) for (u, v) in set(map(frozenset, matching.items()))) ``` ```def is_matching(G, matching): ``` ``` """Decides whether the given set or dictionary represents a valid ``` ``` matching in ``G``. ``` ``` ``` ``` A *matching* in a graph is a set of edges in which no two distinct ``` ``` edges share a common endpoint. ``` ``` ``` ``` Parameters ``` ``` ---------- ``` ``` G : NetworkX graph ``` ``` ``` ``` matching : dict or set ``` ``` A dictionary or set representing a matching. If a dictionary, it ``` ``` must have ``matching[u] == v`` and ``matching[v] == u`` for each ``` ``` edge ``(u, v)`` in the matching. If a set, it must have elements ``` ``` of the form ``(u, v)``, where ``(u, v)`` is an edge in the ``` ``` matching. ``` ``` ``` ``` Returns ``` ``` ------- ``` ``` bool ``` ``` Whether the given set or dictionary represents a valid matching ``` ``` in the graph. ``` ``` ``` ``` """ ``` ``` if isinstance(matching, dict): ``` ``` matching = matching_dict_to_set(matching) ``` ``` # TODO This is parallelizable. ``` ``` return all(len(set(e1) & set(e2)) == 0 ``` ``` for e1, e2 in combinations(matching, 2)) ``` ```def is_maximal_matching(G, matching): ``` ``` """Decides whether the given set or dictionary represents a valid ``` ``` maximal matching in ``G``. ``` ``` ``` ``` A *maximal matching* in a graph is a matching in which adding any ``` ``` edge would cause the set to no longer be a valid matching. ``` ``` ``` ``` Parameters ``` ``` ---------- ``` ``` G : NetworkX graph ``` ``` ``` ``` matching : dict or set ``` ``` A dictionary or set representing a matching. If a dictionary, it ``` ``` must have ``matching[u] == v`` and ``matching[v] == u`` for each ``` ``` edge ``(u, v)`` in the matching. If a set, it must have elements ``` ``` of the form ``(u, v)``, where ``(u, v)`` is an edge in the ``` ``` matching. ``` ``` ``` ``` Returns ``` ``` ------- ``` ``` bool ``` ``` Whether the given set or dictionary represents a valid maximal ``` ``` matching in the graph. ``` ``` ``` ``` """ ``` ``` if isinstance(matching, dict): ``` ``` matching = matching_dict_to_set(matching) ``` ``` # If the given set is not a matching, then it is not a maximal matching. ``` ``` if not is_matching(G, matching): ``` ``` return False ``` ``` # A matching is maximal if adding any unmatched edge to it causes ``` ``` # the resulting set to *not* be a valid matching. ``` ``` # ``` ``` # HACK Since the `matching_dict_to_set` function returns a set of ``` ``` # sets, we need to convert the list of edges to a set of sets in ``` ``` # order for the set difference function to work. Ideally, we would ``` ``` # just be able to do `set(G.edges()) - matching`. ``` ``` all_edges = set(map(frozenset, G.edges())) ``` ``` matched_edges = set(map(frozenset, matching)) ``` ``` unmatched_edges = all_edges - matched_edges ``` ``` # TODO This is parallelizable. ``` ``` return all(not is_matching(G, matching | {e}) for e in unmatched_edges) ``` ```def is_perfect_matching(G, matching): ``` ``` """Decides whether the given set represents a valid perfect matching in ``` ``` ``G``. ``` ``` ``` ``` A *perfect matching* in a graph is a matching in which exactly one edge ``` ``` is incident upon each vertex. ``` ``` ``` ``` Parameters ``` ``` ---------- ``` ``` G : NetworkX graph ``` ``` ``` ``` matching : dict or set ``` ``` A dictionary or set representing a matching. If a dictionary, it ``` ``` must have ``matching[u] == v`` and ``matching[v] == u`` for each ``` ``` edge ``(u, v)`` in the matching. If a set, it must have elements ``` ``` of the form ``(u, v)``, where ``(u, v)`` is an edge in the ``` ``` matching. ``` ``` ``` ``` Returns ``` ``` ------- ``` ``` bool ``` ``` Whether the given set or dictionary represents a valid perfect ``` ``` matching in the graph. ``` ``` ``` ``` """ ``` ``` if isinstance(matching, dict): ``` ``` matching = matching_dict_to_set(matching) ``` ``` if not is_matching(G, matching): ``` ``` return False ``` ``` counts = Counter(sum(matching, ())) ``` ``` return all(counts[v] == 1 for v in G) ``` ```def max_weight_matching(G, maxcardinality=False, weight='weight'): ``` ``` """Compute a maximum-weighted matching of G. ``` ``` ``` ``` A matching is a subset of edges in which no node occurs more than once. ``` ``` The weight of a matching is the sum of the weights of its edges. ``` ``` A maximal matching cannot add more edges and still be a matching. ``` ``` The cardinality of a matching is the number of matched edges. ``` ``` ``` ``` Parameters ``` ``` ---------- ``` ``` G : NetworkX graph ``` ``` Undirected graph ``` ``` ``` ``` maxcardinality: bool, optional (default=False) ``` ``` If maxcardinality is True, compute the maximum-cardinality matching ``` ``` with maximum weight among all maximum-cardinality matchings. ``` ``` ``` ``` weight: string, optional (default='weight') ``` ``` Edge data key corresponding to the edge weight. ``` ``` If key not found, uses 1 as weight. ``` ``` ``` ``` ``` ``` Returns ``` ``` ------- ``` ``` matching : set ``` ``` A maximal matching of the graph. ``` ``` ``` ``` Notes ``` ``` ----- ``` ``` If G has edges with weight attributes the edge data are used as ``` ``` weight values else the weights are assumed to be 1. ``` ``` ``` ``` This function takes time O(number_of_nodes ** 3). ``` ``` ``` ``` If all edge weights are integers, the algorithm uses only integer ``` ``` computations. If floating point weights are used, the algorithm ``` ``` could return a slightly suboptimal matching due to numeric ``` ``` precision errors. ``` ``` ``` ``` This method is based on the "blossom" method for finding augmenting ``` ``` paths and the "primal-dual" method for finding a matching of maximum ``` ``` weight, both methods invented by Jack Edmonds _. ``` ``` ``` ``` Bipartite graphs can also be matched using the functions present in ``` ``` :mod:`networkx.algorithms.bipartite.matching`. ``` ``` ``` ``` References ``` ``` ---------- ``` ``` ..  "Efficient Algorithms for Finding Maximum Matching in Graphs", ``` ``` Zvi Galil, ACM Computing Surveys, 1986. ``` ``` """ ``` ``` # ``` ``` # The algorithm is taken from "Efficient Algorithms for Finding Maximum ``` ``` # Matching in Graphs" by Zvi Galil, ACM Computing Surveys, 1986. ``` ``` # It is based on the "blossom" method for finding augmenting paths and ``` ``` # the "primal-dual" method for finding a matching of maximum weight, both ``` ``` # methods invented by Jack Edmonds. ``` ``` # ``` ``` # A C program for maximum weight matching by Ed Rothberg was used ``` ``` # extensively to validate this new code. ``` ``` # ``` ``` # Many terms used in the code comments are explained in the paper ``` ``` # by Galil. You will probably need the paper to make sense of this code. ``` ``` # ``` ``` class NoNode: ``` ``` """Dummy value which is different from any node.""" ``` ``` pass ``` ``` class Blossom: ``` ``` """Representation of a non-trivial blossom or sub-blossom.""" ``` ``` __slots__ = ['childs', 'edges', 'mybestedges'] ``` ``` # b.childs is an ordered list of b's sub-blossoms, starting with ``` ``` # the base and going round the blossom. ``` ``` # b.edges is the list of b's connecting edges, such that ``` ``` # b.edges[i] = (v, w) where v is a vertex in b.childs[i] ``` ``` # and w is a vertex in b.childs[wrap(i+1)]. ``` ``` # If b is a top-level S-blossom, ``` ``` # b.mybestedges is a list of least-slack edges to neighbouring ``` ``` # S-blossoms, or None if no such list has been computed yet. ``` ``` # This is used for efficient computation of delta3. ``` ``` # Generate the blossom's leaf vertices. ``` ``` def leaves(self): ``` ``` for t in self.childs: ``` ``` if isinstance(t, Blossom): ``` ``` for v in t.leaves(): ``` ``` yield v ``` ``` else: ``` ``` yield t ``` ``` # Get a list of vertices. ``` ``` gnodes = list(G) ``` ``` if not gnodes: ``` ``` return set() # don't bother with empty graphs ``` ``` # Find the maximum edge weight. ``` ``` maxweight = 0 ``` ``` allinteger = True ``` ``` for i, j, d in G.edges(data=True): ``` ``` wt = d.get(weight, 1) ``` ``` if i != j and wt > maxweight: ``` ``` maxweight = wt ``` ``` allinteger = allinteger and (str(type(wt)).split("'") ``` ``` in ('int', 'long')) ``` ``` # If v is a matched vertex, mate[v] is its partner vertex. ``` ``` # If v is a single vertex, v does not occur as a key in mate. ``` ``` # Initially all vertices are single; updated during augmentation. ``` ``` mate = {} ``` ``` # If b is a top-level blossom, ``` ``` # label.get(b) is None if b is unlabeled (free), ``` ``` # 1 if b is an S-blossom, ``` ``` # 2 if b is a T-blossom. ``` ``` # The label of a vertex is found by looking at the label of its top-level ``` ``` # containing blossom. ``` ``` # If v is a vertex inside a T-blossom, label[v] is 2 iff v is reachable ``` ``` # from an S-vertex outside the blossom. ``` ``` # Labels are assigned during a stage and reset after each augmentation. ``` ``` label = {} ``` ``` # If b is a labeled top-level blossom, ``` ``` # labeledge[b] = (v, w) is the edge through which b obtained its label ``` ``` # such that w is a vertex in b, or None if b's base vertex is single. ``` ``` # If w is a vertex inside a T-blossom and label[w] == 2, ``` ``` # labeledge[w] = (v, w) is an edge through which w is reachable from ``` ``` # outside the blossom. ``` ``` labeledge = {} ``` ``` # If v is a vertex, inblossom[v] is the top-level blossom to which v ``` ``` # belongs. ``` ``` # If v is a top-level vertex, inblossom[v] == v since v is itself ``` ``` # a (trivial) top-level blossom. ``` ``` # Initially all vertices are top-level trivial blossoms. ``` ``` inblossom = dict(zip(gnodes, gnodes)) ``` ``` # If b is a sub-blossom, ``` ``` # blossomparent[b] is its immediate parent (sub-)blossom. ``` ``` # If b is a top-level blossom, blossomparent[b] is None. ``` ``` blossomparent = dict(zip(gnodes, repeat(None))) ``` ``` # If b is a (sub-)blossom, ``` ``` # blossombase[b] is its base VERTEX (i.e. recursive sub-blossom). ``` ``` blossombase = dict(zip(gnodes, gnodes)) ``` ``` # If w is a free vertex (or an unreached vertex inside a T-blossom), ``` ``` # bestedge[w] = (v, w) is the least-slack edge from an S-vertex, ``` ``` # or None if there is no such edge. ``` ``` # If b is a (possibly trivial) top-level S-blossom, ``` ``` # bestedge[b] = (v, w) is the least-slack edge to a different S-blossom ``` ``` # (v inside b), or None if there is no such edge. ``` ``` # This is used for efficient computation of delta2 and delta3. ``` ``` bestedge = {} ``` ``` # If v is a vertex, ``` ``` # dualvar[v] = 2 * u(v) where u(v) is the v's variable in the dual ``` ``` # optimization problem (if all edge weights are integers, multiplication ``` ``` # by two ensures that all values remain integers throughout the algorithm). ``` ``` # Initially, u(v) = maxweight / 2. ``` ``` dualvar = dict(zip(gnodes, repeat(maxweight))) ``` ``` # If b is a non-trivial blossom, ``` ``` # blossomdual[b] = z(b) where z(b) is b's variable in the dual ``` ``` # optimization problem. ``` ``` blossomdual = {} ``` ``` # If (v, w) in allowedge or (w, v) in allowedg, then the edge ``` ``` # (v, w) is known to have zero slack in the optimization problem; ``` ``` # otherwise the edge may or may not have zero slack. ``` ``` allowedge = {} ``` ``` # Queue of newly discovered S-vertices. ``` ``` queue = [] ``` ``` # Return 2 * slack of edge (v, w) (does not work inside blossoms). ``` ``` def slack(v, w): ``` ``` return dualvar[v] + dualvar[w] - 2 * G[v][w].get(weight, 1) ``` ``` # Assign label t to the top-level blossom containing vertex w, ``` ``` # coming through an edge from vertex v. ``` ``` def assignLabel(w, t, v): ``` ``` b = inblossom[w] ``` ``` assert label.get(w) is None and label.get(b) is None ``` ``` label[w] = label[b] = t ``` ``` if v is not None: ``` ``` labeledge[w] = labeledge[b] = (v, w) ``` ``` else: ``` ``` labeledge[w] = labeledge[b] = None ``` ``` bestedge[w] = bestedge[b] = None ``` ``` if t == 1: ``` ``` # b became an S-vertex/blossom; add it(s vertices) to the queue. ``` ``` if isinstance(b, Blossom): ``` ``` queue.extend(b.leaves()) ``` ``` else: ``` ``` queue.append(b) ``` ``` elif t == 2: ``` ``` # b became a T-vertex/blossom; assign label S to its mate. ``` ``` # (If b is a non-trivial blossom, its base is the only vertex ``` ``` # with an external mate.) ``` ``` base = blossombase[b] ``` ``` assignLabel(mate[base], 1, base) ``` ``` # Trace back from vertices v and w to discover either a new blossom ``` ``` # or an augmenting path. Return the base vertex of the new blossom, ``` ``` # or NoNode if an augmenting path was found. ``` ``` def scanBlossom(v, w): ``` ``` # Trace back from v and w, placing breadcrumbs as we go. ``` ``` path = [] ``` ``` base = NoNode ``` ``` while v is not NoNode: ``` ``` # Look for a breadcrumb in v's blossom or put a new breadcrumb. ``` ``` b = inblossom[v] ``` ``` if label[b] & 4: ``` ``` base = blossombase[b] ``` ``` break ``` ``` assert label[b] == 1 ``` ``` path.append(b) ``` ``` label[b] = 5 ``` ``` # Trace one step back. ``` ``` if labeledge[b] is None: ``` ``` # The base of blossom b is single; stop tracing this path. ``` ``` assert blossombase[b] not in mate ``` ``` v = NoNode ``` ``` else: ``` ``` assert labeledge[b] == mate[blossombase[b]] ``` ``` v = labeledge[b] ``` ``` b = inblossom[v] ``` ``` assert label[b] == 2 ``` ``` # b is a T-blossom; trace one more step back. ``` ``` v = labeledge[b] ``` ``` # Swap v and w so that we alternate between both paths. ``` ``` if w is not NoNode: ``` ``` v, w = w, v ``` ``` # Remove breadcrumbs. ``` ``` for b in path: ``` ``` label[b] = 1 ``` ``` # Return base vertex, if we found one. ``` ``` return base ``` ``` # Construct a new blossom with given base, through S-vertices v and w. ``` ``` # Label the new blossom as S; set its dual variable to zero; ``` ``` # relabel its T-vertices to S and add them to the queue. ``` ``` def addBlossom(base, v, w): ``` ``` bb = inblossom[base] ``` ``` bv = inblossom[v] ``` ``` bw = inblossom[w] ``` ``` # Create blossom. ``` ``` b = Blossom() ``` ``` blossombase[b] = base ``` ``` blossomparent[b] = None ``` ``` blossomparent[bb] = b ``` ``` # Make list of sub-blossoms and their interconnecting edge endpoints. ``` ``` b.childs = path = [] ``` ``` b.edges = edgs = [(v, w)] ``` ``` # Trace back from v to base. ``` ``` while bv != bb: ``` ``` # Add bv to the new blossom. ``` ``` blossomparent[bv] = b ``` ``` path.append(bv) ``` ``` edgs.append(labeledge[bv]) ``` ``` assert label[bv] == 2 or (label[bv] == 1 and labeledge[ ``` ``` bv] == mate[blossombase[bv]]) ``` ``` # Trace one step back. ``` ``` v = labeledge[bv] ``` ``` bv = inblossom[v] ``` ``` # Add base sub-blossom; reverse lists. ``` ``` path.append(bb) ``` ``` path.reverse() ``` ``` edgs.reverse() ``` ``` # Trace back from w to base. ``` ``` while bw != bb: ``` ``` # Add bw to the new blossom. ``` ``` blossomparent[bw] = b ``` ``` path.append(bw) ``` ``` edgs.append((labeledge[bw], labeledge[bw])) ``` ``` assert label[bw] == 2 or (label[bw] == 1 and labeledge[ ``` ``` bw] == mate[blossombase[bw]]) ``` ``` # Trace one step back. ``` ``` w = labeledge[bw] ``` ``` bw = inblossom[w] ``` ``` # Set label to S. ``` ``` assert label[bb] == 1 ``` ``` label[b] = 1 ``` ``` labeledge[b] = labeledge[bb] ``` ``` # Set dual variable to zero. ``` ``` blossomdual[b] = 0 ``` ``` # Relabel vertices. ``` ``` for v in b.leaves(): ``` ``` if label[inblossom[v]] == 2: ``` ``` # This T-vertex now turns into an S-vertex because it becomes ``` ``` # part of an S-blossom; add it to the queue. ``` ``` queue.append(v) ``` ``` inblossom[v] = b ``` ``` # Compute b.mybestedges. ``` ``` bestedgeto = {} ``` ``` for bv in path: ``` ``` if isinstance(bv, Blossom): ``` ``` if bv.mybestedges is not None: ``` ``` # Walk this subblossom's least-slack edges. ``` ``` nblist = bv.mybestedges ``` ``` # The sub-blossom won't need this data again. ``` ``` bv.mybestedges = None ``` ``` else: ``` ``` # This subblossom does not have a list of least-slack ``` ``` # edges; get the information from the vertices. ``` ``` nblist = [(v, w) ``` ``` for v in bv.leaves() ``` ``` for w in G.neighbors(v) ``` ``` if v != w] ``` ``` else: ``` ``` nblist = [(bv, w) ``` ``` for w in G.neighbors(bv) ``` ``` if bv != w] ``` ``` for k in nblist: ``` ``` (i, j) = k ``` ``` if inblossom[j] == b: ``` ``` i, j = j, i ``` ``` bj = inblossom[j] ``` ``` if (bj != b and label.get(bj) == 1 and ``` ``` ((bj not in bestedgeto) or ``` ``` slack(i, j) < slack(*bestedgeto[bj]))): ``` ``` bestedgeto[bj] = k ``` ``` # Forget about least-slack edge of the subblossom. ``` ``` bestedge[bv] = None ``` ``` b.mybestedges = list(bestedgeto.values()) ``` ``` # Select bestedge[b]. ``` ``` mybestedge = None ``` ``` bestedge[b] = None ``` ``` for k in b.mybestedges: ``` ``` kslack = slack(*k) ``` ``` if mybestedge is None or kslack < mybestslack: ``` ``` mybestedge = k ``` ``` mybestslack = kslack ``` ``` bestedge[b] = mybestedge ``` ``` # Expand the given top-level blossom. ``` ``` def expandBlossom(b, endstage): ``` ``` # Convert sub-blossoms into top-level blossoms. ``` ``` for s in b.childs: ``` ``` blossomparent[s] = None ``` ``` if isinstance(s, Blossom): ``` ``` if endstage and blossomdual[s] == 0: ``` ``` # Recursively expand this sub-blossom. ``` ``` expandBlossom(s, endstage) ``` ``` else: ``` ``` for v in s.leaves(): ``` ``` inblossom[v] = s ``` ``` else: ``` ``` inblossom[s] = s ``` ``` # If we expand a T-blossom during a stage, its sub-blossoms must be ``` ``` # relabeled. ``` ``` if (not endstage) and label.get(b) == 2: ``` ``` # Start at the sub-blossom through which the expanding ``` ``` # blossom obtained its label, and relabel sub-blossoms untili ``` ``` # we reach the base. ``` ``` # Figure out through which sub-blossom the expanding blossom ``` ``` # obtained its label initially. ``` ``` entrychild = inblossom[labeledge[b]] ``` ``` # Decide in which direction we will go round the blossom. ``` ``` j = b.childs.index(entrychild) ``` ``` if j & 1: ``` ``` # Start index is odd; go forward and wrap. ``` ``` j -= len(b.childs) ``` ``` jstep = 1 ``` ``` else: ``` ``` # Start index is even; go backward. ``` ``` jstep = -1 ``` ``` # Move along the blossom until we get to the base. ``` ``` v, w = labeledge[b] ``` ``` while j != 0: ``` ``` # Relabel the T-sub-blossom. ``` ``` if jstep == 1: ``` ``` p, q = b.edges[j] ``` ``` else: ``` ``` q, p = b.edges[j - 1] ``` ``` label[w] = None ``` ``` label[q] = None ``` ``` assignLabel(w, 2, v) ``` ``` # Step to the next S-sub-blossom and note its forward edge. ``` ``` allowedge[(p, q)] = allowedge[(q, p)] = True ``` ``` j += jstep ``` ``` if jstep == 1: ``` ``` v, w = b.edges[j] ``` ``` else: ``` ``` w, v = b.edges[j - 1] ``` ``` # Step to the next T-sub-blossom. ``` ``` allowedge[(v, w)] = allowedge[(w, v)] = True ``` ``` j += jstep ``` ``` # Relabel the base T-sub-blossom WITHOUT stepping through to ``` ``` # its mate (so don't call assignLabel). ``` ``` bw = b.childs[j] ``` ``` label[w] = label[bw] = 2 ``` ``` labeledge[w] = labeledge[bw] = (v, w) ``` ``` bestedge[bw] = None ``` ``` # Continue along the blossom until we get back to entrychild. ``` ``` j += jstep ``` ``` while b.childs[j] != entrychild: ``` ``` # Examine the vertices of the sub-blossom to see whether ``` ``` # it is reachable from a neighbouring S-vertex outside the ``` ``` # expanding blossom. ``` ``` bv = b.childs[j] ``` ``` if label.get(bv) == 1: ``` ``` # This sub-blossom just got label S through one of its ``` ``` # neighbours; leave it be. ``` ``` j += jstep ``` ``` continue ``` ``` if isinstance(bv, Blossom): ``` ``` for v in bv.leaves(): ``` ``` if label.get(v): ``` ``` break ``` ``` else: ``` ``` v = bv ``` ``` # If the sub-blossom contains a reachable vertex, assign ``` ``` # label T to the sub-blossom. ``` ``` if label.get(v): ``` ``` assert label[v] == 2 ``` ``` assert inblossom[v] == bv ``` ``` label[v] = None ``` ``` label[mate[blossombase[bv]]] = None ``` ``` assignLabel(v, 2, labeledge[v]) ``` ``` j += jstep ``` ``` # Remove the expanded blossom entirely. ``` ``` label.pop(b, None) ``` ``` labeledge.pop(b, None) ``` ``` bestedge.pop(b, None) ``` ``` del blossomparent[b] ``` ``` del blossombase[b] ``` ``` del blossomdual[b] ``` ``` # Swap matched/unmatched edges over an alternating path through blossom b ``` ``` # between vertex v and the base vertex. Keep blossom bookkeeping ``` ``` # consistent. ``` ``` def augmentBlossom(b, v): ``` ``` # Bubble up through the blossom tree from vertex v to an immediate ``` ``` # sub-blossom of b. ``` ``` t = v ``` ``` while blossomparent[t] != b: ``` ``` t = blossomparent[t] ``` ``` # Recursively deal with the first sub-blossom. ``` ``` if isinstance(t, Blossom): ``` ``` augmentBlossom(t, v) ``` ``` # Decide in which direction we will go round the blossom. ``` ``` i = j = b.childs.index(t) ``` ``` if i & 1: ``` ``` # Start index is odd; go forward and wrap. ``` ``` j -= len(b.childs) ``` ``` jstep = 1 ``` ``` else: ``` ``` # Start index is even; go backward. ``` ``` jstep = -1 ``` ``` # Move along the blossom until we get to the base. ``` ``` while j != 0: ``` ``` # Step to the next sub-blossom and augment it recursively. ``` ``` j += jstep ``` ``` t = b.childs[j] ``` ``` if jstep == 1: ``` ``` w, x = b.edges[j] ``` ``` else: ``` ``` x, w = b.edges[j - 1] ``` ``` if isinstance(t, Blossom): ``` ``` augmentBlossom(t, w) ``` ``` # Step to the next sub-blossom and augment it recursively. ``` ``` j += jstep ``` ``` t = b.childs[j] ``` ``` if isinstance(t, Blossom): ``` ``` augmentBlossom(t, x) ``` ``` # Match the edge connecting those sub-blossoms. ``` ``` mate[w] = x ``` ``` mate[x] = w ``` ``` # Rotate the list of sub-blossoms to put the new base at the front. ``` ``` b.childs = b.childs[i:] + b.childs[:i] ``` ``` b.edges = b.edges[i:] + b.edges[:i] ``` ``` blossombase[b] = blossombase[b.childs] ``` ``` assert blossombase[b] == v ``` ``` # Swap matched/unmatched edges over an alternating path between two ``` ``` # single vertices. The augmenting path runs through S-vertices v and w. ``` ``` def augmentMatching(v, w): ``` ``` for (s, j) in ((v, w), (w, v)): ``` ``` # Match vertex s to vertex j. Then trace back from s ``` ``` # until we find a single vertex, swapping matched and unmatched ``` ``` # edges as we go. ``` ``` while 1: ``` ``` bs = inblossom[s] ``` ``` assert label[bs] == 1 ``` ``` assert ( ``` ``` labeledge[bs] is None and blossombase[bs] not in mate)\ ``` ``` or (labeledge[bs] == mate[blossombase[bs]]) ``` ``` # Augment through the S-blossom from s to base. ``` ``` if isinstance(bs, Blossom): ``` ``` augmentBlossom(bs, s) ``` ``` # Update mate[s] ``` ``` mate[s] = j ``` ``` # Trace one step back. ``` ``` if labeledge[bs] is None: ``` ``` # Reached single vertex; stop. ``` ``` break ``` ``` t = labeledge[bs] ``` ``` bt = inblossom[t] ``` ``` assert label[bt] == 2 ``` ``` # Trace one more step back. ``` ``` s, j = labeledge[bt] ``` ``` # Augment through the T-blossom from j to base. ``` ``` assert blossombase[bt] == t ``` ``` if isinstance(bt, Blossom): ``` ``` augmentBlossom(bt, j) ``` ``` # Update mate[j] ``` ``` mate[j] = s ``` ``` # Verify that the optimum solution has been reached. ``` ``` def verifyOptimum(): ``` ``` if maxcardinality: ``` ``` # Vertices may have negative dual; ``` ``` # find a constant non-negative number to add to all vertex duals. ``` ``` vdualoffset = max(0, -min(dualvar.values())) ``` ``` else: ``` ``` vdualoffset = 0 ``` ``` # 0. all dual variables are non-negative ``` ``` assert min(dualvar.values()) + vdualoffset >= 0 ``` ``` assert len(blossomdual) == 0 or min(blossomdual.values()) >= 0 ``` ``` # 0. all edges have non-negative slack and ``` ``` # 1. all matched edges have zero slack; ``` ``` for i, j, d in G.edges(data=True): ``` ``` wt = d.get(weight, 1) ``` ``` if i == j: ``` ``` continue # ignore self-loops ``` ``` s = dualvar[i] + dualvar[j] - 2 * wt ``` ``` iblossoms = [i] ``` ``` jblossoms = [j] ``` ``` while blossomparent[iblossoms[-1]] is not None: ``` ``` iblossoms.append(blossomparent[iblossoms[-1]]) ``` ``` while blossomparent[jblossoms[-1]] is not None: ``` ``` jblossoms.append(blossomparent[jblossoms[-1]]) ``` ``` iblossoms.reverse() ``` ``` jblossoms.reverse() ``` ``` for (bi, bj) in zip(iblossoms, jblossoms): ``` ``` if bi != bj: ``` ``` break ``` ``` s += 2 * blossomdual[bi] ``` ``` assert s >= 0 ``` ``` if mate.get(i) == j or mate.get(j) == i: ``` ``` assert mate[i] == j and mate[j] == i ``` ``` assert s == 0 ``` ``` # 2. all single vertices have zero dual value; ``` ``` for v in gnodes: ``` ``` assert (v in mate) or dualvar[v] + vdualoffset == 0 ``` ``` # 3. all blossoms with positive dual value are full. ``` ``` for b in blossomdual: ``` ``` if blossomdual[b] > 0: ``` ``` assert len(b.edges) % 2 == 1 ``` ``` for (i, j) in b.edges[1::2]: ``` ``` assert mate[i] == j and mate[j] == i ``` ``` # Ok. ``` ``` # Main loop: continue until no further improvement is possible. ``` ``` while 1: ``` ``` # Each iteration of this loop is a "stage". ``` ``` # A stage finds an augmenting path and uses that to improve ``` ``` # the matching. ``` ``` # Remove labels from top-level blossoms/vertices. ``` ``` label.clear() ``` ``` labeledge.clear() ``` ``` # Forget all about least-slack edges. ``` ``` bestedge.clear() ``` ``` for b in blossomdual: ``` ``` b.mybestedges = None ``` ``` # Loss of labeling means that we can not be sure that currently ``` ``` # allowable edges remain allowable throughout this stage. ``` ``` allowedge.clear() ``` ``` # Make queue empty. ``` ``` queue[:] = [] ``` ``` # Label single blossoms/vertices with S and put them in the queue. ``` ``` for v in gnodes: ``` ``` if (v not in mate) and label.get(inblossom[v]) is None: ``` ``` assignLabel(v, 1, None) ``` ``` # Loop until we succeed in augmenting the matching. ``` ``` augmented = 0 ``` ``` while 1: ``` ``` # Each iteration of this loop is a "substage". ``` ``` # A substage tries to find an augmenting path; ``` ``` # if found, the path is used to improve the matching and ``` ``` # the stage ends. If there is no augmenting path, the ``` ``` # primal-dual method is used to pump some slack out of ``` ``` # the dual variables. ``` ``` # Continue labeling until all vertices which are reachable ``` ``` # through an alternating path have got a label. ``` ``` while queue and not augmented: ``` ``` # Take an S vertex from the queue. ``` ``` v = queue.pop() ``` ``` assert label[inblossom[v]] == 1 ``` ``` # Scan its neighbours: ``` ``` for w in G.neighbors(v): ``` ``` if w == v: ``` ``` continue # ignore self-loops ``` ``` # w is a neighbour to v ``` ``` bv = inblossom[v] ``` ``` bw = inblossom[w] ``` ``` if bv == bw: ``` ``` # this edge is internal to a blossom; ignore it ``` ``` continue ``` ``` if (v, w) not in allowedge: ``` ``` kslack = slack(v, w) ``` ``` if kslack <= 0: ``` ``` # edge k has zero slack => it is allowable ``` ``` allowedge[(v, w)] = allowedge[(w, v)] = True ``` ``` if (v, w) in allowedge: ``` ``` if label.get(bw) is None: ``` ``` # (C1) w is a free vertex; ``` ``` # label w with T and label its mate with S (R12). ``` ``` assignLabel(w, 2, v) ``` ``` elif label.get(bw) == 1: ``` ``` # (C2) w is an S-vertex (not in the same blossom); ``` ``` # follow back-links to discover either an ``` ``` # augmenting path or a new blossom. ``` ``` base = scanBlossom(v, w) ``` ``` if base is not NoNode: ``` ``` # Found a new blossom; add it to the blossom ``` ``` # bookkeeping and turn it into an S-blossom. ``` ``` addBlossom(base, v, w) ``` ``` else: ``` ``` # Found an augmenting path; augment the ``` ``` # matching and end this stage. ``` ``` augmentMatching(v, w) ``` ``` augmented = 1 ``` ``` break ``` ``` elif label.get(w) is None: ``` ``` # w is inside a T-blossom, but w itself has not ``` ``` # yet been reached from outside the blossom; ``` ``` # mark it as reached (we need this to relabel ``` ``` # during T-blossom expansion). ``` ``` assert label[bw] == 2 ``` ``` label[w] = 2 ``` ``` labeledge[w] = (v, w) ``` ``` elif label.get(bw) == 1: ``` ``` # keep track of the least-slack non-allowable edge to ``` ``` # a different S-blossom. ``` ``` if bestedge.get(bv) is None or \ ``` ``` kslack < slack(*bestedge[bv]): ``` ``` bestedge[bv] = (v, w) ``` ``` elif label.get(w) is None: ``` ``` # w is a free vertex (or an unreached vertex inside ``` ``` # a T-blossom) but we can not reach it yet; ``` ``` # keep track of the least-slack edge that reaches w. ``` ``` if bestedge.get(w) is None or \ ``` ``` kslack < slack(*bestedge[w]): ``` ``` bestedge[w] = (v, w) ``` ``` if augmented: ``` ``` break ``` ``` # There is no augmenting path under these constraints; ``` ``` # compute delta and reduce slack in the optimization problem. ``` ``` # (Note that our vertex dual variables, edge slacks and delta's ``` ``` # are pre-multiplied by two.) ``` ``` deltatype = -1 ``` ``` delta = deltaedge = deltablossom = None ``` ``` # Compute delta1: the minimum value of any vertex dual. ``` ``` if not maxcardinality: ``` ``` deltatype = 1 ``` ``` delta = min(dualvar.values()) ``` ``` # Compute delta2: the minimum slack on any edge between ``` ``` # an S-vertex and a free vertex. ``` ``` for v in G.nodes(): ``` ``` if label.get(inblossom[v]) is None and \ ``` ``` bestedge.get(v) is not None: ``` ``` d = slack(*bestedge[v]) ``` ``` if deltatype == -1 or d < delta: ``` ``` delta = d ``` ``` deltatype = 2 ``` ``` deltaedge = bestedge[v] ``` ``` # Compute delta3: half the minimum slack on any edge between ``` ``` # a pair of S-blossoms. ``` ``` for b in blossomparent: ``` ``` if (blossomparent[b] is None and label.get(b) == 1 and ``` ``` bestedge.get(b) is not None): ``` ``` kslack = slack(*bestedge[b]) ``` ``` if allinteger: ``` ``` assert (kslack % 2) == 0 ``` ``` d = kslack // 2 ``` ``` else: ``` ``` d = kslack / 2.0 ``` ``` if deltatype == -1 or d < delta: ``` ``` delta = d ``` ``` deltatype = 3 ``` ``` deltaedge = bestedge[b] ``` ``` # Compute delta4: minimum z variable of any T-blossom. ``` ``` for b in blossomdual: ``` ``` if (blossomparent[b] is None and label.get(b) == 2 and ``` ``` (deltatype == -1 or blossomdual[b] < delta)): ``` ``` delta = blossomdual[b] ``` ``` deltatype = 4 ``` ``` deltablossom = b ``` ``` if deltatype == -1: ``` ``` # No further improvement possible; max-cardinality optimum ``` ``` # reached. Do a final delta update to make the optimum ``` ``` # verifyable. ``` ``` assert maxcardinality ``` ``` deltatype = 1 ``` ``` delta = max(0, min(dualvar.values())) ``` ``` # Update dual variables according to delta. ``` ``` for v in gnodes: ``` ``` if label.get(inblossom[v]) == 1: ``` ``` # S-vertex: 2*u = 2*u - 2*delta ``` ``` dualvar[v] -= delta ``` ``` elif label.get(inblossom[v]) == 2: ``` ``` # T-vertex: 2*u = 2*u + 2*delta ``` ``` dualvar[v] += delta ``` ``` for b in blossomdual: ``` ``` if blossomparent[b] is None: ``` ``` if label.get(b) == 1: ``` ``` # top-level S-blossom: z = z + 2*delta ``` ``` blossomdual[b] += delta ``` ``` elif label.get(b) == 2: ``` ``` # top-level T-blossom: z = z - 2*delta ``` ``` blossomdual[b] -= delta ``` ``` # Take action at the point where minimum delta occurred. ``` ``` if deltatype == 1: ``` ``` # No further improvement possible; optimum reached. ``` ``` break ``` ``` elif deltatype == 2: ``` ``` # Use the least-slack edge to continue the search. ``` ``` (v, w) = deltaedge ``` ``` assert label[inblossom[v]] == 1 ``` ``` allowedge[(v, w)] = allowedge[(w, v)] = True ``` ``` queue.append(v) ``` ``` elif deltatype == 3: ``` ``` # Use the least-slack edge to continue the search. ``` ``` (v, w) = deltaedge ``` ``` allowedge[(v, w)] = allowedge[(w, v)] = True ``` ``` assert label[inblossom[v]] == 1 ``` ``` queue.append(v) ``` ``` elif deltatype == 4: ``` ``` # Expand the least-z blossom. ``` ``` expandBlossom(deltablossom, False) ``` ``` # End of a this substage. ``` ``` # Paranoia check that the matching is symmetric. ``` ``` for v in mate: ``` ``` assert mate[mate[v]] == v ``` ``` # Stop when no more augmenting path can be found. ``` ``` if not augmented: ``` ``` break ``` ``` # End of a stage; expand all S-blossoms which have zero dual. ``` ``` for b in list(blossomdual.keys()): ``` ``` if b not in blossomdual: ``` ``` continue # already expanded ``` ``` if (blossomparent[b] is None and label.get(b) == 1 and ``` ``` blossomdual[b] == 0): ``` ``` expandBlossom(b, True) ``` ``` # Verify that we reached the optimum solution (only for integer weights). ``` ``` if allinteger: ``` ``` verifyOptimum() ``` ``` return matching_dict_to_set(mate) ```