Use priority queues for delta2

cpp/mwmatching.h

@@ -323,7 +323,8 @@ struct Blossom
     /** Concatenable queue containing all vertices in the blossom. */
     ConcatenableQueue<WeightType, Blossom*, EdgeId> vertex_queue;
 
-    // TODO -- delta2_node
+    /** Optional node in the global delta2 queue. */
+    typename PriorityQueue<WeightType, Blossom*>::Node delta2_node;
 
     /**
      * Accumulated pending lazy updates to the dual variables of the vertices
@@ -433,8 +434,7 @@ struct NonTrivialBlossom : public Blossom<WeightType>
     WeightType dual_var;
 
     /** Top-level T-blossoms are elements in the delta4 queue. */
-    typedef PriorityQueue<WeightType, NonTrivialBlossom*> BlossomQueue;
+    typename PriorityQueue<WeightType, NonTrivialBlossom*>::Node delta4_node;
-    typename BlossomQueue::Node delta4_node;
 
     /** Initialize a non-trivial blossom. */
     NonTrivialBlossom(
@@ -539,6 +539,11 @@ public:
     using EdgeT = Edge<WeightType>;
     using BlossomT = Blossom<WeightType>;
     using NonTrivialBlossomT = NonTrivialBlossom<WeightType>;
+    using EdgeQueueT = PriorityQueue<WeightType, EdgeId>;
+
+    /** Marker value for invalid edge slack, larger than any valid slack. */
+    static constexpr WeightType INVALID_SLACK =
+        std::numeric_limits<WeightType>::max();
 
     /** Specification of a delta step. */
     struct DeltaStep
@@ -618,38 +623,39 @@ public:
     WeightType delta_sum;
 
     /** For each vertex, a node in its top-level blossom. */
-    typedef ConcatenableQueue<WeightType, BlossomT*, EdgeId> VertexQueue;
+    typedef ConcatenableQueue<WeightType, BlossomT*, EdgeId> VertexQueueT;
-    std::vector<typename VertexQueue::Node> vertex_queue_node;
+    std::vector<typename VertexQueueT::Node> vertex_queue_node;
 
-    // TODO -- delta2_queue
+    /**
+     * For each T-vertex or unlabeled vertex, a queue of its edges to
+     * any S-vertex. Priority of an edge is its modified slack.
+     */
+    std::vector<EdgeQueueT> vertex_sedge_queue;
+
+    /** For each edge, a node in an S-edge queue. */
+    std::vector<typename EdgeQueueT::Node> vertex_sedge_node;
+
+    /**
+     * Queue of unlabeled top-level blossoms that have at least one edge
+     * to an S-blossom. The priority of a blossom is its least slack to
+     * an S-blossom plus delta_sum.
+     */
+    PriorityQueue<WeightType, BlossomT*> delta2_queue;
 
     /**
      * Queue of edges between S-vertices in different top-level blossoms.
-     * Priority of an edge is its slack plus 2 * delta_sum.
+     * Priority of an edge is its modified slack.
      */
-    typedef PriorityQueue<WeightType, EdgeId> EdgeQueue;
+    EdgeQueueT delta3_queue;
-    EdgeQueue delta3_queue;
 
     /** For each edge, a node in delta3_queue. */
-    std::vector<typename EdgeQueue::Node> delta3_node;
+    std::vector<typename EdgeQueueT::Node> delta3_node;
 
     /**
      * Queue of top-level non-trivial T-blossoms.
      * Priority of a blossom is its modified dual variable.
      */
-    typedef PriorityQueue<WeightType, NonTrivialBlossomT*> BlossomQueue;
+    PriorityQueue<WeightType, NonTrivialBlossomT*> delta4_queue;
-    BlossomQueue delta4_queue;
-
-    // TODO -- vertex_sedge_queue
-    // TODO -- vertex_sedge_node
-
-// TODO -- delete
-    /**
-     * In case of T-vertex or unlabeled vertex "x",
-     * "vertex_best_edge[x]" is the least-slack edge to any S-vertex,
-     * or "nullptr" if no such edge has been found.
-     */
-    std::vector<const EdgeT*> vertex_best_edge;
 
     /** Queue of S-vertices to be scanned. */
     std::vector<VertexId> scan_queue;
@@ -669,6 +675,8 @@ public:
       : graph(edges_in),
         trivial_blossom(graph.num_vertex),
         vertex_queue_node(graph.num_vertex),
+        vertex_sedge_queue(graph.num_vertex),
+        vertex_sedge_node(edges_in.size()),
         delta3_node(edges_in.size())
     {
         // Initially all vertices are unmatched.
@@ -682,7 +690,7 @@ public:
         // Insert each vertex as the only element in its own blossom.
         for (VertexId x = 0; x < graph.num_vertex; ++x) {
             trivial_blossom[x].vertex_queue.insert(
-                &vertex_queue_node[x], 0, x);
+                &vertex_queue_node[x], INVALID_SLACK, x);
         }
 
         // Vertex duals are initialized to half the maximum edge weight.
@@ -695,9 +703,6 @@ public:
 
         delta_sum = 0;
 
-        // Initialize "vertex_best_edge".
-        vertex_best_edge.resize(graph.num_vertex, nullptr);
-
         // Allocate temporary arrays for path tracing.
         vertex_marker.resize(graph.num_vertex);
         marked_vertex.reserve(graph.num_vertex);
@@ -735,49 +740,6 @@ public:
         return vertex_dual[x] + vertex_dual[y] - weight_factor * edge.weight;
     }
 
-// TODO -- delete
-    WeightType edge_slack(const EdgeT& edge) const
-    {
-        VertexId x = edge.vt.first;
-        VertexId y = edge.vt.second;
-        BlossomT* bx = top_level_blossom(x);
-        BlossomT* by = top_level_blossom(y);
-
-        WeightType ux = vertex_dual[x];
-        if (bx->label == LABEL_S) {
-            ux -= delta_sum;
-        } else if (bx->label == LABEL_T) {
-            ux += delta_sum + bx->vertex_dual_offset;
-        } else {
-            ux += bx->vertex_dual_offset;
-        }
-
-        WeightType uy = vertex_dual[y];
-        if (by->label == LABEL_S) {
-            uy -= delta_sum;
-        } else if (by->label == LABEL_T) {
-            uy += delta_sum + by->vertex_dual_offset;
-        } else {
-            uy += by->vertex_dual_offset;
-        }
-
-        return ux + uy - weight_factor * edge.weight;
-    }
-
-    /**
-     * Reset least-slack edge tracking.
-     *
-     * This function takes time O(n + m).
-     */
-    void lset_reset()
-    {
-        for (VertexId x = 0; x < graph.num_vertex; ++x) {
-            vertex_best_edge[x] = nullptr;
-        }
-
-        delta3_queue.clear();
-    }
-
     /**
      * Add edge "e" for delta2 tracking.
      *
@@ -787,48 +749,109 @@ public:
      */
     void delta2_add_edge(EdgeId e, VertexId y, BlossomT* by)
     {
-        // TODO
+        WeightType prio = edge_pseudo_slack(e);
-        const EdgeT* edge = &graph.edges[e];
+
-        WeightType slack = edge_slack(*edge);
+        // Insert the edge in the S-edge queue of vertex "y".
-        const EdgeT* cur_best_edge = vertex_best_edge[y];
+        // Check whether this improves the min-prio of the vertex.
-        if ((cur_best_edge == nullptr)
+        assert(! vertex_sedge_node[e].valid());
-                || (slack < edge_slack(*cur_best_edge))) {
+        bool improved = vertex_sedge_queue[y].empty()
-            vertex_best_edge[y] = edge;
+                        || (prio < vertex_sedge_queue[y].min_prio());
+        vertex_sedge_queue[y].insert(&vertex_sedge_node[e], prio, e);
+
+        if (improved) {
+            // Update the priority of vertex "y" in its ConcatenableQueue.
+            vertex_queue_node[y].set_prio(prio);
+
+            // Update the global delta2 queue if this blossom is unlabeled
+            // and the new edge becomes its least-slack S-edge.
+            if (by->label == LABEL_NONE) {
+                prio += by->vertex_dual_offset;
+                if (by->delta2_node.valid()) {
+                    if (prio < by->delta2_node.prio()) {
+                        delta2_queue.set_prio(&by->delta2_node, prio);
+                    }
+                } else {
+                    delta2_queue.insert(&by->delta2_node, prio, by);
+                }
+            }
         }
     }
 
+    /**
+     * Enable delta2 tracking for "blossom".
+     *
+     * This function is called when a blossom becomes an unlabeled top-level
+     * blossom. If the blossom has at least one edge to an S-vertex,
+     * the blossom is inserted in the global delta2 queue.
+     *
+     * This function takes time O(log(n)).
+     */
+    void delta2_enable_blossom(BlossomT* blossom)
+    {
+        assert(! blossom->delta2_node.valid());
+        WeightType prio = blossom->vertex_queue.min_prio();
+        if (prio < INVALID_SLACK) {
+            prio += blossom->vertex_dual_offset;
+            delta2_queue.insert(&blossom->delta2_node, prio, blossom);
+        }
+    }
+
+    /**
+     * Disable delta2 tracking for "blossom".
+     *
+     * The blossom is removed from the global delta2 queue.
+     * This function is called when a blossom stops being an unlabeled
+     * top-level blossom.
+     *
+     * This function takes time O(log(n)).
+     */
+    void delta2_disable_blossom(BlossomT* blossom)
+    {
+        if (blossom->delta2_node.valid()) {
+            delta2_queue.remove(&blossom->delta2_node);
+        }
+    }
+
+    /**
+     * Clear delta2 tracking for vertex "x".
+     *
+     * This function is called when "x" becomes an S-vertex.
+     * It is assumed that the blossom containing "x" has already been
+     * disabled for delta2 tracking.
+     *
+     * This function takes time O(k + log(n)),
+     * where "k" is the number of edges incident on "x".
+     */
+    void delta2_clear_vertex(int x)
+    {
+        vertex_sedge_queue[x].clear();
+        vertex_queue_node[x].set_prio(INVALID_SLACK);
+    }
+
     /**
      * Find the least-slack edge between any S-vertex and any unlabeled vertex.
      *
-     * This function takes time O(n).
+     * This function takes time O(1).
-     * TODO -- rework: This function takes time O(log(n)).
      *
      * @return Tuple (edge_index, slack) if there is an S-to-unlabeled edge,
      *         or (NO_EDGE, 0) if there is no such edge.
      */
     std::tuple<EdgeId, WeightType> delta2_get_min_edge()
     {
-        const EdgeT* best_edge = nullptr;
+        if (delta2_queue.empty()) {
-        WeightType best_slack = 0;
-
-        for (VertexId x = 0; x < graph.num_vertex; ++x) {
-            if (top_level_blossom(x)->label == LABEL_NONE) {
-                const EdgeT* edge = vertex_best_edge[x];
-                if (edge != nullptr) {
-                    WeightType slack = edge_slack(*edge);
-                    if ((best_edge == nullptr) || (slack < best_slack)) {
-                        best_edge = edge;
-                        best_slack = slack;
-                    }
-                }
-            }
-        }
-
-        if (best_edge) {
-            return std::make_tuple(best_edge - graph.edges.data(), best_slack);
-        } else {
             return std::make_tuple(NO_EDGE, 0);
         }
+
+        BlossomT* blossom = delta2_queue.min_elem();
+        WeightType prio = delta2_queue.min_prio();
+        assert(! blossom->parent);
+        assert(blossom->label == LABEL_NONE);
+
+        WeightType slack = prio - delta_sum;
+        VertexId x = blossom->vertex_queue.min_elem();
+        EdgeId e = vertex_sedge_queue[x].min_elem();
+
+        return std::make_tuple(e, slack);
     }
 
     /**
@@ -917,6 +940,9 @@ public:
 
         blossom->label = LABEL_S;
 
+        // Labeled blossoms must not be in the delta2 queue.
+        delta2_disable_blossom(blossom);
+
         // Unlabeled blossoms and S-blossoms use different rules
         // for modified blossom duals. Adjust the modified dual variable
         // to preserve the true blossom dual while switching labels.
@@ -945,6 +971,9 @@ public:
                 // Apply adjustment to modified dual variable.
                 vertex_dual[x] += dual_fixup;
 
+                // S-vertices do not keep track of potential delta2 edges.
+                delta2_clear_vertex(x);
+
                 // Add new S-vertices to the scan queue.
                 scan_queue.push_back(x);
             });
@@ -962,6 +991,9 @@ public:
 
         blossom->label = LABEL_T;
 
+        // Labeled blossoms must not be in the delta2 queue.
+        delta2_disable_blossom(blossom);
+
         NonTrivialBlossomT* ntb = blossom->nontrivial();
         if (ntb) {
             // Unlabeled blossoms and T-blossoms use different rules
@@ -1042,6 +1074,9 @@ public:
 
         // Adjust dual offset to preserve true vertex duals.
         blossom->vertex_dual_offset += delta_sum;
+
+        // Enable unlabeled top-level blossom for delta2 tracking.
+        delta2_enable_blossom(blossom);
     }
 
     /**
@@ -1233,7 +1268,7 @@ public:
 
         // Merge concatenable queues.
 // TODO -- avoid temporary array
-        std::vector<VertexQueue*> subqueues;
+        std::vector<VertexQueueT*> subqueues;
         for (BlossomT* sub : subblossoms) {
             subqueues.push_back(&sub->vertex_queue);
         }
@@ -1267,6 +1302,9 @@ public:
         assert(blossom->parent == nullptr);
         assert(blossom->label == LABEL_NONE);
 
+        // Remove blossom from the delta2 queue.
+        delta2_disable_blossom(blossom);
+
         // Split concatenable queue, thus reconstructing the separate
         // concatenable queues of the sub-blossoms.
         blossom->vertex_queue.split();
@@ -1285,6 +1323,9 @@ public:
             // Push pending delta updates to sub-blossom.
             assert(sub_blossom->vertex_dual_offset == 0);
             sub_blossom->vertex_dual_offset = vertex_dual_offset;
+
+            // Add unlabeled blossom to the delta2 queue.
+            delta2_enable_blossom(sub_blossom);
         }
     }
 
@@ -1730,9 +1771,14 @@ public:
      * and the type of delta which obtain the minimum, and the edge or
      * blossom that produces the minimum delta, if applicable.
      *
-     * This function takes time O(n).
+     * This function takes time O((1 + k) * log(n)),
+     * where "k" is the number of intra-blossom edges removed from
+     * the delta3 queue.
      *
-     * @pre There is at least one S-vertex.
+     * At most O(n) delta steps can occur during a stage.
+     * Each edge can be inserted into the delta3 queue at most once per stage.
+     * Therefore, this function takes total time O((n + m) * log(n))
+     * per stage.
      *
      * @return Tuple (delta_type, delta_value, delta_edge, delta_blossom)
      */
@@ -1768,7 +1814,6 @@ public:
         }
 
         // Compute delta4: half minimum dual of a top-level T-blossom.
-        // This takes time O(1).
         if (! delta4_queue.empty()) {
             NonTrivialBlossomT* blossom = delta4_queue.min_elem();
             assert(! blossom->parent);
@@ -1806,7 +1851,16 @@ public:
         }
 
         // Reset least-slack edge tracking.
-        lset_reset();
+        for (BlossomT& blossom : trivial_blossom) {
+            delta2_disable_blossom(&blossom);
+        }
+        for (BlossomT& blossom : nontrivial_blossom) {
+            delta2_disable_blossom(&blossom);
+        }
+        for (VertexId x = 0; x < graph.num_vertex; ++x) {
+            delta2_clear_vertex(x);
+        }
+        delta3_queue.clear();
     }
 
     /**
@@ -1936,7 +1990,7 @@ public:
             cleanup_blossom(&blossom);
         }
 
-        // TODO -- check delta2_queue empty
+        assert(delta2_queue.empty());
         assert(delta3_queue.empty());
         assert(delta4_queue.empty());
     }