1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
|
/* This file is part of the Gudhi Library. The Gudhi library
* (Geometric Understanding in Higher Dimensions) is a generic C++
* library for computational topology.
*
* Author: Francois Godi
*
* Copyright (C) 2015 Inria
*
* Modifications:
* - 2019/06 Vincent Rouvreau : Fix doxygen warning.
*
* This program is free software: you can redistribute it and/or modify
* it under the terms of the GNU General Public License as published by
* the Free Software Foundation, either version 3 of the License, or
* (at your option) any later version.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program. If not, see <http://www.gnu.org/licenses/>.
*/
#ifndef BOTTLENECK_H_
#define BOTTLENECK_H_
#include <gudhi/Graph_matching.h>
#include <vector>
#include <algorithm> // for max
#include <limits> // for numeric_limits
#include <cmath>
#include <cfloat> // FLT_EVAL_METHOD
namespace Gudhi {
namespace persistence_diagram {
inline double bottleneck_distance_approx(Persistence_graph& g, double e) {
double b_lower_bound = 0.;
double b_upper_bound = g.diameter_bound();
const double alpha = std::pow(g.size(), 1. / 5.);
Graph_matching m(g);
Graph_matching biggest_unperfect(g);
while (b_upper_bound - b_lower_bound > 2 * e) {
double step = b_lower_bound + (b_upper_bound - b_lower_bound) / alpha;
#if !defined FLT_EVAL_METHOD || FLT_EVAL_METHOD < 0 || FLT_EVAL_METHOD > 1
// On platforms where double computation is done with excess precision,
// we force it to its true precision so the following test is reliable.
volatile double drop_excess_precision = step;
step = drop_excess_precision;
// Alternative: step = CGAL::IA_force_to_double(step);
#endif
if (step <= b_lower_bound || step >= b_upper_bound) // Avoid precision problem
break;
m.set_r(step);
while (m.multi_augment()) {} // compute a maximum matching (in the graph corresponding to the current r)
if (m.perfect()) {
m = biggest_unperfect;
b_upper_bound = step;
} else {
biggest_unperfect = m;
b_lower_bound = step;
}
}
return (b_lower_bound + b_upper_bound) / 2.;
}
inline double bottleneck_distance_exact(Persistence_graph& g) {
std::vector<double> sd = g.sorted_distances();
long lower_bound_i = 0;
long upper_bound_i = sd.size() - 1;
const double alpha = std::pow(g.size(), 1. / 5.);
Graph_matching m(g);
Graph_matching biggest_unperfect(g);
while (lower_bound_i != upper_bound_i) {
long step = lower_bound_i + static_cast<long> ((upper_bound_i - lower_bound_i - 1) / alpha);
m.set_r(sd.at(step));
while (m.multi_augment()) {} // compute a maximum matching (in the graph corresponding to the current r)
if (m.perfect()) {
m = biggest_unperfect;
upper_bound_i = step;
} else {
biggest_unperfect = m;
lower_bound_i = step + 1;
}
}
return sd.at(lower_bound_i);
}
/** \brief Function to compute the Bottleneck distance between two persistence diagrams.
*
* \tparam Persistence_diagram1,Persistence_diagram2
* models of the concept `PersistenceDiagram`.
*
* \param[in] diag1 The first persistence diagram.
* \param[in] diag2 The second persistence diagram.
*
* \param[in] e
* \parblock
* If `e` is 0, this uses an expensive algorithm to compute the exact distance.
*
* If `e` is not 0, it asks for an additive `e`-approximation, and currently
* also allows a small multiplicative error (the last 2 or 3 bits of the
* mantissa may be wrong). This version of the algorithm takes advantage of the
* limited precision of `double` and is usually a lot faster to compute,
* whatever the value of `e`.
*
* Thus, by default, `e` is the smallest positive double.
* \endparblock
*
* \ingroup bottleneck_distance
*/
template<typename Persistence_diagram1, typename Persistence_diagram2>
double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2,
double e = (std::numeric_limits<double>::min)()) {
Persistence_graph g(diag1, diag2, e);
if (g.bottleneck_alive() == std::numeric_limits<double>::infinity())
return std::numeric_limits<double>::infinity();
return (std::max)(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e));
}
} // namespace persistence_diagram
} // namespace Gudhi
#endif // BOTTLENECK_H_
|
