16 #include <gudhi/Graph_matching.h>
18 #include <CGAL/version.h>
28 #if CGAL_VERSION_NR < 1041101000
29 # error bottleneck_distance is only available for CGAL >= 4.11
34 namespace persistence_diagram {
36 inline double bottleneck_distance_approx(Persistence_graph& g,
double e) {
37 double b_lower_bound = 0.;
38 double b_upper_bound = g.diameter_bound();
39 const double alpha = std::pow(g.size(), 1. / 5.);
41 Graph_matching biggest_unperfect(g);
42 while (b_upper_bound - b_lower_bound > 2 * e) {
43 double step = b_lower_bound + (b_upper_bound - b_lower_bound) / alpha;
44 #if !defined FLT_EVAL_METHOD || FLT_EVAL_METHOD < 0 || FLT_EVAL_METHOD > 1
47 volatile double drop_excess_precision = step;
48 step = drop_excess_precision;
51 if (step <= b_lower_bound || step >= b_upper_bound)
54 while (m.multi_augment()) {}
56 m = biggest_unperfect;
59 biggest_unperfect = m;
63 return (b_lower_bound + b_upper_bound) / 2.;
66 inline double bottleneck_distance_exact(Persistence_graph& g) {
67 std::vector<double> sd = g.sorted_distances();
68 long lower_bound_i = 0;
69 long upper_bound_i = sd.size() - 1;
70 const double alpha = std::pow(g.size(), 1. / 5.);
72 Graph_matching biggest_unperfect(g);
73 while (lower_bound_i != upper_bound_i) {
74 long step = lower_bound_i +
static_cast<long> ((upper_bound_i - lower_bound_i - 1) / alpha);
76 while (m.multi_augment()) {}
78 m = biggest_unperfect;
81 biggest_unperfect = m;
82 lower_bound_i = step + 1;
85 return sd.at(lower_bound_i);
111 template<
typename Persistence_diagram1,
typename Persistence_diagram2>
113 double e = (std::numeric_limits<double>::min)()) {
114 Persistence_graph g(diag1, diag2, e);
115 if (g.bottleneck_alive() == std::numeric_limits<double>::infinity())
116 return std::numeric_limits<double>::infinity();
117 return (std::max)(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e));
124 #endif // BOTTLENECK_H_