diff --git a/cpp/mwmatching.h b/cpp/mwmatching.h index 0b1af14..6c32195 100644 --- a/cpp/mwmatching.h +++ b/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 - * that links to a different S-blossom, or "nullptr" if no such edge - * has been found. + * Accumulated pending lazy updates to the dual variables of the vertices + * inside the blossom. */ - const Edge* 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 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 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 + 2 * delta_sum; + return vertex_dual[x] + vertex_dual[y] - weight_factor * edge.weight; } // 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. - 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; - }); - } + // Expand the unlabeled blossom. + expand_unlabeled_blossom_core(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::max(); - 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]); - } - } + delta.value = init_vertex_dual - delta_sum; // 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(); } };