diff --git a/python/mwmatching.py b/python/mwmatching.py index f05f453..fe52b00 100644 --- a/python/mwmatching.py +++ b/python/mwmatching.py @@ -10,6 +10,8 @@ import collections from collections.abc import Sequence from typing import NamedTuple, Optional +from datastruct import UnionFindQueue, PriorityQueue + def maximum_weight_matching( edges: Sequence[tuple[int, int, float]] @@ -32,7 +34,7 @@ def maximum_weight_matching( no effect on the maximum-weight matching. Edges with negative weight are ignored. - This function takes time O(n**3), where "n" is the number of vertices. + This function takes time O(n**3 + n*m*log(n)), where "n" is the number of vertices. This function uses O(n + m) memory, where "m" is the number of edges. Parameters: @@ -380,11 +382,6 @@ class _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()". @@ -450,11 +447,6 @@ class _NonTrivialBlossom(_Blossom): # New blossoms start with dual variable 0. self.dual_var: 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.""" @@ -530,6 +522,15 @@ class _MatchingContext: # Initially all vertices are trivial top-level blossoms. self.vertex_top_blossom: list[_Blossom] = self.trivial_blossom.copy() + # Running sum of applied delta steps times 2. + self.delta_sum_2x: float = 0 + + # Queue containing edges between S-vertices in different top-level + # blossoms. The priority of each edge is equal to its slack + # plus two times the running sum of delta updates. + self.delta3_queue = PriorityQueue() + self.delta3_set: set[int] = set() + # Every vertex has a variable in the dual LPP. # # "vertex_dual_2x[x]" is 2 times the dual variable of vertex "x". @@ -563,6 +564,23 @@ class _MatchingContext: 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 edge_slack(self, e: int) -> float: + """Return the slack of the edge with index "e". + + This function must only be used for edges between two vertices + in different top-level S-blossoms. If all edge weights are + integers, the result is also an integer. + """ + (x, y, w) = self.graph.edges[e] + assert self.vertex_top_blossom[x] is not self.vertex_top_blossom[y] + dual_2x = self.vertex_dual_2x[x] + self.vertex_dual_2x[y] + if self.graph.integer_weights: + assert dual_2x % 2 == 0 + dual = dual_2x // 2 + else: + dual = dual_2x / 2 + return dual - w + # # Least-slack edge tracking: # @@ -576,16 +594,6 @@ class _MatchingContext: # 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 @@ -604,12 +612,8 @@ class _MatchingContext: 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 + self.delta3_queue.clear() + self.delta3_set.clear() def lset_add_vertex_edge(self, y: int, e: int, slack: float) -> None: """Add edge "e" from an S-vertex to unlabeled vertex or T-vertex "y". @@ -650,155 +654,6 @@ class _MatchingContext: return (best_index, best_slack) - @staticmethod - def lset_new_blossom(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: 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 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. - best_edge_to_blossom: list[int] = num_vertex * [-1] - zero_slack: float = 0 - best_slack_to_blossom: list[float] = num_vertex * [zero_slack] - - # And find the overall least-slack edge to any other S-blossom. - best_edge = -1 - best_slack: 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, 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: float = 0 - - 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: # @@ -954,9 +809,6 @@ class _MatchingContext: if sub.label == _LABEL_T: self.queue.extend(sub.vertices()) - # Merge least-slack edges for the S-sub-blossoms. - self.lset_merge_blossoms(blossom) - @staticmethod def find_path_through_blossom( blossom: _NonTrivialBlossom, @@ -1251,9 +1103,6 @@ class _MatchingContext: # Attach the blossom that contains "x" to the alternating tree. bx.tree_edge = (y, x) - # Start least-slack edge tracking for the S-blossom. - self.lset_new_blossom(bx) - # Add all vertices inside the newly labeled S-blossom to the queue. self.queue.extend(bx.vertices()) @@ -1378,7 +1227,10 @@ class _MatchingContext: elif ylabel == _LABEL_S: # Update tracking of least-slack edges between S-blossoms. - self.lset_add_blossom_edge(bx, e, slack) + if e not in self.delta3_set: + prio = self.edge_slack(e) + self.delta_sum_2x + self.delta3_set.add(e) + self.delta3_queue.insert(prio, e) if ylabel != _LABEL_S: # Update tracking of least-slack edges from vertex "y" to @@ -1435,20 +1287,28 @@ class _MatchingContext: # Compute delta3: half minimum slack of any edge between two top-level # S-blossoms. - (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 - # of any edge between two S blossoms is also an even - # integer. Therefore the delta is an integer. - assert slack % 2 == 0 - slack = slack // 2 - else: - slack = slack / 2 - if slack <= delta_2x: - delta_type = 3 - delta_2x = slack - delta_edge = e + while not self.delta3_queue.empty(): + delta3_node = self.delta3_queue.find_min() + e = delta3_node.data + (x, y, _w) = self.graph.edges[e] + bx = self.vertex_top_blossom[x] + by = self.vertex_top_blossom[y] + assert (bx.label == _LABEL_S) and (by.label == _LABEL_S) + if bx is not by: + # Found edge between different top-level S-blossoms. + slack = self.edge_slack(e) + if slack <= delta_2x: + delta_type = 3 + delta_2x = slack + delta_edge = e + break + + # Reject edges between vertices within the same top-level blossom. + # Although intra-blossom edges are never inserted into the queue, + # existing edges in the queue may become intra-blossom when + # a new blossom is formed. + self.delta3_queue.delete(delta3_node) + self.delta3_set.remove(e) # Compute delta4: half minimum dual variable of a top-level T-blossom. for blossom in self.nontrivial_blossom: @@ -1465,6 +1325,8 @@ class _MatchingContext: num_vertex = self.graph.num_vertex + self.delta_sum_2x += delta_2x + # Apply delta to dual variables of all vertices. for x in range(num_vertex): xlabel = self.vertex_top_blossom[x].label @@ -1498,7 +1360,7 @@ class _MatchingContext: thereby increasing the number of matched edges by 1. If no such path is found, the matching must already be optimal. - This function takes time O(n**2). + This function takes time O(n**2 + m * log(n)). Returns: True if the matching was successfully augmented.