diff options
Diffstat (limited to 'src/Bottleneck_distance/include/gudhi')
| -rw-r--r-- | src/Bottleneck_distance/include/gudhi/Bottleneck.h | 124 | ||||
| -rw-r--r-- | src/Bottleneck_distance/include/gudhi/Graph_matching.h | 162 | ||||
| -rw-r--r-- | src/Bottleneck_distance/include/gudhi/Internal_point.h | 79 | ||||
| -rw-r--r-- | src/Bottleneck_distance/include/gudhi/Neighbors_finder.h | 182 | ||||
| -rw-r--r-- | src/Bottleneck_distance/include/gudhi/Persistence_graph.h | 176 |
5 files changed, 723 insertions, 0 deletions
diff --git a/src/Bottleneck_distance/include/gudhi/Bottleneck.h b/src/Bottleneck_distance/include/gudhi/Bottleneck.h new file mode 100644 index 00000000..82ba9f68 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Bottleneck.h @@ -0,0 +1,124 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: Francois Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - 2019/06 Vincent Rouvreau : Fix doxygen warning. + * - 2019/08 Vincent Rouvreau: Fix issue #10 for CGAL + * - YYYY/MM Author: Description of the modification + */ + +#ifndef BOTTLENECK_H_ +#define BOTTLENECK_H_ + +#include <gudhi/Graph_matching.h> + +#include <CGAL/version.h> // for CGAL_VERSION_NR + +#include <vector> +#include <algorithm> // for max +#include <limits> // for numeric_limits + +#include <cmath> +#include <cfloat> // FLT_EVAL_METHOD + +// Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 +#if CGAL_VERSION_NR < 1041101000 +# error Alpha_complex_3d is only available for CGAL >= 4.11 +#endif + +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_ diff --git a/src/Bottleneck_distance/include/gudhi/Graph_matching.h b/src/Bottleneck_distance/include/gudhi/Graph_matching.h new file mode 100644 index 00000000..7b7623ce --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Graph_matching.h @@ -0,0 +1,162 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: Francois Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef GRAPH_MATCHING_H_ +#define GRAPH_MATCHING_H_ + +#include <gudhi/Neighbors_finder.h> + +#include <vector> +#include <unordered_set> +#include <algorithm> + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief Structure representing a graph matching. The graph is a Persistence_diagrams_graph. + * + * \ingroup bottleneck_distance + */ +class Graph_matching { + public: + /** \internal \brief Constructor constructing an empty matching. */ + explicit Graph_matching(Persistence_graph &g); + /** \internal \brief Is the matching perfect ? */ + bool perfect() const; + /** \internal \brief Augments the matching with a maximal set of edge-disjoint shortest augmenting paths. */ + bool multi_augment(); + /** \internal \brief Sets the maximum length of the edges allowed to be added in the matching, 0 initially. */ + void set_r(double r); + + private: + Persistence_graph* gp; + double r; + /** \internal \brief Given a point from V, provides its matched point in U, null_point_index() if there isn't. */ + std::vector<int> v_to_u; + /** \internal \brief All the unmatched points in U. */ + std::unordered_set<int> unmatched_in_u; + + /** \internal \brief Provides a Layered_neighbors_finder dividing the graph in layers. Basically a BFS. */ + Layered_neighbors_finder layering() const; + /** \internal \brief Augments the matching with a simple path no longer than max_depth. Basically a DFS. */ + bool augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth); + /** \internal \brief Update the matching with the simple augmenting path given as parameter. */ + void update(std::vector<int> & path); +}; + +inline Graph_matching::Graph_matching(Persistence_graph& g) + : gp(&g), r(0.), v_to_u(g.size(), null_point_index()), unmatched_in_u(g.size()) { + for (int u_point_index = 0; u_point_index < g.size(); ++u_point_index) + unmatched_in_u.insert(u_point_index); +} + +inline bool Graph_matching::perfect() const { + return unmatched_in_u.empty(); +} + +inline bool Graph_matching::multi_augment() { + if (perfect()) + return false; + Layered_neighbors_finder layered_nf(layering()); + int max_depth = layered_nf.vlayers_number()*2 - 1; + double rn = sqrt(gp->size()); + // verification of a necessary criterion in order to shortcut if possible + if (max_depth < 0 || (unmatched_in_u.size() > rn && max_depth >= rn)) + return false; + bool successful = false; + std::vector<int> tries(unmatched_in_u.cbegin(), unmatched_in_u.cend()); + for (auto it = tries.cbegin(); it != tries.cend(); it++) + // 'augment' has side-effects which have to be always executed, don't change order + successful = augment(layered_nf, *it, max_depth) || successful; + return successful; +} + +inline void Graph_matching::set_r(double r) { + this->r = r; +} + +inline bool Graph_matching::augment(Layered_neighbors_finder & layered_nf, int u_start_index, int max_depth) { + // V vertices have at most one successor, thus when we backtrack from U we can directly pop_back 2 vertices. + std::vector<int> path; + path.emplace_back(u_start_index); + do { + if (static_cast<int> (path.size()) > max_depth) { + path.pop_back(); + path.pop_back(); + } + if (path.empty()) + return false; + path.emplace_back(layered_nf.pull_near(path.back(), static_cast<int> (path.size()) / 2)); + while (path.back() == null_point_index()) { + path.pop_back(); + path.pop_back(); + if (path.empty()) + return false; + path.pop_back(); + path.emplace_back(layered_nf.pull_near(path.back(), path.size() / 2)); + } + path.emplace_back(v_to_u.at(path.back())); + } while (path.back() != null_point_index()); + // if v_to_u.at(path.back()) has no successor, path.back() is an exposed vertex + path.pop_back(); + update(path); + return true; +} + +inline Layered_neighbors_finder Graph_matching::layering() const { + std::vector<int> u_vertices(unmatched_in_u.cbegin(), unmatched_in_u.cend()); + std::vector<int> v_vertices; + Neighbors_finder nf(*gp, r); + for (int v_point_index = 0; v_point_index < gp->size(); ++v_point_index) + nf.add(v_point_index); + Layered_neighbors_finder layered_nf(*gp, r); + for (int layer = 0; !u_vertices.empty(); layer++) { + // one layer is one step in the BFS + for (auto it1 = u_vertices.cbegin(); it1 != u_vertices.cend(); ++it1) { + std::vector<int> u_succ(nf.pull_all_near(*it1)); + for (auto it2 = u_succ.begin(); it2 != u_succ.end(); ++it2) { + layered_nf.add(*it2, layer); + v_vertices.emplace_back(*it2); + } + } + // When the above for finishes, we have progress of one half-step (from U to V) in the BFS + u_vertices.clear(); + bool end = false; + for (auto it = v_vertices.cbegin(); it != v_vertices.cend(); it++) + if (v_to_u.at(*it) == null_point_index()) + // we stop when a nearest exposed V vertex (from U exposed vertices) has been found + end = true; + else + u_vertices.emplace_back(v_to_u.at(*it)); + // When the above for finishes, we have progress of one half-step (from V to U) in the BFS + if (end) + return layered_nf; + v_vertices.clear(); + } + return layered_nf; +} + +inline void Graph_matching::update(std::vector<int>& path) { + // Must return 1. + unmatched_in_u.erase(path.front()); + for (auto it = path.cbegin(); it != path.cend(); ++it) { + // Be careful, the iterator is incremented twice each time + int tmp = *it; + v_to_u[*(++it)] = tmp; + } +} + + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // GRAPH_MATCHING_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Internal_point.h b/src/Bottleneck_distance/include/gudhi/Internal_point.h new file mode 100644 index 00000000..f829943e --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Internal_point.h @@ -0,0 +1,79 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: Francois Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef INTERNAL_POINT_H_ +#define INTERNAL_POINT_H_ + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief Returns the used index for encoding none of the points */ +int null_point_index(); + +/** \internal \typedef \brief Internal_point is the internal points representation, indexes used outside. */ +struct Internal_point { + double vec[2]; + int point_index; + + Internal_point() { } + + Internal_point(double x, double y, int p_i) { + vec[0] = x; + vec[1] = y; + point_index = p_i; + } + + double x() const { + return vec[ 0 ]; + } + + double y() const { + return vec[ 1 ]; + } + + double& x() { + return vec[ 0 ]; + } + + double& y() { + return vec[ 1 ]; + } + + bool operator==(const Internal_point& p) const { + return point_index == p.point_index; + } + + bool operator!=(const Internal_point& p) const { + return !(*this == p); + } +}; + +inline int null_point_index() { + return -1; +} + +struct Construct_coord_iterator { + typedef const double* result_type; + + const double* operator()(const Internal_point& p) const { + return p.vec; + } + + const double* operator()(const Internal_point& p, int) const { + return p.vec + 2; + } +}; + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // INTERNAL_POINT_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h new file mode 100644 index 00000000..c65e6082 --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Neighbors_finder.h @@ -0,0 +1,182 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: Francois Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef NEIGHBORS_FINDER_H_ +#define NEIGHBORS_FINDER_H_ + +// Inclusion order is important for CGAL patch +#include <CGAL/Kd_tree.h> +#include <CGAL/Search_traits.h> + +#include <gudhi/Persistence_graph.h> +#include <gudhi/Internal_point.h> + +#include <unordered_set> +#include <vector> +#include <algorithm> // for std::max +#include <cmath> // for std::abs + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief Variant of CGAL::Fuzzy_iso_box to ensure that the box ic closed. It isn't clear how necessary that is. + */ +struct Square_query { + typedef CGAL::Dimension_tag<2> D; + typedef Internal_point Point_d; + typedef double FT; + bool contains(Point_d p) const { + return (std::max)(std::abs(p.x()-c.x()), std::abs(p.y()-c.y())) <= size; + } + bool inner_range_intersects(CGAL::Kd_tree_rectangle<FT, D> const&r) const { + return + r.max_coord(0) >= c.x() - size && + r.min_coord(0) <= c.x() + size && + r.max_coord(1) >= c.y() - size && + r.min_coord(1) <= c.y() + size; + } + bool outer_range_contains(CGAL::Kd_tree_rectangle<FT, D> const&r) const { + return + r.min_coord(0) >= c.x() - size && + r.max_coord(0) <= c.x() + size && + r.min_coord(1) >= c.y() - size && + r.max_coord(1) <= c.y() + size; + } + Point_d c; + FT size; +}; + +/** \internal \brief data structure used to find any point (including projections) in V near to a query point from U + * (which can be a projection). + * + * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically + * removed. + * + * \ingroup bottleneck_distance + */ +class Neighbors_finder { + typedef CGAL::Dimension_tag<2> D; + typedef CGAL::Search_traits<double, Internal_point, const double*, Construct_coord_iterator, D> Traits; + typedef CGAL::Kd_tree<Traits> Kd_tree; + + public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Neighbors_finder(const Persistence_graph& g, double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index); + /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if + * there isn't such a point. */ + int pull_near(int u_point_index); + /** \internal \brief Returns and remove all the V points near to the U point given as parameter. */ + std::vector<int> pull_all_near(int u_point_index); + + private: + const Persistence_graph& g; + const double r; + Kd_tree kd_t; + std::unordered_set<int> projections_f; +}; + +/** \internal \brief data structure used to find any point (including projections) in V near to a query point from U + * (which can be a projection) in a layered graph layer given as parmeter. + * + * V points have to be added manually using their index and before the first pull. A neighbor pulled is automatically + * removed. + * + * \ingroup bottleneck_distance + */ +class Layered_neighbors_finder { + public: + /** \internal \brief Constructor taking the near distance definition as parameter. */ + Layered_neighbors_finder(const Persistence_graph& g, double r); + /** \internal \brief A point added will be possibly pulled. */ + void add(int v_point_index, int vlayer); + /** \internal \brief Returns and remove a V point near to the U point given as parameter, null_point_index() if + * there isn't such a point. */ + int pull_near(int u_point_index, int vlayer); + /** \internal \brief Returns the number of layers. */ + int vlayers_number() const; + + private: + const Persistence_graph& g; + const double r; + std::vector<std::unique_ptr<Neighbors_finder>> neighbors_finder; +}; + +inline Neighbors_finder::Neighbors_finder(const Persistence_graph& g, double r) : + g(g), r(r), kd_t(), projections_f() { } + +inline void Neighbors_finder::add(int v_point_index) { + if (g.on_the_v_diagonal(v_point_index)) + projections_f.emplace(v_point_index); + else + kd_t.insert(g.get_v_point(v_point_index)); +} + +inline int Neighbors_finder::pull_near(int u_point_index) { + int tmp; + int c = g.corresponding_point_in_v(u_point_index); + if (g.on_the_u_diagonal(u_point_index) && !projections_f.empty()) { + // Any pair of projection is at distance 0 + tmp = *projections_f.cbegin(); + projections_f.erase(tmp); + } else if (projections_f.count(c) && (g.distance(u_point_index, c) <= r)) { + // Is the query point near to its projection ? + tmp = c; + projections_f.erase(tmp); + } else { + // Is the query point near to a V point in the plane ? + Internal_point u_point = g.get_u_point(u_point_index); + auto neighbor = kd_t.search_any_point(Square_query{u_point, r}); + if (!neighbor) + return null_point_index(); + tmp = neighbor->point_index; + auto point = g.get_v_point(tmp); + int idx = point.point_index; + kd_t.remove(point, [idx](Internal_point const&p){return p.point_index == idx;}); + } + return tmp; +} + +inline std::vector<int> Neighbors_finder::pull_all_near(int u_point_index) { + std::vector<int> all_pull; + int last_pull = pull_near(u_point_index); + while (last_pull != null_point_index()) { + all_pull.emplace_back(last_pull); + last_pull = pull_near(u_point_index); + } + return all_pull; +} + +inline Layered_neighbors_finder::Layered_neighbors_finder(const Persistence_graph& g, double r) : + g(g), r(r), neighbors_finder() { } + +inline void Layered_neighbors_finder::add(int v_point_index, int vlayer) { + for (int l = neighbors_finder.size(); l <= vlayer; l++) + neighbors_finder.emplace_back(std::unique_ptr<Neighbors_finder>(new Neighbors_finder(g, r))); + neighbors_finder.at(vlayer)->add(v_point_index); +} + +inline int Layered_neighbors_finder::pull_near(int u_point_index, int vlayer) { + if (static_cast<int> (neighbors_finder.size()) <= vlayer) + return null_point_index(); + return neighbors_finder.at(vlayer)->pull_near(u_point_index); +} + +inline int Layered_neighbors_finder::vlayers_number() const { + return static_cast<int> (neighbors_finder.size()); +} + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // NEIGHBORS_FINDER_H_ diff --git a/src/Bottleneck_distance/include/gudhi/Persistence_graph.h b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h new file mode 100644 index 00000000..f791e37c --- /dev/null +++ b/src/Bottleneck_distance/include/gudhi/Persistence_graph.h @@ -0,0 +1,176 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author: Francois Godi + * + * Copyright (C) 2015 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef PERSISTENCE_GRAPH_H_ +#define PERSISTENCE_GRAPH_H_ + +#include <gudhi/Internal_point.h> + +#ifdef GUDHI_USE_TBB +#include <tbb/parallel_sort.h> +#endif + +#include <vector> +#include <algorithm> +#include <limits> // for numeric_limits + +namespace Gudhi { + +namespace persistence_diagram { + +/** \internal \brief Structure representing an euclidean bipartite graph containing + * the points from the two persistence diagrams (including the projections). + * + * \ingroup bottleneck_distance + */ +class Persistence_graph { + public: + /** \internal \brief Constructor taking 2 PersistenceDiagrams (concept) as parameters. */ + template<typename Persistence_diagram1, typename Persistence_diagram2> + Persistence_graph(const Persistence_diagram1& diag1, const Persistence_diagram2& diag2, double e); + /** \internal \brief Is the given point from U the projection of a point in V ? */ + bool on_the_u_diagonal(int u_point_index) const; + /** \internal \brief Is the given point from V the projection of a point in U ? */ + bool on_the_v_diagonal(int v_point_index) const; + /** \internal \brief Given a point from V, returns the corresponding (projection or projector) point in U. */ + int corresponding_point_in_u(int v_point_index) const; + /** \internal \brief Given a point from U, returns the corresponding (projection or projector) point in V. */ + int corresponding_point_in_v(int u_point_index) const; + /** \internal \brief Given a point from U and a point from V, returns the distance between those points. */ + double distance(int u_point_index, int v_point_index) const; + /** \internal \brief Returns size = |U| = |V|. */ + int size() const; + /** \internal \brief Is there as many infinite points (alive components) in both diagrams ? */ + double bottleneck_alive() const; + /** \internal \brief Returns the O(n^2) sorted distances between the points. */ + std::vector<double> sorted_distances() const; + /** \internal \brief Returns an upper bound for the diameter of the convex hull of all non infinite points */ + double diameter_bound() const; + /** \internal \brief Returns the corresponding internal point */ + Internal_point get_u_point(int u_point_index) const; + /** \internal \brief Returns the corresponding internal point */ + Internal_point get_v_point(int v_point_index) const; + + private: + std::vector<Internal_point> u; + std::vector<Internal_point> v; + double b_alive; +}; + +template<typename Persistence_diagram1, typename Persistence_diagram2> +Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1, + const Persistence_diagram2 &diag2, double e) + : u(), v(), b_alive(0.) { + std::vector<double> u_alive; + std::vector<double> v_alive; + for (auto it = std::begin(diag1); it != std::end(diag1); ++it) { + if (std::get<1>(*it) == std::numeric_limits<double>::infinity()) + u_alive.push_back(std::get<0>(*it)); + else if (std::get<1>(*it) - std::get<0>(*it) > e) + u.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), u.size())); + } + for (auto it = std::begin(diag2); it != std::end(diag2); ++it) { + if (std::get<1>(*it) == std::numeric_limits<double>::infinity()) + v_alive.push_back(std::get<0>(*it)); + else if (std::get<1>(*it) - std::get<0>(*it) > e) + v.push_back(Internal_point(std::get<0>(*it), std::get<1>(*it), v.size())); + } + if (u.size() < v.size()) + swap(u, v); + std::sort(u_alive.begin(), u_alive.end()); + std::sort(v_alive.begin(), v_alive.end()); + if (u_alive.size() != v_alive.size()) { + b_alive = std::numeric_limits<double>::infinity(); + } else { + for (auto it_u = u_alive.cbegin(), it_v = v_alive.cbegin(); it_u != u_alive.cend(); ++it_u, ++it_v) + b_alive = (std::max)(b_alive, std::fabs(*it_u - *it_v)); + } +} + +inline bool Persistence_graph::on_the_u_diagonal(int u_point_index) const { + return u_point_index >= static_cast<int> (u.size()); +} + +inline bool Persistence_graph::on_the_v_diagonal(int v_point_index) const { + return v_point_index >= static_cast<int> (v.size()); +} + +inline int Persistence_graph::corresponding_point_in_u(int v_point_index) const { + return on_the_v_diagonal(v_point_index) ? + v_point_index - static_cast<int> (v.size()) : v_point_index + static_cast<int> (u.size()); +} + +inline int Persistence_graph::corresponding_point_in_v(int u_point_index) const { + return on_the_u_diagonal(u_point_index) ? + u_point_index - static_cast<int> (u.size()) : u_point_index + static_cast<int> (v.size()); +} + +inline double Persistence_graph::distance(int u_point_index, int v_point_index) const { + if (on_the_u_diagonal(u_point_index) && on_the_v_diagonal(v_point_index)) + return 0.; + Internal_point p_u = get_u_point(u_point_index); + Internal_point p_v = get_v_point(v_point_index); + return (std::max)(std::fabs(p_u.x() - p_v.x()), std::fabs(p_u.y() - p_v.y())); +} + +inline int Persistence_graph::size() const { + return static_cast<int> (u.size() + v.size()); +} + +inline double Persistence_graph::bottleneck_alive() const { + return b_alive; +} + +inline std::vector<double> Persistence_graph::sorted_distances() const { + std::vector<double> distances; + distances.push_back(0.); // for empty diagrams + for (int u_point_index = 0; u_point_index < size(); ++u_point_index) { + distances.push_back(distance(u_point_index, corresponding_point_in_v(u_point_index))); + for (int v_point_index = 0; v_point_index < size(); ++v_point_index) + distances.push_back(distance(u_point_index, v_point_index)); + } +#ifdef GUDHI_USE_TBB + tbb::parallel_sort(distances.begin(), distances.end()); +#else + std::sort(distances.begin(), distances.end()); +#endif + return distances; +} + +inline Internal_point Persistence_graph::get_u_point(int u_point_index) const { + if (!on_the_u_diagonal(u_point_index)) + return u.at(u_point_index); + Internal_point projector = v.at(corresponding_point_in_v(u_point_index)); + double m = (projector.x() + projector.y()) / 2.; + return Internal_point(m, m, u_point_index); +} + +inline Internal_point Persistence_graph::get_v_point(int v_point_index) const { + if (!on_the_v_diagonal(v_point_index)) + return v.at(v_point_index); + Internal_point projector = u.at(corresponding_point_in_u(v_point_index)); + double m = (projector.x() + projector.y()) / 2.; + return Internal_point(m, m, v_point_index); +} + +inline double Persistence_graph::diameter_bound() const { + double max = 0.; + for (auto it = u.cbegin(); it != u.cend(); it++) + max = (std::max)(max, it->y()); + for (auto it = v.cbegin(); it != v.cend(); it++) + max = (std::max)(max, it->y()); + return max; +} + +} // namespace persistence_diagram + +} // namespace Gudhi + +#endif // PERSISTENCE_GRAPH_H_ |