Lazy updates of vertex duals

cpp/mwmatching.h

@@ -324,15 +324,11 @@ struct Blossom
 
     // TODO -- delta2_node
 
-    // TODO -- vertex_dual_offset
-
-// TODO -- remove
     /**
-     * In case of a top-level S-blossom, "best_edge" is the least-slack edge
+     * Accumulated pending lazy updates to the dual variables of the vertices
-     * that links to a different S-blossom, or "nullptr" if no such edge
+     * inside the blossom.
-     * has been found.
      */
-    const Edge<WeightType>* best_edge;
+    WeightType vertex_dual_offset;
 
 protected:
     /** Initialize base class. */
@@ -341,7 +337,7 @@ protected:
         base_vertex(base_vertex),
         label(LABEL_NONE),
         is_nontrivial_blossom(is_nontrivial_blossom),
-        best_edge(nullptr)
+        vertex_dual_offset(0)
     { }
 
 public:
@@ -582,16 +578,34 @@ public:
      */
     std::vector<BlossomT*> vertex_top_blossom;
 
-    // TODO -- start_vertex_dual
-
-// TODO -- description
     /**
-     * Every vertex has a variable in the dual LPP.
+     * Modified dual variable of each vertex.
      *
-     * "vertex_dual[x]" is the dual variable of vertex "x".
+     * Every vertex has a variable in the dual LPP. The true value of the dual
+     * variable changes through delta steps, but the modified dual variables
+     * are invariant under delta steps.
+     *
+     * For an S-vertex "x":
+     *   vertex_dual[x] = u(x) + delta_sum
+     *
+     * For a T-vertex "x":
+     *   vertex_dual[x] = u(x) - delta_sum - B(x).vertex_dual_offset
+     *
+     * For an unlabeled vertex:
+     *   vertex_dual[x] = u(x) - B(x).vertex_dual_offset
+     *
+     * where u(x) is the true dual variable of vertex "x"
+     * and B(x) is the top-level blossom that contains vertex "x".
      */
     std::vector<WeightType> vertex_dual;
 
+    /**
+     * Initial value of all vertex dual variables.
+     *
+     * This is equal to half of the maximum edge weight.
+     */
+    WeightType init_vertex_dual;
+
     /** Running sum of applied delta steps. */
     WeightType delta_sum;
 
@@ -658,7 +672,7 @@ public:
         for (const EdgeT& edge : graph.edges) {
             max_weight = std::max(max_weight, edge.weight);
         }
-        WeightType init_vertex_dual = max_weight * (weight_factor / 2);
+        init_vertex_dual = max_weight * (weight_factor / 2);
         vertex_dual.resize(graph.num_vertex, init_vertex_dual);
 
         delta_sum = 0;
@@ -682,7 +696,7 @@ public:
      *
      * This function takes time O(log(n)).
      */
-    BlossomT* top_level_blossom(VertexId x)
+    BlossomT* top_level_blossom(VertexId x) const
     {
         // TODO
         return vertex_top_blossom[x];
@@ -701,8 +715,7 @@ public:
         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;
-        return vertex_dual[x] + vertex_dual[y] - weight_factor * edge.weight + 2 * delta_sum;
     }
 
 // TODO -- delete
@@ -710,7 +723,28 @@ public:
     {
         VertexId x = edge.vt.first;
         VertexId y = edge.vt.second;
-        return vertex_dual[x] + vertex_dual[y] - weight_factor * edge.weight;
+        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;
     }
 
     /**
@@ -724,14 +758,6 @@ public:
             vertex_best_edge[x] = nullptr;
         }
 
-        for (BlossomT& blossom : trivial_blossom) {
-            blossom.best_edge = nullptr;
-        }
-
-        for (NonTrivialBlossomT& blossom : nontrivial_blossom) {
-            blossom.best_edge = nullptr;
-        }
-
         delta3_queue.clear();
     }
 
@@ -874,9 +900,27 @@ public:
 
         blossom->label = LABEL_S;
 
-        // Add new S-vertices to the scan queue.
+        // Unlabeled vertices and S-vertices use different rules for
+        // modified vertex duals. Calculate the adjustment that must be
+        // applied to modified vertex duals to preserve the true vertex duals
+        // while switching labels.
+        //
+        // Unlabeled vertex: vertex_dual[x] = u(x) - B(x).vertex_dual_offset
+        // S-vertex:         vertex_dual[x] = u(x) + delta_sum
+        //
+        // For S-blossoms, "vertex_dual_offset" is always 0.
+        //
+        WeightType dual_fixup = delta_sum + blossom->vertex_dual_offset;
+        blossom->vertex_dual_offset = 0;
+
+        // Loop over newly labeled S-vertices.
         for_vertices_in_blossom(blossom,
-            [this](VertexId x) {
+            [this,dual_fixup](VertexId x) {
+
+                // Apply adjustment to modified dual variable.
+                vertex_dual[x] += dual_fixup;
+
+                // Add new S-vertices to the scan queue.
                 scan_queue.push_back(x);
             });
     }
@@ -892,6 +936,45 @@ public:
         assert(blossom->label == LABEL_NONE);
 
         blossom->label = LABEL_T;
+
+        // Unlabeled vertices and T-vertices use different rules for
+        // modified vertex duals. Adjust the dual offset to preserve the
+        // true vertex duals while switching labels.
+        //
+        // Unlabeled vertex:
+        //   vertex_dual[x] = u(x) - B(x).vertex_dual_offset
+        //
+        // T-vertex:
+        //   vertex_dual[x] = u(x) - delta_sum - B(x).vertex_dual_offset
+        //
+        blossom->vertex_dual_offset -= delta_sum;
+    }
+
+    /**
+     * Change a top-level S-blossom into an unlabeled blossom.
+     *
+     * For a blossom with "j" vertices and "k" incident edges,
+     * this function takes time O((j + k) * log(n)).
+     *
+     * This function is called at most once per blossom per stage.
+     * It therefore takes total time O((n + m) * log(n)) per stage.
+     */
+    void remove_blossom_label_s(BlossomT* blossom)
+    {
+        assert(! blossom->parent);
+        assert(blossom->label == LABEL_S);
+
+        blossom->label = LABEL_NONE;
+
+        // Unlabeled vertices and S-vertices use different rules for
+        // modified vertex duals. Adjust the modified vertex duals
+        // match the true vertex duals.
+        assert(blossom->vertex_dual_offset == 0);
+        WeightType dual_fixup = -delta_sum;
+        for_vertices_in_blossom(blossom,
+            [this,dual_fixup](VertexId x) {
+                vertex_dual[x] += dual_fixup;
+            });
     }
 
     /**
@@ -905,6 +988,23 @@ public:
         assert(blossom->label == LABEL_T);
 
         blossom->label = LABEL_NONE;
+
+        // Unlabeled vertices and T-vertices use different rules for
+        // modified vertex duals. Adjust the dual offset to preserve the
+        // true vertex duals while switching labels.
+        blossom->vertex_dual_offset += delta_sum;
+    }
+
+    /** Remove blossom label. */
+    void reset_blossom_label(BlossomT* blossom)
+    {
+        if (! blossom->parent) {
+            if (blossom->label == LABEL_S) {
+                remove_blossom_label_s(blossom);
+            } else if (blossom->label == LABEL_T) {
+                remove_blossom_label_t(blossom);
+            }
+        }
     }
 
     /**
@@ -1096,6 +1196,48 @@ public:
         nontrivial_blossom.erase(blossom_it);
     }
 
+    /**
+     * Expand the specified unlabeled blossom but do not yet delete it.
+     *
+     * This function takes time O(n).
+     */
+    void expand_unlabeled_blossom_core(NonTrivialBlossomT* blossom)
+    {
+        assert(blossom->parent == nullptr);
+        assert(blossom->label == LABEL_NONE);
+
+        // Prepare to push pending delta updates down to the sub-blossoms.
+        WeightType vertex_dual_offset = blossom->vertex_dual_offset;
+        blossom->vertex_dual_offset = 0;
+
+        // Convert sub-blossoms into top-level blossoms.
+        for (const auto& sub : blossom->subblossoms) {
+            BlossomT* sub_blossom = sub.blossom;
+            assert(sub_blossom->parent == blossom);
+            assert(sub_blossom->label == LABEL_NONE);
+            sub_blossom->parent = nullptr;
+            for_vertices_in_blossom(sub_blossom,
+                [this,sub_blossom](VertexId x) {
+                    vertex_top_blossom[x] = sub_blossom;
+                });
+
+            // Push pending delta updates to sub-blossom.
+            assert(sub_blossom->vertex_dual_offset == 0);
+            sub_blossom->vertex_dual_offset = vertex_dual_offset;
+        }
+    }
+
+    /**
+     * Expand and delete the specified unlabeled blossom.
+     *
+     * This function takes time O(n).
+     */
+    void expand_unlabeled_blossom(NonTrivialBlossomT* blossom)
+    {
+        expand_unlabeled_blossom_core(blossom);
+        erase_nontrivial_blossom(blossom);
+    }
+
     /**
      * Expand the specified T-blossom.
      *
@@ -1109,17 +1251,8 @@ public:
         // Remove label from blossom.
         remove_blossom_label_t(blossom);
 
-        // Convert sub-blossoms into top-level blossoms.
+        // Expand the unlabeled blossom.
-        for (const auto& sub : blossom->subblossoms) {
+        expand_unlabeled_blossom_core(blossom);
-            BlossomT* sub_blossom = sub.blossom;
-            assert(sub_blossom->parent == blossom);
-            assert(sub_blossom->label == LABEL_NONE);
-            sub_blossom->parent = nullptr;
-            for_vertices_in_blossom(sub_blossom,
-                [this,sub_blossom](VertexId x) {
-                    vertex_top_blossom[x] = sub_blossom;
-                });
-        }
 
         // The expanded blossom was part of an alternating tree.
         // We must now reconstruct the part of the alternating tree
@@ -1185,32 +1318,6 @@ public:
         erase_nontrivial_blossom(blossom);
     }
 
-    /**
-     * Expand the specified unlabeled blossom.
-     *
-     * This function takes time O(n).
-     */
-    void expand_unlabeled_blossom(NonTrivialBlossomT* blossom)
-    {
-        assert(blossom->parent == nullptr);
-        assert(blossom->label == LABEL_NONE);
-
-        // Convert sub-blossoms into top-level blossoms.
-        for (const auto& sub : blossom->subblossoms) {
-            BlossomT* sub_blossom = sub.blossom;
-            assert(sub_blossom->parent == blossom);
-            assert(sub_blossom->label == LABEL_NONE);
-            sub_blossom->parent = nullptr;
-            for_vertices_in_blossom(sub_blossom,
-                [this,sub_blossom](VertexId x) {
-                    vertex_top_blossom[x] = sub_blossom;
-                });
-        }
-
-        // Delete the expanded blossom.
-        erase_nontrivial_blossom(blossom);
-    }
-
     /* **********  Augmenting:  ********** */
 
     /**
@@ -1573,13 +1680,10 @@ public:
         delta.blossom = nullptr;
 
         // 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.kind = 1;
-        delta.value = std::numeric_limits<WeightType>::max();
+        delta.value = init_vertex_dual - delta_sum;
-        for (VertexId x = 0; x < graph.num_vertex; ++x) {
-            if (vertex_top_blossom[x]->label == LABEL_S) {
-                delta.value = std::min(delta.value, vertex_dual[x]);
-            }
-        }
 
         // Compute delta2: minimum slack of any edge between an S-vertex and
         // an unlabeled vertex.
@@ -1618,18 +1722,6 @@ public:
     /** Apply a delta step to the dual LPP variables. */
     void substage_apply_delta_step(WeightType delta)
     {
-        // Apply delta to dual variables of all vertices.
-        for (VertexId x = 0; x < graph.num_vertex; ++x) {
-            BlossomLabel xlabel = vertex_top_blossom[x]->label;
-            if (xlabel == LABEL_S) {
-                // S-vertex: subtract delta from dual variable.
-                vertex_dual[x] -= delta;
-            } else if (xlabel == LABEL_T) {
-                // T-vertex: add delta to dual variable.
-                vertex_dual[x] += delta;
-            }
-        }
-
         // Apply delta to dual variables of top-level non-trivial blossoms.
         for (NonTrivialBlossomT& blossom : nontrivial_blossom) {
             if (blossom.parent == nullptr) {
@@ -1658,10 +1750,10 @@ public:
 
         // Remove blossom labels.
         for (BlossomT& blossom : trivial_blossom) {
-            blossom.label = LABEL_NONE;
+            reset_blossom_label(&blossom);
         }
         for (BlossomT& blossom : nontrivial_blossom) {
-            blossom.label = LABEL_NONE;
+            reset_blossom_label(&blossom);
         }
 
         // Reset least-slack edge tracking.
@@ -1766,6 +1858,42 @@ public:
         return augmented;
     }
 
+    /**
+     * Remove alternating trees and apply lazy updates to dual variables.
+     *
+     * This function takes time O((n + m) * log(n)).
+     * It is called once, at the end of the algorithm.
+     */
+    void cleanup()
+    {
+        assert(scan_queue.empty());
+
+        auto cleanup_blossom = [this](BlossomT* blossom) {
+            assert(blossom->label == LABEL_NONE);
+
+            // Unwind lazy delta updates to vertex dual variables.
+            if (blossom->vertex_dual_offset != 0) {
+                WeightType dual_fixup = blossom->vertex_dual_offset;
+                blossom->vertex_dual_offset = 0;
+                for_vertices_in_blossom(blossom,
+                    [this,dual_fixup](VertexId x) {
+                        vertex_dual[x] += dual_fixup;
+                    });
+            }
+        };
+
+        for (BlossomT& blossom : trivial_blossom) {
+            cleanup_blossom(&blossom);
+        }
+        for (BlossomT& blossom : nontrivial_blossom) {
+            cleanup_blossom(&blossom);
+        }
+
+        // TODO -- check delta2_queue empty
+        assert(delta3_queue.empty());
+        // TODO -- check delta4_queue empty
+    }
+
     /** Run the matching algorithm. */
     void run()
     {
@@ -1777,6 +1905,9 @@ public:
         // This loop runs through at most (n/2 + 1) iterations.
         // Each iteration takes time O(n**2).
         while (run_stage()) ;
+
+        // Clean up and unwind lazy updates to dual variables.
+        cleanup();
     }
 };