diff --git a/python/max_weight_matching.py b/python/max_weight_matching.py index aa0bab5..e829cac 100644 --- a/python/max_weight_matching.py +++ b/python/max_weight_matching.py @@ -299,8 +299,8 @@ class _GraphInfo: # Each vertex is incident to zero or more edges. # - # "adjacent_edges[v]" is the list of edge indices of edges incident - # to the vertex with index "v". + # "adjacent_edges[x]" is the list of edge indices of edges incident + # to the vertex with index "x". # # These data remain unchanged while the algorithm runs. self.adjacent_edges: list[list[int]] = [ @@ -322,7 +322,7 @@ _LABEL_T = 2 class _Blossom: - """Represents a non-trivial blossom in a (partially) matched graph. + """Represents a blossom in a partially matched graph. A blossom is an odd-length alternating cycle over sub-blossoms. An alternating path consists of alternating matched and unmatched edges. @@ -330,30 +330,92 @@ class _Blossom: the same sub-blossom. A single vertex by itself is also a blossom: a "trivial blossom". - We use the term "non-trivial blossom" to refer to a blossom that - contains at least 3 sub-blossoms. + + An instance of this class represents either a trivial blossom, + or a non-trivial blossom which consists of multiple sub-blossoms. Blossoms are recursive structures: A non-trivial blossoms contains sub-blossoms, which may themselves contain sub-blossoms etc. Each blossom contains exactly one vertex that is not matched to another vertex in the same blossom. This is the "base vertex" of the blossom. + """ - Blossoms are created and destroyed by the matching algorithm. + def __init__(self, base_vertex: int) -> None: + """Initialize a new blossom.""" + + # If this is not a top-level blossom, + # "parent" is the blossom in which this blossom is a sub-blossom. + # + # If this is a top-level blossom, + # "parent = None". + self.parent: Optional[_NonTrivialBlossom] = None + + # "base_vertex" is the vertex index of the base of the blossom. + # This is the unique vertex which is contained in the blossom + # but not matched to another vertex in the same blossom. + # + # For trivial blossoms, the base vertex is simply the only + # vertex in the blossom. + self.base_vertex: int = base_vertex + + # A top-level blossom that are part of an alternating tree, + # has labels S or T. Unlabeled top-level blossoms are not (yet) + # part of any alternating tree. + self.label: int = _LABEL_NONE + + # Labeled top-level blossoms keep track of the edge through which + # they are attached to an alternating tree. + # + # "tree_edge = (x, y)" if the blossom is attached to an alternating + # tree via edge "(x, y)" and vertex "y" is contained in the blossom. + # + # "tree_edge = None" if the blossom is the root of an alternating tree. + self.tree_edge: Optional[tuple[int, int]] = None + + # For a top-level S-blossom, + # "best_edge" is the edge index of the least-slack edge to + # a different S-blossom, or -1 if no such edge has been found. + self.best_edge: int = -1 + + # "marker" is a temporary variable used to discover common + # ancestors in the blossom tree. It is normally False, except + # when used by "trace_alternating_paths()". + self.marker: bool = False + + def vertices(self) -> list[int]: + """Return a list of vertex indices contained in the blossom.""" + return [self.base_vertex] + + +class _NonTrivialBlossom(_Blossom): + """Represents a non-trivial blossom in a partially matched graph. + + A non-trivial blossom is a blossom that contains multiple sub-blossoms + (at least 3 sub-blossoms, since all blossoms have odd length). + + Non-trivial blossoms maintain a list of their sub-blossoms and the edges + between their subblossoms. + + Unlike trivial blossoms, each non-trivial blossom is associated with + a variable in the dual LPP problem. + + Non-trivial blossoms are created and destroyed by the matching algorithm. This implies that not every odd-length alternating cycle is a blossom; it only becomes a blossom through an explicit action of the algorithm. - An existing blossom may be changed when the matching is augmented - along a path that runs through the blossom. + An existing blossom may change when the matching is augmented along + a path that runs through the blossom. """ def __init__( self, - subblossoms: list[int], - edges: list[tuple[int, int]], - base_vertex: int + subblossoms: list[_Blossom], + edges: list[tuple[int, int]] ) -> None: """Initialize a new blossom.""" + super().__init__(subblossoms[0].base_vertex) + # Sanity check. n = len(subblossoms) assert len(edges) == n @@ -365,10 +427,7 @@ class _Blossom: # # "subblossoms[0]" is the start and end of the alternating cycle. # "subblossoms[0]" contains the base vertex of the blossom. - # - # "subblossoms[i]" is blossom index (either a non-trivial blossom - # index, or a vertex index if the i-th sub-blossom is trivial). - self.subblossoms: list[int] = subblossoms + self.subblossoms: list[_Blossom] = subblossoms # "edges" is a list of edges linking the sub-blossoms. # Each edge is represented as an ordered pair "(x, y)" where "x" @@ -378,323 +437,46 @@ class _Blossom: # adjacent to vertex "y" in "subblossoms[1]", etc. self.edges: list[tuple[int, int]] = edges - # "base_vertex" is the vertex index of the base of the blossom. - # This is the unique vertex which is contained in the blossom - # but not matched to another vertex in the same blossom. - self.base_vertex: int = base_vertex - - # Every blossom has a variable in the dual LPP. + # Every non-trivial blossom has a variable in the dual LPP. # # "dual_var" is the current value of the dual variable. # New blossoms start with dual variable 0. self.dual_var: int|float = 0 + # For a non-trivial, top-level S-blossom, + # "best_edge_set" is a list of least-slack edges between this blossom + # and other S-blossoms. + self.best_edge_set: Optional[list[int]] = None + + def vertices(self) -> list[int]: + """Return a list of vertex indices contained in the blossom.""" + + # Use an explicit stack to avoid deep recursion. + stack: list[_NonTrivialBlossom] = [self] + nodes: list[int] = [] + + while stack: + b = stack.pop() + for sub in b.subblossoms: + if isinstance(sub, _NonTrivialBlossom): + stack.append(sub) + else: + nodes.append(sub.base_vertex) + + return nodes + class _AlternatingPath(NamedTuple): - """Represents an alternating path or an alternating cycle.""" + """Represents a list of edges forming an alternating path or an + alternating cycle.""" edges: list[tuple[int, int]] -class _VertexBestEdgeTracking: - """Keep track of the least-slack edge from an unlabeled or T-labeled - vertex to any S-vertex. - - The least-slack edge between any S-vertex and unlabeled vertex is needed - to calculate a delta step in the matching algorithm. - - For each unlabeled vertex and each T-vertex, this class keeps track - of the least-slack edge to any S-vertex. Tracking for unlabeled vertices - is done to provide the least-slack edge for the delta step. - Tracking for T-vertices is done because such vertices can turn into - unlabeled vertices if they are part of a T-blossom that gets expanded. - - Note: For a given unlabeled vertex or T-vertex, the identity of the - least-slack edge to any S-blossom remains unchanged during a delta step. - Although the delta step changes edge slack, it changes the slack of - every edge to an S-vertex by the same amount. Therefore the edge that has - least-slack before the delta step, will still have least-slack after - the delta step. - """ - - def __init__(self, num_vertex: int) -> None: - """Initialize least-slack edge tracking.""" - - self.num_vertex = num_vertex - - # "vertex_best_edge[x]" is the edge index of the least-slack edge - # between "x" and any S-vertex, or -1 if no such edge has been found. - self.vertex_best_edge: list[int] = num_vertex * [-1] - - def reset(self) -> None: - """Remove all tracked edges. - - This function takes time O(n). - """ - for x in range(self.num_vertex): - self.vertex_best_edge[x] = -1 - - def add_edge(self, - ctx: _MatchingContext, - y: int, - e: int, - slack: int|float - ) -> None: - """Add edge "e" from an S-vertex to unlabeled vertex or T-vertex "y". - - This function takes time O(1) per call. - This function is called O(m) times per stage. - """ - best_edge = self.vertex_best_edge[y] - if best_edge == -1: - self.vertex_best_edge[y] = e - else: - best_slack = ctx.edge_slack_2x(best_edge) - if slack < best_slack: - self.vertex_best_edge[y] = e - - def get_least_slack_edge( - self, - ctx: _MatchingContext - ) -> tuple[int, int|float]: - """Return the index and slack of the least-slack edge between - any S-vertex and unlabeled vertex. - - This function takes time O(n) per call. - This function takes total time O(n**2) per stage. - - Returns: - Tuple (edge_index, slack_2x) if there is a least-slack edge, - or (-1, 0) if there is no suitable edge. - """ - best_index = -1 - best_slack: int|float = 0 - - for x in range(self.num_vertex): - bx = ctx.vertex_blossom[x] - if ctx.blossom_label[bx] == _LABEL_NONE: - e = self.vertex_best_edge[x] - if e != -1: - slack = ctx.edge_slack_2x(e) - if (best_index == -1) or (slack < best_slack): - best_index = e - best_slack = slack - - return (best_index, best_slack) - - -class _BlossomBestEdgeTracking: - """Keep track of the least-slack edge between top-level S-blossoms. - - The least-slack edge between any two top-level S-blossoms is needed - to calculate a delta step in the matching algorithm. - - For each top-level S-blossom, this class keeps track of the least-slack - edge to any S-vertex not in the same blossom. - - Furthermore, for each top-level S-blossom, this class keeps a list of - least-slack edges to other top-level S-blossoms. For any pair of top-level - S-blossoms, the least-slack edge between them is contained in the edge - list of at least one of the blossoms. An edge list may contain multiple - edges to the same S-blossom. Such redundant edges are pruned during - blossom formation to limit the number of tracked edges. - - Note: For a given S-blossom, the identity of the least-slack edge to other - S-blossoms remains unchanged during a delta step. Although the delta step - changes edge slack, it changes the slack of every edge between S-blossoms - by the same amount. Therefore the edge that has least-slack before the - delta step, will still have least-slack after the delta step. - """ - - def __init__(self, num_vertex: int) -> None: - """Initialize least-slack edge tracking.""" - - self.num_vertex = num_vertex - - # For non-trivial top-level S-blossom "b", - # "blossom_best_edge_set[b]" is a list of edges between blossom "b" - # and other S-blossoms. - self.blossom_best_edge_set: list[Optional[list[int]]] = [ - None for b in range(2 * num_vertex)] - - # For every top-level S-blossom "b", - # "blossom_best_edge[b]" is the edge index of the least-slack edge - # to a different S-blossom, or -1 if no such edge has been found. - self.blossom_best_edge: list[int] = (2 * num_vertex) * [-1] - - def reset(self) -> None: - """Remove all tracked edges. - - This function takes time O(n). - """ - for b in range(2 * self.num_vertex): - self.blossom_best_edge_set[b] = None - self.blossom_best_edge[b] = -1 - - def new_blossom(self, b: int) -> None: - """Start tracking edges for the new blossom "b".""" - assert self.blossom_best_edge_set[b] is None - self.blossom_best_edge_set[b] = [] - - def add_edge( - self, - ctx: _MatchingContext, - b: int, - e: int, - slack: int|float - ) -> None: - """Add edge "e" from S-blossom "b" to another S-blossom. - - This function takes time O(1) per call. - This function is called O(m) times per stage. - """ - # Update the least-slack edge from blossom "b" to any other - # S-blossom. - best_edge = self.blossom_best_edge[b] - if best_edge == -1: - self.blossom_best_edge[b] = e - else: - best_slack = ctx.edge_slack_2x(best_edge) - if slack < best_slack: - self.blossom_best_edge[b] = e - - # Regardless of whether this edge is currently the least-slack - # edge from blossom "b" to any other S-blossom, this edge may - # later become the least-slack edge if other edges become - # unavailable during a merge of S-blossoms. - # - # We therefore add the edge to a list of potential future least-slack - # edges for blossom "b". We do this only for non-trivial blossoms. - if b >= self.num_vertex: - best_edge_set = self.blossom_best_edge_set[b] - assert best_edge_set is not None - best_edge_set.append(e) - - def merge_blossoms( - self, - ctx: _MatchingContext, - b: int, - subblossoms: list[int] - ) -> None: - """Merge one or more S-sub-blossoms to form a new S-blossom "b". - - This function takes time O(n) per call. - This function takes total time O(n**2) per stage. - """ - num_vertex = self.num_vertex - - # Calculate the set of least-slack edges to other S-blossoms. - # We basically merge the edge lists from all sub-blossoms, but reject - # edges that are internal to this blossom, and trim the set such that - # there is at most one edge to each external S-blossom. - # - # Sub-blossoms that were formerly labeled T can be ignored; their - # vertices are in the queue and will discover neighbouring S-blossoms - # via the edge scan process. - # - # Build a temporary array holding the least-slack edge index to - # each top-level S-blossom. - best_edge_to_blossom: list[int] = (2 * num_vertex) * [-1] - zero_slack: int|float = 0 - best_slack_to_blossom: list[int|float] = ( - (2 * num_vertex) * [zero_slack]) - - # And find the overall least-slack edge to any other S-blossom. - best_edge = -1 - best_slack: int|float = 0 - - # Add the least-slack edges of every S-sub-blossom. - for sub in subblossoms: - if sub < num_vertex: - # Trivial blossoms don't have a list of least-slack edges, - # so we just look at all adjacent edges. This happens at most - # once per vertex per stage. - # It adds up to O(m) time per stage. - sub_edge_set = ctx.graph.adjacent_edges[sub] - else: - # Use the edge list of the sub-blossom. - sub_edge_set_opt = self.blossom_best_edge_set[sub] - assert sub_edge_set_opt is not None - sub_edge_set = sub_edge_set_opt - # Delete the edge list from the sub-blossom. - self.blossom_best_edge_set[sub] = None - - # Add edges to the temporary array. - for e in sub_edge_set: - (x, y, _w) = ctx.graph.edges[e] - bx = ctx.vertex_blossom[x] - by = ctx.vertex_blossom[y] - assert (bx == b) or (by == b) - - # Reject internal edges in this blossom. - if bx == by: - continue - - # Set bx = other blossom which is reachable through this edge. -# TODO : generalize over this pattern - bx = by if (bx == b) else bx - - # Reject edges that don't link to an S-blossom. - if ctx.blossom_label[bx] != _LABEL_S: - continue - - # Keep only the least-slack edge to "bx". - slack = ctx.edge_slack_2x(e) - if ((best_edge_to_blossom[bx] == -1) - or (slack < best_slack_to_blossom[bx])): - best_edge_to_blossom[bx] = e - best_slack_to_blossom[bx] = slack - - # Update the overall least-slack edge to any S-blossom. - if (best_edge == -1) or (slack < best_slack): - best_edge = e - best_slack = slack - - # Extract a compact list of least-slack edge indices. - # We can not keep the temporary array because that would blow up - # memory use to O(n**2). - best_edge_set = [e for e in best_edge_to_blossom if e != -1] - self.blossom_best_edge_set[b] = best_edge_set - - # Keep the overall least-slack edge. - self.blossom_best_edge[b] = best_edge - - def get_least_slack_edge( - self, - ctx: _MatchingContext - ) -> tuple[int, int|float]: - """Return the index and slack of the least-slack edge between - any pair of top-level S-blossoms. - - This function takes time O(n) per call. - This function takes total time O(n**2) per stage. - - Returns: - Tuple (edge_index, slack_2x) if there is a least-slack edge, - or (-1, 0) if there is no suitable edge. - """ - best_index = -1 - best_slack: int|float = 0 - - for b in range(2 * self.num_vertex): - if ((ctx.blossom_label[b] == _LABEL_S) - and (ctx.blossom_parent[b] == -1)): - e = self.blossom_best_edge[b] - if e != -1: - slack = ctx.edge_slack_2x(e) - if (best_index == -1) or (slack < best_slack): - best_index = e - best_slack = slack - - return (best_index, best_slack) - - class _MatchingContext: """Holds all data used by the matching algorithm. - It contains a partial solution of the matching problem, - as well as several auxiliary data structures. - - These data change while the algorithm runs. + It contains a partial solution of the matching problem and several + auxiliary data structures. """ def __init__(self, graph: _GraphInfo) -> None: @@ -711,48 +493,35 @@ class _MatchingContext: # # If vertex "x" is matched to vertex "y", # "vertex_mate[x] == y" and "vertex_mate[y] == x". + # # If vertex "x" is unmatched, "vertex_mate[x] == -1". # # Initially all vertices are unmatched. self.vertex_mate: list[int] = num_vertex * [-1] - # Blossoms are indexed by integers in range 0 .. 2*n-1. + # Each vertex is associated with a trivial blossom. + # In addition, non-trivial blossoms may be created and destroyed + # during the course of the matching algorithm. # - # Blossom indices in range 0 .. n-1 refer to the trivial blossoms - # that consist of a single vertex. In this case the blossom index - # is simply equal to the vertex index. - # - # Blossom indices in range n .. 2*n-1 refer to non-trivial blossoms, - # represented by instances of the _Blossom class. - # - # "blossom[b]" (for n <= b < 2*n-1) is an instance of the _Blossom - # class that describes the non-trivial blossom with index "b". - # - # Blossoms are created and destroyed while the algorithm runs. - # Initially there are no blossoms. - self.blossom: list[Optional[_Blossom]] = [ - None for b in range(2 * num_vertex)] + # "trivial_blossom[x]" is the trivial blossom that contains only + # vertex "x". + self.trivial_blossom: list[_Blossom] = [_Blossom(x) + for x in range(num_vertex)] - # List of currently unused blossom indices. - self.unused_blossoms: list[int] = list( - reversed(range(num_vertex, 2 * num_vertex))) - - # Every vertex is part of exactly one top-level blossom, - # possibly a trivial blossom consisting of just that vertex. + # Non-trivial blossoms may be created and destroyed during + # the course of the algorithm. # - # "vertex_blossom[x]" is the index of the top-level blossom that - # contains vertex "x". - # "vertex_blossom[x] == x" if the "x" is a trivial top-level blossom. - # - # Initially all vertices are top-level trivial blossoms. - self.vertex_blossom: list[int] = list(range(num_vertex)) + # Initially there are no non-trivial blossoms. + self.nontrivial_blossom: list[_NonTrivialBlossom] = [] - # "blossom_parent[b]" is the index of the smallest blossom that - # is a strict superset of blossom "b", or - # "blossom_parent[b] == -1" if blossom "b" is a top-level blossom. + # Every vertex is contained in exactly one top-level blossom + # (possibly the trivial blossom that contains just that vertex). + # + # "vertex_top_blossom[x]" is the top-level blossom that contains + # vertex "x". # # Initially all vertices are trivial top-level blossoms. - self.blossom_parent: list[int] = (2 * num_vertex) * [-1] + self.vertex_top_blossom: list[_Blossom] = self.trivial_blossom.copy() # Every vertex has a variable in the dual LPP. # @@ -764,40 +533,15 @@ class _MatchingContext: max_weight = max(w for (_x, _y, w) in graph.edges) self.vertex_dual_2x: list[int|float] = num_vertex * [max_weight] - # Top-level blossoms that are part of an alternating tree are - # labeled S or T. Unlabeled top-level blossoms are not (yet) - # part of any alternating tree. - # - # "blossom_label[b]" is the label of blossom "b". - # - # At the beginning of a stage, all blossoms are unlabeled. - self.blossom_label: list[int] = (2 * num_vertex) * [_LABEL_NONE] + # For each T-vertex or unlabeled vertex "x", + # "vertex_best_edge[x]" is the edge index of the least-slack edge + # between "x" and any S-vertex, or -1 if no such edge has been found. + self.vertex_best_edge: list[int] = num_vertex * [-1] - # For each labeled blossom, we keep track of the edge that attaches - # it to its alternating tree. - # - # "blossom_link[b] = (x, y)" denotes the edge through which - # blossom "b" is attached to the alternating tree, where "x" and "y" - # are vertex indices and vertex "y" is contained in blossom "b". - # - # "blossom_link[b] = None" if "b" is the root of an alternating tree, - # or if "b" is not a labeled, top-level blossom. - self.blossom_link: list[Optional[tuple[int, int]]] = [ - None for b in range(2 * num_vertex)] - - # Track least-slack edges between S-vertex and unlabeled vertex. - self.vertex_best_edges = _VertexBestEdgeTracking(num_vertex) - - # Track least-slack edges between S-blossoms. - self.blossom_best_edges = _BlossomBestEdgeTracking(num_vertex) - - # A list of S-vertices to be scanned. + # "queue" is a list of S-vertices that must be scanned. # We call it a queue, but it is actually a stack. self.queue: list[int] = [] - # Temporary array used to construct new blossoms. - self.blossom_marker: list[bool] = (2 * num_vertex) * [False] - def edge_slack_2x(self, e: int) -> int|float: """Return 2 times the slack of the edge with index "e". @@ -810,32 +554,250 @@ class _MatchingContext: This function is called O(m) times per stage. """ (x, y, w) = self.graph.edges[e] - assert self.vertex_blossom[x] != self.vertex_blossom[y] + assert self.vertex_top_blossom[x] is not self.vertex_top_blossom[y] return self.vertex_dual_2x[x] + self.vertex_dual_2x[y] - 2 * w - def get_blossom(self, b: int) -> _Blossom: - """Return the Blossom instance for blossom index "b".""" - blossom = self.blossom[b] - assert blossom is not None - return blossom + # + # Least-slack edge tracking: + # + # To calculate delta steps, the matching algorithm needs to find + # - the least-slack edge between any S-vertex and an unlabeled vertex; + # - the least-slack edge between any pair of top-level S-blossoms. + # + # For each unlabeled vertex and each T-vertex, we keep track of the + # least-slack edge to any S-vertex. Tracking for unlabeled vertices + # serves to provide the least-slack edge for the delta step. + # Tracking for T-vertices is done because such vertices can turn into + # unlabeled vertices if they are part of a T-blossom that gets expanded. + # + # For each top-level S-blossom, we keep track of the least-slack edge + # to any S-vertex not in the same blossom. + # + # Furthermore, for each top-level S-blossom, we keep a list of least-slack + # edges to other top-level S-blossoms. For any pair of top-level + # S-blossoms, the least-slack edge between them is contained in the edge + # list of at least one of the blossoms. An edge list may contain multiple + # edges to the same S-blossom. Such redundant edges are pruned during + # blossom merging to limit the number of tracked edges. + # + # Note: For a given vertex or blossom, the identity of the least-slack + # edge to any S-blossom remains unchanged during a delta step. + # Although the delta step changes edge slacks, it changes the slack + # of every edge to an S-vertex by the same amount. Therefore the edge + # that had least slack before the delta step, will still have least slack + # after the delta step. + # - def blossom_vertices(self, b: int) -> list[int]: - """Return a list of vertex indices contained in blossom "b".""" + def lset_reset(self) -> None: + """Reset least-slack edge tracking. + + This function takes time O(n). + """ num_vertex = self.graph.num_vertex - if b < num_vertex: - return [b] + + for x in range(num_vertex): + self.vertex_best_edge[x] = -1 + + for blossom in self.trivial_blossom: + blossom.best_edge = -1 + + for blossom in self.nontrivial_blossom: + blossom.best_edge = -1 + blossom.best_edge_set = None + + def lset_add_vertex_edge(self, y: int, e: int, slack: int|float) -> None: + """Add edge "e" from an S-vertex to unlabeled vertex or T-vertex "y". + + This function takes time O(1) per call. + This function is called O(m) times per stage. + """ + best_edge = self.vertex_best_edge[y] + if best_edge == -1: + self.vertex_best_edge[y] = e else: - # Use an explicit stack to avoid deep recursion. - stack: list[int] = [b] - nodes: list[int] = [] - while stack: - b = stack.pop() - for sub in self.get_blossom(b).subblossoms: - if sub < num_vertex: - nodes.append(sub) - else: - stack.append(sub) - return nodes + best_slack = self.edge_slack_2x(best_edge) + if slack < best_slack: + self.vertex_best_edge[y] = e + + def lset_get_best_vertex_edge(self) -> tuple[int, int|float]: + """Return the index and slack of the least-slack edge between + any S-vertex and unlabeled vertex. + + This function takes time O(n) per call. + This function takes total time O(n**2) per stage. + + Returns: + Tuple (edge_index, slack_2x) if there is a least-slack edge, + or (-1, 0) if there is no suitable edge. + """ + best_index = -1 + best_slack: int|float = 0 + + for x in range(self.graph.num_vertex): + if self.vertex_top_blossom[x].label == _LABEL_NONE: + e = self.vertex_best_edge[x] + if e != -1: + slack = self.edge_slack_2x(e) + if (best_index == -1) or (slack < best_slack): + best_index = e + best_slack = slack + + return (best_index, best_slack) + + def lset_new_blossom(self, blossom: _Blossom) -> None: + """Start tracking edges for a new S-blossom.""" + assert blossom.best_edge == -1 + if isinstance(blossom, _NonTrivialBlossom): + assert blossom.best_edge_set is None + blossom.best_edge_set = [] + + def lset_add_blossom_edge( + self, + blossom: _Blossom, + e: int, + slack: int|float + ) -> None: + """Add edge "e" between the specified S-blossom and another S-blossom. + + This function takes time O(1) per call. + This function is called O(m) times per stage. + """ + # Track least-slack edge between this blossom and any other S-blossom. + if blossom.best_edge == -1: + blossom.best_edge = e + else: + best_slack = self.edge_slack_2x(blossom.best_edge) + if slack < best_slack: + blossom.best_edge = e + + # Regardless of whether this edge is currently the least-slack edge, + # this edge may later become the least-slack edge if other edges + # become unavailable when S-blossoms are merged. + # + # We therefore add the edge to a list of potential future least-slack + # edges for this blossom. We do this only for non-trivial blossoms. + if isinstance(blossom, _NonTrivialBlossom): + assert blossom.best_edge_set is not None + blossom.best_edge_set.append(e) + + def lset_merge_blossoms(self, blossom: _NonTrivialBlossom) -> None: + """Update least-slack edge tracking after merging sub-blossoms + into a new S-blossom. + + This function takes time O(n) per call. + This function takes total time O(n**2) per stage. + """ + num_vertex = self.graph.num_vertex + + # Calculate the set of least-slack edges to other S-blossoms. + # We basically merge the edge lists from all sub-blossoms, but reject + # edges that are internal to this blossom, and trim the set such that + # there is at most one edge to each external S-blossom. + # + # Sub-blossoms that were formerly labeled T can be ignored; their + # vertices are in the queue and will discover neighbouring S-blossoms + # via the edge scan process. + # + # Build a temporary array holding the least-slack edge index to + # each top-level S-blossom. This array is indexed by the base vertex + # of the blossoms. +# TODO : consider using pre-allocated arrays for this purpose + best_edge_to_blossom: list[int] = num_vertex * [-1] + zero_slack: int|float = 0 + best_slack_to_blossom: list[int|float] = num_vertex * [zero_slack] + + # And find the overall least-slack edge to any other S-blossom. + best_edge = -1 + best_slack: int|float = 0 + + # Add the least-slack edges of every S-sub-blossom. + for sub in blossom.subblossoms: + + if sub.label != _LABEL_S: + continue + + if isinstance(sub, _NonTrivialBlossom): + # Pull the edge list from the sub-blossom. + assert sub.best_edge_set is not None + sub_edge_set = sub.best_edge_set + # Delete the edge list from the sub-blossom. + sub.best_edge_set = None + else: + # Trivial blossoms don't have a list of least-slack edges, + # so we just look at all adjacent edges. This happens at most + # once per vertex per stage. + # It adds up to O(m) time per stage. + sub_edge_set = self.graph.adjacent_edges[sub.base_vertex] + + # Add edges to the temporary array. + for e in sub_edge_set: + (x, y, _w) = self.graph.edges[e] + bx = self.vertex_top_blossom[x] + by = self.vertex_top_blossom[y] + assert (bx is blossom) or (by is blossom) + + # Reject internal edges in this blossom. + if bx is by: + continue + + # Set bx = blossom at the other end of this edge. + bx = by if (bx is blossom) else bx + + # Reject edges that don't link to an S-blossom. + if bx.label != _LABEL_S: + continue + + # Keep only the least-slack edge to "bx". + slack = self.edge_slack_2x(e) + bx_base = bx.base_vertex + if ((best_edge_to_blossom[bx_base] == -1) + or (slack < best_slack_to_blossom[bx_base])): + best_edge_to_blossom[bx_base] = e + best_slack_to_blossom[bx_base] = slack + + # Update the overall least-slack edge to any S-blossom. + if (best_edge == -1) or (slack < best_slack): + best_edge = e + best_slack = slack + + # Extract a compact list of least-slack edge indices. + # We can not keep the temporary array because that would blow up + # memory use to O(n**2). + best_edge_set = [e for e in best_edge_to_blossom if e != -1] + blossom.best_edge_set = best_edge_set + + # Keep the overall least-slack edge. + blossom.best_edge = best_edge + + def lset_get_best_blossom_edge(self) -> tuple[int, int|float]: + """Return the index and slack of the least-slack edge between + any pair of top-level S-blossoms. + + This function takes time O(n) per call. + This function takes total time O(n**2) per stage. + + Returns: + Tuple (edge_index, slack_2x) if there is a least-slack edge, + or (-1, 0) if there is no suitable edge. + """ + best_index = -1 + best_slack: int|float = 0 + +# TODO : do we really want to split trivial/nontrivial blossoms ? + for blossom in self.trivial_blossom + self.nontrivial_blossom: + if (blossom.label == _LABEL_S) and (blossom.parent is None): + e = blossom.best_edge + if e != -1: + slack = self.edge_slack_2x(e) + if (best_index == -1) or (slack < best_slack): + best_index = e + best_slack = slack + + return (best_index, best_slack) + + # + # General support routines: + # def reset_stage(self) -> None: """Reset data which are only valid during a stage. @@ -844,19 +806,16 @@ class _MatchingContext: and resets tracking of least-slack edges. """ - num_vertex = self.graph.num_vertex - # Remove blossom labels. - for b in range(2 * num_vertex): - self.blossom_label[b] = _LABEL_NONE - self.blossom_link[b] = None + for blossom in self.trivial_blossom + self.nontrivial_blossom: + blossom.label = _LABEL_NONE + blossom.tree_edge = None # Clear the queue. self.queue.clear() # Reset least-slack edge tracking. - self.vertex_best_edges.reset() - self.blossom_best_edges.reset() + self.lset_reset() def trace_alternating_paths(self, x: int, y: int) -> _AlternatingPath: """Trace back through the alternating trees from vertices "x" and "y". @@ -877,8 +836,7 @@ class _MatchingContext: blossoms. """ - blossom_marker = self.blossom_marker - marked_blossoms: list[int] = [] + marked_blossoms: list[_Blossom] = [] # "xedges" is a list of edges used while tracing from "x". # "yedges" is a list of edges used while tracing from "y". @@ -893,27 +851,26 @@ class _MatchingContext: # This ensures that the search time is bounded by the size of the # newly found blossom. # TODO : this code is a bit shady; maybe reconsider the swapping trick - first_common = -1 + first_common: Optional[_Blossom] = None while x != -1 or y != -1: # Check if we found a common ancestor. - bx = self.vertex_blossom[x] - if blossom_marker[bx]: + bx = self.vertex_top_blossom[x] + if bx.marker: first_common = bx break # Mark blossom as a potential common ancestor. - blossom_marker[bx] = True + bx.marker = True marked_blossoms.append(bx) # Track back through the link in the alternating tree. - link = self.blossom_link[bx] - if link is None: + if bx.tree_edge is None: # Reached the root of this alternating tree. x = -1 else: - xedges.append(link) - x = link[0] + xedges.append(bx.tree_edge) + x = bx.tree_edge[0] # Swap "x" and "y" to alternate between paths. if y != -1: @@ -922,14 +879,14 @@ class _MatchingContext: # Remove all markers we placed. for b in marked_blossoms: - blossom_marker[b] = False + b.marker = False # If we found a common ancestor, trim the paths so they end there. # TODO : also this is just plain ugly - try to rework - if first_common != -1: - assert self.vertex_blossom[xedges[-1][0]] == first_common - while (yedges - and (self.vertex_blossom[yedges[-1][0]] != first_common)): + if first_common is not None: + assert self.vertex_top_blossom[xedges[-1][0]] is first_common + while yedges and (self.vertex_top_blossom[yedges[-1][0]] + is not first_common): yedges.pop() # Fuse the two paths. @@ -942,6 +899,10 @@ class _MatchingContext: return _AlternatingPath(path_edges) + # + # Merge and expand blossoms: + # + def make_blossom(self, path: _AlternatingPath) -> None: """Create a new blossom from an alternating cycle. @@ -952,86 +913,67 @@ class _MatchingContext: This function takes total time O(n**2) per stage. """ - num_vertex = self.graph.num_vertex - # Check that the path is odd-length. assert len(path.edges) % 2 == 1 assert len(path.edges) >= 3 # Construct the list of sub-blossoms (current top-level blossoms). - subblossoms = [self.vertex_blossom[x] for (x, y) in path.edges] + subblossoms = [self.vertex_top_blossom[x] for (x, y) in path.edges] # Check that the path is cyclic. # Note the path may not start and end with the same _vertex_, # but it must start and end in the same _blossom_. - subblossoms_next = [self.vertex_blossom[y] for (x, y) in path.edges] + subblossoms_next = [self.vertex_top_blossom[y] + for (x, y) in path.edges] assert subblossoms[0] == subblossoms_next[-1] assert subblossoms[1:] == subblossoms_next[:-1] - # Determine the base vertex of the new blossom. - base_blossom = subblossoms[0] - if base_blossom >= num_vertex: - base_vertex = self.get_blossom(base_blossom).base_vertex - else: - base_vertex = base_blossom - # Create the new blossom object. - blossom = _Blossom(subblossoms, path.edges, base_vertex) + blossom = _NonTrivialBlossom(subblossoms, path.edges) - # Allocate a new blossom index and create the blossom object. - b = self.unused_blossoms.pop() - self.blossom[b] = blossom - self.blossom_parent[b] = -1 + # Insert into the blossom array. +# TODO : rework this + self.nontrivial_blossom.append(blossom) # Link the subblossoms to the their new parent. for sub in subblossoms: - self.blossom_parent[sub] = b + sub.parent = blossom # Update blossom-membership of all vertices in the new blossom. # NOTE: This step takes O(n) time per blossom formation, and adds up # to O(n**2) total time per stage. # This could be improved through a union-find datastructure, or # by re-using the blossom index of the largest sub-blossom. - for x in self.blossom_vertices(b): - self.vertex_blossom[x] = b + for x in blossom.vertices(): + self.vertex_top_blossom[x] = blossom # Assign label S to the new blossom. - assert self.blossom_label[base_blossom] == _LABEL_S - self.blossom_label[b] = _LABEL_S - self.blossom_link[b] = self.blossom_link[base_blossom] + assert subblossoms[0].label == _LABEL_S + blossom.label = _LABEL_S + blossom.tree_edge = subblossoms[0].tree_edge # Former T-vertices which are part of this blossom now become # S-vertices. Add them to the queue. for sub in subblossoms: - if self.blossom_label[sub] == _LABEL_T: - if sub < num_vertex: - self.queue.append(sub) - else: - self.queue.extend(self.blossom_vertices(sub)) + if sub.label == _LABEL_T: + self.queue.extend(sub.vertices()) - # Merge least-slack edges for the S-sub-blossoms and transfer them - # to the new blossom index. - subblossoms_s = [ - sub for sub in subblossoms if self.blossom_label[sub] == _LABEL_S] - self.blossom_best_edges.merge_blossoms(self, b, subblossoms_s) + # Merge least-slack edges for the S-sub-blossoms. + self.lset_merge_blossoms(blossom) def find_path_through_blossom( self, - b: int, - sub: int - ) -> tuple[list[int], list[tuple[int, int]]]: - """Construct a path through blossom "b" from sub-blossom "sub" - to the base of the blossom. + blossom: _NonTrivialBlossom, + sub: _Blossom + ) -> tuple[list[_Blossom], list[tuple[int, int]]]: + """Construct a path through the specified blossom, + from sub-blossom "sub" to the base of the blossom. Return: Tuple (nodes, edges). """ - blossom = self.blossom[b] - assert blossom is not None -# TODO : consider whether we can do without the explicit list of nodes (if not, that's fine too) - - nodes: list[int] = [sub] + nodes: list[_Blossom] = [sub] edges: list[tuple[int, int]] = [] # Walk around the blossom from "sub" to its base. @@ -1068,28 +1010,24 @@ class _MatchingContext: return (nodes, edges) - def expand_t_blossom(self, b: int) -> None: + def expand_t_blossom(self, blossom: _NonTrivialBlossom) -> None: """Expand the specified T-blossom. This function takes time O(n). """ - num_vertex = self.graph.num_vertex - - blossom = self.blossom[b] - assert blossom is not None - assert self.blossom_parent[b] == -1 - assert self.blossom_label[b] == _LABEL_T + assert blossom.parent is None + assert blossom.label == _LABEL_T # Convert sub-blossoms into top-level blossoms. for sub in blossom.subblossoms: - self.blossom_parent[sub] = -1 - if sub < num_vertex: - self.vertex_blossom[sub] = sub + assert sub.label == _LABEL_NONE + sub.parent = None + if isinstance(sub, _NonTrivialBlossom): + for x in sub.vertices(): + self.vertex_top_blossom[x] = sub else: - for x in self.blossom_vertices(sub): - self.vertex_blossom[x] = sub - assert self.blossom_label[sub] == _LABEL_NONE + self.vertex_top_blossom[sub.base_vertex] = sub # The expanding blossom was part of an alternating tree, linked to # a parent node in the tree via one of its subblossoms, and linked to @@ -1100,19 +1038,20 @@ class _MatchingContext: # Find the sub-blossom that is attached to the parent node in # the alternating tree. # TODO : uglyness with the assertion - entry_link = self.blossom_link[b] + entry_link = blossom.tree_edge assert entry_link is not None (x, y) = entry_link - sub = self.vertex_blossom[y] + sub = self.vertex_top_blossom[y] # Assign label T to that sub-blossom. - self.blossom_label[sub] = _LABEL_T - self.blossom_link[sub] = (x, y) + sub.label = _LABEL_T + sub.tree_edge = blossom.tree_edge # Walk through the expanded blossom from "sub" to the base vertex. # Assign alternating S and T labels to the sub-blossoms and attach # them to the alternating tree. - (path_nodes, path_edges) = self.find_path_through_blossom(b, sub) + (path_nodes, path_edges) = self.find_path_through_blossom(blossom, + sub) for p in range(0, len(path_edges), 2): # @@ -1126,54 +1065,53 @@ class _MatchingContext: (y, x) = path_edges[p] self.assign_label_s(x) - # Assign label T to path_nodes[i+2] and attach it to path_nodes[p+1]. + # Assign label T to path_nodes[i+2] and attach it + # to path_nodes[p+1]. sub = path_nodes[p+2] - self.blossom_label[sub] = _LABEL_T - self.blossom_link[sub] = path_edges[p+1] + sub.label = _LABEL_T + sub.tree_edge = path_edges[p+1] - # Unlabel and delete the expanded blossom. Recycle its blossom index. - self.blossom_label[b] = _LABEL_NONE - self.blossom_link[b] = None - self.blossom[b] = None - self.unused_blossoms.append(b) + # Delete the expanded blossom. +# TODO : avoid O(n) here if possible + self.nontrivial_blossom.remove(blossom) - def expand_blossom_rec(self, b: int, stack: list[int]) -> None: - """Expand blossom "b" and recursively expand any sub-blossoms that - have dual variable zero. + def expand_blossom_rec( + self, + blossom: _NonTrivialBlossom, + stack: list[_NonTrivialBlossom] + ) -> None: + """Expand the specified blossom and recursively expand any + sub-blossoms that have dual variable zero. Use the stack object instead of making direct recursive calls. """ - num_vertex = self.graph.num_vertex + assert blossom.parent is None - blossom = self.blossom[b] - assert blossom is not None - assert self.blossom_parent[b] == -1 - - # Examine sub-blossoms of "b". + # Examine sub-blossoms. for sub in blossom.subblossoms: # Mark the sub-blossom as a top-level blossom. - self.blossom_parent[sub] = -1 + sub.parent = None - if sub < num_vertex: - # Trivial sub-blossom. Mark it as top-level vertex. - self.vertex_blossom[sub] = sub - else: + if isinstance(sub, _NonTrivialBlossom): # Non-trivial sub-blossom. # If its dual variable is zero, expand it recursively. - if self.get_blossom(sub).dual_var == 0: + if sub.dual_var == 0: stack.append(sub) else: # This sub-blossom will not be expanded; # it now becomes top-level. Update its vertices # to point to this sub-blossom. - for x in self.blossom_vertices(sub): - self.vertex_blossom[x] = sub + for x in sub.vertices(): + self.vertex_top_blossom[x] = sub + else: + # Trivial sub-blossom. Mark it as top-level vertex. + self.vertex_top_blossom[sub.base_vertex] = sub - # Delete the expanded blossom. Recycle its blossom index. - self.blossom[b] = None - self.unused_blossoms.append(b) + # Delete the expanded blossom. +# TODO : avoid O(n) here if possible + self.nontrivial_blossom.remove(blossom) def expand_zero_dual_blossoms(self) -> None: """Expand all blossoms with zero dual variable (recursively). @@ -1184,78 +1122,76 @@ class _MatchingContext: This function takes time O(n). """ - num_vertex = self.graph.num_vertex - # Use an explicit stack to avoid deep recursion. # The stack contains a list of blossoms to be expanded. - stack: list[int] = [] + stack: list[_NonTrivialBlossom] = [] # Find top-level blossoms with zero slack. - for b in range(num_vertex, 2 * num_vertex): - blossom = self.blossom[b] - if (blossom is not None) and (self.blossom_parent[b] == -1): + for blossom in self.nontrivial_blossom: + if blossom.parent is None: # We typically expand only S-blossoms that were created after # the most recent delta step. Those blossoms have _exactly_ # zero dual. So this comparison is reliable, even in case # of floating point edge weights. if blossom.dual_var == 0: - stack.append(b) + stack.append(blossom) # Expand blossoms. while stack: - b = stack.pop() - self.expand_blossom_rec(b, stack) + blossom = stack.pop() + self.expand_blossom_rec(blossom, stack) + + # + # Augmenting: + # def augment_blossom_rec( self, - b: int, - sub: int, - stack: list[tuple[int, int]] + blossom: _NonTrivialBlossom, + sub: _Blossom, + stack: list[tuple[_NonTrivialBlossom, _Blossom]] ) -> None: - """Augment along an alternating path through blossom "b", + """Augment along an alternating path through the specified blossom, from sub-blossom "sub" to the base vertex of the blossom. Modify the blossom to reflect that sub-blossom "sub" contains the base vertex after augmenting. - Recursively augment any sub-blossoms on the alternating path, - except for sub-blossom "sub" which has already been augmented. - Use the stack object instead of making direct recursive calls. + Mark any sub-blossoms on the alternating path for recursive + augmentation, except for sub-blossom "sub" which has already been + augmented. Use the stack instead of making direct recursive calls. """ - num_vertex = self.graph.num_vertex - blossom = self.blossom[b] - assert blossom is not None - - # Walk through the expanded blossom from "sub" to the base vertex. - (path_nodes, path_edges) = self.find_path_through_blossom(b, sub) + # Walk through the blossom from "sub" to the base vertex. + (path_nodes, path_edges) = self.find_path_through_blossom(blossom, + sub) for p in range(0, len(path_edges), 2): # Before augmentation: # path_nodes[p] is matched to path_nodes[p+1] # - # (p) ===== (p+1) ---(i,j)--- (p+2) + # (p) ===== (p+1) ---(x,y)--- (p+2) # # After augmentation: # path_nodes[p+1] matched to path_nodes[p+2] via edge (i,j) # - # (p) ----- (p+1) ===(i,j)=== (p+2) + # (p) ----- (p+1) ===(x,y)=== (p+2) # - # Pull the edge (i, j) into the matching. + # Pull the edge (x, y) into the matching. (x, y) = path_edges[p+1] self.vertex_mate[x] = y self.vertex_mate[y] = x - # Augment through the subblossoms touching the edge (i, j). + # Augment through the subblossoms touching the edge (x, y). # Nothing needs to be done for trivial subblossoms. bx = path_nodes[p+1] - if bx >= num_vertex: - stack.append((bx, x)) + if isinstance(bx, _NonTrivialBlossom): + stack.append((bx, self.trivial_blossom[x])) by = path_nodes[p+2] - if by >= num_vertex: - stack.append((by, y)) + if isinstance(by, _NonTrivialBlossom): + stack.append((by, self.trivial_blossom[y])) # Rotate the subblossom list so the new base ends up in position 0. p = blossom.subblossoms.index(sub) @@ -1266,40 +1202,45 @@ class _MatchingContext: # Update the base vertex. # We can pull this from the sub-blossom where we started since # its augmentation has already finished. - if sub < num_vertex: - blossom.base_vertex = sub - else: - blossom.base_vertex = self.get_blossom(sub).base_vertex + blossom.base_vertex = sub.base_vertex - def augment_blossom(self, b: int, sub: int) -> None: - """Augment along an alternating path through blossom "b", + def augment_blossom( + self, + blossom: _NonTrivialBlossom, + sub: _Blossom + ) -> None: + """Augment along an alternating path through the specified blossom, from sub-blossom "sub" to the base vertex of the blossom. + Recursively augment any sub-blossoms on the alternating path. + This function takes time O(n). """ # Use an explicit stack to avoid deep recursion. - stack = [(b, sub)] + stack = [(blossom, sub)] while stack: (outer_blossom, sub) = stack.pop() - b = self.blossom_parent[sub] + assert sub.parent is not None + blossom = sub.parent - if b != outer_blossom: + if blossom != outer_blossom: # Sub-blossom "sub" is an indirect (nested) child of - # the blossom we are supposed to be augmenting. + # the "outer_blossom" we are supposed to be augmenting. # - # Blossom "b" is the direct parent of "sub". - # Let's first augment "b" from "sub" to its base vertex. - # Then continue bay augmenting the parent of "b" from - # "b" to its base vertex, and so on until we get to the - # outer blossom. + # "blossom" is the direct parent of "sub". + # Let's first augment "blossom" from "sub" to its base vertex. + # Then continue by augmenting the parent of "blossom", + # from "blossom" to its base vertex, and so on until we + # get to the "outer_blossom". # - # Set up to continue augmenting through the parent of "b". - stack.append((outer_blossom, b)) + # Set up to continue augmenting through the parent of + # "blossom". + stack.append((outer_blossom, blossom)) - # Augment blossom "b" from "sub" to the base vertex of "b". - self.augment_blossom_rec(b, sub, stack) + # Augment "blossom" from "sub" to the base vertex. + self.augment_blossom_rec(blossom, sub, stack) def augment_matching(self, path: _AlternatingPath) -> None: """Augment the matching through the specified augmenting path. @@ -1311,10 +1252,8 @@ class _MatchingContext: # an unmatched vertex or a blossom with unmatched base. assert len(path.edges) % 2 == 1 for x in (path.edges[0][0], path.edges[-1][1]): - b = self.vertex_blossom[x] - if b != x: - x = self.get_blossom(b).base_vertex - assert self.vertex_mate[x] == -1 + b = self.vertex_top_blossom[x] + assert self.vertex_mate[b.base_vertex] == -1 # The augmenting path looks like this: # @@ -1331,18 +1270,22 @@ class _MatchingContext: # Augment the non-trivial blossoms on either side of this edge. # No action is necessary for trivial blossoms. - bx = self.vertex_blossom[x] - if bx != x: - self.augment_blossom(bx, x) + bx = self.vertex_top_blossom[x] + if isinstance(bx, _NonTrivialBlossom): + self.augment_blossom(bx, self.trivial_blossom[x]) - by = self.vertex_blossom[y] - if by != y: - self.augment_blossom(by, y) + by = self.vertex_top_blossom[y] + if isinstance(by, _NonTrivialBlossom): + self.augment_blossom(by, self.trivial_blossom[y]) # Pull the edge into the matching. self.vertex_mate[x] = y self.vertex_mate[y] = x + # + # Labeling and alternating tree expansion: + # + def assign_label_s(self, x: int) -> None: """Assign label S to the unlabeled blossom that contains vertex "x". @@ -1357,37 +1300,34 @@ class _MatchingContext: a T-vertex via a tight edge. """ - # Assign label S to the blossom that contains vertex "v". - bx = self.vertex_blossom[x] - assert self.blossom_label[bx] == _LABEL_NONE - self.blossom_label[bx] = _LABEL_S + # Assign label S to the blossom that contains vertex "x". + bx = self.vertex_top_blossom[x] + assert bx.label == _LABEL_NONE + bx.label = _LABEL_S y = self.vertex_mate[x] if y == -1: # Vertex "x" is unmatched. # It must be either a top-level vertex or the base vertex of # a top-level blossom. - assert (bx == x) or (self.get_blossom(bx).base_vertex == x) + assert bx.base_vertex == x - # Mark the blossom that contains "v" as root of an alternating tree. - self.blossom_link[bx] = None + # Mark the blossom as root of an alternating tree. + bx.tree_edge = None else: # Vertex "x" is matched to T-vertex "y". - by = self.vertex_blossom[y] - assert self.blossom_label[by] == _LABEL_T + by = self.vertex_top_blossom[y] + assert by.label == _LABEL_T # Attach the blossom that contains "x" to the alternating tree. - self.blossom_link[bx] = (y, x) + bx.tree_edge = (y, x) # Start least-slack edge tracking for the S-blossom. - self.blossom_best_edges.new_blossom(bx) + self.lset_new_blossom(bx) # Add all vertices inside the newly labeled S-blossom to the queue. - if bx == x: - self.queue.append(x) - else: - self.queue.extend(self.blossom_vertices(bx)) + self.queue.extend(bx.vertices()) def assign_label_t(self, x: int, y: int) -> None: """Assign label T to the unlabeled blossom that contains vertex "y". @@ -1400,17 +1340,16 @@ class _MatchingContext: - "y" is an unlabeled, matched vertex. - There is a tight edge between vertices "x" and "y". """ - assert self.blossom_label[self.vertex_blossom[x]] == _LABEL_S + assert self.vertex_top_blossom[x].label == _LABEL_S # Assign label T to the unlabeled blossom. - by = self.vertex_blossom[y] - assert self.blossom_label[by] == _LABEL_NONE - self.blossom_label[by] = _LABEL_T - self.blossom_link[by] = (x, y) + by = self.vertex_top_blossom[y] + assert by.label == _LABEL_NONE + by.label = _LABEL_T + by.tree_edge = (x, y) # Assign label S to the blossom that contains the mate of vertex "y". - ybase = y if by == y else self.get_blossom(by).base_vertex - z = self.vertex_mate[ybase] + z = self.vertex_mate[by.base_vertex] assert z != -1 self.assign_label_s(z) @@ -1437,7 +1376,7 @@ class _MatchingContext: # but not necessarily in the same vertex within that blossom. p = path.edges[0][0] q = path.edges[-1][1] - if self.vertex_blossom[p] == self.vertex_blossom[q]: + if self.vertex_top_blossom[p] is self.vertex_top_blossom[q]: self.make_blossom(path) return None else: @@ -1466,8 +1405,8 @@ class _MatchingContext: x = self.queue.pop() # Double-check that "x" is an S-vertex. - bx = self.vertex_blossom[x] - assert self.blossom_label[bx] == _LABEL_S + bx = self.vertex_top_blossom[x] + assert bx.label == _LABEL_S # Scan the edges that are incident on "x". # This loop runs through O(m) iterations per stage. @@ -1480,14 +1419,14 @@ class _MatchingContext: # Note: blossom index of vertex "x" may change during # this loop, so we need to refresh it here. - bx = self.vertex_blossom[x] - by = self.vertex_blossom[y] + bx = self.vertex_top_blossom[x] + by = self.vertex_top_blossom[y] # Ignore edges that are internal to a blossom. - if bx == by: + if bx is by: continue - ylabel = self.blossom_label[by] + ylabel = by.label # Check whether this edge is tight (has zero slack). # Only tight edges may be part of an alternating tree. @@ -1505,19 +1444,25 @@ class _MatchingContext: elif ylabel == _LABEL_S: # Update tracking of least-slack edges between S-blossoms. - self.blossom_best_edges.add_edge(self, bx, e, slack) + self.lset_add_blossom_edge(bx, e, slack) if ylabel != _LABEL_S: # Update tracking of least-slack edges from vertex "y" to # any S-vertex. We do this for T-vertices and unlabeled # vertices. Edges which already have zero slack are still # tracked. - self.vertex_best_edges.add_edge(self, y, e, slack) + self.lset_add_vertex_edge(y, e, slack) # No further S vertices to scan, and no augmenting path found. return None - def substage_calc_dual_delta(self) -> tuple[int, float|int, int, int]: + # + # Delta steps: + # + + def substage_calc_dual_delta( + self + ) -> tuple[int, float|int, int, Optional[_NonTrivialBlossom]]: """Calculate a delta step in the dual LPP problem. This function returns the minimum of the 4 types of delta values, @@ -1537,18 +1482,18 @@ class _MatchingContext: num_vertex = self.graph.num_vertex delta_edge = -1 - delta_blossom = -1 + delta_blossom: Optional[_NonTrivialBlossom] = None # Compute delta1: minimum dual variable of any S-vertex. delta_type = 1 delta_2x = min( self.vertex_dual_2x[x] for x in range(num_vertex) - if self.blossom_label[self.vertex_blossom[x]] == _LABEL_S) + if self.vertex_top_blossom[x].label == _LABEL_S) # Compute delta2: minimum slack of any edge between an S-vertex and # an unlabeled vertex. - (e, slack) = self.vertex_best_edges.get_least_slack_edge(self) + (e, slack) = self.lset_get_best_vertex_edge() if (e != -1) and (slack <= delta_2x): delta_type = 2 delta_2x = slack @@ -1556,7 +1501,7 @@ class _MatchingContext: # Compute delta3: half minimum slack of any edge between two top-level # S-blossoms. - (e, slack) = self.blossom_best_edges.get_least_slack_edge(self) + (e, slack) = self.lset_get_best_blossom_edge() if e != -1: if self.graph.integer_weights: # If all edge weights are even integers, the slack @@ -1566,20 +1511,18 @@ class _MatchingContext: slack = slack // 2 else: slack = slack / 2 - if slack < delta_2x: + if slack <= delta_2x: delta_type = 3 delta_2x = slack delta_edge = e # Compute delta4: half minimum dual variable of a top-level T-blossom. - for b in range(num_vertex, 2 * num_vertex): - if (self.blossom_label[b] == _LABEL_T - and self.blossom_parent[b] == -1): - slack = self.get_blossom(b).dual_var - if slack < delta_2x: + for blossom in self.nontrivial_blossom: + if (blossom.label == _LABEL_T) and (blossom.parent is None): + if blossom.dual_var <= delta_2x: delta_type = 4 - delta_2x = slack - delta_blossom = b + delta_2x = blossom.dual_var + delta_blossom = blossom return (delta_type, delta_2x, delta_edge, delta_blossom) @@ -1590,7 +1533,7 @@ class _MatchingContext: # Apply delta to dual variables of all vertices. for x in range(num_vertex): - xlabel = self.blossom_label[self.vertex_blossom[x]] + xlabel = self.vertex_top_blossom[x].label if xlabel == _LABEL_S: # S-vertex: subtract delta from dual variable. self.vertex_dual_2x[x] -= delta_2x @@ -1599,15 +1542,19 @@ class _MatchingContext: self.vertex_dual_2x[x] += delta_2x # Apply delta to dual variables of top-level non-trivial blossoms. - for b in range(num_vertex, 2 * num_vertex): - if self.blossom_parent[b] == -1: - blabel = self.blossom_label[b] + for blossom in self.nontrivial_blossom: + if blossom.parent is None: + blabel = blossom.label if blabel == _LABEL_S: # S-blossom: add 2*delta to dual variable. - self.get_blossom(b).dual_var += delta_2x + blossom.dual_var += delta_2x elif blabel == _LABEL_T: # T-blossom: subtract 2*delta from dual variable. - self.get_blossom(b).dual_var -= delta_2x + blossom.dual_var -= delta_2x + + # + # Main stage function: + # def run_stage(self) -> bool: """Run one stage of the matching algorithm. @@ -1664,12 +1611,12 @@ class _MatchingContext: # Use the edge from S-vertex to unlabeled vertex that got # unlocked through the delta update. (x, y, _w) = self.graph.edges[delta_edge] - if self.blossom_label[self.vertex_blossom[x]] != _LABEL_S: + if self.vertex_top_blossom[x].label != _LABEL_S: (x, y) = (y, x) self.assign_label_t(x, y) elif delta_type == 3: - # Use the S-to-S edge that got unlocked through the delta update. + # Use the S-to-S edge that got unlocked by the delta update. # This may reveal an augmenting path. (x, y, _w) = self.graph.edges[delta_edge] augmenting_path = self.add_s_to_s_edge(x, y) @@ -1679,6 +1626,7 @@ class _MatchingContext: elif delta_type == 4: # Expand the T-blossom that reached dual value 0 through # the delta update. + assert delta_blossom is not None self.expand_t_blossom(delta_blossom) else: @@ -1686,9 +1634,6 @@ class _MatchingContext: assert delta_type == 1 break - # Remove all labels, clear queue. - self.reset_stage() - # Augment the matching if an augmenting path was found. if augmenting_path is not None: self.augment_matching(augmenting_path) @@ -1698,14 +1643,20 @@ class _MatchingContext: # get expanded as soon as their dual variable hits zero. self.expand_zero_dual_blossoms() + # Remove all labels, clear queue. + self.reset_stage() + # Return True if the matching was augmented. return (augmenting_path is not None) +# TODO : clean up this whole mess + def _verify_optimum(ctx: _MatchingContext) -> None: """Verify that the optimum solution has been found. This function takes time O(m * n). +TODO : really ?? Raises: AssertionError: If the solution is not optimal. @@ -1716,11 +1667,6 @@ def _verify_optimum(ctx: _MatchingContext) -> None: vertex_mate = ctx.vertex_mate vertex_dual_var_2x = ctx.vertex_dual_2x - # Extract dual values of blossoms - blossom_dual_var = [ - (blossom.dual_var if blossom is not None else 0) - for blossom in ctx.blossom] - # Double-check that each matching edge actually exists in the graph. num_matched_vertex = 0 for x in range(num_vertex): @@ -1736,50 +1682,51 @@ def _verify_optimum(ctx: _MatchingContext) -> None: # Check that all dual variables are non-negative. assert min(vertex_dual_var_2x) >= 0 - assert min(blossom_dual_var) >= 0 + for blossom in ctx.nontrivial_blossom: + assert blossom.dual_var >= 0 # Count the number of vertices in each blossom. - blossom_nvertex = (2 * num_vertex) * [0] + blossom_nvertex = {id(blossom): 0 for blossom in ctx.nontrivial_blossom} for x in range(num_vertex): - b = ctx.blossom_parent[x] - while b != -1: - blossom_nvertex[b] += 1 - b = ctx.blossom_parent[b] + b = ctx.trivial_blossom[x] + while b.parent is not None: + b = b.parent + blossom_nvertex[id(b)] += 1 # Calculate slack of each edge. # Also count the number of matched edges in each blossom. - blossom_nmatched = (2 * num_vertex) * [0] + blossom_nmatched = {id(blossom): 0 for blossom in ctx.nontrivial_blossom} for (x, y, w) in ctx.graph.edges: # List blossoms that contain vertex "x". xblossoms = [] - bx = ctx.blossom_parent[x] - while bx != -1: + bx = ctx.trivial_blossom[x] + while bx.parent is not None: + bx = bx.parent xblossoms.append(bx) - bx = ctx.blossom_parent[bx] # List blossoms that contain vertex "y". yblossoms = [] - by = ctx.blossom_parent[y] - while by != -1: + by = ctx.trivial_blossom[y] + while by.parent is not None: + by = by.parent yblossoms.append(by) - by = ctx.blossom_parent[by] # List blossoms that contain the edge (x, y). - edge_blossoms = [] + edge_blossoms: list[_NonTrivialBlossom] = [] for (bx, by) in zip(reversed(xblossoms), reversed(yblossoms)): - if bx != by: + if bx is not by: break edge_blossoms.append(bx) # Calculate edge slack = # dual[x] + dual[y] - weight - # + sum(dual[b] for blossoms "b" containing the edge) + # + sum(blossom.dual_var for "blossom" containing the edge) # # Multiply weights by 2 to ensure integer values. slack = vertex_dual_var_2x[x] + vertex_dual_var_2x[y] - 2 * w - slack += 2 * sum(blossom_dual_var[b] for b in edge_blossoms) + slack += 2 * sum(blossom.dual_var for blossom in edge_blossoms) # Check that all edges have non-negative slack. assert slack >= 0 @@ -1791,7 +1738,7 @@ def _verify_optimum(ctx: _MatchingContext) -> None: # Update number of matched edges in each blossom. if vertex_mate[x] == y: for b in edge_blossoms: - blossom_nmatched[b] += 1 + blossom_nmatched[id(b)] += 1 # Check that all unmatched vertices have zero dual. for x in range(num_vertex): @@ -1801,8 +1748,8 @@ def _verify_optimum(ctx: _MatchingContext) -> None: # Check that all blossoms with positive dual are "full". # A blossom is full if all except one of its vertices are matched # to another vertex in the same blossom. - for b in range(num_vertex, 2 * num_vertex): - if blossom_dual_var[b] > 0: - assert blossom_nvertex[b] == 2 * blossom_nmatched[b] + 1 + for blossom in ctx.nontrivial_blossom: + if blossom.dual_var > 0: + assert blossom_nvertex[id(blossom)] == 2 * blossom_nmatched[id(blossom)] + 1 # Optimum solution confirmed.