Implement heap-based tracking for delta3

cpp/mwmatching.h

@@ -18,6 +18,8 @@
 #include <utility>
 #include <vector>
 
+#include "datastruct.h"
+
 
 namespace mwmatching {
 
@@ -67,6 +69,11 @@ namespace impl {
 /** Value used to mark an invalid or undefined vertex. */
 constexpr VertexId NO_VERTEX = std::numeric_limits<VertexId>::max();
 
+/** Type representing an index in the edge list. */
+using EdgeIndex = unsigned int;
+
+/** Value used to mark an invalid or undefined vertex. */
+constexpr EdgeIndex NO_EDGE = std::numeric_limits<EdgeIndex>::max();
 
 /** Top-level blossoms may be labeled "S" or "T" or unlabeled. */
 enum BlossomLabel { LABEL_NONE = 0, LABEL_S = 1, LABEL_T = 2 };
@@ -99,6 +106,10 @@ void check_input_graph(const std::vector<Edge<WeightType>>& edges)
 {
     const VertexId max_num_vertex = std::numeric_limits<VertexId>::max();
 
+    if (edges.size() >= NO_EDGE) {
+        throw std::invalid_argument("Too many edges");
+    }
+
     for (const Edge<WeightType>& edge : edges) {
 
         // Check that vertex IDs are valid.
@@ -167,7 +178,7 @@ struct Graph
     const VertexId num_vertex;
 
     /** For each vertex, a vector of pointers to its incident edges. */
-    const std::vector<std::vector<const EdgeT*>> adjacent_edges;
+    const std::vector<std::vector<EdgeIndex>> adjacent_edges;
 
     /**
      * Initialize the graph representation and prepare adjacent edge lists.
@@ -233,7 +244,7 @@ struct Graph
      * @param  edges       List of edges in the graph.
      * @return Vector of incident edges for each vertex.
      */
-    static std::vector<std::vector<const EdgeT*>> build_adjacent_edges(
+    static std::vector<std::vector<EdgeIndex>> build_adjacent_edges(
         const std::vector<EdgeT>& edges,
         VertexId num_vertex)
     {
@@ -247,13 +258,14 @@ struct Graph
         }
 
         // Build adjacency lists.
-        std::vector<std::vector<const EdgeT*>> adjacent_edges(num_vertex);
+        std::vector<std::vector<EdgeIndex>> adjacent_edges(num_vertex);
         for (VertexId i = 0; i < num_vertex; ++i) {
             adjacent_edges[i].reserve(vertex_degree[i]);
         }
-        for (const EdgeT& edge : edges) {
+        for (EdgeIndex e = 0; e < edges.size(); e++) {
-            adjacent_edges[edge.vt.first].push_back(&edge);
+            const EdgeT& edge = edges[e];
-            adjacent_edges[edge.vt.second].push_back(&edge);
+            adjacent_edges[edge.vt.first].push_back(e);
+            adjacent_edges[edge.vt.second].push_back(e);
         }
 
         return adjacent_edges;
@@ -306,6 +318,17 @@ struct Blossom
     /** Optional edge that attaches this blossom to the alternating tree. */
     VertexPair tree_edge;
 
+    // TODO -- tree_blossoms
+
+    // TOOD -- vertex_queue
+
+    // TODO -- delta2_node
+
+    // TODO -- vertex_dual_offset
+
+    // TODO -- marker
+
+// TODO -- remove
     /**
      * In case of a top-level S-blossom, "best_edge" is the least-slack edge
      * that links to a different S-blossom, or "nullptr" if no such edge
@@ -393,14 +416,11 @@ struct NonTrivialBlossom : public Blossom<WeightType>
      */
     std::list<SubBlossom> subblossoms;
 
+// TODO -- description
     /** Dual LPP variable for this blossom. */
     WeightType dual_var;
 
-    /**
+    // TODO -- delta4_node
-     * In case of a top-level S-blossom, "best_edge_set" is a list of
-     * least-slack edges between this blossom and other S-blossoms.
-     */
-    std::list<const Edge<WeightType>*> best_edge_set;
 
     /** Initialize a non-trivial blossom. */
     NonTrivialBlossom(
@@ -552,6 +572,9 @@ public:
 // NOTE - this MUST be a list, because we delete items from it while keeping pointers to other items
     std::list<NonTrivialBlossomT> nontrivial_blossom;
 
+    // TODO -- vertex_queue_node
+
+// TODO -- remove
     /**
      * Every vertex is contained in exactly one top-level blossom
      * (possibly the trivial blossom that contains just that vertex).
@@ -561,6 +584,9 @@ public:
      */
     std::vector<BlossomT*> vertex_top_blossom;
 
+    // TODO -- start_vertex_dual
+
+// TODO -- description
     /**
      * Every vertex has a variable in the dual LPP.
      *
@@ -568,6 +594,24 @@ public:
      */
     std::vector<WeightType> vertex_dual;
 
+    /** Running sum of applied delta steps. */
+    WeightType delta_sum;
+
+    // TODO -- delta2_queue
+
+    /** Queue of edges between S-vertices in different top-level blossoms. */
+    typedef PriorityQueue<WeightType, EdgeIndex> EdgeQueue;
+    EdgeQueue delta3_queue;
+
+    /** For each edge, a node in delta3_queue. */
+    std::vector<typename EdgeQueue::Node> delta3_node;
+
+    // TODO -- 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,
@@ -575,12 +619,15 @@ public:
      */
     std::vector<const EdgeT*> vertex_best_edge;
 
+// TODO -- rename to scan_queue
     /** Queue of S-vertices to be scanned. */
     std::deque<VertexId> queue;
 
+// TODO -- delete
     /** Markers placed while tracing an alternating path. */
     std::vector<bool> vertex_marker;
 
+// TODO -- delete
     /** Vertices marked while tracing an alternating path. */
     std::vector<VertexId> marked_vertex;
 
@@ -590,7 +637,8 @@ public:
      * This function takes time O(n + m).
      */
     explicit MatchingContext(const std::vector<EdgeT>& edges_in)
-      : graph(edges_in)
+      : graph(edges_in),
+        delta3_node(edges_in.size())
     {
         // Initially all vertices are unmatched.
         vertex_mate.resize(graph.num_vertex, NO_VERTEX);
@@ -615,6 +663,8 @@ public:
         WeightType init_vertex_dual = max_weight * (weight_factor / 2);
         vertex_dual.resize(graph.num_vertex, init_vertex_dual);
 
+        delta_sum = 0;
+
         // Initialize "vertex_best_edge".
         vertex_best_edge.resize(graph.num_vertex, nullptr);
 
@@ -627,48 +677,44 @@ public:
     MatchingContext(const MatchingContext&) = delete;
     MatchingContext& operator=(const MatchingContext&) = delete;
 
-    /* **********  General support routines:  ********** */
+    /* **********  Find top-level blossom:  ********** */
 
-    /** Calculate edge slack. */
+    /**
+     * Find the top-level blossom that contains vertex "x".
+     *
+     * This function takes time O(log(n)).
+     */
+    BlossomT* top_level_blossom(VertexId x)
+    {
+        // TODO
+        return vertex_top_blossom[x];
+    }
+
+    /* **********  Least slack edge tracking:  ********** */
+
+    /**
+     * Return the modified edge slack of the specified edge.
+     *
+     * Modified edge slack is related to true edge slack, but adjusted
+     * to make it invariant under delta steps.
+     */
+    WeightType edge_pseudo_slack(EdgeIndex e) const
+    {
+        const EdgeT& edge = graph.edges[e];
+        VertexId x = edge.vt.first;
+        VertexId y = edge.vt.second;
+        // TODO -- remove delta_sum here
+        return vertex_dual[x] + vertex_dual[y] - weight_factor * edge.weight + 2 * delta_sum;
+    }
+
+// TODO -- delete
     WeightType edge_slack(const EdgeT& edge) const
     {
         VertexId x = edge.vt.first;
         VertexId y = edge.vt.second;
-        assert(vertex_top_blossom[x] != vertex_top_blossom[y]);
         return vertex_dual[x] + vertex_dual[y] - weight_factor * edge.weight;
     }
 
-    /*
-     * Least-slack edge tracking:
-     *
-     * To calculate delta steps, the matching algorithm needs to find
-     *  - the least-slack edge between any S-vertex and an unlabeled vertex;
-     *  - the least-slack edge between any pair of top-level S-blossoms.
-     *
-     * For each unlabeled vertex and each T-vertex, we keep track of the
-     * least-slack edge to any S-vertex. Tracking for unlabeled vertices
-     * serves to provide the least-slack edge for the delta step.
-     * 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
-     * of every edge to an S-vertex by the same amount. Therefore the edge
-     * that had least slack before the delta step, will still have least slack
-     * after the delta step.
-     */
-
     /**
      * Reset least-slack edge tracking.
      *
@@ -686,8 +732,9 @@ public:
 
         for (NonTrivialBlossomT& blossom : nontrivial_blossom) {
             blossom.best_edge = nullptr;
-            blossom.best_edge_set.clear();
         }
+
+        delta3_queue.clear();
     }
 
     /**
@@ -736,162 +783,72 @@ public:
         return std::make_pair(best_edge, best_slack);
     }
 
-    /** Start tracking edges for a new S-blossom. */
-    void lset_new_blossom(BlossomT* blossom)
-    {
-        assert(blossom->best_edge == nullptr);
-        assert((blossom->nontrivial() == nullptr)
-               || blossom->nontrivial()->best_edge_set.empty());
-    }
-
     /**
-     * Add edge "e" between the specified S-blossom and another S-blossom.
+     * Add specified edge for delta3 tracking.
      *
-     * This function takes time O(1) per call.
+     * This function is called if a vertex becomes an S-vertex and edge "e"
-     * This function is called O(m) times per stage.
+     * connects it to an S-vertex in a different top-level blossom.
+     *
+     * This function takes time O(log(n)).
      */
-    void lset_add_blossom_edge(BlossomT* blossom,
+    void delta3_add_edge(EdgeIndex e)
-                               const EdgeT* edge,
-                               WeightType slack)
     {
-        const EdgeT* cur_best_edge = blossom->best_edge;
+        // The edge may already be in the delta3 queue, if it was previously
-        if ((cur_best_edge == nullptr)
+        // iscovered in the opposite direction. In that case do nothing.
-                || (slack < edge_slack(*cur_best_edge))) {
+        if (! delta3_node[e].valid()) {
-            blossom->best_edge = edge;
+            // Insert edge. Use modified edge slack as priority.
-        }
+            WeightType prio = edge_pseudo_slack(e);
-
+            delta3_queue.insert(&delta3_node[e], prio, e);
-        NonTrivialBlossomT* ntb = blossom->nontrivial();
-        if (ntb) {
-            ntb->best_edge_set.push_back(edge);
         }
     }
 
     /**
-     * Update least-slack edge tracking after merging sub-blossoms
+     * Remove specified edge from delta3 tracking.
-     * into a new S-blossom.
      *
-     * This function takes total time O(n**2) per stage.
+     * This function is called if a former S-vertex becomes unlabeled
+     * and edge "e" connects it to another S-vertex.
+     *
+     * This function takes time O(log(n)).
      */
-    void lset_merge_blossoms(NonTrivialBlossomT* blossom)
+    void delta3_remove_edge(EdgeIndex e)
     {
-        assert(blossom->best_edge == nullptr);
+        if (delta3_node[e].valid()) {
-        assert(blossom->best_edge_set.empty());
+            delta3_queue.remove(&delta3_node[e]);
-
-        // Collect edges from the sub-blossoms that used to be S-blossoms.
-        for (auto& subblossom_item : blossom->subblossoms) {
-            BlossomT* sub = subblossom_item.blossom;
-            if (sub->label == LABEL_S) {
-                NonTrivialBlossomT* ntb = sub->nontrivial();
-                if (ntb) {
-                    // Take least-slack edge set from this subblossom.
-                    blossom->best_edge_set.splice(
-                        blossom->best_edge_set.end(),
-                        ntb->best_edge_set);
-                } else {
-                    // Trivial blossoms don't maintain a least-slack edge set.
-                    // Just consider all incident edges.
-                    for (const EdgeT* edge : graph.adjacent_edges[sub->base_vertex]) {
-                        // Only take edges between different S-blossoms.
-                        VertexId x = edge->vt.first;
-                        VertexId y = edge->vt.second;
-                        BlossomT* bx = vertex_top_blossom[x];
-                        BlossomT* by = vertex_top_blossom[y];
-                        if ((bx != by)
-                                && (bx->label == LABEL_S)
-                                && (by->label == LABEL_S)) {
-                            blossom->best_edge_set.push_back(edge);
-                        }
-                    }
-                }
-            }
         }
+    }
 
-        // Build a temporary array holding the least-slack edge index to
+    /**
-        // each top-level S-blossom. This array is indexed by the base vertex
+     * Find the least-slack edge between any pair of S-vertices
-        // of the blossoms.
+     * in different top-level blossoms.
-        std::vector<std::pair<const EdgeT*, WeightType>>
+     *
-            best_edge_to_blossom(graph.num_vertex);
+     * This function takes time O(1 + k * log(n)),
-
+     * where "k" is the number of intra-blossom edges removed from the queue.
-        for (const EdgeT* edge : blossom->best_edge_set) {
+     *
-            BlossomT* bx = vertex_top_blossom[edge->vt.first];
+     * @return Tuple (edge_index, slack) if there is an S-to-S edge,
-            BlossomT* by = vertex_top_blossom[edge->vt.second];
+     *         or (NO_EDGE, 0) if there is no suitable edge.
-            assert(bx == blossom || by == blossom);
+     */
-
+    std::tuple<EdgeIndex, WeightType> delta3_get_min_edge()
-            // Ignore internal edges.
+    {
+        while (! delta3_queue.empty()) {
+            EdgeIndex e = delta3_queue.min_elem();
+            const EdgeT& edge = graph.edges[e];
+            BlossomT* bx = top_level_blossom(edge.vt.first);
+            BlossomT* by = top_level_blossom(edge.vt.second);
+            assert(bx->label == LABEL_S && by->label == LABEL_S);
             if (bx != by) {
-                bx = (bx == blossom) ? by : bx;
+                WeightType slack = delta3_node[e].prio() - 2 * delta_sum;
-
+                return std::make_tuple(e, slack);
-                // Only consider edges to S-blossoms.
-                if (bx->label == LABEL_S) {
-
-                    // Keep only the least-slack edge to blossom "bx".
-                    WeightType slack = edge_slack(*edge);
-                    VertexId bx_base = bx->base_vertex;
-                    auto& best_edge_item = best_edge_to_blossom[bx_base];
-                    if ((best_edge_item.first == nullptr)
-                            || (slack < best_edge_item.second)) {
-                        best_edge_item.first = edge;
-                        best_edge_item.second = slack;
-                    }
-                }
             }
-        };
 
-        // Rebuild a compact list of least-slack edges.
+            // Reject edges between vertices in the same top-level blossom.
-        // Also find the overall least-slack edge to any other S-blossom.
+            // Although such edges are never inserted into the queue,
-        WeightType best_slack = 0;
+            // existing edges in the queue may become intra-blossom when
-        blossom->best_edge_set.clear();
+            // blossoms are merged.
-        for (auto& best_edge_item : best_edge_to_blossom) {
+            delta3_queue.remove(&delta3_node[e]);
-            const EdgeT* edge = best_edge_item.first;
-            if (edge != nullptr) {
-                blossom->best_edge_set.push_back(edge);
-                WeightType slack = best_edge_item.second;
-                if ((blossom->best_edge == nullptr) || (slack < best_slack)) {
-                    blossom->best_edge = edge;
-                    best_slack = slack;
-                }
-            }
-        }
-    }
-
-    /**
-     * 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.
-     *
-     * @return Pair (edge, slack) if there is a least-slack edge,
-     *         or (nullptr, 0) if there is no suitable edge.
-     */
-    std::pair<const EdgeT*, WeightType> lset_get_best_blossom_edge()
-    {
-        const EdgeT* best_edge = nullptr;
-        WeightType best_slack = 0;
-
-        auto consider_blossom =
-            [this,&best_edge,&best_slack](BlossomT* blossom)
-        {
-            if ((blossom->parent == nullptr) && (blossom->label == LABEL_S)) {
-                const EdgeT* edge = blossom->best_edge;
-                if (edge != nullptr) {
-                    WeightType slack = edge_slack(*edge);
-                    if ((best_edge == nullptr) || (slack < best_slack)) {
-                        best_edge = edge;
-                        best_slack = slack;
-                    }
-                }
-            }
-        };
-
-        for (BlossomT& blossom : trivial_blossom) {
-            consider_blossom(&blossom);
-        }
-        for (BlossomT& blossom : nontrivial_blossom) {
-            consider_blossom(&blossom);
         }
 
-        return std::make_pair(best_edge, best_slack);
+        // Queue empty; no suitable edge exists.
+        return std::make_tuple(NO_EDGE, 0);
     }
 
     /* **********  Creating and expanding blossoms:  ********** */
@@ -1052,9 +1009,6 @@ public:
                     });
             }
         }
-
-        // Merge least-slack edges for the S-sub-blossoms.
-        lset_merge_blossoms(blossom);
     }
 
     /** Erase the specified non-trivial blossom. */
@@ -1386,9 +1340,6 @@ public:
             bx->tree_edge = std::make_pair(y, x);
         }
 
-        // Start least-slack edge tracking for the S-blossom.
-        lset_new_blossom(bx);
-
         // Add all vertices inside the newly labeled S-blossom to the queue.
         for_vertices_in_blossom(bx,
             [this](VertexId v) {
@@ -1495,9 +1446,10 @@ public:
 
             // Scan the edges that are incident on "x".
             // This loop runs through O(m) iterations per stage.
-            for (const EdgeT* edge : graph.adjacent_edges[x]) {
+            for (EdgeIndex e : graph.adjacent_edges[x]) {
-                VertexId y = (edge->vt.first != x) ? edge->vt.first
+                const EdgeT& edge = graph.edges[e];
-                                                   : edge->vt.second;
+                VertexId y = (edge.vt.first != x) ? edge.vt.first
+                                                  : edge.vt.second;
 
                 // Note: The top-level blossom of vertex "x" may change
                 // during this loop, so we need to refresh it in each pass.
@@ -1512,7 +1464,7 @@ public:
 
                 // Check whether this edge is tight (has zero slack).
                 // Only tight edges may be part of an alternating tree.
-                WeightType slack = edge_slack(*edge);
+                WeightType slack = edge_slack(edge);
                 if (slack <= 0) {
                     if (ylabel == LABEL_NONE) {
                         // Found a tight edge to an unlabeled blossom.
@@ -1529,14 +1481,14 @@ public:
                 } else if (ylabel == LABEL_S) {
                     // Found a non-tight edge between two S-blossoms.
                     // Pass it to the least-slack edge tracker.
-                    lset_add_blossom_edge(bx, edge, slack);
+                    delta3_add_edge(e);
                 }
 
                 if (ylabel != LABEL_S) {
                     // Found an to a T-vertex or unlabeled vertex "y".
                     // Pass it to the least-slack edge tracker.
                     // Tight edges must also be tracked in this way.
-                    lset_add_vertex_edge(y, edge, slack);
+                    lset_add_vertex_edge(y, &edge, slack);
                 }
             }
         }
@@ -1587,11 +1539,12 @@ public:
 
         // Compute delta3: half minimum slack of any edge between two
         // top-level S-blossoms.
-        std::tie(edge, slack) = lset_get_best_blossom_edge();
+        EdgeIndex e;
-        if ((edge != nullptr) && (slack / 2 <= delta.value)) {
+        std::tie(e, slack) = delta3_get_min_edge();
+        if ((e != NO_EDGE) && (slack / 2 <= delta.value)) {
             delta.kind = 3;
             delta.value = slack / 2;
-            delta.edge = edge->vt;
+            delta.edge = graph.edges[e].vt;
         }
 
         // Compute delta4: half minimum dual of a top-level T-blossom.
@@ -1635,6 +1588,8 @@ public:
                 }
             }
         }
+
+        delta_sum += delta;
     }
 
     /* **********  Main algorithm:  ********** */
@@ -1962,7 +1917,8 @@ private:
                     // For each incident edge, find the smallest blossom
                     // that contains it.
                     VertexId x = sub->base_vertex;
-                    for (const EdgeT* edge : graph.adjacent_edges[x]) {
+                    for (EdgeIndex e : graph.adjacent_edges[x]) {
+                        const EdgeT* edge = &graph.edges[e];
                         // Only consider edges pointing out from "x".
                         if (edge->vt.first == x) {
                             VertexId y = edge->vt.second;