Implement heap-based edge tracking for delta3

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:
+        self.delta3_queue.clear()
-            blossom.best_edge = -1
+        self.delta3_set.clear()
-
-        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: 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()
+        while not self.delta3_queue.empty():
-        if e != -1:
+            delta3_node = self.delta3_queue.find_min()
-            if self.graph.integer_weights:
+            e = delta3_node.data
-                # If all edge weights are even integers, the slack
+            (x, y, _w) = self.graph.edges[e]
-                # of any edge between two S blossoms is also an even
+            bx = self.vertex_top_blossom[x]
-                # integer. Therefore the delta is an integer.
+            by = self.vertex_top_blossom[y]
-                assert slack % 2 == 0
+            assert (bx.label == _LABEL_S) and (by.label == _LABEL_S)
-                slack = slack // 2
+            if bx is not by:
-            else:
+                # Found edge between different top-level S-blossoms.
-                slack = slack / 2
+                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.