From f98edae416607751b907f58ef47d6c2117c23b76 Mon Sep 17 00:00:00 2001 From: Joris van Rantwijk Date: Thu, 28 Nov 2024 22:38:15 +0100 Subject: [PATCH 01/18] weighted_matching_test check input file Without this check, the test program declares all tests passed if it fails to open the input file. --- test/weighted_matching_test.cpp | 1 + 1 file changed, 1 insertion(+) diff --git a/test/weighted_matching_test.cpp b/test/weighted_matching_test.cpp index 2baa0adcf..30042c553 100644 --- a/test/weighted_matching_test.cpp +++ b/test/weighted_matching_test.cpp @@ -139,6 +139,7 @@ Graph make_graph(typename graph_traits< Graph >::vertices_size_type num_v, int main(int, char*[]) { std::ifstream in_file("weighted_matching.dat"); + BOOST_TEST(in_file.good()); std::string line; while (std::getline(in_file, line)) { From f30ba13a24360ab040afae09eb5df640e35289ea Mon Sep 17 00:00:00 2001 From: Joris van Rantwijk Date: Sun, 1 Dec 2024 12:18:18 +0100 Subject: [PATCH 02/18] Add more tests for maximum_weighted_matching - Hand-picked graphs to explore basic functionality. - Graphs that are known to trigger bugs in the old implementation of maximum_weighted_matching(). - Random small graphs. --- test/Jamfile.v2 | 1 + test/weighted_matching_test2.cpp | 522 +++++++++++++++++++++++++++++++ 2 files changed, 523 insertions(+) create mode 100644 test/weighted_matching_test2.cpp diff --git a/test/Jamfile.v2 b/test/Jamfile.v2 index b2656eb07..b2e29dd02 100644 --- a/test/Jamfile.v2 +++ b/test/Jamfile.v2 @@ -114,6 +114,7 @@ alias graph_test_regular : [ run king_ordering.cpp ] [ run matching_test.cpp ] [ run weighted_matching_test.cpp ] + [ run weighted_matching_test2.cpp ] [ run max_flow_test.cpp ] [ run boykov_kolmogorov_max_flow_test.cpp ] [ run cycle_ratio_tests.cpp ../build//boost_graph : $(CYCLE_RATIO_INPUT_FILE) ] diff --git a/test/weighted_matching_test2.cpp b/test/weighted_matching_test2.cpp new file mode 100644 index 000000000..3db775f0b --- /dev/null +++ b/test/weighted_matching_test2.cpp @@ -0,0 +1,522 @@ +//======================================================================= +// Copyright (c) 2024 Joris van Rantwijk +// +// Distributed under the Boost Software License, Version 1.0. +// (See accompanying file LICENSE_1_0.txt or copy at +// http://www.boost.org/LICENSE_1_0.txt) +// +//======================================================================= + +#include + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +using namespace boost; + + +template < typename WeightType > +using adj_vec_test_graph = adjacency_list< + vecS, // OutEdgeList + vecS, // VertexList + undirectedS, // DirectedS + no_property, // VertexProperty (vertex_index is implicit) + property< edge_weight_t, WeightType > >; // EdgeProperty + +template < typename WeightType > +using adj_list_test_graph = adjacency_list< + listS, // OutEdgeList + listS, // VertexList + undirectedS, // DirectedS + property< vertex_index_t, unsigned int >, // VertexProperty + property< edge_weight_t, WeightType > >; // EdgeProperty + +template < typename WeightType > +using adj_matrix_test_graph = adjacency_matrix< + undirectedS, // Directed + no_property, // VertexProperty (vertex_index is implicit) + property< edge_weight_t, WeightType > >; // EdgeProperty + + +template < typename WeightType > +struct edge_info +{ + unsigned int x, y; + WeightType w; +}; + + +// Initialize vertex_index if necessary. +template < typename Graph > +struct vertex_index_installer +{ + static void install(Graph&) {} +}; + +template < typename WeightType > +struct vertex_index_installer< adj_list_test_graph< WeightType > > +{ + static void install(adj_list_test_graph< WeightType >& g) + { + auto vrange = vertices(g); + unsigned int i = 0; + for (auto it = vrange.first; it != vrange.second; ++it, ++i) + put(vertex_index, g, *it, i); + } +}; + + +template < typename Graph, typename WeightType > +Graph make_graph( + unsigned int nv, + const std::vector< edge_info< WeightType > >& edges_info) +{ + typedef property< edge_weight_t, WeightType > EdgeProperty; + Graph g(nv); + vertex_index_installer< Graph >::install(g); + for (const auto& e : edges_info) + add_edge(vertex(e.x, g), vertex(e.y, g), EdgeProperty(e.w), g); + return g; +} + + +template < typename Graph > +using graph_edge_weight_type = typename property_traits< + typename property_map< Graph, edge_weight_t >::type >::value_type; + + +template < typename Graph, typename MateMap > +std::pair< graph_edge_weight_type< Graph >, bool > +check_matching(const Graph& g, const MateMap& mate) +{ + graph_edge_weight_type< Graph > matching_weight = 0; + unsigned int n_matched_edges = 0; + + auto edge_range = edges(g); + for (auto it = edge_range.first; it != edge_range.second; ++it) + { + if (mate[source(*it, g)] == target(*it, g)) + { + ++n_matched_edges; + matching_weight += get(edge_weight, g, *it); + } + } + + unsigned int n_matched_vertices = 0; + auto vertex_range = vertices(g); + for (auto it = vertex_range.first; it != vertex_range.second; ++it) + { + if (mate[*it] != graph_traits< Graph >::null_vertex()) + { + if (mate[mate[*it]] != *it) + return std::make_pair(0, false); + ++n_matched_vertices; + } + } + + if (n_matched_vertices != 2 * n_matched_edges) + return std::make_pair(0, false); + + return std::make_pair(matching_weight, true); +} + + +template < typename Graph > +std::string show_graph(const Graph& g) +{ + std::ostringstream ss; + auto edge_range = edges(g); + for (auto it = edge_range.first; it != edge_range.second; ++it) + { + unsigned int x = get(vertex_index, g, source(*it, g)); + unsigned int y = get(vertex_index, g, target(*it, g)); + auto w = get(edge_weight, g, *it); + if (ss.tellp() > 0) + ss << " "; + ss << "{" << x << "," << y << "," << w << "}"; + } + return ss.str(); +} + + +template < typename Graph, typename WeightType > +void test_matching(const Graph& g, WeightType answer) +{ + typedef typename property_map< Graph, vertex_index_t >::type + vertex_index_map_t; + typedef vector_property_map< + typename graph_traits< Graph>::vertex_descriptor, vertex_index_map_t > + mate_t; + + mate_t mate(num_vertices(g)); + maximum_weighted_matching(g, mate); + + WeightType weight; + bool consistent_matching; + std::tie(weight, consistent_matching) = check_matching(g, mate); + BOOST_TEST(consistent_matching); + if (! consistent_matching) + { + std::cout << std::endl + << "Inconsistent answer for graph" << std::endl + << " " << show_graph(g) << std::endl; + } + else + { + bool same_answer; + if (std::numeric_limits< WeightType >::is_integer) + { + same_answer = (weight == answer); + } + else + { + WeightType max_error = std::max(1e-14, answer * 1e-14); + same_answer = (std::abs(weight - answer) <= max_error); + } + + BOOST_TEST(same_answer); + if (! same_answer) + { + std::cout << std::endl + << "Wrong answer for graph" << std::endl + << " " << show_graph(g) << std::endl + << " found weight " << weight + << " while expecting " << answer << std::endl; + } + } +} + + +template < typename WeightType = long > +void run_test_graph( + unsigned int nv, + const std::vector< edge_info< WeightType > >& edges_info, + WeightType answer) +{ + typedef adj_vec_test_graph< WeightType > vec_graph_t; + test_matching( + make_graph< vec_graph_t, WeightType >(nv, edges_info), + answer); + + typedef adj_list_test_graph< WeightType > list_graph_t; + test_matching( + make_graph< list_graph_t, WeightType >(nv, edges_info), + answer); + + typedef adj_matrix_test_graph< WeightType > matrix_graph_t; + test_matching( + make_graph< matrix_graph_t, WeightType >(nv, edges_info), + answer); +} + + +template < typename WeightType > +WeightType find_brute_force_answer( + unsigned int nv, + const std::vector< edge_info< WeightType > >& edges_info) +{ + typedef adj_vec_test_graph< WeightType > vec_graph_t; + typedef typename graph_traits< vec_graph_t >::vertex_descriptor vertex_t; + + vec_graph_t g = make_graph< vec_graph_t, WeightType >(nv, edges_info); + vector_property_map< vertex_t > mate(nv); + brute_force_maximum_weighted_matching(g, mate); + + WeightType weight; + bool ok; + std::tie(weight, ok) = check_matching(g, mate); + BOOST_TEST(ok); + return weight; +} + + +template < typename WeightType > +void run_random_graph(unsigned int n, unsigned int m, std::mt19937& rng) +{ + using edge_weight_dist_type = typename std::conditional< + std::is_integral< WeightType >::value, + std::uniform_int_distribution< WeightType >, + std::uniform_real_distribution< WeightType > >::type; + + std::vector< edge_info< WeightType > > edges_info; + for (unsigned int i = 0; i < m; ++i) + { + while (true) + { + std::uniform_int_distribution< unsigned int > xdist{0, n - 1}; + unsigned int x = xdist(rng); + std::uniform_int_distribution< unsigned int > ydist{0, n - 2}; + unsigned int y = ydist(rng); + if (y >= x) + y += 1; + bool ok = true; + for (const auto& e : edges_info) + { + if (((e.x == x) && (e.y == y)) || ((e.x == y) && (e.y == x))) + { + ok = false; + break; + } + } + if (ok) + { + edge_weight_dist_type edge_weight_dist(0, 1000); + WeightType w = edge_weight_dist(rng); + edges_info.push_back(edge_info< WeightType >{x, y, w}); + break; + } + } + } + WeightType answer = find_brute_force_answer(n, edges_info); + run_test_graph< WeightType >(n, edges_info, answer); +} + + +int main(int, char*[]) +{ + // Simple test cases: + run_test_graph(0, {}, 0); + run_test_graph(1, {}, 0); + run_test_graph(2, {}, 0); + run_test_graph(2, {{0, 1, 1}}, 1); + run_test_graph(3, {{0, 1, 10}, {1, 2, 11}}, 11); + run_test_graph(4, {{0, 1, 5}, {1, 2, 11}, {2, 3, 5}}, 11); + run_test_graph(4, {{0, 1, 5}, {1, 2, 11}, {2, 3, 7}}, 12); + + // Floating point edge weights: + run_test_graph(4, { + {0, 1, 3.1415}, {1, 2, 2.7183}, {0, 2, 3.0}, {0, 3, 1.4142}}, + 4.1325); + + // Negative edge weights: + run_test_graph(4, { + {0, 1, 2}, {0, 2, -2}, {1, 2, 1}, {1, 3, -1}, {2, 3, -6}}, + 2); + + // Blossoms: + run_test_graph(4, {{0, 1, 8}, {0, 2, 9}, {1, 2, 10}, {2, 3, 7}}, 15); + + run_test_graph(6, { + {0, 1, 8}, {0, 2, 9}, {1, 2, 10}, {2, 3, 7}, {0, 5, 5}, {3, 4, 6}}, + 21); + + run_test_graph(6, { + {0, 1, 9}, {0, 2, 8}, {1, 2, 10}, {0, 3, 5}, {3, 4, 4}, {0, 5, 3}}, + 17); + + run_test_graph(6, { + {0, 1, 9}, {0, 2, 8}, {1, 2, 10}, {0, 3, 5}, {3, 4, 3}, {0, 5, 4}}, + 17); + + run_test_graph(6, { + {0, 1, 9}, {0, 2, 8}, {1, 2, 10}, {0, 3, 5}, {3, 4, 3}, {2, 5, 4}}, + 16); + + run_test_graph(6, { + {0, 1, 9}, {0, 2, 9}, {1, 2, 10}, {1, 3, 8}, {2, 4, 8}, + {3, 4, 10}, {4, 5, 6}}, + 23); + + run_test_graph(8, { + {0, 1, 10}, {0, 6, 10}, {1, 2, 12}, {2, 3, 20}, {2, 4, 20}, + {3, 4, 25}, {4, 5, 10}, {5, 6, 10}, {6, 7, 8}}, + 48); + + run_test_graph(8, { + {0, 1, 8}, {0, 2, 8}, {1, 2, 10}, {1, 3, 12}, {2, 4, 12}, + {3, 4, 14}, {3, 5, 12}, {4, 6, 12}, {5, 6, 14}, {6, 7, 12}}, + 44); + + run_test_graph(8, { + {0, 1, 23}, {0, 4, 22}, {0, 5, 15}, {1, 2, 25}, {2, 3, 22}, + {3, 4, 25}, {3, 7, 14}, {4, 6, 13}}, + 67); + + run_test_graph(8, { + {0, 1, 19}, {0, 2, 20}, {0, 7, 8}, {1, 2, 25}, {1, 3, 18}, + {2, 4, 18}, {3, 4, 13}, {3, 6, 7}, {4, 5, 7}}, + 47); + + // Somewhat tricky test cases: + run_test_graph(10, { + {0, 1, 45}, {0, 4, 45}, {1, 2, 50}, {2, 3, 45}, {3, 4, 50}, + {0, 5, 30}, {2, 8, 35}, {3, 7, 35}, {4, 6, 26}, {8, 9, 5}}, + 146); + + run_test_graph(10, { + {0, 1, 45}, {0, 4, 45}, {1, 2, 50}, {2, 3, 45}, {3, 4, 50}, + {0, 5, 30}, {2, 8, 35}, {3, 7, 26}, {4, 6, 40}, {8, 9, 5}}, + 151); + + run_test_graph(10, { + {0, 1, 45}, {0, 4, 45}, {1, 2, 50}, {2, 3, 45}, {3, 4, 50}, + {0, 5, 30}, {2, 8, 35}, {3, 7, 28}, {4, 6, 26}, {8, 9, 5}}, + 139); + + run_test_graph(12, { + {0, 1, 45}, {0, 6, 45}, {1, 2, 50}, {2, 3, 45}, {3, 4, 95}, + {3, 5, 94}, {4, 5, 94}, {5, 6, 50}, {0, 7, 30}, {2, 10, 35}, + {4, 8, 36}, {6, 9, 26}, {10, 11, 5}}, + 241); + + run_test_graph(10, { + {0, 1, 40}, {0, 2, 40}, {1, 2, 60}, {1, 3, 55}, {2, 4, 55}, + {3, 4, 50}, {0, 7, 15}, {4, 6, 30}, {6, 5, 10}, {7, 9, 10}, + {3, 8, 30}}, + 145); + + run_test_graph(6, { + {0, 1, 2}, {0, 4, 3}, {1, 2, 7}, {1, 5, 2}, {2, 3, 9}, + {2, 5, 4}, {3, 4, 8}, {3, 5, 4}}, + 15); + + run_test_graph(6, { + {0, 1, 8}, {0, 2, 8}, {1, 2, 9}, {1, 3, 6}, {2, 4, 7}, + {3, 4, 8}, {3, 5, 5}}, + 20); + + run_test_graph(7, { + {0, 1, 7}, {0, 2, 7}, {1, 2, 9}, {0, 3, 7}, {0, 4, 7}, + {3, 4, 9}, {5, 6, 2}}, + 20); + + run_test_graph(7, { + {0, 1, 7}, {0, 4, 6}, {1, 2, 9}, {2, 3, 8}, {3, 4, 9}, + {3, 5, 1}, {4, 5, 1}}, + 18); + + run_test_graph(5, { + {0, 2, 7}, {0, 3, 3}, {1, 3, 1}, {2, 3, 5}, {2, 4, 5}}, + 8); + + run_test_graph(7, { + {0, 1, 15}, {0, 2, 10}, {0, 3, 11}, {0, 5, 17}, {1, 4, 12}, + {1, 5, 8}, {4, 6, 15}, {5, 6, 7}}, + 34); + + run_test_graph(7, { + {0, 1, 12}, {0, 3, 11}, {0, 4, 11}, {1, 3, 12}, {1, 5, 11}, + {2, 3, 14}, {2, 6, 14}, {3, 6, 11}}, + 37); + + run_test_graph(7, { + {0, 1, 19}, {0, 3, 17}, {0, 4, 19}, {1, 2, 15}, {1, 4, 21}, + {3, 5, 18}, {3, 6, 11}, {4, 5, 19}}, + 52); + + run_test_graph(6, { + {0, 1, 19}, {0, 2, 19}, {0, 3, 15}, {0, 4, 17}, {1, 2, 21}, + {3, 4, 16}, {4, 5, 10}}, + 46); + + run_test_graph(6, { + {0, 1, 48}, {0, 3, 42}, {0, 4, 57}, {1, 3, 51}, {1, 5, 36}, + {2, 3, 23}, {4, 5, 46}}, + 117); + + run_test_graph(6, { + {0, 2, 7}, {0, 5, 20}, {1, 2, 50}, {1, 4, 46}, {2, 3, 35}, + {2, 4, 8}, {2, 5, 25}, {3, 5, 47}}, + 101); + + run_test_graph(5, { + {0, 1, 48}, {0, 2, 44}, {0, 4, 48}, {1, 4, 36}, {3, 4, 31}}, + 80); + + run_test_graph(6, { + {3, 5, 71}, {2, 4, 36}, {3, 4, 27}, {5, 1, 29}, {2, 0, 48}, + {1, 0, 60}, {1, 3, 11}, {4, 5, 54}}, + 167); + + run_test_graph(6, { + {0, 2, 92}, {5, 4, 50}, {5, 0, 86}, {2, 1, 87}, {1, 3, 39}, + {4, 0, 2}, {4, 1, 83}, {0, 1, 56}}, + 181); + + // Triangles, unit weight: + run_test_graph(9, { + {0, 1, 1}, {0, 2, 1}, {1, 2, 1}, {0, 3, 1}, + {3, 4, 1}, {3, 5, 1}, {4, 5, 1}, {4, 7, 1}, + {6, 7, 1}, {6, 8, 1}, {7, 8, 1}}, + 4); + + // Trigger known bugs in a previous version of maximum_weighted_matching: + + // wrong answer + run_test_graph(8, { + {4, 7, 453}, {0, 4, 627}, {4, 6, 853}, {6, 7, 344}, {5, 6, 906}, + {4, 5, 689}, {2, 3, 741}, {2, 7, 746}, {3, 6, 647}, {1, 5, 385}, + {5, 7, 215}, {3, 7, 640}}, + 2405); + + // wrong answer + run_test_graph(12, { + {0, 3, 448}, {0, 4, 919}, {0, 5, 918}, {0, 9, 72}, {1, 3, 830}, + {1, 4, 687}, {1, 5, 559}, {1, 6, 580}, {1, 7, 679}, {1, 8, 627}, + {2, 4, 585}, {2, 5, 835}, {2, 6, 822}, {2, 8, 462}, {2, 11, 380}, + {3, 4, 328}, {3, 5, 860}, {3, 7, 297}, {3, 8, 590}, {3, 11, 235}, + {4, 6, 692}, {4, 10, 227}, {4, 11, 354}, {5, 6, 160}, {5, 9, 400}, + {6, 8, 410}, {6, 9, 420}, {7, 11, 924}, {9, 11, 825}, {10, 11, 790}}, + 4233); + + // assertion + run_test_graph(9, { + {4, 6, 796}, {3, 4, 553}, {0, 1, 792}, {1, 7, 1000}, {4, 5, 360}, + {5, 7, 183}, {4, 8, 694}, {0, 4, 741}, {0, 2, 483}, {5, 8, 228}, + {7, 8, 644}, {1, 3, 236}, {6, 7, 895}, {1, 6, 913}, {1, 8, 617}}, + 2593); + + // assertion + run_test_graph(10, { + {0, 6, 892}, {2, 4, 517}, {4, 7, 560}, {1, 5, 828}, {1, 7, 831}, + {6, 7, 397}, {4, 5, 43}, {5, 6, 944}, {6, 9, 215}, {5, 7, 753}, + {4, 6, 901}, {1, 2, 530}, {2, 8, 384}, {3, 4, 499}, {5, 8, 190}}, + 2674); + + // segmentation fault + run_test_graph(7, { + {5, 6, 799}, {0, 4, 601}, {0, 1, 578}, {0, 3, 373}, {4, 6, 675}, + {4, 5, 925}, {0, 5, 697}, {3, 4, 260}, {1, 3, 464}, {1, 5, 845}, + {0, 2, 176}, {2, 5, 685}}, + 1938); + + // segmentation fault + run_test_graph(8, { + {2, 7, 420}, {6, 7, 414}, {2, 3, 421}, {2, 6, 854}, {1, 6, 997}, + {4, 5, 46}, {1, 2, 467}, {2, 4, 230}, {3, 7, 555}, {0, 3, 334}, + {0, 7, 341}, {0, 6, 634}}, + 1805); + + // hang + run_test_graph(7, { + {2, 3, 837}, {1, 4, 458}, {3, 4, 291}, {5, 6, 601}, {1, 2, 202}, + {0, 4, 491}, {4, 5, 910}, {0, 1, 159}, {3, 5, 684}, {4, 6, 139}, + {0, 2, 792}, {1, 3, 232}}, + 1934); + + // hang + run_test_graph(12, { + {0, 1, 856}, {0, 4, 462}, {0, 6, 874}, {0, 7, 406}, {0, 10, 294}, + {1, 3, 936}, {1, 5, 852}, {1, 7, 501}, {1, 8, 555}, {1, 10, 41}, + {2, 3, 325}, {2, 6, 748}, {2, 9, 808}, {3, 5, 870}, {3, 7, 25}, + {3, 8, 663}, {3, 9, 897}, {4, 7, 617}, {4, 9, 435}, {4, 10, 818}, + {4, 11, 933}, {5, 6, 608}, {5, 8, 636}, {6, 9, 274}, {6, 10, 279}, + {7, 9, 705}, {7, 10, 114}, {8, 11, 602}, {9, 10, 764}, {10, 11, 413}}, + 4599); + + // Random graphs: + std::mt19937 rng(123456); + for (int k = 0; k < 20000; ++k) + run_random_graph< long >(10, 15, rng); + for (int k = 0; k < 10000; ++k) + run_random_graph< double >(10, 15, rng); + + return boost::report_errors(); +} From fb66ac492fd6cd725a38969af55394d5b71598bb Mon Sep 17 00:00:00 2001 From: Joris van Rantwijk Date: Sat, 11 Jan 2025 18:28:40 +0100 Subject: [PATCH 03/18] Replace maximum_weighted_matching implementation The new code runs in O(V^3). It also solves a number of known issues. --- .../boost/graph/maximum_weighted_matching.hpp | 2434 ++++++++++------- 1 file changed, 1428 insertions(+), 1006 deletions(-) diff --git a/include/boost/graph/maximum_weighted_matching.hpp b/include/boost/graph/maximum_weighted_matching.hpp index 098a9c2aa..cd3ded6f1 100644 --- a/include/boost/graph/maximum_weighted_matching.hpp +++ b/include/boost/graph/maximum_weighted_matching.hpp @@ -1,5 +1,6 @@ //======================================================================= // Copyright (c) 2018 Yi Ji +// Copyright (c) 2025 Joris van Rantwijk // // Distributed under the Boost Software License, Version 1.0. // (See accompanying file LICENSE_1_0.txt or copy at @@ -10,13 +11,27 @@ #ifndef BOOST_GRAPH_MAXIMUM_WEIGHTED_MATCHING_HPP #define BOOST_GRAPH_MAXIMUM_WEIGHTED_MATCHING_HPP -#include // for std::iter_swap -#include -#include -#include +#include +#include +#include +#include +#include +#include +#include // for std::tie +#include // for std::pair, std::swap +#include + +#include +#include +#include +#include +#include +#include +#include // for empty_matching namespace boost { + template < typename Graph, typename MateMap, typename VertexIndexMap > typename property_traits< typename property_map< Graph, edge_weight_t >::type >::value_type @@ -49,1265 +64,1672 @@ matching_weight_sum(const Graph& g, MateMap mate) return matching_weight_sum(g, mate, get(vertex_index, g)); } +// brute-force matcher searches all possible combinations of matched edges to +// get the maximum weighted matching which can be used for testing on small +// graphs (within dozens vertices) template < typename Graph, typename MateMap, typename VertexIndexMap > -class weighted_augmenting_path_finder +class brute_force_matching { public: - template < typename T > struct map_vertex_to_ - { - typedef boost::iterator_property_map< - typename std::vector< T >::iterator, VertexIndexMap > - type; - }; - typedef typename graph::detail::VERTEX_STATE vertex_state_t; - typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t; typedef typename graph_traits< Graph >::vertex_descriptor vertex_descriptor_t; - typedef typename std::vector< vertex_descriptor_t >::const_iterator - vertex_vec_iter_t; + typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t; typedef - typename graph_traits< Graph >::out_edge_iterator out_edge_iterator_t; - typedef typename graph_traits< Graph >::edge_descriptor edge_descriptor_t; + typename std::vector< vertex_descriptor_t >::iterator vertex_vec_iter_t; typedef typename graph_traits< Graph >::edge_iterator edge_iterator_t; - typedef typename property_traits< typename property_map< Graph, - edge_weight_t >::type >::value_type edge_property_t; - typedef std::deque< vertex_descriptor_t > vertex_list_t; - typedef std::vector< edge_descriptor_t > edge_list_t; - typedef typename map_vertex_to_< vertex_descriptor_t >::type + typedef boost::iterator_property_map< vertex_vec_iter_t, VertexIndexMap > vertex_to_vertex_map_t; - typedef - typename map_vertex_to_< edge_property_t >::type vertex_to_weight_map_t; - typedef typename map_vertex_to_< bool >::type vertex_to_bool_map_t; - typedef typename map_vertex_to_< std::pair< vertex_descriptor_t, - vertex_descriptor_t > >::type vertex_to_pair_map_t; - typedef - typename map_vertex_to_< std::pair< edge_descriptor_t, bool > >::type - vertex_to_edge_map_t; - typedef typename map_vertex_to_< vertex_to_edge_map_t >::type - vertex_pair_to_edge_map_t; - class blossom + brute_force_matching( + const Graph& arg_g, MateMap arg_mate, VertexIndexMap arg_vm) + : g(arg_g) + , vm(arg_vm) + , mate_vector(num_vertices(g)) + , best_mate_vector(num_vertices(g)) + , mate(mate_vector.begin(), vm) + , best_mate(best_mate_vector.begin(), vm) { - public: - typedef boost::shared_ptr< blossom > blossom_ptr_t; - std::vector< blossom_ptr_t > sub_blossoms; - edge_property_t dual_var; - blossom_ptr_t father; - - blossom() : dual_var(0), father(blossom_ptr_t()) {} - - // get the base vertex of a blossom by recursively getting - // its base sub-blossom, which is always the first one in - // sub_blossoms because of how we create and maintain blossoms - virtual vertex_descriptor_t get_base() const - { - const blossom* b = this; - while (!b->sub_blossoms.empty()) - b = b->sub_blossoms[0].get(); - return b->get_base(); - } + vertex_iterator_t vi, vi_end; + for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) + best_mate[*vi] = mate[*vi] = get(arg_mate, *vi); + } - // set a sub-blossom as a blossom's base by exchanging it - // with its first sub-blossom - void set_base(const blossom_ptr_t& sub) - { - for (blossom_iterator_t bi = sub_blossoms.begin(); - bi != sub_blossoms.end(); ++bi) - { - if (sub.get() == bi->get()) - { - std::iter_swap(sub_blossoms.begin(), bi); - break; - } - } - } + template < typename PropertyMap > void find_matching(PropertyMap pm) + { + edge_iterator_t ei; + boost::tie(ei, ei_end) = edges(g); + select_edge(ei); + + vertex_iterator_t vi, vi_end; + for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) + put(pm, *vi, best_mate[*vi]); + } + +private: + const Graph& g; + VertexIndexMap vm; + std::vector< vertex_descriptor_t > mate_vector, best_mate_vector; + vertex_to_vertex_map_t mate, best_mate; + edge_iterator_t ei_end; - // get all vertices inside recursively - virtual std::vector< vertex_descriptor_t > vertices() const + void select_edge(edge_iterator_t ei) + { + if (ei == ei_end) { - std::vector< vertex_descriptor_t > all_vertices; - for (typename std::vector< blossom_ptr_t >::const_iterator bi - = sub_blossoms.begin(); - bi != sub_blossoms.end(); ++bi) + if (matching_weight_sum(g, mate) + > matching_weight_sum(g, best_mate)) { - std::vector< vertex_descriptor_t > some_vertices - = (*bi)->vertices(); - all_vertices.insert(all_vertices.end(), some_vertices.begin(), - some_vertices.end()); + vertex_iterator_t vi, vi_end; + for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) + best_mate[*vi] = mate[*vi]; } - return all_vertices; + return; } - }; - // a trivial_blossom only has one vertex and no sub-blossom; - // for each vertex v, in_blossom[v] is the trivial_blossom that contains it - // directly - class trivial_blossom : public blossom - { - public: - trivial_blossom(vertex_descriptor_t v) : trivial_vertex(v) {} - virtual vertex_descriptor_t get_base() const { return trivial_vertex; } + vertex_descriptor_t v, w; + v = source(*ei, g); + w = target(*ei, g); + + select_edge(++ei); - virtual std::vector< vertex_descriptor_t > vertices() const + if (mate[v] == graph_traits< Graph >::null_vertex() + && mate[w] == graph_traits< Graph >::null_vertex()) { - std::vector< vertex_descriptor_t > all_vertices; - all_vertices.push_back(trivial_vertex); - return all_vertices; + mate[v] = w; + mate[w] = v; + select_edge(ei); + mate[v] = mate[w] = graph_traits< Graph >::null_vertex(); } + } +}; - private: - vertex_descriptor_t trivial_vertex; - }; - - typedef boost::shared_ptr< blossom > blossom_ptr_t; - typedef typename std::vector< blossom_ptr_t >::iterator blossom_iterator_t; - typedef - typename map_vertex_to_< blossom_ptr_t >::type vertex_to_blossom_map_t; - - weighted_augmenting_path_finder( - const Graph& arg_g, MateMap arg_mate, VertexIndexMap arg_vm) - : g(arg_g) - , vm(arg_vm) - , null_edge(std::pair< edge_descriptor_t, bool >( - num_edges(g) == 0 ? edge_descriptor_t() : *edges(g).first, false)) - , mate_vector(num_vertices(g)) - , label_S_vector(num_vertices(g), graph_traits< Graph >::null_vertex()) - , label_T_vector(num_vertices(g), graph_traits< Graph >::null_vertex()) - , outlet_vector(num_vertices(g), graph_traits< Graph >::null_vertex()) - , tau_idx_vector(num_vertices(g), graph_traits< Graph >::null_vertex()) - , dual_var_vector(std::vector< edge_property_t >( - num_vertices(g), std::numeric_limits< edge_property_t >::min())) - , pi_vector(std::vector< edge_property_t >( - num_vertices(g), std::numeric_limits< edge_property_t >::max())) - , gamma_vector(std::vector< edge_property_t >( - num_vertices(g), std::numeric_limits< edge_property_t >::max())) - , tau_vector(std::vector< edge_property_t >( - num_vertices(g), std::numeric_limits< edge_property_t >::max())) - , in_blossom_vector(num_vertices(g)) - , old_label_vector(num_vertices(g)) - , critical_edge_vectors(num_vertices(g), - std::vector< std::pair< edge_descriptor_t, bool > >( - num_vertices(g), null_edge)) - , - - mate(mate_vector.begin(), vm) - , label_S(label_S_vector.begin(), vm) - , label_T(label_T_vector.begin(), vm) - , outlet(outlet_vector.begin(), vm) - , tau_idx(tau_idx_vector.begin(), vm) - , dual_var(dual_var_vector.begin(), vm) - , pi(pi_vector.begin(), vm) - , gamma(gamma_vector.begin(), vm) - , tau(tau_vector.begin(), vm) - , in_blossom(in_blossom_vector.begin(), vm) - , old_label(old_label_vector.begin(), vm) - { - vertex_iterator_t vi, vi_end; - edge_iterator_t ei, ei_end; - - edge_property_t max_weight - = std::numeric_limits< edge_property_t >::min(); - for (boost::tie(ei, ei_end) = edges(g); ei != ei_end; ++ei) - max_weight = std::max(max_weight, get(edge_weight, g, *ei)); - - typename std::vector< - std::vector< std::pair< edge_descriptor_t, bool > > >::iterator vei; +template < typename Graph, typename MateMap, typename VertexIndexMap > +void brute_force_maximum_weighted_matching( + const Graph& g, MateMap mate, VertexIndexMap vm) +{ + empty_matching< Graph, MateMap >::find_matching(g, mate); + brute_force_matching< Graph, MateMap, VertexIndexMap > brute_force_matcher( + g, mate, vm); + brute_force_matcher.find_matching(mate); +} - for (boost::tie(vi, vi_end) = vertices(g), - vei = critical_edge_vectors.begin(); - vi != vi_end; ++vi, ++vei) - { - vertex_descriptor_t u = *vi; - mate[u] = get(arg_mate, u); - dual_var[u] = 2 * max_weight; - in_blossom[u] = boost::make_shared< trivial_blossom >(u); - outlet[u] = u; - critical_edge_vector.push_back( - vertex_to_edge_map_t(vei->begin(), vm)); - } +template < typename Graph, typename MateMap > +inline void brute_force_maximum_weighted_matching(const Graph& g, MateMap mate) +{ + brute_force_maximum_weighted_matching(g, mate, get(vertex_index, g)); +} - critical_edge - = vertex_pair_to_edge_map_t(critical_edge_vector.begin(), vm); +namespace graph +{ +namespace detail +{ - init(); +/** Check that vertex indices are unique and in range [0, V). */ +template < typename Graph, typename VertexIndexMap > +void check_vertex_index_range(const Graph& g, VertexIndexMap vm) +{ + typedef typename property_traits< VertexIndexMap >::value_type index_t; + typedef typename std::make_unsigned::type unsigned_index_t; + auto nv = num_vertices(g); + std::vector got_vertex(nv); + for (const auto& x : make_iterator_range(vertices(g))) + { + index_t i = get(vm, x); + if ((i < 0) || (static_cast(i) >= nv)) + throw bad_graph("Invalid vertex index."); + if (got_vertex[i]) + throw bad_graph("Duplicate vertex index."); + got_vertex[i] = true; } +} - // return the top blossom where v is contained inside - blossom_ptr_t in_top_blossom(vertex_descriptor_t v) const +/** Check that edge weights are valid. */ +template < typename Graph, typename EdgeWeightMap > +void check_maximum_weighted_matching_edge_weights( + const Graph& g, EdgeWeightMap edge_weights) +{ + for (const auto& e : make_iterator_range(edges(g))) { - blossom_ptr_t b = in_blossom[v]; - while (b->father != blossom_ptr_t()) - b = b->father; - return b; + auto w = get(edge_weights, e); + auto max_weight = (std::numeric_limits< decltype(w) >::max)() / 4; + if (! (w <= max_weight)) // inverted logic to catch NaN + throw bad_graph("Edge weight exceeds maximum supported value."); } +} - // check if vertex v is in blossom b - bool is_in_blossom(blossom_ptr_t b, vertex_descriptor_t v) const +/** Implementation of the matching algorithm. */ +template < typename Graph, typename VertexIndexMap, typename EdgeWeightMap > +struct maximum_weighted_matching_context +{ + typedef typename graph_traits< Graph >::vertex_descriptor vertex_t; + typedef typename graph_traits< Graph >::edge_descriptor edge_t; + typedef typename graph_traits< Graph >::vertices_size_type vertices_size_t; + typedef typename graph_traits< Graph >::edges_size_type edges_size_t; + typedef typename property_traits< EdgeWeightMap >::value_type weight_t; + + /** Ordered pair of vertices. */ + typedef std::pair< vertex_t, vertex_t > vertex_pair_t; + + /** + * List of edges forming an alternating path or alternating cycle. + * + * The path is defined over top-level blossoms; it skips parts of the path + * that are internal to blossoms. Vertex pairs are oriented to match the + * direction of the path. + */ + typedef std::deque< vertex_pair_t > alternating_path_t; + + /** Top-level blossoms may be labeled "S" or "T" or unlabeled. */ + enum blossom_label_t { LABEL_NONE = 0, LABEL_S = 1, LABEL_T = 2 }; + + struct nontrivial_blossom_t; // forward declaration + + /** + * A blossom is either a single vertex, or an odd-length alternating + * cycle over sub-blossoms. + */ + struct blossom_t { - if (v == graph_traits< Graph >::null_vertex()) - return false; - blossom_ptr_t vb = in_blossom[v]->father; - while (vb != blossom_ptr_t()) + /** Parent of this blossom, or "nullptr" for top-level blossom. */ + nontrivial_blossom_t* parent; + + /** + * Base vertex of this blossom. + * + * This is the unique vertex inside the blossom which is not matched + * to another vertex in the same blossom. + */ + vertex_t base_vertex; + + /** Label S or T or NONE, only valid for top-level blossoms. */ + blossom_label_t label; + + /** True if this is an instance of nontrivial_blossom. */ + const bool is_nontrivial_blossom; + + /** Edge that attaches this blossom to the alternating tree. */ + vertex_pair_t tree_edge; + + /** Base vertex of the blossom at the root of the alternating tree. */ + vertex_t tree_root; + + /** Least-slack edge to a different S-blossom. */ + optional< edge_t > best_edge; + + /** Initialize a blossom instance. */ + blossom_t(vertex_t base_vertex = graph_traits< Graph >::null_vertex(), + bool is_nontrivial_blossom = false) + : parent(nullptr) + , base_vertex(base_vertex) + , label(LABEL_NONE) + , is_nontrivial_blossom(is_nontrivial_blossom) + { } + + /** + * Cast this blossom_t instance to nontrivial_blossom_t if possible, + * otherwise return "nullptr". + */ + nontrivial_blossom_t* nontrivial() { - if (vb.get() == b.get()) - return true; - vb = vb->father; + // This is much faster than dynamic_cast. + return (is_nontrivial_blossom ? + static_cast< nontrivial_blossom_t* >(this) : nullptr); } - return false; - } - // return the base vertex of the top blossom that contains v - inline vertex_descriptor_t base_vertex(vertex_descriptor_t v) const - { - return in_top_blossom(v)->get_base(); - } - - // add an existed top blossom of base vertex v into new top - // blossom b as its sub-blossom - void add_sub_blossom(blossom_ptr_t b, vertex_descriptor_t v) - { - blossom_ptr_t sub = in_top_blossom(v); - sub->father = b; - b->sub_blossoms.push_back(sub); - if (sub->sub_blossoms.empty()) - return; - for (blossom_iterator_t bi = top_blossoms.begin(); - bi != top_blossoms.end(); ++bi) + const nontrivial_blossom_t* nontrivial() const { - if (bi->get() == sub.get()) - { - top_blossoms.erase(bi); - break; - } + return (is_nontrivial_blossom ? + static_cast< const nontrivial_blossom_t* >(this) : nullptr); } - } + }; - // when a top blossom is created or its base vertex getting an S-label, - // add all edges incident to this blossom into even_edges - void bloom(blossom_ptr_t b) + /** + * A non-trivial blossom is an odd-length alternating cycle over + * 3 or more sub-blossoms. + */ + struct nontrivial_blossom_t : public blossom_t { - std::vector< vertex_descriptor_t > vertices_of_b = b->vertices(); - vertex_vec_iter_t vi; - for (vi = vertices_of_b.begin(); vi != vertices_of_b.end(); ++vi) + struct sub_blossom_t { + /** Pointer to sub-blossom. */ + blossom_t* blossom; + + /** Vertex pair (x, y) where "x" is a vertex in "blossom" and + "y" is a vertex in the next sub-blossom. */ + vertex_pair_t edge; + }; + + /** List of sub-blossoms, ordered along the alternating cycle. */ + std::list< sub_blossom_t > subblossoms; + + /** Dual LPP variable for this blossom. */ + weight_t dual_var; + + /** Least-slack edges to other S-blossoms. */ + std::list< edge_t > best_edge_set; + + /** Initialize a non-trivial blossom. */ + nontrivial_blossom_t( + const std::vector< blossom_t* >& blossoms, + const std::deque< vertex_pair_t >& edges) + : blossom_t(blossoms.front()->base_vertex, true) + , dual_var(0) { - out_edge_iterator_t oei, oei_end; - for (boost::tie(oei, oei_end) = out_edges(*vi, g); oei != oei_end; - ++oei) - { - if (target(*oei, g) != *vi && mate[*vi] != target(*oei, g)) - even_edges.push_back(*oei); + BOOST_ASSERT(blossoms.size() == edges.size()); + BOOST_ASSERT(blossoms.size() % 2 == 1); + BOOST_ASSERT(blossoms.size() >= 3); + + auto blossom_it = blossoms.begin(); + auto blossom_end = blossoms.end(); + auto edge_it = edges.begin(); + while (blossom_it != blossom_end) { + subblossoms.push_back({*blossom_it, *edge_it}); + ++blossom_it; + ++edge_it; } } - } - - // assigning a T-label to a non S-vertex, along with outlet and updating pi - // value if updated pi[v] equals zero, augment the matching from its mate - // vertex - void put_T_label(vertex_descriptor_t v, vertex_descriptor_t T_label, - vertex_descriptor_t outlet_v, edge_property_t pi_v) - { - if (label_S[v] != graph_traits< Graph >::null_vertex()) - return; - label_T[v] = T_label; - outlet[v] = outlet_v; - pi[v] = pi_v; - - vertex_descriptor_t v_mate = mate[v]; - if (pi[v] == 0) + /** Find the position of the specified subblossom. */ + std::pair< vertices_size_t, + typename std::list< sub_blossom_t >::iterator > + find_subblossom(blossom_t* child) { - label_T[v_mate] = graph_traits< Graph >::null_vertex(); - label_S[v_mate] = v; - bloom(in_top_blossom(v_mate)); + vertices_size_t pos = 0; + auto it = subblossoms.begin(); + while (it->blossom != child) { + ++it; + ++pos; + BOOST_ASSERT(it != subblossoms.end()); + } + return std::make_pair(pos, it); } - } + }; - // get the missing T-label for a to-be-expanded base vertex - // the missing T-label is the last vertex of the path from outlet[v] to v - std::pair< vertex_descriptor_t, vertex_descriptor_t > missing_label( - vertex_descriptor_t b_base) + /** Specification of a delta step. */ + struct delta_step_t { - vertex_descriptor_t missing_outlet = outlet[b_base]; + /** Type of delta step: 1, 2, 3 or 4. */ + int kind; - if (outlet[b_base] == b_base) - return std::make_pair( - graph_traits< Graph >::null_vertex(), missing_outlet); + /** Delta value. */ + weight_t value; - vertex_iterator_t vi, vi_end; - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) - old_label[*vi] = std::make_pair(label_T[*vi], outlet[*vi]); + /** Edge on which the minimum delta occurs (for delta type 2 or 3). */ + edge_t edge; - std::pair< vertex_descriptor_t, vertex_state_t > child( - outlet[b_base], graph::detail::V_EVEN); - blossom_ptr_t b = in_blossom[child.first]; - for (; b->father->father != blossom_ptr_t(); b = b->father) - ; - child.first = b->get_base(); - - if (child.first == b_base) - return std::make_pair( - graph_traits< Graph >::null_vertex(), missing_outlet); + /** Blossom on which the minimum delta occurs (for delta type 4). */ + nontrivial_blossom_t* blossom; + }; - while (true) - { - std::pair< vertex_descriptor_t, vertex_state_t > child_parent - = parent(child, true); + /** Similar to vector_property_map, but uses a fixed-size vector. */ + template < typename T > + struct vertex_map + { + typedef typename property_traits::key_type key_type; + std::vector< T > vec; + VertexIndexMap vm; - for (b = in_blossom[child_parent.first]; - b->father->father != blossom_ptr_t(); b = b->father) - ; - missing_outlet = child_parent.first; - child_parent.first = b->get_base(); + vertex_map(vertices_size_t arg_size, VertexIndexMap arg_vm) + : vec(arg_size) + , vm(arg_vm) + { } - if (child_parent.first == b_base) - break; - else - child = child_parent; + const T& operator[](const key_type& v) const + { + return vec[get(vm, v)]; } - return std::make_pair(child.first, missing_outlet); - } - // expand a top blossom, put all its non-trivial sub-blossoms into - // top_blossoms - blossom_iterator_t expand_blossom( - blossom_iterator_t bi, std::vector< blossom_ptr_t >& new_ones) - { - blossom_ptr_t b = *bi; - for (blossom_iterator_t i = b->sub_blossoms.begin(); - i != b->sub_blossoms.end(); ++i) + T& operator[](const key_type& v) { - blossom_ptr_t sub_blossom = *i; - vertex_descriptor_t sub_base = sub_blossom->get_base(); - label_S[sub_base] = label_T[sub_base] - = graph_traits< Graph >::null_vertex(); - outlet[sub_base] = sub_base; - sub_blossom->father = blossom_ptr_t(); - // new top blossoms cannot be pushed back into top_blossoms - // immediately, because push_back() may cause reallocation and then - // invalid iterators - if (!sub_blossom->sub_blossoms.empty()) - new_ones.push_back(sub_blossom); + return vec[get(vm, v)]; } - return top_blossoms.erase(bi); - } + }; - // when expanding a T-blossom with base v, it requires more operations: - // supply the missing T-labels for new base vertices by picking the minimum - // tau from vertices of each corresponding new top-blossoms; when label_T[v] - // is null or we have a smaller tau from missing_label(v), replace T-label - // and outlet of v (but don't bloom v) - blossom_iterator_t expand_T_blossom( - blossom_iterator_t bi, std::vector< blossom_ptr_t >& new_ones) + /** Keep track of the least-slack edge out of a bunch of edges. */ + struct least_slack_edge_t { - blossom_ptr_t b = *bi; - - vertex_descriptor_t b_base = b->get_base(); - std::pair< vertex_descriptor_t, vertex_descriptor_t > T_and_outlet - = missing_label(b_base); + optional< edge_t > edge; + weight_t slack; - blossom_iterator_t next_bi = expand_blossom(bi, new_ones); + least_slack_edge_t() : slack(0) {} - for (blossom_iterator_t i = b->sub_blossoms.begin(); - i != b->sub_blossoms.end(); ++i) + void update(const edge_t& e, weight_t s) { - blossom_ptr_t sub_blossom = *i; - vertex_descriptor_t sub_base = sub_blossom->get_base(); - vertex_descriptor_t min_tau_v - = graph_traits< Graph >::null_vertex(); - edge_property_t min_tau - = std::numeric_limits< edge_property_t >::max(); - - std::vector< vertex_descriptor_t > sub_vertices - = sub_blossom->vertices(); - for (vertex_vec_iter_t v = sub_vertices.begin(); - v != sub_vertices.end(); ++v) + if ((! edge.has_value()) || (s < slack)) { - if (tau[*v] < min_tau) - { - min_tau = tau[*v]; - min_tau_v = *v; - } + edge = e; + slack = s; } - - if (min_tau < std::numeric_limits< edge_property_t >::max()) - put_T_label( - sub_base, tau_idx[min_tau_v], min_tau_v, tau[min_tau_v]); } + }; - if (label_T[b_base] == graph_traits< Graph >::null_vertex() - || tau[old_label[b_base].second] < pi[b_base]) - boost::tie(label_T[b_base], outlet[b_base]) = T_and_outlet; - - return next_bi; - } + /** Scale integer edge weights to enable integer-only calculations. */ + static constexpr weight_t weight_factor = + std::numeric_limits< weight_t >::is_integer ? 2 : 1; - // when vertices v and w are matched to each other by augmenting, - // we must set v/w as base vertex of any blossom who contains v/w and - // is a sub-blossom of their lowest (smallest) common blossom - void adjust_blossom(vertex_descriptor_t v, vertex_descriptor_t w) + /** Input graph. */ + const Graph& g; + const VertexIndexMap vm; + const EdgeWeightMap edge_weights; + + const vertex_t null_vertex; + + /** Current mate of each vertex. */ + vertex_map< vertex_t > vertex_mate; + + /** + * For each vertex, the trivial blossom that contains it. + * + * Pointers to blossoms must remain valid for the life time of + * this data structure, therefore the underlying vector must + * have a fixed size. + */ + vertex_map< blossom_t > trivial_blossom; + + /** + * List of non-trivial blossoms. + * + * This must be a linked list to ensure that elements can be added + * and removed without invalidating pointers to other elements. + */ + std::list< nontrivial_blossom_t > nontrivial_blossom; + + /** For each vertex, the unique top-level blossom that contains it. */ + vertex_map< blossom_t* > vertex_top_blossom; + + /** For each vertex, a variable in the dual LPP. */ + vertex_map< weight_t > vertex_dual; + + /** For T-vertex or unlabeled vertex, least-slack edge to any S-vertex. */ + vertex_map< optional< edge_t > > vertex_best_edge; + + /** Queue of S-vertices to be scanned. */ + std::deque< vertex_t > scan_queue; + + /** Initialize the matching algorithm. */ + explicit maximum_weighted_matching_context( + const Graph& arg_g, VertexIndexMap arg_vm, EdgeWeightMap weights) + : g(arg_g) + , vm(arg_vm) + , edge_weights(weights) + , null_vertex(graph_traits< Graph >::null_vertex()) + , vertex_mate(num_vertices(g), arg_vm) + , trivial_blossom(num_vertices(g), arg_vm) + , vertex_top_blossom(num_vertices(g), arg_vm) + , vertex_dual(num_vertices(g), arg_vm) + , vertex_best_edge(num_vertices(g), arg_vm) { - blossom_ptr_t vb = in_blossom[v], wb = in_blossom[w], - lowest_common_blossom; - std::vector< blossom_ptr_t > v_ancestors, w_ancestors; + // Vertex duals are initialized to half the maximum edge weight. + weight_t max_weight = 0; + for (const edge_t& e : make_iterator_range(edges(g))) + max_weight = (std::max)(max_weight, get(weights, e)); + weight_t init_vertex_dual = max_weight * (weight_factor / 2); - while (vb->father != blossom_ptr_t()) - { - v_ancestors.push_back(vb->father); - vb = vb->father; - } - while (wb->father != blossom_ptr_t()) + for (const vertex_t& x : make_iterator_range(vertices(g))) { - w_ancestors.push_back(wb->father); - wb = wb->father; + vertex_mate[x] = null_vertex; + trivial_blossom[x].base_vertex = x; + vertex_top_blossom[x] = &trivial_blossom[x]; + vertex_dual[x] = init_vertex_dual; } + } - typename std::vector< blossom_ptr_t >::reverse_iterator i, j; - i = v_ancestors.rbegin(); - j = w_ancestors.rbegin(); - while (i != v_ancestors.rend() && j != w_ancestors.rend() - && i->get() == j->get()) - { - lowest_common_blossom = *i; - ++i; - ++j; + /** Call a function for every vertex inside the specified blossom. */ + template < typename Func > + static void for_vertices_in_blossom(const blossom_t* blossom, Func func) + { + const nontrivial_blossom_t* ntb = blossom->nontrivial(); + if (ntb) { + // Visit all vertices in the non-trivial blossom. + // Use an explicit stack to avoid deep call chains. + std::vector< const nontrivial_blossom_t* > stack; + stack.push_back(ntb); + while (! stack.empty()) { + const nontrivial_blossom_t* cur = stack.back(); + stack.pop_back(); + for (const auto& sub : cur->subblossoms) { + ntb = sub.blossom->nontrivial(); + if (ntb) { + stack.push_back(ntb); + } else { + func(sub.blossom->base_vertex); + } + } + } + } else { + // A trivial blossom contains just one vertex. + func(blossom->base_vertex); } + } - vb = in_blossom[v]; - wb = in_blossom[w]; - while (vb->father != lowest_common_blossom) - { - vb->father->set_base(vb); - vb = vb->father; - } - while (wb->father != lowest_common_blossom) - { - wb->father->set_base(wb); - wb = wb->father; - } + /** Mark blossom as the top-level blossom of its vertices. */ + void update_top_level_blossom(blossom_t* blossom) + { + BOOST_ASSERT(! blossom->parent); + for_vertices_in_blossom(blossom, + [this,blossom](vertex_t x) + { + vertex_top_blossom[x] = blossom; + }); } - // every edge weight is multiplied by 4 to ensure integer weights - // throughout the algorithm if all input weights are integers - inline edge_property_t slack(const edge_descriptor_t& e) const + /** + * Calculate edge slack. + * The two vertices must be in different top-level blossoms. + */ + weight_t edge_slack(const edge_t& e) const { - vertex_descriptor_t v, w; - v = source(e, g); - w = target(e, g); - return dual_var[v] + dual_var[w] - 4 * get(edge_weight, g, e); + vertex_t x = source(e, g); + vertex_t y = target(e, g); + weight_t w = get(edge_weights, e); + BOOST_ASSERT(vertex_top_blossom[x] != vertex_top_blossom[y]); + return vertex_dual[x] + vertex_dual[y] - weight_factor * w; } - // backtrace one step on vertex v along the augmenting path - // by its labels and its vertex state; - // boolean parameter "use_old" means whether we are updating labels, - // if we are, then we use old labels to backtrace and also we - // don't jump to its base vertex when we reach an odd vertex - std::pair< vertex_descriptor_t, vertex_state_t > parent( - std::pair< vertex_descriptor_t, vertex_state_t > v, - bool use_old = false) const + /** Clear least-slack edge tracking. */ + void clear_best_edge() { - if (v.second == graph::detail::V_EVEN) + for (vertex_t x : make_iterator_range(vertices(g))) { - // a paranoid check: label_S shoule be the same as mate in - // backtracing - if (label_S[v.first] == graph_traits< Graph >::null_vertex()) - label_S[v.first] = mate[v.first]; - return std::make_pair(label_S[v.first], graph::detail::V_ODD); + vertex_best_edge[x].reset(); + trivial_blossom[x].best_edge.reset(); } - else if (v.second == graph::detail::V_ODD) + + for (nontrivial_blossom_t& b : nontrivial_blossom) { - vertex_descriptor_t w = use_old ? old_label[v.first].first - : base_vertex(label_T[v.first]); - return std::make_pair(w, graph::detail::V_EVEN); + b.best_edge.reset(); + b.best_edge_set.clear(); } - return std::make_pair(v.first, graph::detail::V_UNREACHED); } - // backtrace from vertices v and w to their free (unmatched) ancesters, - // return the nearest common ancestor (null_vertex if none) of v and w - vertex_descriptor_t nearest_common_ancestor(vertex_descriptor_t v, - vertex_descriptor_t w, vertex_descriptor_t& v_free_ancestor, - vertex_descriptor_t& w_free_ancestor) const + /** Add edge from unlabeled verter or T-vertex "y" to an S-vertex. */ + void add_delta2_edge(vertex_t y, const edge_t& e, weight_t slack) { - std::pair< vertex_descriptor_t, vertex_state_t > v_up( - v, graph::detail::V_EVEN); - std::pair< vertex_descriptor_t, vertex_state_t > w_up( - w, graph::detail::V_EVEN); - vertex_descriptor_t nca; - nca = w_free_ancestor = v_free_ancestor - = graph_traits< Graph >::null_vertex(); - - std::vector< bool > ancestor_of_w_vector(num_vertices(g), false); - std::vector< bool > ancestor_of_v_vector(num_vertices(g), false); - vertex_to_bool_map_t ancestor_of_w(ancestor_of_w_vector.begin(), vm); - vertex_to_bool_map_t ancestor_of_v(ancestor_of_v_vector.begin(), vm); - - while (nca == graph_traits< Graph >::null_vertex() - && (v_free_ancestor == graph_traits< Graph >::null_vertex() - || w_free_ancestor == graph_traits< Graph >::null_vertex())) + auto& best_edge = vertex_best_edge[y]; + if ((! best_edge.has_value()) || slack < edge_slack(*best_edge)) + best_edge = e; + } + + /** Return least-slack edge between any unlabeled vertex and S-vertex. */ + least_slack_edge_t get_least_slack_delta2_edge() + { + least_slack_edge_t best_edge; + for (vertex_t x : make_iterator_range(vertices(g))) { - ancestor_of_w[w_up.first] = true; - ancestor_of_v[v_up.first] = true; - - if (w_free_ancestor == graph_traits< Graph >::null_vertex()) - w_up = parent(w_up); - if (v_free_ancestor == graph_traits< Graph >::null_vertex()) - v_up = parent(v_up); - - if (mate[v_up.first] == graph_traits< Graph >::null_vertex()) - v_free_ancestor = v_up.first; - if (mate[w_up.first] == graph_traits< Graph >::null_vertex()) - w_free_ancestor = w_up.first; - - if (ancestor_of_w[v_up.first] == true || v_up.first == w_up.first) - nca = v_up.first; - else if (ancestor_of_v[w_up.first] == true) - nca = w_up.first; - else if (v_free_ancestor == w_free_ancestor - && v_free_ancestor != graph_traits< Graph >::null_vertex()) - nca = v_up.first; + if (vertex_top_blossom[x]->label == LABEL_NONE) + { + const auto& vertex_edge = vertex_best_edge[x]; + if (vertex_edge.has_value()) + best_edge.update(*vertex_edge, edge_slack(*vertex_edge)); + } } + return best_edge; + } - return nca; + /** Add edge between different S-blossoms. */ + void add_delta3_edge(blossom_t* b, const edge_t& e, weight_t slack) + { + auto& best_edge = b->best_edge; + if ((! best_edge.has_value()) || slack < edge_slack(*best_edge)) + best_edge = e; + + nontrivial_blossom_t* ntb = b->nontrivial(); + if (ntb) + ntb->best_edge_set.push_back(e); } - // when a new top blossom b is created by connecting (v, w), we add - // sub-blossoms into b along backtracing from v_prime and w_prime to - // stop_vertex (the base vertex); also, we set labels and outlet for each - // base vertex we pass by - void make_blossom(blossom_ptr_t b, vertex_descriptor_t w_prime, - vertex_descriptor_t v_prime, vertex_descriptor_t stop_vertex) + /** Update least-slack edge tracking after merging blossoms. */ + void merge_delta3_blossoms(nontrivial_blossom_t* blossom) { - std::pair< vertex_descriptor_t, vertex_state_t > u( - v_prime, graph::detail::V_ODD); - std::pair< vertex_descriptor_t, vertex_state_t > u_up( - w_prime, graph::detail::V_EVEN); + // Build a temporary array holding the least-slack edges to + // other S-blossoms. The array is indexed by base vertex. + std::vector< least_slack_edge_t > tmp_best_edge(num_vertices(g)); - for (; u_up.first != stop_vertex; u = u_up, u_up = parent(u)) + // Collect edges from sub-blossoms that were S-blossoms. + for (auto& sub : blossom->subblossoms) { - if (u_up.second == graph::detail::V_EVEN) + blossom_t* b = sub.blossom; + if (b->label == LABEL_S) { - if (!in_top_blossom(u_up.first)->sub_blossoms.empty()) - outlet[u_up.first] = label_T[u.first]; - label_T[u_up.first] = outlet[u.first]; + b->best_edge.reset(); + nontrivial_blossom_t* ntb = b->nontrivial(); + if (ntb) + { + // Use least-slack edges from subblossom. + for (const edge_t& e : ntb->best_edge_set) + { + blossom_t* bx = vertex_top_blossom[source(e, g)]; + blossom_t* by = vertex_top_blossom[target(e, g)]; + BOOST_ASSERT(bx == blossom); + // Only use edges between top-level blossoms. + if (bx != by) + { + BOOST_ASSERT(by->label == LABEL_S); + tmp_best_edge[get(vm, by->base_vertex)].update( + e, edge_slack(e)); + } + } + ntb->best_edge_set.clear(); + } + else + { + // Trivial blossoms don't maintain a least-slack edge set. + // Consider all incident edges. + for (const edge_t& e : + make_iterator_range(out_edges(b->base_vertex, g))) + { + BOOST_ASSERT(source(e, g) == b->base_vertex); + vertex_t y = target(e, g); + blossom_t* by = vertex_top_blossom[y]; + // Only use edges to different S-blossoms. + // Ignore edges with negative weight. + if ((by != blossom) + && (by->label == LABEL_S) + && (get(edge_weights, e) >= 0)) + { + tmp_best_edge[get(vm, by->base_vertex)].update( + e, edge_slack(e)); + } + } + } } - else if (u_up.second == graph::detail::V_ODD) - label_S[u_up.first] = u.first; + } - add_sub_blossom(b, u_up.first); + // Build a compact list of edges from the temporary array. + // Also find the overall least-slack edge to any other S-blossom. + BOOST_ASSERT(blossom->best_edge_set.empty()); + BOOST_ASSERT(! blossom->best_edge.has_value()); + least_slack_edge_t best_edge; + for (const least_slack_edge_t& item : tmp_best_edge) + { + if (item.edge.has_value()) + { + blossom->best_edge_set.push_back(*item.edge); + best_edge.update(*item.edge, item.slack); + } } + blossom->best_edge = best_edge.edge; } - // the design of recursively expanding augmenting path in - // (reversed_)retrieve_augmenting_path functions is inspired by same - // functions in max_cardinality_matching.hpp; except that in weighted - // matching, we use "outlet" vertices instead of "bridge" vertex pairs: if - // blossom b is the smallest non-trivial blossom that contains its base - // vertex v, then v and outlet[v] are where augmenting path enters and - // leaves b - void retrieve_augmenting_path( - vertex_descriptor_t v, vertex_descriptor_t w, vertex_state_t v_state) + /** Return least-slack edge between any pair of S-blossoms. */ + least_slack_edge_t get_least_slack_delta3_edge() { - if (v == w) - aug_path.push_back(v); - else if (v_state == graph::detail::V_EVEN) - { - aug_path.push_back(v); - retrieve_augmenting_path(label_S[v], w, graph::detail::V_ODD); - } - else if (v_state == graph::detail::V_ODD) + least_slack_edge_t best_edge; + + for (vertex_t x : make_iterator_range(vertices(g))) { - if (outlet[v] == v) - aug_path.push_back(v); - else - reversed_retrieve_augmenting_path( - outlet[v], v, graph::detail::V_EVEN); - retrieve_augmenting_path(label_T[v], w, graph::detail::V_EVEN); + blossom_t* b = vertex_top_blossom[x]; + BOOST_ASSERT(! b->parent); + if ((b->label == LABEL_S) && (b->best_edge.has_value())) + best_edge.update(*b->best_edge, edge_slack(*b->best_edge)); } + + return best_edge; + } + + /** Add the vertices in a blossom to the scan queue. */ + void add_vertices_to_scan_queue(blossom_t* blossom) + { + for_vertices_in_blossom(blossom, + [this](vertex_t x) + { + scan_queue.push_back(x); + }); } - void reversed_retrieve_augmenting_path( - vertex_descriptor_t v, vertex_descriptor_t w, vertex_state_t v_state) + /** Trace back through the alternating trees from vertices "x" and "y". */ + alternating_path_t trace_alternating_paths(vertex_t x, vertex_t y) { - if (v == w) - aug_path.push_back(v); - else if (v_state == graph::detail::V_EVEN) + // Initialize a path containing only the edge (x, y). + alternating_path_t path; + path.emplace_back(x, y); + + // Trace alternating path from "x" to the root of the tree. + while (x != null_vertex) { - reversed_retrieve_augmenting_path( - label_S[v], w, graph::detail::V_ODD); - aug_path.push_back(v); + blossom_t* bx = vertex_top_blossom[x]; + x = bx->tree_edge.first; + if (x != null_vertex) + path.push_front(bx->tree_edge); } - else if (v_state == graph::detail::V_ODD) + + // Trace alternating path from "y" to the root of the tree. + while (y != null_vertex) { - reversed_retrieve_augmenting_path( - label_T[v], w, graph::detail::V_EVEN); - if (outlet[v] != v) - retrieve_augmenting_path(outlet[v], v, graph::detail::V_EVEN); - else - aug_path.push_back(v); + blossom_t* by = vertex_top_blossom[y]; + y = by->tree_edge.first; + if (y != null_vertex) + path.emplace_back(by->tree_edge.second, y); } + + // If we found a common ancestor, trim the paths so they end there. + while (path.front().second == path.back().first) + { + BOOST_ASSERT(path.size() > 2); + path.pop_front(); + path.pop_back(); + } + + // Any alternating path between S-blossoms must have odd length. + BOOST_ASSERT(path.size() % 2 == 1); + + return path; } - // correct labels for vertices in the augmenting path - void relabel(vertex_descriptor_t v) + /** Create a new S-blossom from an alternating cycle. */ + void make_blossom(const alternating_path_t& path) { - blossom_ptr_t b = in_blossom[v]->father; - - if (!is_in_blossom(b, mate[v])) - { // if v is a new base vertex - std::pair< vertex_descriptor_t, vertex_state_t > u( - v, graph::detail::V_EVEN); - while (label_S[u.first] != u.first - && is_in_blossom(b, label_S[u.first])) - u = parent(u, true); - - vertex_descriptor_t old_base = u.first; - if (label_S[old_base] != old_base) - { // if old base is not exposed - label_T[v] = label_S[old_base]; - outlet[v] = old_base; - } - else - { // if old base is exposed then new label_T[v] is not in b, - // we must (i) make b2 the smallest blossom containing v but not - // as base vertex (ii) backtrace from b2's new base vertex to b - label_T[v] = graph_traits< Graph >::null_vertex(); - for (b = b->father; b != blossom_ptr_t() && b->get_base() == v; - b = b->father) - ; - if (b != blossom_ptr_t()) - { - u = std::make_pair(b->get_base(), graph::detail::V_ODD); - while (!is_in_blossom( - in_blossom[v]->father, old_label[u.first].first)) - u = parent(u, true); - label_T[v] = u.first; - outlet[v] = old_label[u.first].first; - } - } + BOOST_ASSERT(path.size() % 2 == 1); + BOOST_ASSERT(path.size() >= 3); + + // Collect pointers to sub-blossoms. + std::vector< blossom_t* > subblossoms; + subblossoms.reserve(path.size()); + for (const vertex_pair_t& edge : path) + subblossoms.push_back(vertex_top_blossom[edge.first]); + + // Check that the path is cyclic. + vertices_size_t pos = 0; + for (const vertex_pair_t& edge : path) + { + pos = (pos + 1) % subblossoms.size(); + BOOST_ASSERT(vertex_top_blossom[edge.second] == subblossoms[pos]); } - else if (label_S[v] == v || !is_in_blossom(b, label_S[v])) - { // if v is an old base vertex - // let u be the new base vertex; backtrace from u's old T-label - std::pair< vertex_descriptor_t, vertex_state_t > u( - b->get_base(), graph::detail::V_ODD); - while ( - old_label[u.first].first != graph_traits< Graph >::null_vertex() - && old_label[u.first].first != v) - u = parent(u, true); - label_T[v] = old_label[u.first].second; - outlet[v] = v; + + // Create the new blossom. + nontrivial_blossom.emplace_back(subblossoms, path); + nontrivial_blossom_t* blossom = &nontrivial_blossom.back(); + + // Link the subblossoms to the new parent. + // Insert former T-vertices into the scan queue. + for (blossom_t* b : subblossoms) + { + BOOST_ASSERT(! b->parent); + b->parent = blossom; + if (b->label == LABEL_T) + add_vertices_to_scan_queue(b); } - else // if v is neither a new nor an old base vertex - label_T[v] = label_S[v]; + + // Mark vertices as belonging to the new blossom. + update_top_level_blossom(blossom); + + // Assign label S to the new blossom and link to the alternating tree. + BOOST_ASSERT(subblossoms.front()->label == LABEL_S); + blossom->label = LABEL_S; + blossom->tree_edge = subblossoms.front()->tree_edge; + blossom->tree_root = subblossoms.front()->tree_root; + + // Merge least-slack edges for the S-sub-blossoms. + merge_delta3_blossoms(blossom); } - void augmenting(vertex_descriptor_t v, vertex_descriptor_t v_free_ancestor, - vertex_descriptor_t w, vertex_descriptor_t w_free_ancestor) + /** Expand a T-blossom. */ + void expand_t_blossom(nontrivial_blossom_t* blossom) { - vertex_iterator_t vi, vi_end; + BOOST_ASSERT(! blossom->parent); + BOOST_ASSERT(blossom->label == LABEL_T); - // retrieve the augmenting path and put it in aug_path - reversed_retrieve_augmenting_path( - v, v_free_ancestor, graph::detail::V_EVEN); - retrieve_augmenting_path(w, w_free_ancestor, graph::detail::V_EVEN); - - // augment the matching along aug_path - vertex_descriptor_t a, b; - vertex_list_t reversed_aug_path; - while (!aug_path.empty()) + // Convert sub-blossoms into top-level blossoms. + for (const auto& sub : blossom->subblossoms) { - a = aug_path.front(); - aug_path.pop_front(); - reversed_aug_path.push_back(a); - b = aug_path.front(); - aug_path.pop_front(); - reversed_aug_path.push_back(b); - - mate[a] = b; - mate[b] = a; - - // reset base vertex for every blossom in augment path - adjust_blossom(a, b); + blossom_t* b = sub.blossom; + BOOST_ASSERT(b->parent == blossom); + b->parent = nullptr; + b->label = LABEL_NONE; + update_top_level_blossom(b); } - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) - old_label[*vi] = std::make_pair(label_T[*vi], outlet[*vi]); + // Reconstruct the alternating tree via the sub-blossoms. + // Find the sub-blossom that attaches the expanding blossom to + // the alternating tree. + blossom_t* entry = vertex_top_blossom[blossom->tree_edge.second]; - // correct labels for in-blossom vertices along aug_path - while (!reversed_aug_path.empty()) - { - a = reversed_aug_path.front(); - reversed_aug_path.pop_front(); + auto subblossom_loc = blossom->find_subblossom(entry); + auto entry_pos = subblossom_loc.first; + auto entry_it = subblossom_loc.second; - if (in_blossom[a]->father != blossom_ptr_t()) - relabel(a); - } + // Get the edge that attached this blossom to the alternating tree. + vertex_t x, y; + std::tie(x, y) = blossom->tree_edge; - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) + // Reconstruct the tree in an even number of steps from entry to base. + auto sub_it = entry_it; + if (entry_pos % 2 == 0) { - vertex_descriptor_t u = *vi; - if (mate[u] != graph_traits< Graph >::null_vertex()) - label_S[u] = mate[u]; + // Walk backward around the blossom. + auto sub_begin = blossom->subblossoms.begin(); + while (sub_it != sub_begin) + { + extend_tree_s_to_t(x, y); + --sub_it; + BOOST_ASSERT(sub_it != sub_begin); + --sub_it; + std::tie(y, x) = sub_it->edge; + } } - - // expand blossoms with zero dual variables - std::vector< blossom_ptr_t > new_top_blossoms; - for (blossom_iterator_t bi = top_blossoms.begin(); - bi != top_blossoms.end();) + else { - if ((*bi)->dual_var <= 0) - bi = expand_blossom(bi, new_top_blossoms); - else - ++bi; + // Walk forward around the blossom. + auto sub_end = blossom->subblossoms.end(); + while (sub_it != sub_end) { + extend_tree_s_to_t(x, y); + ++sub_it; + BOOST_ASSERT(sub_it != sub_end); + std::tie(x, y) = sub_it->edge; + ++sub_it; + } } - top_blossoms.insert(top_blossoms.end(), new_top_blossoms.begin(), - new_top_blossoms.end()); - init(); + + // Assign label T to the base sub-blossom. + blossom_t* base = blossom->subblossoms.front().blossom; + base->label = LABEL_T; + base->tree_edge = std::make_pair(x, y); + base->tree_root = blossom->tree_root; + + // Delete the expanded blossom. + auto blossom_it = std::find_if( + nontrivial_blossom.begin(), + nontrivial_blossom.end(), + [blossom](const nontrivial_blossom_t& b) + { + return (&b == blossom); + }); + BOOST_ASSERT(blossom_it != nontrivial_blossom.end()); + nontrivial_blossom.erase(blossom_it); } - // create a new blossom and set labels for vertices inside - void blossoming(vertex_descriptor_t v, vertex_descriptor_t v_prime, - vertex_descriptor_t w, vertex_descriptor_t w_prime, - vertex_descriptor_t nca) + void augment_blossom_rec( + nontrivial_blossom_t* blossom, blossom_t* entry, + std::stack< std::pair< nontrivial_blossom_t*, blossom_t* > >& stack) { - vertex_iterator_t vi, vi_end; + auto subblossom_loc = blossom->find_subblossom(entry); + auto entry_pos = subblossom_loc.first; + auto entry_it = subblossom_loc.second; + + // Walk from entry to the base in an even number of steps. + auto sub_begin = blossom->subblossoms.begin(); + auto sub_end = blossom->subblossoms.end(); + auto sub_it = entry_it; + while ((sub_it != sub_begin) && (sub_it != sub_end)) { + vertex_t x, y; + blossom_t* bx; + blossom_t* by; + + if (entry_pos % 2 == 0) + { + // Walk backward around the blossom. + --sub_it; + by = sub_it->blossom; + BOOST_ASSERT(sub_it != sub_begin); + --sub_it; + bx = sub_it->blossom; + std::tie(x, y) = sub_it->edge; + } + else + { + // Walk forward around the blossom. + ++sub_it; + BOOST_ASSERT(sub_it != sub_end); + std::tie(x, y) = sub_it->edge; + bx = sub_it->blossom; + ++sub_it; + by = (sub_it == sub_end) ? + blossom->subblossoms.front().blossom : + sub_it->blossom; + } - std::vector< bool > is_old_base_vector(num_vertices(g)); - vertex_to_bool_map_t is_old_base(is_old_base_vector.begin(), vm); - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) - { - if (*vi == base_vertex(*vi)) - is_old_base[*vi] = true; + // Pull this edge into the matching. + vertex_mate[x] = y; + vertex_mate[y] = x; + + // Augment through any non-trivial subblossoms touching this edge. + nontrivial_blossom_t* bx_ntb = bx->nontrivial(); + if (bx_ntb) + stack.emplace(bx_ntb, &trivial_blossom[x]); + + nontrivial_blossom_t* by_ntb = by->nontrivial(); + if (by_ntb) + stack.emplace(by_ntb, &trivial_blossom[y]); } - blossom_ptr_t b = boost::make_shared< blossom >(); - add_sub_blossom(b, nca); + // Re-orient the blossom. + if (entry_it != sub_begin) + blossom->subblossoms.splice( + sub_begin, blossom->subblossoms, entry_it, sub_end); + blossom->base_vertex = entry->base_vertex; + } - label_T[w_prime] = v; - label_T[v_prime] = w; - outlet[w_prime] = w; - outlet[v_prime] = v; + void augment_blossom(nontrivial_blossom_t* blossom, blossom_t* entry) + { + // Use an explicit stack to avoid deep call chains. + std::stack< std::pair< nontrivial_blossom_t*, blossom_t* > > stack; + stack.emplace(blossom, entry); - make_blossom(b, w_prime, v_prime, nca); - make_blossom(b, v_prime, w_prime, nca); + while (! stack.empty()) { + nontrivial_blossom_t* outer_blossom; + blossom_t* inner_entry; + std::tie(outer_blossom, inner_entry) = stack.top(); - label_T[nca] = graph_traits< Graph >::null_vertex(); - outlet[nca] = nca; + nontrivial_blossom_t* inner_blossom = inner_entry->parent; + BOOST_ASSERT(inner_blossom); - top_blossoms.push_back(b); - bloom(b); + if (inner_blossom != outer_blossom) + stack.top() = std::make_pair(outer_blossom, inner_blossom); + else + stack.pop(); - // set gamma[b_base] = min_slack{critical_edge(b_base, other_base)} - // where each critical edge is updated before, by - // argmin{slack(old_bases_in_b, other_base)}; - vertex_vec_iter_t i, j; - std::vector< vertex_descriptor_t > b_vertices = b->vertices(), - old_base_in_b, other_base; - vertex_descriptor_t b_base = b->get_base(); - for (i = b_vertices.begin(); i != b_vertices.end(); ++i) + augment_blossom_rec(inner_blossom, inner_entry, stack); + } + } + + /** Augment the matching through the specified augmenting path. */ + void augment_matching(const alternating_path_t& path) + { + BOOST_ASSERT(path.size() % 2 == 1); + + // Process the unmatched edges on the augmenting path. + auto path_it = path.begin(); + auto path_end = path.end(); + while (path_it != path_end) { - if (is_old_base[*i]) - old_base_in_b.push_back(*i); + vertex_t x = path_it->first; + vertex_t y = path_it->second; + + // Augment any non-trivial blossoms that touch this edge. + blossom_t* bx = vertex_top_blossom[x]; + nontrivial_blossom_t* bx_ntb = bx->nontrivial(); + if (bx_ntb) + augment_blossom(bx_ntb, &trivial_blossom[x]); + + blossom_t* by = vertex_top_blossom[y]; + nontrivial_blossom_t* by_ntb = by->nontrivial(); + if (by_ntb) + augment_blossom(by_ntb, &trivial_blossom[y]); + + // Pull this edge into the matching. + vertex_mate[x] = y; + vertex_mate[y] = x; + + // Go to the next unmatched edge on the path. + ++path_it; + if (path_it == path_end) + break; + ++path_it; } - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) + } + + /** + * Remove non-S-vertices from the scan queue. + * This is necessary after removing labels from S-blossoms. + */ + void refresh_scan_queue() + { + std::deque< vertex_t > new_scan_queue; + for (const vertex_t& x : scan_queue) { - if (*vi != b_base && *vi == base_vertex(*vi)) - other_base.push_back(*vi); + if (vertex_top_blossom[x]->label == LABEL_S) + new_scan_queue.push_back(x); } - for (i = other_base.begin(); i != other_base.end(); ++i) + scan_queue = std::move(new_scan_queue); + } + + /** Remove edges to non-S-vertices from delta3 edge tracking. */ + void refresh_delta3_blossom(blossom_t* b) + { + // Do nothing if this blossom was not tracking any delta3 edge. + if (! b->best_edge.has_value()) + return; + + nontrivial_blossom_t* ntb = b->nontrivial(); + if (ntb) { - edge_property_t min_slack - = std::numeric_limits< edge_property_t >::max(); - std::pair< edge_descriptor_t, bool > b_vi = null_edge; - for (j = old_base_in_b.begin(); j != old_base_in_b.end(); ++j) + // Remove bad edges from best_edge_set and refresh best_edge. + least_slack_edge_t best_edge; + auto it = ntb->best_edge_set.begin(); + auto it_end = ntb->best_edge_set.end(); + while (it != it_end) { - if (critical_edge[*j][*i] != null_edge - && min_slack > slack(critical_edge[*j][*i].first)) + vertex_t y = target(*it, g); + blossom_t* by = vertex_top_blossom[y]; + BOOST_ASSERT(by != b); + if (by->label == LABEL_S) + { + best_edge.update(*it, edge_slack(*it)); + ++it; + } + else { - min_slack = slack(critical_edge[*j][*i].first); - b_vi = critical_edge[*j][*i]; + // Remove edge to non-S-blossom. + it = ntb->best_edge_set.erase(it); } } - critical_edge[b_base][*i] = critical_edge[*i][b_base] = b_vi; + b->best_edge = best_edge.edge; } - gamma[b_base] = std::numeric_limits< edge_property_t >::max(); - for (i = other_base.begin(); i != other_base.end(); ++i) + else { - if (critical_edge[b_base][*i] != null_edge) - gamma[b_base] = std::min( - gamma[b_base], slack(critical_edge[b_base][*i].first)); + // Trivial blossom does not maintain best_edge_set. + // If its current best_edge is invalid, recompute it. + vertex_t x = b->base_vertex; + vertex_t y = target(*b->best_edge, g); + blossom_t* by = vertex_top_blossom[y]; + BOOST_ASSERT(by != b); + if (by->label != LABEL_S) + { + // Consider all incident edges to recompute best_edge. + least_slack_edge_t best_edge; + for (const edge_t& e : make_iterator_range(out_edges(x, g))) + { + BOOST_ASSERT(source(e, g) == x); + y = target(e, g); + by = vertex_top_blossom[y]; + if ((by != b) && (by->label == LABEL_S)) + best_edge.update(e, edge_slack(e)); + } + b->best_edge = best_edge.edge; + } } } - void init() + /** Recompute vertex_best_edge for the specified vertex. */ + void recompute_vertex_best_edge(vertex_t x) { - even_edges.clear(); + least_slack_edge_t best_edge; - vertex_iterator_t vi, vi_end; - typename std::vector< - std::vector< std::pair< edge_descriptor_t, bool > > >::iterator vei; + for (const edge_t& e : make_iterator_range(out_edges(x, g))) + { + BOOST_ASSERT(source(e, g) == x); + vertex_t y = target(e, g); + blossom_t* by = vertex_top_blossom[y]; + if (by->label == LABEL_S) + best_edge.update(e, edge_slack(e)); + } + vertex_best_edge[x] = best_edge.edge; + } - for (boost::tie(vi, vi_end) = vertices(g), - vei = critical_edge_vectors.begin(); - vi != vi_end; ++vi, ++vei) + /** Remove the alternating trees with specified root vertices. */ + void remove_alternating_tree(vertex_t r1, vertex_t r2) + { + // Find blossoms that are part of the specified alternating trees. + std::vector< blossom_t* > former_s_blossoms; + for (vertex_t x : make_iterator_range(vertices(g))) { - vertex_descriptor_t u = *vi; - out_edge_iterator_t ei, ei_end; + blossom_t* b = vertex_top_blossom[x]; + if ((! b->parent) + && (b->label != LABEL_NONE) + && (b->tree_root == r1 || b->tree_root == r2)) + { + // Build list of former S-blossoms. + if (b->label == LABEL_S) + former_s_blossoms.push_back(b); - gamma[u] = tau[u] = pi[u] - = std::numeric_limits< edge_property_t >::max(); - std::fill(vei->begin(), vei->end(), null_edge); + // Remove label. + b->label = LABEL_NONE; + } + } - if (base_vertex(u) != u) - continue; + vertex_map< char > blossom_recompute_best_edge(num_vertices(g), vm); + vertex_map< char > vertex_recompute_best_edge(num_vertices(g), vm); - label_S[u] = label_T[u] = graph_traits< Graph >::null_vertex(); - outlet[u] = u; + // Visit former S-blossoms. + for (blossom_t* b : former_s_blossoms) + { + // Clear best_edge tracking. + b->best_edge.reset(); + nontrivial_blossom_t* ntb = b->nontrivial(); + if (ntb) + ntb->best_edge_set.clear(); + + // Visit edges that touch this blossom. + for_vertices_in_blossom(b, + [&](vertex_t x) + { + // Mark former S-vertices. + vertex_recompute_best_edge[x] = 1; - if (mate[u] == graph_traits< Graph >::null_vertex()) - { - label_S[u] = u; - bloom(in_top_blossom(u)); - } + for (const edge_t& e : make_iterator_range(out_edges(x, g))) + { + // Mark S-blossoms adjacent to the former S-blossom. + BOOST_ASSERT(source(e, g) == x); + vertex_t y = target(e, g); + blossom_t* by = vertex_top_blossom[y]; + if (by->label == LABEL_S) + blossom_recompute_best_edge[by->base_vertex] = 1; + + // Mark non-S-vertices with least-slack edge to + // former S-blossom. + if (by->label != LABEL_S) + { + const auto& best_edge = vertex_best_edge[y]; + if (best_edge.has_value() && (*best_edge == e)) + vertex_recompute_best_edge[y] = 1; + } + } + }); } + + // Recompute delta3 tracking of affected S-blossoms. + for (vertex_t x : make_iterator_range(vertices(g))) + { + if (blossom_recompute_best_edge[x]) + refresh_delta3_blossom(vertex_top_blossom[x]); + } + + // Recompute vertex_best_edge of affected non-S-vertices. + for (vertex_t x : make_iterator_range(vertices(g))) + { + if (vertex_recompute_best_edge[x]) + recompute_vertex_best_edge(x); + } + + refresh_scan_queue(); } - bool augment_matching() + /** + * Extend the alternating tree via the edge from S-vertex "x" + * to unlabeled vertex "y". + * + * Assign label T to the blossom that contains "y", then assign + * label S to the blossom matched to the newly labeled T-blossom. + */ + void extend_tree_s_to_t(vertex_t x, vertex_t y) { - vertex_descriptor_t v, w, w_free_ancestor, v_free_ancestor; - v = w = w_free_ancestor = v_free_ancestor - = graph_traits< Graph >::null_vertex(); - bool found_alternating_path = false; + blossom_t* bx = vertex_top_blossom[x]; + blossom_t* by = vertex_top_blossom[y]; + + BOOST_ASSERT(bx->label == LABEL_S); + BOOST_ASSERT(by->label == LABEL_NONE); + by->label = LABEL_T; + by->tree_edge = std::make_pair(x, y); + by->tree_root = bx->tree_root; + + vertex_t y2 = by->base_vertex; + vertex_t z = vertex_mate[y2]; + BOOST_ASSERT(z != null_vertex); + + blossom_t* bz = vertex_top_blossom[z]; + BOOST_ASSERT(bz->label == LABEL_NONE); + BOOST_ASSERT(! bz->best_edge.has_value()); + bz->label = LABEL_S; + bz->tree_edge = std::make_pair(y2, z); + bz->tree_root = by->tree_root; + add_vertices_to_scan_queue(bz); + } - // note that we only use edges of zero slack value for augmenting - while (!even_edges.empty() && !found_alternating_path) + /** + * Add the edge between S-vertices "x" and "y". + * + * If the edge connects blossoms that are part of the same alternating + * tree, a new S-blossom is created. + * + * If the edge connects two different alternating trees, an augmenting + * path has been discovered. In this case the matching is augmented. + * + * @return true if the matching was augmented; otherwise false. + */ + bool add_s_to_s_edge(vertex_t x, vertex_t y) + { + blossom_t* bx = vertex_top_blossom[x]; + blossom_t* by = vertex_top_blossom[y]; + BOOST_ASSERT(bx->label == LABEL_S); + BOOST_ASSERT(by->label == LABEL_S); + BOOST_ASSERT(bx != by); + + alternating_path_t path = trace_alternating_paths(x, y); + + // Check whether the path starts and ends in the same blossom. + blossom_t* b1 = vertex_top_blossom[path.front().first]; + blossom_t* b2 = vertex_top_blossom[path.back().second]; + if (b1 == b2) + { + make_blossom(path); + return false; + } + else { - // search for augmenting paths depth-first - edge_descriptor_t current_edge = even_edges.back(); - even_edges.pop_back(); + // Remove labels the two alternating trees on the augmenting path. + remove_alternating_tree(bx->tree_root, by->tree_root); - v = source(current_edge, g); - w = target(current_edge, g); + // Augment matching. + augment_matching(path); + return true; + } + } - vertex_descriptor_t v_prime = base_vertex(v); - vertex_descriptor_t w_prime = base_vertex(w); + /** + * Scan incident edges of newly labeled S-vertices. + * + * @return true if the matching was augmented; otherwise false. + */ + bool scan_new_s_vertices() + { + while (! scan_queue.empty()) + { + vertex_t x = scan_queue.front(); + scan_queue.pop_front(); - // w_prime == v_prime implies that we get an edge that has been - // shrunk into a blossom - if (v_prime == w_prime) - continue; + BOOST_ASSERT(vertex_top_blossom[x]->label == LABEL_S); - // a paranoid check - if (label_S[v_prime] == graph_traits< Graph >::null_vertex()) + for (const edge_t& e : make_iterator_range(out_edges(x, g))) { - std::swap(v_prime, w_prime); - std::swap(v, w); - } + BOOST_ASSERT(source(e, g) == x); + vertex_t y = target(e, g); - // w_prime may be unlabeled or have a T-label; replace the existed - // T-label if the edge slack is smaller than current pi[w_prime] and - // update it. Note that a T-label is "deserved" only when pi equals - // zero. also update tau and tau_idx so that tau_idx becomes T-label - // when a T-blossom is expanded - if (label_S[w_prime] == graph_traits< Graph >::null_vertex()) - { - if (slack(current_edge) < pi[w_prime]) - put_T_label(w_prime, v, w, slack(current_edge)); - if (slack(current_edge) < tau[w]) - { - if (in_blossom[w]->father == blossom_ptr_t() - || label_T[w_prime] == v - || label_T[w_prime] - == graph_traits< Graph >::null_vertex() - || nearest_common_ancestor(v_prime, label_T[w_prime], - v_free_ancestor, w_free_ancestor) - == graph_traits< Graph >::null_vertex()) - { - tau[w] = slack(current_edge); - tau_idx[w] = v; - } - } - } + // Note: top level blossom of x may change during this loop, + // so we must look it up on each pass. + blossom_t* bx = vertex_top_blossom[x]; + blossom_t* by = vertex_top_blossom[y]; - else - { - if (slack(current_edge) > 0) + // Ignore internal edges in blossom. + if (bx == by) + continue; + + // Ignore edges with negative slack to prevent numeric overflow. + if (get(edge_weights, e) < 0) + continue; + + weight_t slack = edge_slack(e); + if (slack <= 0) { - // update gamma and critical_edges when we have a smaller - // edge slack - gamma[v_prime] - = std::min(gamma[v_prime], slack(current_edge)); - gamma[w_prime] - = std::min(gamma[w_prime], slack(current_edge)); - if (critical_edge[v_prime][w_prime] == null_edge - || slack(critical_edge[v_prime][w_prime].first) - > slack(current_edge)) + // Tight edge. + if (by->label == LABEL_NONE) + extend_tree_s_to_t(x, y); + else if (by->label == LABEL_S) { - critical_edge[v_prime][w_prime] - = std::pair< edge_descriptor_t, bool >( - current_edge, true); - critical_edge[w_prime][v_prime] - = std::pair< edge_descriptor_t, bool >( - current_edge, true); + bool augmented = add_s_to_s_edge(x, y); + if (augmented) + return true; } - continue; } - else if (slack(current_edge) == 0) + else { - // if nca is null_vertex then we have an augmenting path; - // otherwise we have a new top blossom with nca as its base - // vertex - vertex_descriptor_t nca = nearest_common_ancestor( - v_prime, w_prime, v_free_ancestor, w_free_ancestor); - - if (nca == graph_traits< Graph >::null_vertex()) - found_alternating_path - = true; // to break out of the loop - else - blossoming(v, v_prime, w, w_prime, nca); + // Track non-tight edges between S-blossoms. + if (by->label == LABEL_S) + add_delta3_edge(bx, e, slack); } + + // Track all edges between S-blossom and non-S-blossom. + if (by->label != LABEL_S) + add_delta2_edge(y, e, slack); } } - if (!found_alternating_path) - return false; - - augmenting(v, v_free_ancestor, w, w_free_ancestor); - return true; + return false; } - // slack the vertex and blossom dual variables when there is no augmenting - // path found according to the primal-dual method - bool adjust_dual() + /** Calculate a delta step in the dual LPP problem. */ + delta_step_t calc_dual_delta_step() { - edge_property_t delta1, delta2, delta3, delta4, delta; - delta1 = delta2 = delta3 = delta4 - = std::numeric_limits< edge_property_t >::max(); - - vertex_iterator_t vi, vi_end; + delta_step_t delta{}; - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) + // Compute delta1: minimum dual variable of any S-vertex. + delta.kind = 1; + delta.value = (std::numeric_limits< weight_t >::max)(); + for (vertex_t x : make_iterator_range(vertices(g))) { - delta1 = std::min(delta1, dual_var[*vi]); - delta4 = pi[*vi] > 0 ? std::min(delta4, pi[*vi]) : delta4; - if (*vi == base_vertex(*vi)) - delta3 = std::min(delta3, gamma[*vi] / 2); + if (vertex_top_blossom[x]->label == LABEL_S) + delta.value = (std::min)(delta.value, vertex_dual[x]); } - for (blossom_iterator_t bi = top_blossoms.begin(); - bi != top_blossoms.end(); ++bi) + // Compute delta2: minimum slack of edge from S-vertex to unlabeled. + least_slack_edge_t best_edge = get_least_slack_delta2_edge(); + if (best_edge.edge.has_value() && (best_edge.slack <= delta.value)) { - vertex_descriptor_t b_base = (*bi)->get_base(); - if (label_T[b_base] != graph_traits< Graph >::null_vertex() - && pi[b_base] == 0) - delta2 = std::min(delta2, (*bi)->dual_var / 2); + delta.kind = 2; + delta.edge = *best_edge.edge; + delta.value = best_edge.slack; } - delta = std::min(std::min(delta1, delta2), std::min(delta3, delta4)); - - // start updating dual variables, note that the order is important - - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) + // Compute delta3: half minimum slack of edge between S-blossoms. + best_edge = get_least_slack_delta3_edge(); + if (best_edge.edge.has_value() && (best_edge.slack / 2 <= delta.value)) { - vertex_descriptor_t v = *vi, v_prime = base_vertex(v); - - if (label_S[v_prime] != graph_traits< Graph >::null_vertex()) - dual_var[v] -= delta; - else if (label_T[v_prime] != graph_traits< Graph >::null_vertex() - && pi[v_prime] == 0) - dual_var[v] += delta; - - if (v == v_prime) - gamma[v] -= 2 * delta; + delta.kind = 3; + delta.edge = *best_edge.edge; + delta.value = best_edge.slack / 2; } - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) + // Compute delta4: half minimum dual of a top-level T-blossom. + for (nontrivial_blossom_t& blossom : nontrivial_blossom) { - vertex_descriptor_t v_prime = base_vertex(*vi); - if (pi[v_prime] > 0) - tau[*vi] -= delta; + if ((! blossom.parent) + && (blossom.label == LABEL_T) + && (blossom.dual_var / 2 <= delta.value)) + { + delta.kind = 4; + delta.blossom = &blossom; + delta.value = blossom.dual_var / 2; + } } - for (blossom_iterator_t bi = top_blossoms.begin(); - bi != top_blossoms.end(); ++bi) - { - vertex_descriptor_t b_base = (*bi)->get_base(); - if (label_T[b_base] != graph_traits< Graph >::null_vertex() - && pi[b_base] == 0) - (*bi)->dual_var -= 2 * delta; - if (label_S[b_base] != graph_traits< Graph >::null_vertex()) - (*bi)->dual_var += 2 * delta; - } + return delta; + } - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) + /** Apply a delta step to the dual LPP variables. */ + void apply_delta_step(weight_t delta) + { + for (vertex_t x : make_iterator_range(vertices(g))) { - vertex_descriptor_t v = *vi; - if (pi[v] > 0) - pi[v] -= delta; - - // when some T-vertices have zero pi value, bloom their mates so - // that matching can be further augmented - if (label_T[v] != graph_traits< Graph >::null_vertex() - && pi[v] == 0) - put_T_label(v, label_T[v], outlet[v], pi[v]); + blossom_t* b = vertex_top_blossom[x]; + if (b->label == LABEL_S) + vertex_dual[x] -= delta; + else if (b->label == LABEL_T) + vertex_dual[x] += delta; } - // optimal solution reached, halt - if (delta == delta1) - return false; - - // expand odd blossoms with zero dual variables and zero pi value of - // their base vertices - if (delta == delta2 && delta != delta3) + for (nontrivial_blossom_t& blossom : nontrivial_blossom) { - std::vector< blossom_ptr_t > new_top_blossoms; - for (blossom_iterator_t bi = top_blossoms.begin(); - bi != top_blossoms.end();) + if (! blossom.parent) { - const blossom_ptr_t b = *bi; - vertex_descriptor_t b_base = b->get_base(); - if (b->dual_var == 0 - && label_T[b_base] != graph_traits< Graph >::null_vertex() - && pi[b_base] == 0) - bi = expand_T_blossom(bi, new_top_blossoms); - else - ++bi; + if (blossom.label == LABEL_S) + blossom.dual_var += 2 * delta; + else if (blossom.label == LABEL_T) + blossom.dual_var -= 2 * delta; } - top_blossoms.insert(top_blossoms.end(), new_top_blossoms.begin(), - new_top_blossoms.end()); } + } - while (true) - { - // find a zero-slack critical edge (v, w) of zero gamma values - std::pair< edge_descriptor_t, bool > best_edge = null_edge; - std::vector< vertex_descriptor_t > base_nodes; - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) + /** + * Run one stage of the matching algorithm. + * + * @return True if the matching was successfully augmented; + * false if no further improvement is possible. + */ + bool run_stage() + { + // Run substages. + while (true) { + bool augmented = scan_new_s_vertices(); + if (augmented) + return true; + + delta_step_t delta = calc_dual_delta_step(); + apply_delta_step(delta.value); + + if (delta.kind == 2) { - if (*vi == base_vertex(*vi)) - base_nodes.push_back(*vi); + vertex_t x = source(delta.edge, g); + vertex_t y = target(delta.edge, g); + if (vertex_top_blossom[x]->label != LABEL_S) + std::swap(x, y); + extend_tree_s_to_t(x, y); } - for (vertex_vec_iter_t i = base_nodes.begin(); - i != base_nodes.end(); ++i) + else if (delta.kind == 3) { - if (gamma[*i] == 0) - { - for (vertex_vec_iter_t j = base_nodes.begin(); - j != base_nodes.end(); ++j) - { - if (critical_edge[*i][*j] != null_edge - && slack(critical_edge[*i][*j].first) == 0) - best_edge = critical_edge[*i][*j]; - } - } + vertex_t x = source(delta.edge, g); + vertex_t y = target(delta.edge, g); + augmented = add_s_to_s_edge(x, y); + if (augmented) + return true; } - - // if not found, continue finding other augment matching - if (best_edge == null_edge) + else if (delta.kind == 4) { - bool augmented = augment_matching(); - return augmented || delta != delta1; + BOOST_ASSERT(delta.blossom); + expand_t_blossom(delta.blossom); } - // if found, determine either augmenting or blossoming - vertex_descriptor_t v = source(best_edge.first, g), - w = target(best_edge.first, g); - vertex_descriptor_t v_prime = base_vertex(v), - w_prime = base_vertex(w), v_free_ancestor, - w_free_ancestor; - vertex_descriptor_t nca = nearest_common_ancestor( - v_prime, w_prime, v_free_ancestor, w_free_ancestor); - if (nca == graph_traits< Graph >::null_vertex()) + else { - augmenting(v, v_free_ancestor, w, w_free_ancestor); - return true; + // No further improvement possible. + BOOST_ASSERT(delta.kind == 1); + return false; } - else - blossoming(v, v_prime, w, w_prime, nca); } - - return false; } - template < typename PropertyMap > void get_current_matching(PropertyMap pm) + /** Run the matching algorithm. */ + void run() { - vertex_iterator_t vi, vi_end; - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) - put(pm, *vi, mate[*vi]); - } + // Assign label S to all vertices and put them in the queue. + for (vertex_t x : make_iterator_range(vertices(g))) + { + BOOST_ASSERT(vertex_mate[x] == null_vertex); + blossom_t* bx = vertex_top_blossom[x]; + BOOST_ASSERT(bx->label == LABEL_NONE); + BOOST_ASSERT(bx->base_vertex == x); + bx->label = LABEL_S; + bx->tree_edge = std::make_pair(null_vertex, x); + bx->tree_root = x; + scan_queue.push_back(x); + } -private: - const Graph& g; - VertexIndexMap vm; - const std::pair< edge_descriptor_t, bool > null_edge; - - // storage for the property maps below - std::vector< vertex_descriptor_t > mate_vector; - std::vector< vertex_descriptor_t > label_S_vector, label_T_vector; - std::vector< vertex_descriptor_t > outlet_vector; - std::vector< vertex_descriptor_t > tau_idx_vector; - std::vector< edge_property_t > dual_var_vector; - std::vector< edge_property_t > pi_vector, gamma_vector, tau_vector; - std::vector< blossom_ptr_t > in_blossom_vector; - std::vector< std::pair< vertex_descriptor_t, vertex_descriptor_t > > - old_label_vector; - std::vector< vertex_to_edge_map_t > critical_edge_vector; - std::vector< std::vector< std::pair< edge_descriptor_t, bool > > > - critical_edge_vectors; - - // iterator property maps - vertex_to_vertex_map_t mate; - vertex_to_vertex_map_t label_S; // v has an S-label -> v can be an even - // vertex, label_S[v] is its mate - vertex_to_vertex_map_t - label_T; // v has a T-label -> v can be an odd vertex, label_T[v] is its - // predecessor in aug_path - vertex_to_vertex_map_t outlet; - vertex_to_vertex_map_t tau_idx; - vertex_to_weight_map_t dual_var; - vertex_to_weight_map_t pi, gamma, tau; - vertex_to_blossom_map_t - in_blossom; // map any vertex v to the trivial blossom containing v - vertex_to_pair_map_t old_label; // before - // relabeling or expanding T-blossoms - vertex_pair_to_edge_map_t - critical_edge; // an not matched edge (v, w) is critical if v and w - // belongs to different S-blossoms - - vertex_list_t aug_path; - edge_list_t even_edges; - std::vector< blossom_ptr_t > top_blossoms; -}; + // Improve the solution until no further improvement is possible. + while (run_stage()) ; -template < typename Graph, typename MateMap, typename VertexIndexMap > -void maximum_weighted_matching(const Graph& g, MateMap mate, VertexIndexMap vm) -{ - empty_matching< Graph, MateMap >::find_matching(g, mate); - weighted_augmenting_path_finder< Graph, MateMap, VertexIndexMap > augmentor( - g, mate, vm); + // Clear temporary data structures. + scan_queue.clear(); + clear_best_edge(); + } - // can have |V| times augmenting at most - for (std::size_t t = 0; t < num_vertices(g); ++t) + /** Copy matching to specified map. */ + template < typename MateMap > + void extract_matching(MateMap mate) { - bool augmented = false; - while (!augmented) + for (vertex_t x : make_iterator_range(vertices(g))) + put(mate, x, vertex_mate[x]); + } + + /** Check that the array "vertex_mate" is consistent. */ + bool verify_vertex_mate() + { + vertices_size_t num_matched_vertex = 0; + for (vertex_t x : make_iterator_range(vertices(g))) { - augmented = augmentor.augment_matching(); - if (!augmented) + vertex_t y = vertex_mate[x]; + if (y != null_vertex) { - // halt if adjusting dual variables can't bring potential - // augment - if (!augmentor.adjust_dual()) - break; + ++num_matched_vertex; + if (vertex_mate[y] != x) + return false; } } - if (!augmented) - break; - } - - augmentor.get_current_matching(mate); -} -template < typename Graph, typename MateMap > -inline void maximum_weighted_matching(const Graph& g, MateMap mate) -{ - maximum_weighted_matching(g, mate, get(vertex_index, g)); -} + edges_size_t num_matched_edge = 0; + for (const edge_t& e : make_iterator_range(edges(g))) + { + if (vertex_mate[source(e, g)] == target(e, g)) + ++num_matched_edge; + } -// brute-force matcher searches all possible combinations of matched edges to -// get the maximum weighted matching which can be used for testing on small -// graphs (within dozens vertices) -template < typename Graph, typename MateMap, typename VertexIndexMap > -class brute_force_matching -{ -public: - typedef - typename graph_traits< Graph >::vertex_descriptor vertex_descriptor_t; - typedef typename graph_traits< Graph >::vertex_iterator vertex_iterator_t; - typedef - typename std::vector< vertex_descriptor_t >::iterator vertex_vec_iter_t; - typedef typename graph_traits< Graph >::edge_iterator edge_iterator_t; - typedef boost::iterator_property_map< vertex_vec_iter_t, VertexIndexMap > - vertex_to_vertex_map_t; + return num_matched_vertex == 2 * num_matched_edge; + } - brute_force_matching( - const Graph& arg_g, MateMap arg_mate, VertexIndexMap arg_vm) - : g(arg_g) - , vm(arg_vm) - , mate_vector(num_vertices(g)) - , best_mate_vector(num_vertices(g)) - , mate(mate_vector.begin(), vm) - , best_mate(best_mate_vector.begin(), vm) + /** + * Verify that vertex dual variables are non-negative, + * and all unmatched vertices have zero dual. + */ + bool verify_vertex_duals() { - vertex_iterator_t vi, vi_end; - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) - best_mate[*vi] = mate[*vi] = get(arg_mate, *vi); + for (vertex_t x : make_iterator_range(vertices(g))) + { + if (vertex_dual[x] < 0) + return false; + if ((vertex_mate[x] == null_vertex) && (vertex_dual[x] != 0)) + return false; + } + return true; } - template < typename PropertyMap > void find_matching(PropertyMap pm) + /** Verify that blossom dual variables are non-negative. */ + bool verify_blossom_duals() { - edge_iterator_t ei; - boost::tie(ei, ei_end) = edges(g); - select_edge(ei); - - vertex_iterator_t vi, vi_end; - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) - put(pm, *vi, best_mate[*vi]); + for (const nontrivial_blossom_t& blossom : nontrivial_blossom) + { + if (blossom.dual_var < 0) + return false; + } + return true; } -private: - const Graph& g; - VertexIndexMap vm; - std::vector< vertex_descriptor_t > mate_vector, best_mate_vector; - vertex_to_vertex_map_t mate, best_mate; - edge_iterator_t ei_end; - - void select_edge(edge_iterator_t ei) + /** Verify slack of edges between trivial blossoms. */ + bool verify_trivial_blossom_edges() { - if (ei == ei_end) + for (const edge_t& e : make_iterator_range(edges(g))) { - if (matching_weight_sum(g, mate) - > matching_weight_sum(g, best_mate)) + vertex_t x = source(e, g); + vertex_t y = target(e, g); + if ((! trivial_blossom[x].parent) && (! trivial_blossom[y].parent)) { - vertex_iterator_t vi, vi_end; - for (boost::tie(vi, vi_end) = vertices(g); vi != vi_end; ++vi) - best_mate[*vi] = mate[*vi]; + weight_t w = weight_factor * get(edge_weights, e); + weight_t dual_sum = vertex_dual[x] + vertex_dual[y]; + // Edge slack must be non-negative. + if (w > dual_sum) + return false; + // Matched edges must have slack zero. + if ((vertex_mate[x] == y) && (w != dual_sum)) + return false; } - return; } + return true; + } - vertex_descriptor_t v, w; - v = source(*ei, g); - w = target(*ei, g); + /** Verify one top-level blossom and its edges. */ + bool verify_blossom(const nontrivial_blossom_t* blossom) + { + // This function recursively descends down the blossom structure. + // It checks the number of matched edges in every sub-blossom. + // It also checks the slack of edges contained in this blossom + // and of edges incident on this blossom. - select_edge(++ei); + // For each vertex, deepest level on current recursion path + // that contains the vertex. + vertex_map< vertices_size_t > vertex_depth(num_vertices(g), vm); - if (mate[v] == graph_traits< Graph >::null_vertex() - && mate[w] == graph_traits< Graph >::null_vertex()) + // At each recursion depth, sum of duals of containing blossoms. + std::vector< weight_t > path_sum_dual = {0}; + + // At each recursion depth, number of internally matched vertices. + std::vector< vertices_size_t > path_num_matched = {0}; + + // Use an explicit stack to avoid deep call chains. + using sub_blossom_iterator = typename std::list< + typename nontrivial_blossom_t::sub_blossom_t >::const_iterator; + std::stack< std::pair< + const nontrivial_blossom_t*, sub_blossom_iterator > > stack; + + stack.emplace(blossom, blossom->subblossoms.cbegin()); + while (! stack.empty()) { - mate[v] = w; - mate[w] = v; - select_edge(ei); - mate[v] = mate[w] = graph_traits< Graph >::null_vertex(); + vertices_size_t depth = stack.size(); + blossom = stack.top().first; + auto subblossom_it = stack.top().second; + + if (subblossom_it == blossom->subblossoms.cbegin()) + { + // Entering a new (sub)-blossom. + if (blossom->subblossoms.size() < 3) { + return false; + } + + // Update the depth of vertices in this sub-blossom. + for_vertices_in_blossom(blossom, + [&vertex_depth,depth](vertex_t x) { + vertex_depth[x] = depth; + }); + + // Calculate the sum of blossom duals at the new depth. + path_sum_dual.push_back( + path_sum_dual.back() + blossom->dual_var); + + // Initialize the number of matched edges at the new depth. + path_num_matched.push_back(0); + } + + if (subblossom_it != blossom->subblossoms.cend()) + { + // Update the sub-blossom pointer at the current depth. + ++(stack.top().second); + + // Examine a sub-blossom. + blossom_t* b = subblossom_it->blossom; + BOOST_ASSERT(b->parent == blossom); + nontrivial_blossom_t* ntb = b->nontrivial(); + if (ntb) + { + // Prepare to descend into the sub-blossom. + stack.emplace(ntb, ntb->subblossoms.cbegin()); + } + else + { + // Handle single vertex. + vertex_t x = b->base_vertex; + for (const edge_t& e : make_iterator_range(out_edges(x, g))) + { + BOOST_ASSERT(source(e, g) == x); + vertex_t y = target(e, g); + // Find the smallest blossom that contains this edge. + vertices_size_t edge_depth = vertex_depth[y]; + // Verify edge slack. + weight_t w = weight_factor * get(edge_weights, e); + weight_t dual_sum = vertex_dual[x] + vertex_dual[y] + + path_sum_dual[edge_depth]; + // Edge slack must be non-negative. + if (w > dual_sum) + return false; + // Matched edges must have slack zero. + if ((vertex_mate[x] == y) && (w != dual_sum)) + return false; + // Update number of internally matched vertices. + if ((vertex_mate[x] == y) && (edge_depth > 0)) + path_num_matched[edge_depth] += 1; + } + } + } + else + { + // Leaving the current (sub)-blossom. + // Count the number of vertices inside this blossom. + vertices_size_t blossom_num_vertex = 0; + for_vertices_in_blossom(blossom, + [&blossom_num_vertex](vertex_t) { + ++blossom_num_vertex; + }); + + // Check that the (sub)-blossom is "full". + // A blossom is full if all except one of its vertices + // are matched to another vertex within the blossom. + vertices_size_t blossom_num_matched = path_num_matched[depth]; + if (blossom_num_vertex != blossom_num_matched + 1) + return false; + + // Update the number of matched edges in the parent blossom. + path_num_matched[depth - 1] += path_num_matched[depth]; + + // Restore the depth of vertices. + for_vertices_in_blossom(blossom, + [&vertex_depth,depth](vertex_t x) { + vertex_depth[x] = depth - 1; + }); + + // Remove the blossom from the stack. + path_sum_dual.pop_back(); + path_num_matched.pop_back(); + stack.pop(); + } + } + + return true; + } + + /** Verify blossoms and their edges. */ + bool verify_blossoms() + { + // Recursively verify each non-trivial top-level blossom. + // This takes total time O(V^2). + for (const nontrivial_blossom_t& blossom : nontrivial_blossom) + { + if (! blossom.parent) + { + if (! verify_blossom(&blossom)) + return false; + } } + return true; + } + + /** Verify that the optimum solution has been found. */ + bool verify() + { + return (verify_vertex_mate() + && verify_vertex_duals() + && verify_blossom_duals() + && verify_trivial_blossom_edges() + && verify_blossoms()); + return true; } }; +} // namespace detail +} // namespace graph + +/** + * Compute a maximum-weighted matching in an undirected weighted graph. + * + * This function takes time O(V^3). + * This function uses memory size O(V + E). + * + * @tparam Verify True to verify the matching. + * This is only supported for integer edge weights. + * Verification will always succeed unless there is a bug + * in the matching algorithm. + * + * @param g Input graph. + * The graph type must be a model of VertexListGraph + * and EdgeListGraph and IncidenceGraph. + * The graph must be undirected. + * The graph may not contain parallel edges. + * + * @param mate ReadWritePropertyMap, mapping vertices to vertices. + * This map returns the result of the computation. + * For each vertex "v", "mate[v]" is the vertex that is + * matched to "v", or "null_vertex()" if "v" is not matched. + * + * @param vm ReadablePropertyMap, mapping vertices to indexes + * in range [0, num_vertices(g)). + * + * @param weights ReadablePropertyMap, mapping edges to edge weights. + * Edge weights must be integers or floating point values. + * + * @throw boost::bad_graph If the input graph is not valid. + */ +template < typename Graph, typename MateMap, typename VertexIndexMap, + typename EdgeWeightMap, bool Verify = true > +void maximum_weighted_matching( + const Graph& g, MateMap mate, VertexIndexMap vm, EdgeWeightMap weights) +{ + BOOST_CONCEPT_ASSERT((VertexListGraphConcept< Graph >)); + BOOST_CONCEPT_ASSERT((EdgeListGraphConcept< Graph >)); + BOOST_CONCEPT_ASSERT((IncidenceGraphConcept< Graph >)); + BOOST_STATIC_ASSERT(is_undirected_graph< Graph >::value); + + typedef typename graph_traits< Graph >::vertex_descriptor vertex_t; + typedef typename graph_traits< Graph >::edge_descriptor edge_t; + BOOST_CONCEPT_ASSERT( + (ReadWritePropertyMapConcept< MateMap, vertex_t >)); + BOOST_CONCEPT_ASSERT( + (ReadablePropertyMapConcept< VertexIndexMap, vertex_t >)); + BOOST_CONCEPT_ASSERT( + (ReadablePropertyMapConcept< EdgeWeightMap, edge_t >)); + + graph::detail::check_vertex_index_range(g, vm); + graph::detail::check_maximum_weighted_matching_edge_weights(g, weights); + + graph::detail::maximum_weighted_matching_context< + Graph, VertexIndexMap, EdgeWeightMap > + matching(g, vm, weights); + matching.run(); + + typedef typename property_traits< EdgeWeightMap >::value_type weight_t; + if (Verify && std::numeric_limits< weight_t >::is_integer) + { + if (! matching.verify()) + throw std::logic_error("Incorrect solution."); + } + + matching.extract_matching(mate); +} + +/** + * Compute a maximum-weighted matching in an undirected weighted graph. + * + * This overloaded function obtains edge weights from "get(edge_weight, g)". + */ template < typename Graph, typename MateMap, typename VertexIndexMap > -void brute_force_maximum_weighted_matching( +inline void maximum_weighted_matching( const Graph& g, MateMap mate, VertexIndexMap vm) { - empty_matching< Graph, MateMap >::find_matching(g, mate); - brute_force_matching< Graph, MateMap, VertexIndexMap > brute_force_matcher( - g, mate, vm); - brute_force_matcher.find_matching(mate); + maximum_weighted_matching(g, mate, vm, get(edge_weight, g)); } +/** + * Compute a maximum-weighted matching in an undirected weighted graph. + * + * This overloaded function obtains vertex indices from "get(vertex_index, g)" + * and edge weights from "get(edge_weight, g)". + */ template < typename Graph, typename MateMap > -inline void brute_force_maximum_weighted_matching(const Graph& g, MateMap mate) +inline void maximum_weighted_matching(const Graph& g, MateMap mate) { - brute_force_maximum_weighted_matching(g, mate, get(vertex_index, g)); + maximum_weighted_matching(g, mate, get(vertex_index, g)); } -} +} // namespace boost -#endif +#endif // BOOST_GRAPH_MAXIMUM_WEIGHTED_MATCHING_HPP From 059365ffe86142d9888983371885dec83422f1f2 Mon Sep 17 00:00:00 2001 From: Joris van Rantwijk Date: Sun, 12 Jan 2025 00:52:52 +0100 Subject: [PATCH 04/18] Update maximum_weighted_matching documentation --- doc/bibliography.html | 5 + doc/maximum_weighted_matching.html | 143 +++++++++++++++++++---------- 2 files changed, 97 insertions(+), 51 deletions(-) diff --git a/doc/bibliography.html b/doc/bibliography.html index b541abf62..e79333f52 100644 --- a/doc/bibliography.html +++ b/doc/bibliography.html @@ -453,6 +453,11 @@

Bibliography

Data Structures for Weighted Matching and Nearest Common Ancestors with Linking
Proceedings of the First Annual ACM-SIAM Symposium on Discrete Algorithms, pp. 434-443, 1990. +

77 +
Zvi Galil
+Efficient Algorithms for Finding Maximum Matching in Graphs
+ACM Computing Surveys (CSUR), 18(1), 23-38, 1986. +
diff --git a/doc/maximum_weighted_matching.html b/doc/maximum_weighted_matching.html index ab40f1f45..01b48fdf7 100644 --- a/doc/maximum_weighted_matching.html +++ b/doc/maximum_weighted_matching.html @@ -1,5 +1,6 @@