11#ifndef PERSISTENCE_GRAPH_H_
12#define PERSISTENCE_GRAPH_H_
14#include <gudhi/Internal_point.h>
17#include <tbb/parallel_sort.h>
27namespace persistence_diagram {
34class Persistence_graph {
37 template<
typename Persistence_diagram1,
typename Persistence_diagram2>
38 Persistence_graph(
const Persistence_diagram1& diag1,
const Persistence_diagram2& diag2,
double e);
40 bool on_the_u_diagonal(
int u_point_index)
const;
42 bool on_the_v_diagonal(
int v_point_index)
const;
44 int corresponding_point_in_u(
int v_point_index)
const;
46 int corresponding_point_in_v(
int u_point_index)
const;
48 double distance(
int u_point_index,
int v_point_index)
const;
52 double bottleneck_alive()
const;
54 std::vector<double> sorted_distances()
const;
56 double max_dist_to_diagonal()
const;
58 Internal_point get_u_point(
int u_point_index)
const;
60 Internal_point get_v_point(
int v_point_index)
const;
63 std::vector<Internal_point> u;
64 std::vector<Internal_point> v;
68template<
typename Persistence_diagram1,
typename Persistence_diagram2>
69Persistence_graph::Persistence_graph(
const Persistence_diagram1 &diag1,
70 const Persistence_diagram2 &diag2,
double e)
71 : u(), v(), b_alive(0.) {
72 std::vector<double> u_alive;
73 std::vector<double> v_alive;
74 std::vector<double> u_nalive;
75 std::vector<double> v_nalive;
78 double inf = std::numeric_limits<double>::infinity();
81 for (
auto it = std::begin(diag1); it != std::end(diag1); ++it) {
82 if (std::get<0>(*it) != inf && std::get<1>(*it) != neginf){
83 if (std::get<0>(*it) == neginf && std::get<1>(*it) == inf)
85 else if (std::get<0>(*it) == neginf)
86 u_nalive.push_back(std::get<1>(*it));
87 else if (std::get<1>(*it) == inf)
88 u_alive.push_back(std::get<0>(*it));
89 else if (std::get<1>(*it) - std::get<0>(*it) > e)
90 u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size()));
93 for (
auto it = std::begin(diag2); it != std::end(diag2); ++it) {
94 if (std::get<0>(*it) != inf && std::get<1>(*it) != neginf){
95 if (std::get<0>(*it) == neginf && std::get<1>(*it) == inf)
97 else if (std::get<0>(*it) == neginf)
98 v_nalive.push_back(std::get<1>(*it));
99 else if (std::get<1>(*it) == inf)
100 v_alive.push_back(std::get<0>(*it));
101 else if (std::get<1>(*it) - std::get<0>(*it) > e)
102 v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size()));
105 if (u.size() < v.size())
108 if (u_alive.size() != v_alive.size() || u_nalive.size() != v_nalive.size() || u_inf != v_inf) {
109 b_alive = std::numeric_limits<double>::infinity();
111 std::sort(u_alive.begin(), u_alive.end());
112 std::sort(v_alive.begin(), v_alive.end());
113 std::sort(u_nalive.begin(), u_nalive.end());
114 std::sort(v_nalive.begin(), v_nalive.end());
115 for (
auto it_u = u_alive.cbegin(), it_v = v_alive.cbegin(); it_u != u_alive.cend(); ++it_u, ++it_v)
116 b_alive = (std::max)(b_alive, std::fabs(*it_u - *it_v));
117 for (
auto it_u = u_nalive.cbegin(), it_v = v_nalive.cbegin(); it_u != u_nalive.cend(); ++it_u, ++it_v)
118 b_alive = (std::max)(b_alive, std::fabs(*it_u - *it_v));
122inline bool Persistence_graph::on_the_u_diagonal(
int u_point_index)
const {
123 return u_point_index >=
static_cast<int> (u.size());
126inline bool Persistence_graph::on_the_v_diagonal(
int v_point_index)
const {
127 return v_point_index >=
static_cast<int> (v.size());
130inline int Persistence_graph::corresponding_point_in_u(
int v_point_index)
const {
131 return on_the_v_diagonal(v_point_index) ?
132 v_point_index -
static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size());
135inline int Persistence_graph::corresponding_point_in_v(
int u_point_index)
const {
136 return on_the_u_diagonal(u_point_index) ?
137 u_point_index -
static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size());
140inline double Persistence_graph::distance(
int u_point_index,
int v_point_index)
const {
141 if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index))
143 Internal_point p_u = get_u_point(u_point_index);
144 Internal_point p_v = get_v_point(v_point_index);
145 return (std::max)(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y()));
148inline int Persistence_graph::size()
const {
149 return static_cast<int> (u.size() + v.size());
152inline double Persistence_graph::bottleneck_alive()
const {
156inline std::vector<double> Persistence_graph::sorted_distances()
const {
157 std::vector<double> distances;
158 distances.push_back(0.);
159 for (
int u_point_index = 0; u_point_index < size(); ++u_point_index) {
160 distances.push_back(distance(u_point_index, corresponding_point_in_v(u_point_index)));
161 for (
int v_point_index = 0; v_point_index < size(); ++v_point_index)
162 distances.push_back(distance(u_point_index, v_point_index));
165 tbb::parallel_sort(distances.begin(), distances.end());
167 std::sort(distances.begin(), distances.end());
172inline Internal_point Persistence_graph::get_u_point(
int u_point_index)
const {
173 if (!on_the_u_diagonal(u_point_index))
174 return u.at(u_point_index);
175 Internal_point projector = v.at(corresponding_point_in_v(u_point_index));
176 double m = (projector.x() + projector.y()) / 2.;
177 return Internal_point(m, m, u_point_index);
180inline Internal_point Persistence_graph::get_v_point(
int v_point_index)
const {
181 if (!on_the_v_diagonal(v_point_index))
182 return v.at(v_point_index);
183 Internal_point projector = u.at(corresponding_point_in_u(v_point_index));
184 double m = (projector.x() + projector.y()) / 2.;
185 return Internal_point(m, m, v_point_index);
188inline double Persistence_graph::max_dist_to_diagonal()
const {
191 max = (std::max)(max, p.y() - p.x());
193 max = (std::max)(max, p.y() - p.x());
Gudhi namespace.
Definition SimplicialComplexForAlpha.h:14