diff --git a/cpp/mwmatching.h b/cpp/mwmatching.h index 8bc5d45..270868f 100644 --- a/cpp/mwmatching.h +++ b/cpp/mwmatching.h @@ -148,29 +148,117 @@ void check_input_graph(const std::vector>& edges) } +/* ************************************************** + * ** struct Graph ** + * ************************************************** */ + /** - * Copy edges and remove edges with negative weight. - * - * This function takes time O(m). - * - * @param edges Vector of weighted edges defining the graph. - * @return Vector of edges with non-negative weight. + * Representation of the input graph. */ template -std::vector> remove_negative_weight_edges( - const std::vector>& edges) +struct Graph { - std::vector> real_edges; - real_edges.reserve(edges.size()); + using EdgeT = Edge; - for (const Edge& edge : edges) { - if (edge.weight >= 0) { - real_edges.push_back(edge); + /** Graph edges. */ + const std::vector edges; + + /** Number of vertices. */ + const VertexId num_vertex; + + /** For each vertex, a vector of pointers to its incident edges. */ + const std::vector> adjacent_edges; + + /** + * Initialize the graph representation and prepare adjacent edge lists. + * + * This function takes time O(n + m). + */ + explicit Graph(const std::vector& edges_in) + : edges(remove_negative_weight_edges(edges_in)), + num_vertex(count_num_vertex(edges)), + adjacent_edges(build_adjacent_edges(edges, num_vertex)) + { } + + // Prevent copying. + Graph(const Graph&) = delete; + Graph& operator=(const Graph&) = delete; + + /** + * Copy edges and remove edges with negative weight. + * + * This function takes time O(m). + * + * @param edges Vector of weighted edges defining the graph. + * @return Vector of edges with non-negative weight. + */ + static std::vector remove_negative_weight_edges( + const std::vector& edges_in) + { + std::vector real_edges; + real_edges.reserve(edges_in.size()); + + for (const Edge& edge : edges_in) { + if (edge.weight >= 0) { + real_edges.push_back(edge); + } } + + return real_edges; } - return real_edges; -} + /** Count the number of vertices in the graph. */ + static VertexId count_num_vertex(const std::vector& edges) + { + VertexId num_vertex = 0; + for (const Edge& edge : edges) { + VertexId m = std::max(edge.vt.first, edge.vt.second); + assert(m < std::numeric_limits::max()); + num_vertex = std::max(num_vertex, m + 1); + } + + return num_vertex; + } + + /** + * Build adjacency lists for a graph. + * + * Note that the adjacency lists contain pointers to edge instances, + * and therefore become invalid if the input edge vector is modified + * in any way that invalidates its iterators. + * + * This function takes time O(m). + * + * @param num_vertex Number of vertices in the graph. + * @param edges List of edges in the graph. + * @return Vector of incident edges for each vertex. + */ + static std::vector> build_adjacent_edges( + const std::vector& edges, + VertexId num_vertex) + { + // Count the number of incident edges per vertex. + std::vector vertex_degree(num_vertex); + for (const EdgeT& edge : edges) { + assert(edge.vt.first < num_vertex); + assert(edge.vt.second < num_vertex); + vertex_degree[edge.vt.first]++; + vertex_degree[edge.vt.second]++; + } + + // Build adjacency lists. + std::vector> adjacent_edges(num_vertex); + for (VertexId i = 0; i < num_vertex; ++i) { + adjacent_edges[i].reserve(vertex_degree[i]); + } + for (const EdgeT& edge : edges) { + adjacent_edges[edge.vt.first].push_back(&edge); + adjacent_edges[edge.vt.second].push_back(&edge); + } + + return adjacent_edges; + } +}; /* ************************************************** @@ -410,8 +498,10 @@ struct AlternatingPath * auxiliary data structures. */ template -struct MatchingContext +class MatchingContext { +public: + using EdgeT = Edge; using BlossomT = Blossom; using NonTrivialBlossomT = NonTrivialBlossom; @@ -436,14 +526,8 @@ struct MatchingContext static constexpr WeightType weight_factor = std::numeric_limits::is_integer ? 2 : 1; - /** Graph edges. */ - const std::vector edges; - - /** Number of vertices. */ - const VertexId num_vertex; - - /** For each vertex, a vector of pointers to its incident edges. */ - const std::vector> adjacent_edges; + /** Input graph. */ + const Graph graph; /** * Current state of the matching. @@ -506,99 +590,43 @@ struct MatchingContext * This function takes time O(n + m). */ explicit MatchingContext(const std::vector& edges_in) - : edges(remove_negative_weight_edges(edges_in)), - num_vertex(count_num_vertex(edges)), - adjacent_edges(build_adjacent_edges(edges, num_vertex)) + : graph(edges_in) { // Initially all vertices are unmatched. - vertex_mate.resize(num_vertex, NO_VERTEX); + vertex_mate.resize(graph.num_vertex, NO_VERTEX); // Create a trivial blossom for each vertex. - trivial_blossom.reserve(num_vertex); - for (VertexId x = 0; x < num_vertex; ++x) { + trivial_blossom.reserve(graph.num_vertex); + for (VertexId x = 0; x < graph.num_vertex; ++x) { trivial_blossom.emplace_back(x); } // Initially all vertices are trivial top-level blossoms. - vertex_top_blossom.reserve(num_vertex); - for (VertexId x = 0; x < num_vertex; ++x) { + vertex_top_blossom.reserve(graph.num_vertex); + for (VertexId x = 0; x < graph.num_vertex; ++x) { vertex_top_blossom.push_back(&trivial_blossom[x]); } // Vertex duals are initialized to half the maximum edge weight. WeightType max_weight = 0; - for (const EdgeT& edge : edges) { + for (const EdgeT& edge : graph.edges) { max_weight = std::max(max_weight, edge.weight); } WeightType init_vertex_dual = max_weight * (weight_factor / 2); - vertex_dual.resize(num_vertex, init_vertex_dual); + vertex_dual.resize(graph.num_vertex, init_vertex_dual); // Initialize "vertex_best_edge". - vertex_best_edge.resize(num_vertex, nullptr); + vertex_best_edge.resize(graph.num_vertex, nullptr); // Allocate temporary arrays for path tracing. - vertex_marker.resize(num_vertex); - marked_vertex.reserve(num_vertex); + vertex_marker.resize(graph.num_vertex); + marked_vertex.reserve(graph.num_vertex); } // Prevent copying. MatchingContext(const MatchingContext&) = delete; MatchingContext& operator=(const MatchingContext&) = delete; - /* ********** Static methods used by the constructor: ********** */ - - /** Count the number of vertices in the graph. */ - static VertexId count_num_vertex(const std::vector& edges) - { - VertexId num_vertex = 0; - for (const Edge& edge : edges) { - VertexId m = std::max(edge.vt.first, edge.vt.second); - assert(m < std::numeric_limits::max()); - num_vertex = std::max(num_vertex, m + 1); - } - - return num_vertex; - } - - /** - * Build adjacency lists for a graph. - * - * Note that the adjacency lists contain pointers to edge instances, - * and therefore become invalid if the input edge vector is modified - * in any way that invalidates its iterators. - * - * This function takes time O(m). - * - * @param num_vertex Number of vertices in the graph. - * @param edges List of edges in the graph. - * @return Vector of incident edges for each vertex. - */ - static std::vector> build_adjacent_edges( - const std::vector& edges, - VertexId num_vertex) - { - // Count the number of incident edges per vertex. - std::vector vertex_degree(num_vertex); - for (const EdgeT& edge : edges) { - assert(edge.vt.first < num_vertex); - assert(edge.vt.second < num_vertex); - vertex_degree[edge.vt.first]++; - vertex_degree[edge.vt.second]++; - } - - // Build adjacency lists. - std::vector> adjacent_edges(num_vertex); - for (VertexId i = 0; i < num_vertex; ++i) { - adjacent_edges[i].reserve(vertex_degree[i]); - } - for (const EdgeT& edge : edges) { - adjacent_edges[edge.vt.first].push_back(&edge); - adjacent_edges[edge.vt.second].push_back(&edge); - } - - return adjacent_edges; - } - /* ********** General support routines: ********** */ /** Calculate edge slack. */ @@ -648,7 +676,7 @@ struct MatchingContext */ void lset_reset() { - for (VertexId x = 0; x < num_vertex; ++x) { + for (VertexId x = 0; x < graph.num_vertex; ++x) { vertex_best_edge[x] = nullptr; } @@ -692,7 +720,7 @@ struct MatchingContext const EdgeT* best_edge = nullptr; WeightType best_slack = 0; - for (VertexId x = 0; x < num_vertex; ++x) { + for (VertexId x = 0; x < graph.num_vertex; ++x) { if (vertex_top_blossom[x]->label == LABEL_NONE) { const EdgeT* edge = vertex_best_edge[x]; if (edge != nullptr) { @@ -762,7 +790,7 @@ struct MatchingContext } else { // Trivial blossoms don't maintain a least-slack edge set. // Just consider all incident edges. - for (const EdgeT* edge : adjacent_edges[sub->base_vertex]) { + 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; @@ -782,7 +810,7 @@ struct MatchingContext // each top-level S-blossom. This array is indexed by the base vertex // of the blossoms. std::vector> - best_edge_to_blossom(num_vertex); + best_edge_to_blossom(graph.num_vertex); for (const EdgeT* edge : blossom->best_edge_set) { BlossomT* bx = vertex_top_blossom[edge->vt.first]; @@ -1467,7 +1495,7 @@ struct MatchingContext // Scan the edges that are incident on "x". // This loop runs through O(m) iterations per stage. - for (const EdgeT* edge : adjacent_edges[x]) { + for (const EdgeT* edge : graph.adjacent_edges[x]) { VertexId y = (edge->vt.first != x) ? edge->vt.first : edge->vt.second; @@ -1540,7 +1568,7 @@ struct MatchingContext // Compute delta1: minimum dual variable of any S-vertex. delta.kind = 1; delta.value = std::numeric_limits::max(); - for (VertexId x = 0; x < num_vertex; ++x) { + 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]); } @@ -1584,7 +1612,7 @@ struct MatchingContext void substage_apply_delta_step(WeightType delta) { // Apply delta to dual variables of all vertices. - for (VertexId x = 0; x < num_vertex; ++x) { + 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. @@ -1650,7 +1678,7 @@ struct MatchingContext bool run_stage() { // Assign label S to all unmatched vertices and put them in the queue. - for (VertexId x = 0; x < num_vertex; ++x) { + for (VertexId x = 0; x < graph.num_vertex; ++x) { if (vertex_mate[x] == NO_VERTEX) { assign_label_s(x); } @@ -1747,24 +1775,36 @@ struct MatchingContext /** Helper class to verify that an optimal solution has been found. */ template -struct MatchingVerifier +class MatchingVerifier { +public: + MatchingVerifier(const MatchingContext& ctx) + : ctx(ctx), + graph(ctx.graph), + edge_duals(ctx.graph.edges.size()) + { } + + /** + * Verify that the optimum solution has been found. + * + * This function takes time O(n**2). + * + * @return True if the solution is optimal; otherwise false. + */ + bool verify() + { + return (verify_vertex_mate() + && verify_vertex_duals() + && verify_blossom_duals() + && verify_blossoms_and_calc_edge_duals() + && verify_edge_slack()); + } + +private: using EdgeT = Edge; using BlossomT = Blossom; using NonTrivialBlossomT = NonTrivialBlossom; - // Reference to the MatchingContext instance. - const MatchingContext& ctx; - - // For each edge, the sum of duals of its incident vertices - // and duals of all blossoms that contain the edge. - std::vector edge_duals; - - MatchingVerifier(const MatchingContext& ctx) - : ctx(ctx), - edge_duals(ctx.edges.size()) - { } - static bool checked_add(WeightType& result, WeightType a, WeightType b) { if (a > std::numeric_limits::max() - b) { @@ -1778,7 +1818,7 @@ struct MatchingVerifier /** Convert edge pointer to its index in the vector "edges". */ std::size_t edge_index(const EdgeT* edge) { - return edge - ctx.edges.data(); + return edge - graph.edges.data(); } /** Check that the array "vertex_mate" is consistent. */ @@ -1786,7 +1826,7 @@ struct MatchingVerifier { // Count matched vertices and check symmetry of "vertex_mate". VertexId num_matched_vertex = 0; - for (VertexId x = 0; x < ctx.num_vertex; ++x) { + for (VertexId x = 0; x < graph.num_vertex; ++x) { VertexId y = ctx.vertex_mate[x]; if (y != NO_VERTEX) { ++num_matched_vertex; @@ -1798,7 +1838,7 @@ struct MatchingVerifier // Count matched edges. VertexId num_matched_edge = 0; - for (const EdgeT& edge : ctx.edges) { + for (const EdgeT& edge : graph.edges) { if (ctx.vertex_mate[edge.vt.first] == edge.vt.second) { ++num_matched_edge; } @@ -1814,7 +1854,7 @@ struct MatchingVerifier */ bool verify_vertex_duals() { - for (VertexId x = 0; x < ctx.num_vertex; ++x) { + for (VertexId x = 0; x < graph.num_vertex; ++x) { if (ctx.vertex_dual[x] < 0) { return false; } @@ -1859,7 +1899,7 @@ struct MatchingVerifier // For each vertex "x", // "vertex_depth[x]" is the depth of the smallest blossom on // the current descent path that contains "x". - std::vector vertex_depth(ctx.num_vertex); + std::vector vertex_depth(graph.num_vertex); // At each depth, keep track of the sum of blossom duals // along the current descent path. @@ -1922,7 +1962,7 @@ struct MatchingVerifier // For each incident edge, find the smallest blossom // that contains it. VertexId x = sub->base_vertex; - for (const EdgeT* edge : ctx.adjacent_edges[x]) { + for (const EdgeT* edge : graph.adjacent_edges[x]) { // Only consider edges pointing out from "x". if (edge->vt.first == x) { VertexId y = edge->vt.second; @@ -1991,7 +2031,7 @@ struct MatchingVerifier bool verify_blossoms_and_calc_edge_duals() { // For each edge, calculate the sum of its vertex duals. - for (const EdgeT& edge : ctx.edges) { + for (const EdgeT& edge : graph.edges) { if (checked_add(edge_duals[edge_index(&edge)], ctx.vertex_dual[edge.vt.first], ctx.vertex_dual[edge.vt.second])) { @@ -2022,7 +2062,7 @@ struct MatchingVerifier */ bool verify_edge_slack() { - for (const EdgeT& edge : ctx.edges) { + for (const EdgeT& edge : graph.edges) { WeightType duals = edge_duals[edge_index(&edge)]; WeightType weight = ctx.weight_factor * edge.weight; @@ -2040,21 +2080,18 @@ struct MatchingVerifier return true; } +private: + /** Reference to the MatchingContext instance. */ + const MatchingContext& ctx; + + /** Reference to the input graph. */ + const Graph& graph; + /** - * Verify that the optimum solution has been found. - * - * This function takes time O(n**2). - * - * @return True if the solution is optimal; otherwise false. + * For each edge, the sum of duals of its incident vertices + * and duals of all blossoms that contain the edge. */ - bool verify() - { - return (verify_vertex_mate() - && verify_vertex_duals() - && verify_blossom_duals() - && verify_blossoms_and_calc_edge_duals() - && verify_edge_slack()); - } + std::vector edge_duals; }; @@ -2114,7 +2151,7 @@ std::vector maximum_weight_matching( // Extract the matched edges. std::vector solution; - for (const Edge& edge : matching.edges) { + for (const Edge& edge : matching.graph.edges) { if (matching.vertex_mate[edge.vt.first] == edge.vt.second) { solution.push_back(edge.vt); }