The original graph looks like
+----+ 5 +----+ 1 +----+ 6
| 1 | ---- | 2 | ---- | 6 |--------------------+
+----+ +----+ +----+ |
| | |
+-----------------------------+-----------------------+ |
| | |
| +----+ 1 +----+ 7 +----+ 5 +----+ 4 +----+ 6 +----+ 7 +----+
| | 0 | ---- | 4 | ---- | 5 | ---- | 9 | ---- | 14 | ---- | 15 | ---- | 10 |
| +----+ +----+ +----+ +----+ +----+ +----+ +----+
| | |
| | 2 |
| | |
| 5 +----+ +----+ 4 +----+ 2 +----+ |
| | 3 | | 8 | ---- | 13 | ---- | 12 | |
| +----+ +----+ +----+ +----+ |
| | 4 |
| | 4 +-----------------------------------------------------------+
| | |
| +----+ 5 +----+ 6 +----+ 5 +----+
+-- | 7 | ---- | 11 | ---- | 16 | ---- | 17 |
+----+ +----+ +----+ +----+
Your graph on the other hand looks like
0.1
+--------------------------------------+
| |
+---+ 0.5 +---+ 0.4 +---+ 0.1 +---+
| 0 | ------ | 2 | ------ | 1 | ------ | 3 |
+---+ +---+ +---+ +---+
I note that the bruteforce method works: Live On Coliru
====== brute_force_maximum_weighted_matching ======
Found a matching:
Matching size is 2, total weight is 0.6
The matching is:
{0, 2}
{1, 3}
Indeed as you mentioned, changing the weights to integers did allow it to work:
add_edge(0, 2, EdgeProperty(/*0.*/5), g);
add_edge(0, 3, EdgeProperty(/*0.*/1), g);
add_edge(1, 2, EdgeProperty(/*0.*/4), g);
add_edge(1, 3, EdgeProperty(/*0.*/1), g);
Live On Coliru printing:
In the following graph:
0.1
+--------------------------------------+
| |
+---+ 0.5 +---+ 0.4 +---+ 0.1 +---+
| 0 | ------ | 2 | ------ | 1 | ------ | 3 |
+---+ +---+ +---+ +---+
====== maximum_weighted_matching ======
Found a matching:
Matching size is 2, total weight is 6
The matching is:
{0, 2}
{1, 3}
====== brute_force_maximum_weighted_matching ======
Found a matching:
Matching size is 2, total weight is 6
The matching is:
{0, 2}
{1, 3}
Why, Though
The $1m question.
I noticed:
None of this made any difference.
Then I suspected floating point accuracy issues. So I started scaling the weights by factors of 10 (nEx for X = 1, 2, 3, ... so that the total weights would be 6, 60, 600, ...). E.g.
add_edge(0, 2, EdgeProperty(5e5), g);
add_edge(0, 3, EdgeProperty(1e5), g);
add_edge(1, 2, EdgeProperty(4e5), g);
add_edge(1, 3, EdgeProperty(1e5), g);
In this approach X = -1 is identical to your question weights
This started failing at X = 10, again running indefinitely. At that point I selected double instead of float again, and lo and behold: Live On Coliru:
====== maximum_weighted_matching ======
Found a matching:
Matching size is 2, total weight is 6e+10
The matching is:
{0, 2}
{1, 3}
So... Floating Point Strikes Again?
Before jumping to conclusions, I read the documentation on maximum_weighted_matching. First off
Both maximum_weighted_matching and brute_force_maximum_weighted_matching find a maximum weighted matching in any undirected graph.
This is good, because I was wondering since your graph doesn't "look like" the examples used.
The maximum weighted matching problem was solved by Edmonds in [74]. The implementation of maximum_weighted_matching followed Chapter 6, Section 10 of [20] and was written in a consistent style with edmonds_maximum_cardinality_matching because of their algorithmic similarity. In addition, a brute-force verifier brute_force_maximum_weighted_matching simply searches all possible matchings in any graph and selects one with the maximum weight sum.
Ah. This increases the likelihood that there is an implementation specific bug/undocumented limitation the edmonds_maximum_cardinality_matching implementation, which is not present in the brute_force_maximum_weighted_matching variant.
WORKAROUND #1: Use the brute_force_maximum_weighted_matching algorithm instead
Now since
For maximum_weighted_matching, the management of blossoms is much more involved than in the case of max_cardinality_matching
I thought to also apply max_cardinality_matching itself:
//maximum_weighted_matching(g, &mate[0]);
boost::edmonds_maximum_cardinality_matching(g, mate.data());
Yes, that works without a hitch for weights that hang up maximum_weighted_matching (obviously, since the weights aren't used). So far, so good.
More Docs...
Why is a verification algorithm needed? Edmonds' algorithm is fairly complex, and it's nearly impossible for a human without a few days of spare time to figure out if the matching produced by edmonds_matching on a graph with, say, 100 vertices and 500 edges is indeed a maximum cardinality matching
Oh. Wow. This doesn't immediately inspire confidence.
But nothing in the docs really gives me any indication why it could hang/suffer extreme worst case behaviour on certain weights.
Since floating point accuracy is at fault/involved, let me present secondary workardounds:
- Workaround #2: Use
long double
- Workaround #2: Use decimal floats
Workaround #2: Use long double
Annoyingly, this Just Works™
Since this kind of "guess-work" solution feels bad to me, let's be a little bit more methodical:
Workaround #3: Use decimal floats
Boost Multiprecision has our backs. But just
using Weight = boost::multiprecision::cpp_dec_float_50;
doesn't cut it, because BGL uses std::min to get the minimum of two expressions, and due to expression templates the template argument cannot be deduced¹.
So short of fixing that¹ lets disable expression templates:
using Weight = // boost::multiprecision::cpp_dec_float_50;
boost::multiprecision::number<
boost::multiprecision::cpp_dec_float<50>,
boost::multiprecision::et_off >;
This also does the trick: Live On Coliru.
SUMMARY
These are three workarounds. I suggest workaround #1 because it is reliable. It might not suit your performance needs, however.
In that case I'd consider Workaround #3 while also reporting the current test case as a bug with the library developers.
Listing
Anti-bitrot listing of the Workaround #3:
#include <boost/graph/adjacency_list.hpp>
#include <boost/multiprecision/cpp_dec_float.hpp>
#include <boost/graph/maximum_weighted_matching.hpp>
#include <iostream>
#include <string>
#include <vector>
using Weight = // boost::multiprecision::cpp_dec_float_50;
boost::multiprecision::number<
boost::multiprecision::cpp_dec_float<50>,
boost::multiprecision::et_off >;
using EdgeProperty =
boost::property<boost::edge_weight_t, Weight>;
using my_graph =
boost::adjacency_list<
boost::vecS,
boost::vecS,
boost::undirectedS,
boost::no_property,
EdgeProperty>;
using V = boost::graph_traits<my_graph>::vertex_descriptor;
using E = boost::graph_traits<my_graph>::edge_descriptor;
static auto report(my_graph const& g, std::vector<V> const& mate) {
auto sum = matching_weight_sum(g, &mate[0]);
std::cout << "Found a matching:" << std::endl;
std::cout << "Matching size is " << matching_size(g, &mate[0])
<< ", total weight is " << sum
<< std::endl;
std::cout << std::endl;
std::cout << "The matching is:" << std::endl;
for (V v : boost::make_iterator_range(vertices(g))) {
if (mate[v] != g.null_vertex() && v < mate[v]) {
std::cout << "{" << v << ", " << mate[v] << "}" << std::endl;
}
}
std::cout << std::endl;
return sum;
}
int main() {
// vertices can be refered by integers because my_graph use vector to store
// them
my_graph g(4);
add_edge(0, 2, EdgeProperty(5e-1), g);
add_edge(0, 3, EdgeProperty(1e-1), g);
add_edge(1, 2, EdgeProperty(4e-1), g);
add_edge(1, 3, EdgeProperty(1e-1), g);
// print the ascii graph into terminal (better to use fixed-width font)
std::cout << R"(In the following graph:
0.1
+--------------------------------------+
| |
+---+ 0.5 +---+ 0.4 +---+ 0.1 +---+
| 0 | ------ | 2 | ------ | 1 | ------ | 3 |
+---+ +---+ +---+ +---+
)" << std::endl << std::endl;
Weight sum1 = 0, sum2 = 0;
if (1) {
std::cout << "====== maximum_weighted_matching ======\n";
std::vector<V> mate(num_vertices(g));
maximum_weighted_matching(g, &mate[0]);
sum1 = report(g, mate);
}
// now we check the correctness by compare the weight sum to a brute-force
// matching result note that two matchings may be different because of
// multiple optimal solutions
if (1) {
std::cout << "====== brute_force_maximum_weighted_matching ======\n";
std::vector<V> mate(num_vertices(g));
brute_force_maximum_weighted_matching(g, &mate[0]);
sum2 = report(g, mate);
}
assert(sum1 == sum2);
}
¹ this could easily be fixed in the implementation as well:
//delta3 = std::min(delta3, gamma[*vi] / 2);
delta3 = std::min<edge_property_t>(delta3, gamma[*vi] / 2); // SEHE WAS HERE
//delta2 = std::min(delta2, (*bi)->dual_var / 2);
delta2 = std::min<edge_property_t>(delta2, (*bi)->dual_var / 2); // SEHE WAS HERE
I might makemade a pull request for this one
