23 #ifndef PERSISTENCE_GRAPH_H_ 24 #define PERSISTENCE_GRAPH_H_ 26 #include <gudhi/Internal_point.h> 29 #include <tbb/parallel_sort.h> 38 namespace persistence_diagram {
45 class Persistence_graph {
48 template<
typename Persistence_diagram1,
typename Persistence_diagram2>
49 Persistence_graph(
const Persistence_diagram1& diag1,
const Persistence_diagram2& diag2,
double e);
51 bool on_the_u_diagonal(
int u_point_index)
const;
53 bool on_the_v_diagonal(
int v_point_index)
const;
55 int corresponding_point_in_u(
int v_point_index)
const;
57 int corresponding_point_in_v(
int u_point_index)
const;
59 double distance(
int u_point_index,
int v_point_index)
const;
63 double bottleneck_alive()
const;
65 std::vector<double> sorted_distances()
const;
67 double diameter_bound()
const;
69 Internal_point get_u_point(
int u_point_index)
const;
71 Internal_point get_v_point(
int v_point_index)
const;
74 std::vector<Internal_point> u;
75 std::vector<Internal_point> v;
79 template<
typename Persistence_diagram1,
typename Persistence_diagram2>
80 Persistence_graph::Persistence_graph(
const Persistence_diagram1 &diag1,
81 const Persistence_diagram2 &diag2,
double e)
82 : u(), v(), b_alive(0.) {
83 std::vector<double> u_alive;
84 std::vector<double> v_alive;
85 for (
auto it = std::begin(diag1); it != std::end(diag1); ++it) {
86 if (std::get<1>(*it) == std::numeric_limits<double>::infinity())
87 u_alive.push_back(std::get<0>(*it));
88 else if (std::get<1>(*it) - std::get<0>(*it) > e)
89 u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size()));
91 for (
auto it = std::begin(diag2); it != std::end(diag2); ++it) {
92 if (std::get<1>(*it) == std::numeric_limits<double>::infinity())
93 v_alive.push_back(std::get<0>(*it));
94 else if (std::get<1>(*it) - std::get<0>(*it) > e)
95 v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size()));
97 if (u.size() < v.size())
99 std::sort(u_alive.begin(), u_alive.end());
100 std::sort(v_alive.begin(), v_alive.end());
101 if (u_alive.size() != v_alive.size()) {
102 b_alive = std::numeric_limits<double>::infinity();
104 for (
auto it_u = u_alive.cbegin(), it_v = v_alive.cbegin(); it_u != u_alive.cend(); ++it_u, ++it_v)
105 b_alive = (std::max)(b_alive, std::fabs(*it_u - *it_v));
109 inline bool Persistence_graph::on_the_u_diagonal(
int u_point_index)
const {
110 return u_point_index >=
static_cast<int> (u.size());
113 inline bool Persistence_graph::on_the_v_diagonal(
int v_point_index)
const {
114 return v_point_index >=
static_cast<int> (v.size());
117 inline int Persistence_graph::corresponding_point_in_u(
int v_point_index)
const {
118 return on_the_v_diagonal(v_point_index) ?
119 v_point_index -
static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size());
122 inline int Persistence_graph::corresponding_point_in_v(
int u_point_index)
const {
123 return on_the_u_diagonal(u_point_index) ?
124 u_point_index -
static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size());
127 inline double Persistence_graph::distance(
int u_point_index,
int v_point_index)
const {
128 if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index))
130 Internal_point p_u = get_u_point(u_point_index);
131 Internal_point p_v = get_v_point(v_point_index);
132 return (std::max)(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y()));
135 inline int Persistence_graph::size()
const {
136 return static_cast<int> (u.size() + v.size());
139 inline double Persistence_graph::bottleneck_alive()
const {
143 inline std::vector<double> Persistence_graph::sorted_distances()
const {
144 std::vector<double> distances;
145 distances.push_back(0.);
146 for (
int u_point_index = 0; u_point_index < size(); ++u_point_index) {
147 distances.push_back(distance(u_point_index, corresponding_point_in_v(u_point_index)));
148 for (
int v_point_index = 0; v_point_index < size(); ++v_point_index)
149 distances.push_back(distance(u_point_index, v_point_index));
152 tbb::parallel_sort(distances.begin(), distances.end());
154 std::sort(distances.begin(), distances.end());
159 inline Internal_point Persistence_graph::get_u_point(
int u_point_index)
const {
160 if (!on_the_u_diagonal(u_point_index))
161 return u.at(u_point_index);
162 Internal_point projector = v.at(corresponding_point_in_v(u_point_index));
163 double m = (projector.x() + projector.y()) / 2.;
164 return Internal_point(m, m, u_point_index);
167 inline Internal_point Persistence_graph::get_v_point(
int v_point_index)
const {
168 if (!on_the_v_diagonal(v_point_index))
169 return v.at(v_point_index);
170 Internal_point projector = u.at(corresponding_point_in_u(v_point_index));
171 double m = (projector.x() + projector.y()) / 2.;
172 return Internal_point(m, m, v_point_index);
175 inline double Persistence_graph::diameter_bound()
const {
177 for (
auto it = u.cbegin(); it != u.cend(); it++)
178 max = (std::max)(max, it->y());
179 for (
auto it = v.cbegin(); it != v.cend(); it++)
180 max = (std::max)(max, it->y());
188 #endif // PERSISTENCE_GRAPH_H_ Definition: SimplicialComplexForAlpha.h:26