Implement heap-based tracking for delta4

python/mwmatching.py

@@ -34,8 +34,9 @@ def maximum_weight_matching(
     no effect on the maximum-weight matching.
     Edges with negative weight are ignored.
 
-    This function takes time O(n**3 + n*m*log(n)), where "n" is the number of vertices.
+    This function takes time O(n**3 + n*m*log(n)),
-    This function uses O(n + m) memory, where "m" is the number of edges.
+    where "n" is the number of vertices and "m" is the number of edges.
+    This function uses O(n + m) memory.
 
     Parameters:
         edges: List of edges, each edge specified as a tuple "(x, y, w)"
@@ -447,6 +448,11 @@ class _NonTrivialBlossom(_Blossom):
         # New blossoms start with dual variable 0.
         self.dual_var: float = 0
 
+        # If this is a top-level T-blossom,
+        # "delta4_node" is the corresponding node in the delta4 queue.
+        # Otherwise "delta4_node" is None.
+        self.delta4_node: Optional[PriorityQueue.Node] = None
+
     def vertices(self) -> list[int]:
         """Return a list of vertex indices contained in the blossom."""
 
@@ -522,24 +528,31 @@ 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.
+        # Initial dual value of all vertices 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".
         # Multiplication by 2 ensures that the values are integers
         # if all edge weights are integers.
         #
         # Vertex duals are initialized to half the maximum edge weight.
-        max_weight = max(w for (_x, _y, w) in graph.edges)
+        self.start_vertex_dual_2x = max(w for (_x, _y, w) in graph.edges)
-        self.vertex_dual_2x: list[float] = num_vertex * [max_weight]
+
+        # Every vertex has a variable in the dual LPP.
+        # "vertex_dual_2x[x]" is 2 times the dual variable of vertex "x".
+        self.vertex_dual_2x: list[float]
+        self.vertex_dual_2x = num_vertex * [self.start_vertex_dual_2x]
+
+        # 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 an edge is its slack plus two times the
+        # running sum of delta updates.
+        self.delta3_queue: PriorityQueue[int] = PriorityQueue()
+        self.delta3_set: set[int] = set()
+
+        # Queue containing top-level non-trivial T-blossoms.
+        # The priority of a blossom is its dual plus two times the running
+        # sum of delta updates.
+        self.delta4_queue: PriorityQueue[_NonTrivialBlossom] = PriorityQueue()
 
         # For each T-vertex or unlabeled vertex "x",
         # "vertex_best_edge[x]" is the edge index of the least-slack edge
@@ -612,9 +625,6 @@ class _MatchingContext:
         for x in range(num_vertex):
             self.vertex_best_edge[x] = -1
 
-        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".
 
@@ -663,6 +673,8 @@ class _MatchingContext:
 
         Marks all blossoms as unlabeled, clears the queue,
         and resets tracking of least-slack edges.
+
+        This function takes time O(n).
         """
 
         # Remove blossom labels.
@@ -676,6 +688,14 @@ class _MatchingContext:
         # Reset least-slack edge tracking.
         self.lset_reset()
 
+        # Reset delta queues.
+        self.delta3_queue.clear()
+        self.delta3_set.clear()
+        self.delta4_queue.clear()
+
+        for ntb in self.nontrivial_blossom:
+            ntb.delta4_node = None
+
     def trace_alternating_paths(self, x: int, y: int) -> _AlternatingPath:
         """Trace back through the alternating trees from vertices "x" and "y".
 
@@ -809,6 +829,13 @@ class _MatchingContext:
             if sub.label == _LABEL_T:
                 self.queue.extend(sub.vertices())
 
+        # Remove sub-blossoms from the delta4 queue.
+        for sub in subblossoms:
+            if ((isinstance(sub, _NonTrivialBlossom))
+                    and (sub.delta4_node is not None)):
+                self.delta4_queue.delete(sub.delta4_node)
+                sub.delta4_node = None
+
     @staticmethod
     def find_path_through_blossom(
             blossom: _NonTrivialBlossom,
@@ -845,6 +872,10 @@ class _MatchingContext:
         assert blossom.parent is None
         assert blossom.label == _LABEL_T
 
+        # Remove expanded blossom from the delta4 queue.
+        assert blossom.delta4_node is not None
+        self.delta4_queue.delete(blossom.delta4_node)
+
         # Convert sub-blossoms into top-level blossoms.
         for sub in blossom.subblossoms:
             assert sub.label == _LABEL_NONE
@@ -871,6 +902,12 @@ class _MatchingContext:
         sub.label = _LABEL_T
         sub.tree_edge = blossom.tree_edge
 
+        # Insert T sub-blossom into the delta4 queue.
+        if isinstance(sub, _NonTrivialBlossom):
+            assert sub.delta4_node is None
+            prio = sub.dual_var + self.delta_sum_2x
+            sub.delta4_node = self.delta4_queue.insert(prio, sub)
+
         # 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.
@@ -895,6 +932,12 @@ class _MatchingContext:
             sub.label = _LABEL_T
             sub.tree_edge = path_edges[p+1]
 
+            # Insert T sub-blossom into the delta4 queue.
+            if isinstance(sub, _NonTrivialBlossom):
+                assert sub.delta4_node is None
+                prio = sub.dual_var + self.delta_sum_2x
+                sub.delta4_node = self.delta4_queue.insert(prio, sub)
+
         # Delete the expanded blossom.
         self.nontrivial_blossom.remove(blossom)
 
@@ -1131,6 +1174,12 @@ class _MatchingContext:
         by.label = _LABEL_T
         by.tree_edge = (x, y)
 
+        # Insert T sub-blossom into the delta4 queue.
+        if isinstance(by, _NonTrivialBlossom):
+            assert by.delta4_node is None
+            prio = by.dual_var + self.delta_sum_2x
+            by.delta4_node = self.delta4_queue.insert(prio, by)
+
         # Assign label S to the blossom that is mated to the T-blossom.
         z = self.vertex_mate[by.base_vertex]
         assert z != -1
@@ -1265,17 +1314,14 @@ class _MatchingContext:
         Returns:
             Tuple (delta_type, delta_2x, delta_edge, delta_blossom).
         """
-        num_vertex = self.graph.num_vertex
-
         delta_edge = -1
         delta_blossom: Optional[_NonTrivialBlossom] = None
 
         # Compute delta1: minimum dual variable of any S-vertex.
+        # All unmatched vertices have the same dual value, and this is
+        # the minimum value among all S-vertices.
         delta_type = 1
-        delta_2x = min(
+        delta_2x = self.start_vertex_dual_2x - self.delta_sum_2x
-            self.vertex_dual_2x[x]
-            for x in range(num_vertex)
-            if self.vertex_top_blossom[x].label == _LABEL_S)
 
         # Compute delta2: minimum slack of any edge between an S-vertex and
         # an unlabeled vertex.
@@ -1311,12 +1357,14 @@ class _MatchingContext:
             self.delta3_set.remove(e)
 
         # Compute delta4: half minimum dual variable of a top-level T-blossom.
-        for blossom in self.nontrivial_blossom:
+        if not self.delta4_queue.empty():
-            if (blossom.label == _LABEL_T) and (blossom.parent is None):
+            blossom = self.delta4_queue.find_min().data
-                if blossom.dual_var <= delta_2x:
+            assert blossom.label == _LABEL_T
-                    delta_type = 4
+            assert blossom.parent is None
-                    delta_2x = blossom.dual_var
+            if blossom.dual_var <= delta_2x:
-                    delta_blossom = blossom
+                delta_type = 4
+                delta_2x = blossom.dual_var
+                delta_blossom = blossom
 
         return (delta_type, delta_2x, delta_edge, delta_blossom)