Lazy updates of blossom duals

cpp/mwmatching.h

@@ -410,8 +410,24 @@ struct NonTrivialBlossom : public Blossom<WeightType>
      */
     std::list<SubBlossom> subblossoms;
 
-// TODO -- description
+    /**
-    /** Dual LPP variable for this blossom. */
+     * Modified dual variable of this blossom.
+     *
+     * Every non-trivial blossom has a variable in the dual LPP. The true
+     * value of the dual variable changes through delta steps, but the
+     * modified dual variable is invariant under delta steps.
+     *
+     * For an S-blossom "b":
+     *   b.dual_var = z(b) - 2 * delta_sum
+     *
+     * For a T-blossom "b":
+     *   b.dual_var = z(b) + 2 * delta_sum
+     *
+     * For an unlabeled blossom or non-top-level blossom "b":
+     *   b.dual_var = z(b)
+     *
+     * where z(b) is the true dual variable of blossom "b".
+     */
     WeightType dual_var;
 
     // TODO -- delta4_node
@@ -900,6 +916,14 @@ public:
 
         blossom->label = LABEL_S;
 
+        // 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.
+        NonTrivialBlossomT* ntb = blossom->nontrivial();
+        if (ntb) {
+            ntb->dual_var -= 2 * delta_sum;
+        }
+
         // 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
@@ -937,6 +961,14 @@ public:
 
         blossom->label = LABEL_T;
 
+        // Unlabeled blossoms and T-blossoms use different rules
+        // for modified blossom duals. Adjust the modified dual variable
+        // to preserve the true blossom dual while switching labels.
+        NonTrivialBlossomT* ntb = blossom->nontrivial();
+        if (ntb) {
+            ntb->dual_var += 2 * delta_sum;
+        }
+
         // 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.
@@ -966,9 +998,13 @@ public:
 
         blossom->label = LABEL_NONE;
 
-        // Unlabeled vertices and S-vertices use different rules for
+        // Adjust modified blossom dual to preserve true blossom dual.
-        // modified vertex duals. Adjust the modified vertex duals
+        NonTrivialBlossomT* ntb = blossom->nontrivial();
-        // match the true vertex duals.
+        if (ntb) {
+            ntb->dual_var += 2 * delta_sum;
+        }
+
+        // Adjust modified vertex duals to preserve true vertex dual.
         assert(blossom->vertex_dual_offset == 0);
         WeightType dual_fixup = -delta_sum;
         for_vertices_in_blossom(blossom,
@@ -989,22 +1025,14 @@ public:
 
         blossom->label = LABEL_NONE;
 
-        // Unlabeled vertices and T-vertices use different rules for
+        // Adjust modified blossom dual to preserve true blossom dual.
-        // modified vertex duals. Adjust the dual offset to preserve the
+        NonTrivialBlossomT* ntb = blossom->nontrivial();
-        // true vertex duals while switching labels.
+        if (ntb) {
-        blossom->vertex_dual_offset += delta_sum;
+            ntb->dual_var -= 2 * 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);
-            }
         }
+
+        // Adjust dual offset to preserve true vertex duals.
+        blossom->vertex_dual_offset += delta_sum;
     }
 
     /**
@@ -1018,6 +1046,25 @@ public:
         assert(blossom->label == LABEL_S);
 
         blossom->label = LABEL_NONE;
+
+        // The subblossom becomes unlabeled.
+        // Adjust its modified blossom dual to preserve the true blossom dual.
+        NonTrivialBlossomT* ntb = blossom->nontrivial();
+        if (ntb) {
+            ntb->dual_var += 2 * 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);
+            }
+        }
     }
 
     /* **********  Creating and expanding blossoms:  ********** */
@@ -1168,6 +1215,9 @@ public:
         nontrivial_blossom.emplace_back(subblossoms, path.edges);
         NonTrivialBlossomT* blossom = &nontrivial_blossom.back();
 
+        // Initialize modified blossom dual to set true dual to 0.
+        blossom->dual_var = -2 * delta_sum;
+
         // Link the subblossoms to the their new parent.
         for (BlossomT* sub : subblossoms) {
             sub->parent = blossom;
@@ -1708,9 +1758,10 @@ public:
         // Compute delta4: half minimum dual of a top-level T-blossom.
         for (NonTrivialBlossomT& blossom : nontrivial_blossom) {
             if ((! blossom.parent) && (blossom.label == LABEL_T)) {
-                if (blossom.dual_var / 2 <= delta.value) {
+                WeightType true_dual_var = blossom.dual_var - 2 * delta_sum;
+                if (true_dual_var / 2 <= delta.value) {
                     delta.kind = 4;
-                    delta.value = blossom.dual_var / 2;
+                    delta.value = true_dual_var / 2;
                     delta.blossom = &blossom;
                 }
             }
@@ -1719,23 +1770,6 @@ public:
         return delta;
     }
 
-    /** Apply a delta step to the dual LPP variables. */
-    void substage_apply_delta_step(WeightType delta)
-    {
-        // Apply delta to dual variables of top-level non-trivial blossoms.
-        for (NonTrivialBlossomT& blossom : nontrivial_blossom) {
-            if (blossom.parent == nullptr) {
-                if (blossom.label == LABEL_S) {
-                    // S-blossom: add 2*delta to dual variable.
-                    blossom.dual_var += 2 * delta;
-                } else if (blossom.label == LABEL_T) {
-                    // T-blossom: subtract 2*delta from dual variable.
-                    blossom.dual_var -= 2 * delta;
-                }
-            }
-        }
-    }
-
     /* **********  Main algorithm:  ********** */
 
     /**
@@ -1812,10 +1846,8 @@ public:
             // Calculate delta step in the dual LPP problem.
             DeltaStep delta = calc_dual_delta_step();
 
-            // Apply the delta step to the dual variables.
+            // Update the running sum of delta steps. This implicitly updates
-            substage_apply_delta_step(delta.value);
+            // the dual variables of vertices and blossoms.
-
-            // Update the running sum of delta steps.
             delta_sum += delta.value;
 
             if (delta.kind == 2) {