Lazy delta updates of S-vertex duals

python/mwmatching.py

@@ -536,7 +536,15 @@ class _MatchingContext:
         self.start_vertex_dual_2x = max(w for (_x, _y, w) in graph.edges)
 
         # Every vertex has a variable in the dual LPP.
+        #
+        # For an unlabeled vertex "x",
         # "vertex_dual_2x[x]" is 2 times the dual variable of vertex "x".
+        #
+        # For an S-vertex "x",
+        # "vertex_dual_2x[x]" is 2 times the dual variable of vertex "x"
+        # plus two times the running sum of delta updates.
+        #
+        # For a T-vertex "x", ... TODO
         self.vertex_dual_2x: list[float]
         self.vertex_dual_2x = num_vertex * [self.start_vertex_dual_2x]
 
@@ -560,39 +568,33 @@ class _MatchingContext:
         self.vertex_best_edge: list[int] = num_vertex * [-1]
 
         # Queue of S-vertices to be scanned.
-        self.queue: collections.deque[int] = collections.deque()
+        self.scan_queue: collections.deque[int] = collections.deque()
 
-    def edge_slack_2x(self, e: int) -> float:
+    def edge_slack_2x(
-        """Return 2 times the slack of the edge with index "e".
+            self,
+            x: int,
+            y: int,
+            bx: _Blossom,
+            by: _Blossom,
+            w: float
+            ) -> float:
+        """Calculate edge slack.
 
-        The result is only valid for edges that are not between vertices
+        Return 2 times the slack of an edge from vertex "x" in blossom "bx"
-        that belong to the same top-level blossom.
+        to vertex "y" in blossom "by" with weight "w".
+
+        The two vertices must be in different top-level blossoms.
 
         Multiplication by 2 ensures that the return value is an integer
         if all edge weights are integers.
-
-        This function is called O(n**2) times per stage.
         """
-        (x, y, w) = self.graph.edges[e]
+        assert bx is not by
-        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:
+        if bx.label == _LABEL_S:
-            assert dual_2x % 2 == 0
+            dual_2x -= self.delta_sum_2x
-            dual = dual_2x // 2
+        if by.label == _LABEL_S:
-        else:
+            dual_2x -= self.delta_sum_2x
-            dual = dual_2x / 2
+        return dual_2x - 2 * w
-        return dual - w
 
     #
     # Least-slack edge tracking:
@@ -635,7 +637,10 @@ class _MatchingContext:
         if best_edge == -1:
             self.vertex_best_edge[y] = e
         else:
-            best_slack = self.edge_slack_2x(best_edge)
+            (xx, yy, w) = self.graph.edges[best_edge]
+            bx = self.vertex_top_blossom[xx]
+            by = self.vertex_top_blossom[yy]
+            best_slack = self.edge_slack_2x(xx, yy, bx, by, w)
             if slack < best_slack:
                 self.vertex_best_edge[y] = e
 
@@ -657,7 +662,10 @@ class _MatchingContext:
             if self.vertex_top_blossom[x].label == _LABEL_NONE:
                 e = self.vertex_best_edge[x]
                 if e != -1:
-                    slack = self.edge_slack_2x(e)
+                    (x, y, w) = self.graph.edges[e]
+                    bx = self.vertex_top_blossom[x]
+                    by = self.vertex_top_blossom[y]
+                    slack = self.edge_slack_2x(x, y, bx, by, w)
                     if (best_index == -1) or (slack < best_slack):
                         best_index = e
                         best_slack = slack
@@ -668,6 +676,70 @@ class _MatchingContext:
     # General support routines:
     #
 
+    def assign_blossom_label_s(self, blossom: _Blossom) -> None:
+        """Assign label S to an unlabeled top-level blossom.
+
+        All vertices in the newly labeled blossom are added to the scan queue.
+        Dual variables are adjusted.
+        """
+
+        assert blossom.parent is None
+        assert blossom.label == _LABEL_NONE
+
+        # Assign label S.
+        blossom.label = _LABEL_S
+
+        # Adjust the vertex dual variables of the new S-vertices.
+        vertices = blossom.vertices()
+        for x in vertices:
+            self.vertex_dual_2x[x] += self.delta_sum_2x
+
+        # Add all vertices inside the newly labeled S-blossom to the queue.
+        self.scan_queue.extend(vertices)
+
+    def remove_blossom_label_s(self, blossom: _Blossom) -> None:
+        """Remove label S from a top-level S-blossom."""
+
+        assert blossom.parent is None
+        assert blossom.label == _LABEL_S
+
+        # Remove label.
+        blossom.label = _LABEL_NONE
+
+        # Adjust the vertex dual variables of the former S-vertices.
+        for x in blossom.vertices():
+            self.vertex_dual_2x[x] -= self.delta_sum_2x
+
+    def assign_blossom_label_t(self, blossom: _Blossom) -> None:
+        """Assign label T to an unlabeled top-level blossom."""
+
+        assert blossom.parent is None
+        assert blossom.label == _LABEL_NONE
+
+        # Assign label T.
+        blossom.label = _LABEL_T
+
+        # Insert blossom into the delta4 queue.
+        if isinstance(blossom, _NonTrivialBlossom):
+            assert blossom.delta4_node is None
+            prio = blossom.dual_var + self.delta_sum_2x
+            blossom.delta4_node = self.delta4_queue.insert(prio, blossom)
+
+    def remove_blossom_label_t(self, blossom: _Blossom) -> None:
+        """Remove label T from a top-level T-blossom."""
+
+        assert blossom.parent is None
+        assert blossom.label == _LABEL_T
+
+        # Remove label.
+        blossom.label = _LABEL_NONE
+
+        # Remove blossom from delta4 queue.
+        if isinstance(blossom, _NonTrivialBlossom):
+            assert blossom.delta4_node is not None
+            self.delta4_queue.delete(blossom.delta4_node)
+            blossom.delta4_node = None
+
     def reset_stage(self) -> None:
         """Reset data which are only valid during a stage.
 
@@ -679,11 +751,16 @@ class _MatchingContext:
 
         # Remove blossom labels.
         for blossom in self.trivial_blossom + self.nontrivial_blossom:
+            if blossom.parent is None:
+                if blossom.label == _LABEL_S:
+                    self.remove_blossom_label_s(blossom)
+                elif blossom.label == _LABEL_T:
+                    self.remove_blossom_label_t(blossom)
             blossom.label = _LABEL_NONE
             blossom.tree_edge = None
 
-        # Clear the queue.
+        # Clear the scan queue.
-        self.queue.clear()
+        self.scan_queue.clear()
 
         # Reset least-slack edge tracking.
         self.lset_reset()
@@ -691,10 +768,7 @@ class _MatchingContext:
         # Reset delta queues.
         self.delta3_queue.clear()
         self.delta3_set.clear()
-        self.delta4_queue.clear()
+        assert self.delta4_queue.empty()
-
-        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".
@@ -797,16 +871,30 @@ class _MatchingContext:
         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_,
+        # Note the path will not always start and end with the same _vertex_,
         # but it must start and end in the same _blossom_.
         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]
 
+        # Blossom must start and end with an S-sub-blossom.
+        assert subblossoms[0].label == _LABEL_S
+
+        # Relabel the T-sub-blossoms to S-sub-blossoms.
+        # This also adds new S-vertices to the scan queue.
+        for sub in subblossoms:
+            if sub.label == _LABEL_T:
+                self.remove_blossom_label_t(sub)
+                self.assign_blossom_label_s(sub)
+
         # Create the new blossom object.
         blossom = _NonTrivialBlossom(subblossoms, path.edges)
 
+        # Assign label S to the new blossom and link it to the tree.
+        blossom.label = _LABEL_S
+        blossom.tree_edge = subblossoms[0].tree_edge
+
         # Insert into the blossom array.
         self.nontrivial_blossom.append(blossom)
 
@@ -818,24 +906,6 @@ class _MatchingContext:
         for x in blossom.vertices():
             self.vertex_top_blossom[x] = blossom
 
-        # Assign label S to the new 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 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,
@@ -875,6 +945,7 @@ class _MatchingContext:
         # Remove expanded blossom from the delta4 queue.
         assert blossom.delta4_node is not None
         self.delta4_queue.delete(blossom.delta4_node)
+        blossom.delta4_node = None
 
         # Convert sub-blossoms into top-level blossoms.
         for sub in blossom.subblossoms:
@@ -899,15 +970,9 @@ class _MatchingContext:
         sub = self.vertex_top_blossom[y]
 
         # Assign label T to that sub-blossom.
-        sub.label = _LABEL_T
+        self.assign_blossom_label_t(sub)
         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.
@@ -924,20 +989,14 @@ class _MatchingContext:
 
             # Assign label S to path_nodes[p+1].
             (y, x) = path_edges[p]
-            self.assign_label_s(x)
+            self.extend_tree_s(x)
 
             # Assign label T to path_nodes[i+2] and attach it
             # to path_nodes[p+1].
             sub = path_nodes[p+2]
-            sub.label = _LABEL_T
+            self.assign_blossom_label_t(sub)
             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)
 
@@ -1109,7 +1168,7 @@ class _MatchingContext:
     # Labeling and alternating tree expansion:
     #
 
-    def assign_label_s(self, x: int) -> None:
+    def extend_tree_s(self, x: int) -> None:
         """Assign label S to the unlabeled blossom that contains vertex "x".
 
         If vertex "x" is matched, it is attached to the alternating tree
@@ -1125,8 +1184,7 @@ class _MatchingContext:
 
         # Assign label S to the blossom that contains vertex "x".
         bx = self.vertex_top_blossom[x]
-        assert bx.label == _LABEL_NONE
+        self.assign_blossom_label_s(bx)
-        bx.label = _LABEL_S
 
         y = self.vertex_mate[x]
         if y == -1:
@@ -1146,10 +1204,7 @@ class _MatchingContext:
             # Attach the blossom that contains "x" to the alternating tree.
             bx.tree_edge = (y, x)
 
-        # Add all vertices inside the newly labeled S-blossom to the queue.
+    def extend_tree_t(self, x: int, y: int) -> None:
-        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".
 
         Attach it to the alternating tree via edge (x, y).
@@ -1170,20 +1225,13 @@ class _MatchingContext:
             by = self.vertex_top_blossom[y]
 
         # Assign label T to the unlabeled blossom.
-        assert by.label == _LABEL_NONE
+        self.assign_blossom_label_t(by)
-        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
-        self.assign_label_s(z)
+        self.extend_tree_s(z)
 
     def add_s_to_s_edge(self, x: int, y: int) -> Optional[_AlternatingPath]:
         """Add the edge between S-vertices "x" and "y".
@@ -1231,10 +1279,10 @@ class _MatchingContext:
 
         # Process S-vertices waiting to be scanned.
         # This loop runs through O(n) iterations per stage.
-        while self.queue:
+        while self.scan_queue:
 
             # Take a vertex from the queue.
-            x = self.queue.popleft()
+            x = self.scan_queue.popleft()
 
             # Double-check that "x" is an S-vertex.
             bx = self.vertex_top_blossom[x]
@@ -1243,7 +1291,7 @@ class _MatchingContext:
             # Scan the edges that are incident on "x".
             # This loop runs through O(m) iterations per stage.
             for e in adjacent_edges[x]:
-                (p, q, _w) = edges[e]
+                (p, q, w) = edges[e]
                 y = p if p != x else q
 
                 # Consider the edge between vertices "x" and "y".
@@ -1262,11 +1310,11 @@ class _MatchingContext:
 
                 # Check whether this edge is tight (has zero slack).
                 # Only tight edges may be part of an alternating tree.
-                slack = self.edge_slack_2x(e)
+                slack = self.edge_slack_2x(x, y, bx, by, w)
                 if slack <= 0:
                     if ylabel == _LABEL_NONE:
                         # Assign label T to the blossom that contains "y".
-                        self.assign_label_t(x, y)
+                        self.extend_tree_t(x, y)
                     elif ylabel == _LABEL_S:
                         # This edge connects two S-blossoms. Use it to find
                         # either a new blossom or an augmenting path.
@@ -1276,8 +1324,17 @@ class _MatchingContext:
 
                 elif ylabel == _LABEL_S:
                     # Update tracking of least-slack edges between S-blossoms.
+                    # Priority is edge slack plus 2 times the running sum of
+                    # delta updates.
                     if e not in self.delta3_set:
-                        prio = self.edge_slack(e) + self.delta_sum_2x
+                        if self.graph.integer_weights:
+                            # If all edge weights are integers, the slack of
+                            # any edge between S-vertices is also an integer.
+                            assert slack % 2 == 0
+                            slack_1x = slack // 2
+                        else:
+                            slack_1x = slack / 2
+                        prio = slack_1x + self.delta_sum_2x
                         self.delta3_set.add(e)
                         self.delta3_queue.insert(prio, e)
 
@@ -1342,7 +1399,7 @@ class _MatchingContext:
             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)
+                slack = delta3_node.prio - self.delta_sum_2x
                 if slack <= delta_2x:
                     delta_type = 3
                     delta_2x = slack
@@ -1378,10 +1435,7 @@ class _MatchingContext:
         # Apply delta to dual variables of all vertices.
         for x in range(num_vertex):
             xlabel = self.vertex_top_blossom[x].label
-            if xlabel == _LABEL_S:
+            if xlabel == _LABEL_T:
-                # S-vertex: subtract delta from dual variable.
-                self.vertex_dual_2x[x] -= delta_2x
-            elif xlabel == _LABEL_T:
                 # T-vertex: add delta to dual variable.
                 self.vertex_dual_2x[x] += delta_2x
 
@@ -1420,12 +1474,12 @@ class _MatchingContext:
         # Assign label S to all unmatched vertices and put them in the queue.
         for x in range(num_vertex):
             if self.vertex_mate[x] == -1:
-                self.assign_label_s(x)
+                self.extend_tree_s(x)
 
         # Stop if all vertices are matched.
         # No further improvement is possible in that case.
         # This avoids messy calculations of delta steps without any S-vertex.
-        if not self.queue:
+        if not self.scan_queue:
             return False
 
         # Each pass through the following loop is a "substage".
@@ -1457,7 +1511,7 @@ class _MatchingContext:
                 (x, y, _w) = self.graph.edges[delta_edge]
                 if self.vertex_top_blossom[x].label != _LABEL_S:
                     (x, y) = (y, x)
-                self.assign_label_t(x, y)
+                self.extend_tree_t(x, y)
 
             elif delta_type == 3:
                 # Use the S-to-S edge that got unlocked by the delta update.