Persistence_graph.h
1 /* This file is part of the Gudhi Library. The Gudhi library
2  * (Geometric Understanding in Higher Dimensions) is a generic C++
3  * library for computational topology.
4  *
5  * Author: Francois Godi
6  *
7  * Copyright (C) 2015 Inria
8  *
9  * This program is free software: you can redistribute it and/or modify
10  * it under the terms of the GNU General Public License as published by
11  * the Free Software Foundation, either version 3 of the License, or
12  * (at your option) any later version.
13  *
14  * This program is distributed in the hope that it will be useful,
15  * but WITHOUT ANY WARRANTY; without even the implied warranty of
16  * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
17  * GNU General Public License for more details.
18  *
19  * You should have received a copy of the GNU General Public License
20  * along with this program. If not, see <http://www.gnu.org/licenses/>.
21  */
22 
23 #ifndef PERSISTENCE_GRAPH_H_
24 #define PERSISTENCE_GRAPH_H_
25 
26 #include <gudhi/Internal_point.h>
27 
28 #ifdef GUDHI_USE_TBB
29 #include <tbb/parallel_sort.h>
30 #endif
31 
32 #include <vector>
33 #include <algorithm>
34 #include <limits> // for numeric_limits
35 
36 namespace Gudhi {
37 
38 namespace persistence_diagram {
39 
45 class Persistence_graph {
46  public:
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;
61  int size() 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;
72 
73  private:
74  std::vector<Internal_point> u;
75  std::vector<Internal_point> v;
76  double b_alive;
77 };
78 
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()));
90  }
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()));
96  }
97  if (u.size() < v.size())
98  swap(u, v);
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();
103  } else {
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));
106  }
107 }
108 
109 inline bool Persistence_graph::on_the_u_diagonal(int u_point_index) const {
110  return u_point_index >= static_cast<int> (u.size());
111 }
112 
113 inline bool Persistence_graph::on_the_v_diagonal(int v_point_index) const {
114  return v_point_index >= static_cast<int> (v.size());
115 }
116 
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());
120 }
121 
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());
125 }
126 
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))
129  return 0.;
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()));
133 }
134 
135 inline int Persistence_graph::size() const {
136  return static_cast<int> (u.size() + v.size());
137 }
138 
139 inline double Persistence_graph::bottleneck_alive() const {
140  return b_alive;
141 }
142 
143 inline std::vector<double> Persistence_graph::sorted_distances() const {
144  std::vector<double> distances;
145  distances.push_back(0.); // for empty diagrams
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));
150  }
151 #ifdef GUDHI_USE_TBB
152  tbb::parallel_sort(distances.begin(), distances.end());
153 #else
154  std::sort(distances.begin(), distances.end());
155 #endif
156  return distances;
157 }
158 
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);
165 }
166 
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);
173 }
174 
175 inline double Persistence_graph::diameter_bound() const {
176  double max = 0.;
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());
181  return max;
182 }
183 
184 } // namespace persistence_diagram
185 
186 } // namespace Gudhi
187 
188 #endif // PERSISTENCE_GRAPH_H_
Definition: SimplicialComplexForAlpha.h:26
GUDHI  Version 2.2.0  - C++ library for Topological Data Analysis (TDA) and Higher Dimensional Geometry Understanding.  - Copyright : GPL v3 Generated on Fri Jun 22 2018 08:12:19 for GUDHI by Doxygen 1.8.13