summaryrefslogtreecommitdiff
path: root/include
diff options
context:
space:
mode:
Diffstat (limited to 'include')
-rw-r--r--include/gudhi/Bitmap_cubical_complex.h27
-rw-r--r--include/gudhi/Bottleneck.h4
-rw-r--r--include/gudhi/Clock.h20
-rw-r--r--include/gudhi/Euclidean_strong_witness_complex.h2
-rw-r--r--include/gudhi/Euclidean_witness_complex.h2
-rw-r--r--include/gudhi/Hasse_complex.h11
-rw-r--r--include/gudhi/Kd_tree_search.h81
-rw-r--r--include/gudhi/Neighbors_finder.h8
-rw-r--r--include/gudhi/Persistence_graph.h8
-rw-r--r--include/gudhi/Persistent_cohomology.h13
-rw-r--r--include/gudhi/Simplex_tree.h26
-rw-r--r--include/gudhi/Strong_witness_complex.h4
-rw-r--r--include/gudhi/Tangential_complex.h21
-rw-r--r--include/gudhi/Witness_complex.h4
-rw-r--r--include/gudhi/pick_n_random_points.h4
-rw-r--r--include/gudhi/random_point_generators.h32
-rw-r--r--include/gudhi/reader_utils.h145
-rw-r--r--include/gudhi/sparsify_point_set.h2
18 files changed, 253 insertions, 161 deletions
diff --git a/include/gudhi/Bitmap_cubical_complex.h b/include/gudhi/Bitmap_cubical_complex.h
index 5a87b9b8..f395de65 100644
--- a/include/gudhi/Bitmap_cubical_complex.h
+++ b/include/gudhi/Bitmap_cubical_complex.h
@@ -97,7 +97,7 @@ class Bitmap_cubical_complex : public T {
* with filtration on top dimensional cells.
**/
Bitmap_cubical_complex(const std::vector<unsigned>& dimensions,
- const std::vector<typename T::filtration_type>& top_dimensional_cells) :
+ const std::vector<Filtration_value>& top_dimensional_cells) :
T(dimensions, top_dimensional_cells),
key_associated_to_simplex(this->total_number_of_cells + 1) {
for (size_t i = 0; i != this->total_number_of_cells; ++i) {
@@ -111,13 +111,13 @@ class Bitmap_cubical_complex : public T {
/**
* Constructor that requires vector of elements of type unsigned, which gives number of top dimensional cells
- * in the following directions and vector of element of a type T::filtration_type
+ * in the following directions and vector of element of a type Filtration_value
* with filtration on top dimensional cells. The last parameter of the constructor is a vector of boolean of a length
* equal to the dimension of cubical complex.
* If the position i on this vector is true, then we impose periodic boundary conditions in this direction.
**/
Bitmap_cubical_complex(const std::vector<unsigned>& dimensions,
- const std::vector<typename T::filtration_type>& top_dimensional_cells,
+ const std::vector<Filtration_value>& top_dimensional_cells,
std::vector< bool > directions_in_which_periodic_b_cond_are_to_be_imposed) :
T(dimensions, top_dimensional_cells, directions_in_which_periodic_b_cond_are_to_be_imposed),
key_associated_to_simplex(this->total_number_of_cells + 1) {
@@ -170,20 +170,20 @@ class Bitmap_cubical_complex : public T {
if (globalDbg) {
std::cerr << "unsigned dimension(const Simplex_handle& sh)\n";
}
- if (sh != std::numeric_limits<Simplex_handle>::max()) return this->get_dimension_of_a_cell(sh);
+ if (sh != null_simplex()) return this->get_dimension_of_a_cell(sh);
return -1;
}
/**
* Return the filtration of a cell pointed by the Simplex_handle.
**/
- typename T::filtration_type filtration(Simplex_handle sh) {
+ Filtration_value filtration(Simplex_handle sh) {
if (globalDbg) {
- std::cerr << "T::filtration_type filtration(const Simplex_handle& sh)\n";
+ std::cerr << "Filtration_value filtration(const Simplex_handle& sh)\n";
}
// Returns the filtration value of a simplex.
- if (sh != std::numeric_limits<Simplex_handle>::max()) return this->data[sh];
- return std::numeric_limits<Simplex_handle>::max();
+ if (sh != null_simplex()) return this->data[sh];
+ return std::numeric_limits<Filtration_value>::infinity();
}
/**
@@ -203,7 +203,7 @@ class Bitmap_cubical_complex : public T {
if (globalDbg) {
std::cerr << "Simplex_key key(const Simplex_handle& sh)\n";
}
- if (sh != std::numeric_limits<Simplex_handle>::max()) {
+ if (sh != null_simplex()) {
return this->key_associated_to_simplex[sh];
}
return this->null_key();
@@ -216,7 +216,7 @@ class Bitmap_cubical_complex : public T {
if (globalDbg) {
std::cerr << "Simplex_handle simplex(Simplex_key key)\n";
}
- if (key != std::numeric_limits<Simplex_handle>::max()) {
+ if (key != null_key()) {
return this->simplex_associated_to_key[ key ];
}
return null_simplex();
@@ -229,7 +229,7 @@ class Bitmap_cubical_complex : public T {
if (globalDbg) {
std::cerr << "void assign_key(Simplex_handle& sh, Simplex_key key)\n";
}
- if (key == std::numeric_limits<Simplex_handle>::max()) return;
+ if (key == null_key()) return;
this->key_associated_to_simplex[sh] = key;
this->simplex_associated_to_key[key] = sh;
}
@@ -566,8 +566,9 @@ class is_before_in_filtration {
bool operator()(const typename Bitmap_cubical_complex<T>::Simplex_handle& sh1,
const typename Bitmap_cubical_complex<T>::Simplex_handle& sh2) const {
// Not using st_->filtration(sh1) because it uselessly tests for null_simplex.
- typename T::filtration_type fil1 = CC_->data[sh1];
- typename T::filtration_type fil2 = CC_->data[sh2];
+ typedef typename T::filtration_type Filtration_value;
+ Filtration_value fil1 = CC_->data[sh1];
+ Filtration_value fil2 = CC_->data[sh2];
if (fil1 != fil2) {
return fil1 < fil2;
}
diff --git a/include/gudhi/Bottleneck.h b/include/gudhi/Bottleneck.h
index b90a0ee0..8c97dce9 100644
--- a/include/gudhi/Bottleneck.h
+++ b/include/gudhi/Bottleneck.h
@@ -101,11 +101,11 @@ double bottleneck_distance_exact(Persistence_graph& g) {
*/
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()) {
+ 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));
+ return (std::max)(g.bottleneck_alive(), e == 0. ? bottleneck_distance_exact(g) : bottleneck_distance_approx(g, e));
}
} // namespace persistence_diagram
diff --git a/include/gudhi/Clock.h b/include/gudhi/Clock.h
index 77f196ca..b83de2f5 100644
--- a/include/gudhi/Clock.h
+++ b/include/gudhi/Clock.h
@@ -23,9 +23,9 @@
#ifndef CLOCK_H_
#define CLOCK_H_
-#include <boost/date_time/posix_time/posix_time.hpp>
-
+#include <iostream>
#include <string>
+#include <chrono>
namespace Gudhi {
@@ -33,20 +33,20 @@ class Clock {
public:
// Construct and start the timer
Clock(const std::string& msg_ = std::string())
- : startTime(boost::posix_time::microsec_clock::local_time()),
+ : startTime(std::chrono::system_clock::now()),
end_called(false),
msg(msg_) { }
// Restart the timer
void begin() const {
end_called = false;
- startTime = boost::posix_time::microsec_clock::local_time();
+ startTime = std::chrono::system_clock::now();
}
// Stop the timer
void end() const {
end_called = true;
- endTime = boost::posix_time::microsec_clock::local_time();
+ endTime = std::chrono::system_clock::now();
}
std::string message() const {
@@ -62,7 +62,7 @@ class Clock {
if (!clock.msg.empty())
stream << clock.msg << ": ";
- stream << clock.num_seconds() << "s";
+ stream << clock.num_seconds() << "s\n";
return stream;
}
@@ -71,15 +71,15 @@ class Clock {
// - or now otherwise. In this case, the timer is not stopped.
double num_seconds() const {
if (!end_called) {
- auto end = boost::posix_time::microsec_clock::local_time();
- return (end - startTime).total_milliseconds() / 1000.;
+ auto end = std::chrono::system_clock::now();
+ return std::chrono::duration_cast<std::chrono::milliseconds>(end-startTime).count() / 1000.;
} else {
- return (endTime - startTime).total_milliseconds() / 1000.;
+ return std::chrono::duration_cast<std::chrono::milliseconds>(endTime-startTime).count() / 1000.;
}
}
private:
- mutable boost::posix_time::ptime startTime, endTime;
+ mutable std::chrono::time_point<std::chrono::system_clock> startTime, endTime;
mutable bool end_called;
std::string msg;
};
diff --git a/include/gudhi/Euclidean_strong_witness_complex.h b/include/gudhi/Euclidean_strong_witness_complex.h
index fb669ef8..4f3cef4f 100644
--- a/include/gudhi/Euclidean_strong_witness_complex.h
+++ b/include/gudhi/Euclidean_strong_witness_complex.h
@@ -84,7 +84,7 @@ class Euclidean_strong_witness_complex
: landmarks_(std::begin(landmarks), std::end(landmarks)), landmark_tree_(landmarks_) {
nearest_landmark_table_.reserve(boost::size(witnesses));
for (auto w : witnesses)
- nearest_landmark_table_.push_back(landmark_tree_.query_incremental_nearest_neighbors(w));
+ nearest_landmark_table_.push_back(landmark_tree_.incremental_nearest_neighbors(w));
}
/** \brief Returns the point corresponding to the given vertex.
diff --git a/include/gudhi/Euclidean_witness_complex.h b/include/gudhi/Euclidean_witness_complex.h
index 6afe9a5d..ff8bb139 100644
--- a/include/gudhi/Euclidean_witness_complex.h
+++ b/include/gudhi/Euclidean_witness_complex.h
@@ -86,7 +86,7 @@ class Euclidean_witness_complex
: landmarks_(std::begin(landmarks), std::end(landmarks)), landmark_tree_(landmarks) {
nearest_landmark_table_.reserve(boost::size(witnesses));
for (auto w : witnesses)
- nearest_landmark_table_.push_back(landmark_tree_.query_incremental_nearest_neighbors(w));
+ nearest_landmark_table_.push_back(landmark_tree_.incremental_nearest_neighbors(w));
}
/** \brief Returns the point corresponding to the given vertex.
diff --git a/include/gudhi/Hasse_complex.h b/include/gudhi/Hasse_complex.h
index 8b06b771..e67f7609 100644
--- a/include/gudhi/Hasse_complex.h
+++ b/include/gudhi/Hasse_complex.h
@@ -30,6 +30,7 @@
#include <algorithm>
#include <utility> // for pair
#include <vector>
+#include <limits> // for infinity value
#ifdef GUDHI_USE_TBB
#include <tbb/parallel_for.h>
@@ -104,7 +105,6 @@ class Hasse_complex {
Hasse_complex(Complex_ds & cpx)
: complex_(cpx.num_simplices())
, vertices_()
- , threshold_(cpx.filtration())
, num_vertices_()
, dim_max_(cpx.dimension()) {
int size = complex_.size();
@@ -125,7 +125,6 @@ class Hasse_complex {
Hasse_complex()
: complex_()
, vertices_()
- , threshold_(0)
, num_vertices_(0)
, dim_max_(-1) { }
@@ -157,15 +156,11 @@ class Hasse_complex {
Filtration_value filtration(Simplex_handle sh) {
if (sh == null_simplex()) {
- return filtration();
+ return std::numeric_limits<Filtration_value>::infinity();
}
return complex_[sh].filtration_;
}
- Filtration_value filtration() {
- return threshold_;
- }
-
int dimension(Simplex_handle sh) {
if (complex_[sh].boundary_.empty()) return 0;
return complex_[sh].boundary_.size() - 1;
@@ -206,7 +201,6 @@ class Hasse_complex {
std::vector< Hasse_simp, Gudhi::no_init_allocator<Hasse_simp> > complex_;
std::vector<Simplex_handle> vertices_;
- Filtration_value threshold_;
size_t num_vertices_;
int dim_max_;
};
@@ -245,7 +239,6 @@ std::istream& operator>>(std::istream & is
}
hcpx.dim_max_ = max_dim;
- hcpx.threshold_ = max_fil;
return is;
}
diff --git a/include/gudhi/Kd_tree_search.h b/include/gudhi/Kd_tree_search.h
index 6728d56e..ef428002 100644
--- a/include/gudhi/Kd_tree_search.h
+++ b/include/gudhi/Kd_tree_search.h
@@ -27,6 +27,7 @@
#include <CGAL/Orthogonal_incremental_neighbor_search.h>
#include <CGAL/Search_traits.h>
#include <CGAL/Search_traits_adapter.h>
+#include <CGAL/Fuzzy_sphere.h>
#include <CGAL/property_map.h>
#include <boost/property_map/property_map.hpp>
@@ -41,19 +42,19 @@ namespace spatial_searching {
/**
* \class Kd_tree_search Kd_tree_search.h gudhi/Kd_tree_search.h
- * \brief Spatial tree data structure to perform (approximate) nearest and farthest neighbor search.
+ * \brief Spatial tree data structure to perform (approximate) nearest and furthest neighbor search.
*
* \ingroup spatial_searching
*
* \details
* The class Kd_tree_search is a tree-based data structure, based on
* <a target="_blank" href="http://doc.cgal.org/latest/Spatial_searching/index.html">CGAL dD spatial searching data structures</a>.
- * It provides a simplified API to perform (approximate) nearest and farthest neighbor searches. Contrary to CGAL default behavior, the tree
+ * It provides a simplified API to perform (approximate) nearest and furthest neighbor searches. Contrary to CGAL default behavior, the tree
* does not store the points themselves, but stores indices.
*
- * There are two types of queries: the <i>k-nearest or k-farthest neighbor query</i>, where <i>k</i> is fixed and the <i>k</i> nearest
- * or farthest points are computed right away,
- * and the <i>incremental nearest or farthest neighbor query</i>, where no number of neighbors is provided during the call, as the
+ * There are two types of queries: the <i>k-nearest or k-furthest neighbor query</i>, where <i>k</i> is fixed and the <i>k</i> nearest
+ * or furthest points are computed right away,
+ * and the <i>incremental nearest or furthest neighbor query</i>, where no number of neighbors is provided during the call, as the
* neighbors will be computed incrementally when the iterator on the range is incremented.
*
* \tparam Search_traits must be a model of the <a target="_blank"
@@ -87,11 +88,15 @@ class Kd_tree_search {
std::ptrdiff_t,
Point_property_map,
Traits_base> STraits;
+ typedef CGAL::Distance_adapter<
+ std::ptrdiff_t,
+ Point_property_map,
+ CGAL::Euclidean_distance<Traits_base> > Orthogonal_distance;
typedef CGAL::Orthogonal_k_neighbor_search<STraits> K_neighbor_search;
typedef typename K_neighbor_search::Tree Tree;
typedef typename K_neighbor_search::Distance Distance;
- /// \brief The range returned by a k-nearest or k-farthest neighbor search.
+ /// \brief The range returned by a k-nearest or k-furthest neighbor search.
/// Its value type is `std::pair<std::size_t, FT>` where `first` is the index
/// of a point P and `second` is the squared distance between P and the query point.
typedef K_neighbor_search KNS_range;
@@ -99,11 +104,12 @@ class Kd_tree_search {
typedef CGAL::Orthogonal_incremental_neighbor_search<
STraits, Distance, CGAL::Sliding_midpoint<STraits>, Tree>
Incremental_neighbor_search;
- /// \brief The range returned by an incremental nearest or farthest neighbor search.
+ /// \brief The range returned by an incremental nearest or furthest neighbor search.
/// Its value type is `std::pair<std::size_t, FT>` where `first` is the index
/// of a point P and `second` is the squared distance between P and the query point.
typedef Incremental_neighbor_search INS_range;
+ typedef CGAL::Fuzzy_sphere<STraits> Fuzzy_sphere;
/// \brief Constructor
/// @param[in] points Const reference to the point range. This range
/// is not copied, so it should not be destroyed or modified afterwards.
@@ -164,9 +170,9 @@ class Kd_tree_search {
/// @param[in] k Number of nearest points to search.
/// @param[in] sorted Indicates if the computed sequence of k-nearest neighbors needs to be sorted.
/// @param[in] eps Approximation factor.
- /// @return A range containing the k-nearest neighbors.
- KNS_range query_k_nearest_neighbors(const
- Point &p,
+ /// @return A range (whose `value_type` is `std::size_t`) containing the k-nearest neighbors.
+ KNS_range k_nearest_neighbors(
+ Point const& p,
unsigned int k,
bool sorted = true,
FT eps = FT(0)) const {
@@ -179,8 +185,7 @@ class Kd_tree_search {
k,
eps,
true,
- CGAL::Distance_adapter<std::ptrdiff_t, Point_property_map, CGAL::Euclidean_distance<Traits_base> >(
- std::begin(m_points)), sorted);
+ Orthogonal_distance(std::begin(m_points)), sorted);
return search;
}
@@ -188,10 +193,11 @@ class Kd_tree_search {
/// \brief Search incrementally for the nearest neighbors from a query point.
/// @param[in] p The query point.
/// @param[in] eps Approximation factor.
- /// @return A range containing the neighbors sorted by their distance to p.
+ /// @return A range (whose `value_type` is `std::size_t`) containing the
+ /// neighbors sorted by their distance to p.
/// All the neighbors are not computed by this function, but they will be
/// computed incrementally when the iterator on the range is incremented.
- INS_range query_incremental_nearest_neighbors(const Point &p, FT eps = FT(0)) const {
+ INS_range incremental_nearest_neighbors(Point const& p, FT eps = FT(0)) const {
// Initialize the search structure, and search all N points
// Note that we need to pass the Distance explicitly since it needs to
// know the property map
@@ -200,20 +206,19 @@ class Kd_tree_search {
p,
eps,
true,
- CGAL::Distance_adapter<std::ptrdiff_t, Point_property_map, CGAL::Euclidean_distance<Traits_base> >(
- std::begin(m_points)) );
+ Orthogonal_distance(std::begin(m_points)) );
return search;
}
- /// \brief Search for the k-farthest points from a query point.
+ /// \brief Search for the k-furthest points from a query point.
/// @param[in] p The query point.
- /// @param[in] k Number of farthest points to search.
- /// @param[in] sorted Indicates if the computed sequence of k-farthest neighbors needs to be sorted.
+ /// @param[in] k Number of furthest points to search.
+ /// @param[in] sorted Indicates if the computed sequence of k-furthest neighbors needs to be sorted.
/// @param[in] eps Approximation factor.
- /// @return A range containing the k-farthest neighbors.
- KNS_range query_k_farthest_neighbors(const
- Point &p,
+ /// @return A range (whose `value_type` is `std::size_t`) containing the k-furthest neighbors.
+ KNS_range k_furthest_neighbors(
+ Point const& p,
unsigned int k,
bool sorted = true,
FT eps = FT(0)) const {
@@ -226,19 +231,19 @@ class Kd_tree_search {
k,
eps,
false,
- CGAL::Distance_adapter<std::ptrdiff_t, Point_property_map, CGAL::Euclidean_distance<Traits_base> >(
- std::begin(m_points)), sorted);
+ Orthogonal_distance(std::begin(m_points)), sorted);
return search;
}
- /// \brief Search incrementally for the farthest neighbors from a query point.
+ /// \brief Search incrementally for the furthest neighbors from a query point.
/// @param[in] p The query point.
/// @param[in] eps Approximation factor.
- /// @return A range containing the neighbors sorted by their distance to p.
+ /// @return A range (whose `value_type` is `std::size_t`)
+ /// containing the neighbors sorted by their distance to p.
/// All the neighbors are not computed by this function, but they will be
/// computed incrementally when the iterator on the range is incremented.
- INS_range query_incremental_farthest_neighbors(const Point &p, FT eps = FT(0)) const {
+ INS_range incremental_furthest_neighbors(Point const& p, FT eps = FT(0)) const {
// Initialize the search structure, and search all N points
// Note that we need to pass the Distance explicitly since it needs to
// know the property map
@@ -247,12 +252,30 @@ class Kd_tree_search {
p,
eps,
false,
- CGAL::Distance_adapter<std::ptrdiff_t, Point_property_map, CGAL::Euclidean_distance<Traits_base> >(
- std::begin(m_points)) );
+ Orthogonal_distance(std::begin(m_points)) );
return search;
}
+ /// \brief Search for all the neighbors in a ball.
+ /// @param[in] p The query point.
+ /// @param[in] radius The search radius
+ /// @param[out] it The points that lie inside the sphere of center `p` and radius `radius`.
+ /// Note: `it` is used this way: `*it++ = each_point`.
+ /// @param[in] eps Approximation factor.
+ template <typename OutputIterator>
+ void all_near_neighbors(Point const& p,
+ FT radius,
+ OutputIterator it,
+ FT eps = FT(0)) const {
+ m_tree.search(it, Fuzzy_sphere(p, radius, eps, m_tree.traits()));
+ }
+
+ int tree_depth() const
+ {
+ return m_tree.root()->depth();
+ }
+
private:
Point_range const& m_points;
Tree m_tree;
diff --git a/include/gudhi/Neighbors_finder.h b/include/gudhi/Neighbors_finder.h
index bdc47578..a6b9b021 100644
--- a/include/gudhi/Neighbors_finder.h
+++ b/include/gudhi/Neighbors_finder.h
@@ -44,16 +44,16 @@ struct Square_query {
typedef Internal_point Point_d;
typedef double FT;
bool contains(Point_d p) const {
- return std::abs(p.x()-c.x())<=size && std::abs(p.y()-c.y())<=size;
+ return std::abs(p.x()-c.x()) <= size && std::abs(p.y()-c.y()) <= size;
}
- bool inner_range_intersects(CGAL::Kd_tree_rectangle<FT,D> const&r) const {
+ 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 {
+ 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 &&
@@ -146,7 +146,7 @@ inline int Neighbors_finder::pull_near(int u_point_index) {
// 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)
+ if (!neighbor)
return null_point_index();
tmp = neighbor->point_index;
auto point = g.get_v_point(tmp);
diff --git a/include/gudhi/Persistence_graph.h b/include/gudhi/Persistence_graph.h
index 44f4b827..622b0691 100644
--- a/include/gudhi/Persistence_graph.h
+++ b/include/gudhi/Persistence_graph.h
@@ -102,7 +102,7 @@ Persistence_graph::Persistence_graph(const Persistence_diagram1 &diag1,
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));
+ b_alive = (std::max)(b_alive, std::fabs(*it_u - *it_v));
}
}
@@ -129,7 +129,7 @@ inline double Persistence_graph::distance(int u_point_index, int 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()));
+ 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 {
@@ -175,9 +175,9 @@ inline Internal_point Persistence_graph::get_v_point(int v_point_index) const {
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());
+ max = (std::max)(max, it->y());
for (auto it = v.cbegin(); it != v.cend(); it++)
- max = std::max(max, it->y());
+ max = (std::max)(max, it->y());
return max;
}
diff --git a/include/gudhi/Persistent_cohomology.h b/include/gudhi/Persistent_cohomology.h
index 672fda48..e0a147b3 100644
--- a/include/gudhi/Persistent_cohomology.h
+++ b/include/gudhi/Persistent_cohomology.h
@@ -591,10 +591,17 @@ class Persistent_cohomology {
std::ofstream diagram_out(diagram_name.c_str());
cmp_intervals_by_length cmp(cpx_);
std::sort(std::begin(persistent_pairs_), std::end(persistent_pairs_), cmp);
+ bool has_infinity = std::numeric_limits<Filtration_value>::has_infinity;
for (auto pair : persistent_pairs_) {
- diagram_out << cpx_->dimension(get<0>(pair)) << " "
- << cpx_->filtration(get<0>(pair)) << " "
- << cpx_->filtration(get<1>(pair)) << std::endl;
+ // Special case on windows, inf is "1.#INF"
+ if (has_infinity && cpx_->filtration(get<1>(pair)) == std::numeric_limits<Filtration_value>::infinity()) {
+ diagram_out << cpx_->dimension(get<0>(pair)) << " "
+ << cpx_->filtration(get<0>(pair)) << " inf" << std::endl;
+ } else {
+ diagram_out << cpx_->dimension(get<0>(pair)) << " "
+ << cpx_->filtration(get<0>(pair)) << " "
+ << cpx_->filtration(get<1>(pair)) << std::endl;
+ }
}
}
diff --git a/include/gudhi/Simplex_tree.h b/include/gudhi/Simplex_tree.h
index 317bce23..37b3ea97 100644
--- a/include/gudhi/Simplex_tree.h
+++ b/include/gudhi/Simplex_tree.h
@@ -289,7 +289,6 @@ class Simplex_tree {
/** \brief Constructs an empty simplex tree. */
Simplex_tree()
: null_vertex_(-1),
- threshold_(0),
root_(nullptr, null_vertex_),
filtration_vect_(),
dimension_(-1) { }
@@ -297,7 +296,6 @@ class Simplex_tree {
/** \brief User-defined copy constructor reproduces the whole tree structure. */
Simplex_tree(const Simplex_tree& simplex_source)
: null_vertex_(simplex_source.null_vertex_),
- threshold_(simplex_source.threshold_),
root_(nullptr, null_vertex_ , simplex_source.root_.members_),
filtration_vect_(),
dimension_(simplex_source.dimension_) {
@@ -323,12 +321,10 @@ class Simplex_tree {
/** \brief User-defined move constructor moves the whole tree structure. */
Simplex_tree(Simplex_tree && old)
: null_vertex_(std::move(old.null_vertex_)),
- threshold_(std::move(old.threshold_)),
root_(std::move(old.root_)),
filtration_vect_(std::move(old.filtration_vect_)),
dimension_(std::move(old.dimension_)) {
old.dimension_ = -1;
- old.threshold_ = 0;
old.root_ = Siblings(nullptr, null_vertex_);
}
@@ -356,7 +352,6 @@ class Simplex_tree {
/** \brief Checks if two simplex trees are equal. */
bool operator==(Simplex_tree& st2) {
if ((null_vertex_ != st2.null_vertex_) ||
- (threshold_ != st2.threshold_) ||
(dimension_ != st2.dimension_))
return false;
return rec_equal(&root_, &st2.root_);
@@ -407,14 +402,14 @@ class Simplex_tree {
/** \brief Returns the filtration value of a simplex.
*
- * Called on the null_simplex, returns INFINITY.
+ * Called on the null_simplex, it returns infinity.
* If SimplexTreeOptions::store_filtration is false, returns 0.
*/
static Filtration_value filtration(Simplex_handle sh) {
if (sh != null_simplex()) {
return sh->second.filtration();
} else {
- return INFINITY;
+ return std::numeric_limits<Filtration_value>::infinity();
}
}
@@ -427,11 +422,6 @@ class Simplex_tree {
sh->second.assign_filtration(fv);
}
- /** \brief Returns an upper bound of the filtration values of the simplices. */
- Filtration_value filtration() const {
- return threshold_;
- }
-
/** \brief Returns a Simplex_handle different from all Simplex_handles
* associated to the simplices in the simplicial complex.
*
@@ -757,11 +747,6 @@ class Simplex_tree {
return &root_;
}
- /** Set an upper bound for the filtration values. */
- void set_filtration(Filtration_value fil) {
- threshold_ = fil;
- }
-
/** Set a dimension for the simplicial complex. */
void set_dimension(int dimension) {
dimension_ = dimension;
@@ -1215,8 +1200,6 @@ class Simplex_tree {
private:
Vertex_handle null_vertex_;
- /** \brief Upper bound on the filtration values of the simplices.*/
- Filtration_value threshold_;
/** \brief Total number of simplices in the complex, without the empty simplex.*/
/** \brief Set of simplex tree Nodes representing the vertices.*/
Siblings root_;
@@ -1244,7 +1227,6 @@ std::istream& operator>>(std::istream & is, Simplex_tree<T...> & st) {
typedef Simplex_tree<T...> ST;
std::vector<typename ST::Vertex_handle> simplex;
typename ST::Filtration_value fil;
- typename ST::Filtration_value max_fil = 0;
int max_dim = -1;
while (read_simplex(is, simplex, fil)) {
// read all simplices in the file as a list of vertices
@@ -1253,15 +1235,11 @@ std::istream& operator>>(std::istream & is, Simplex_tree<T...> & st) {
if (max_dim < dim) {
max_dim = dim;
}
- if (max_fil < fil) {
- max_fil = fil;
- }
// insert every simplex in the simplex tree
st.insert_simplex(simplex, fil);
simplex.clear();
}
st.set_dimension(max_dim);
- st.set_filtration(max_fil);
return is;
}
diff --git a/include/gudhi/Strong_witness_complex.h b/include/gudhi/Strong_witness_complex.h
index a973ddb7..6f4bcf60 100644
--- a/include/gudhi/Strong_witness_complex.h
+++ b/include/gudhi/Strong_witness_complex.h
@@ -105,10 +105,6 @@ class Strong_witness_complex {
<< "non-negative.\n";
return false;
}
- if (limit_dimension < 0) {
- std::cerr << "Strong witness complex cannot create complex - limit dimension must be non-negative.\n";
- return false;
- }
for (auto w : nearest_landmark_table_) {
ActiveWitness aw(w);
typeVectorVertex simplex;
diff --git a/include/gudhi/Tangential_complex.h b/include/gudhi/Tangential_complex.h
index cfc82eb1..a5cefd6a 100644
--- a/include/gudhi/Tangential_complex.h
+++ b/include/gudhi/Tangential_complex.h
@@ -1059,7 +1059,6 @@ class Tangential_complex {
Triangulation &triangulation, bool verbose = false) {
int tangent_space_dim = tsb.dimension();
const Tr_traits &local_tr_traits = triangulation.geom_traits();
- Tr_vertex_handle center_vertex;
// Kernel functor & objects
typename K::Squared_distance_d k_sqdist = m_k.squared_distance_d_object();
@@ -1084,7 +1083,7 @@ class Tangential_complex {
proj_wp = project_point_and_compute_weight(wp, tsb, local_tr_traits);
}
- center_vertex = triangulation.insert(proj_wp);
+ Tr_vertex_handle center_vertex = triangulation.insert(proj_wp);
center_vertex->data() = i;
if (verbose)
std::cerr << "* Inserted point #" << i << "\n";
@@ -1094,8 +1093,8 @@ class Tangential_complex {
std::size_t num_inserted_points = 1;
#endif
// const int NUM_NEIGHBORS = 150;
- // KNS_range ins_range = m_points_ds.query_k_nearest_neighbors(center_pt, NUM_NEIGHBORS);
- INS_range ins_range = m_points_ds.query_incremental_nearest_neighbors(center_pt);
+ // KNS_range ins_range = m_points_ds.k_nearest_neighbors(center_pt, NUM_NEIGHBORS);
+ INS_range ins_range = m_points_ds.incremental_nearest_neighbors(center_pt);
// While building the local triangulation, we keep the radius
// of the sphere "star sphere" centered at "center_vertex"
@@ -1204,7 +1203,7 @@ class Tangential_complex {
Point center_point = compute_perturbed_point(i);
// Among updated point, what is the closer from our center point?
std::size_t closest_pt_index =
- updated_pts_ds.query_k_nearest_neighbors(center_point, 1, false).begin()->first;
+ updated_pts_ds.k_nearest_neighbors(center_point, 1, false).begin()->first;
typename K::Construct_weighted_point_d k_constr_wp =
m_k.construct_weighted_point_d_object();
@@ -1316,11 +1315,10 @@ class Tangential_complex {
m_k.compute_coordinate_d_object();
#ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
- KNS_range kns_range = m_points_ds_for_tse.query_k_nearest_neighbors(
- p, num_pts_for_pca, false);
+ KNS_range kns_range = m_points_ds_for_tse.k_nearest_neighbors(p, num_pts_for_pca, false);
const Points &points_for_pca = m_points_for_tse;
#else
- KNS_range kns_range = m_points_ds.query_k_nearest_neighbors(p, num_pts_for_pca, false);
+ KNS_range kns_range = m_points_ds.k_nearest_neighbors(p, num_pts_for_pca, false);
const Points &points_for_pca = m_points;
#endif
@@ -1414,11 +1412,10 @@ class Tangential_complex {
const Point &p = m_points[*it_index];
#ifdef GUDHI_TC_USE_ANOTHER_POINT_SET_FOR_TANGENT_SPACE_ESTIM
- KNS_range kns_range = m_points_ds_for_tse.query_k_nearest_neighbors(
- p, num_pts_for_pca, false);
+ KNS_range kns_range = m_points_ds_for_tse.k_nearest_neighbors(p, num_pts_for_pca, false);
const Points &points_for_pca = m_points_for_tse;
#else
- KNS_range kns_range = m_points_ds.query_k_nearest_neighbors(p, num_pts_for_pca, false);
+ KNS_range kns_range = m_points_ds.k_nearest_neighbors(p, num_pts_for_pca, false);
const Points &points_for_pca = m_points;
#endif
@@ -1823,7 +1820,6 @@ class Tangential_complex {
bool is_inconsistent = false;
Star const& star = m_stars[tr_index];
- Tr_vertex_handle center_vh = m_triangulations[tr_index].center_vertex();
// For each incident simplex
Star::const_iterator it_inc_simplex = star.begin();
@@ -1950,7 +1946,6 @@ class Tangential_complex {
bool is_star_inconsistent = false;
Triangulation const& tr = it_tr->tr();
- Tr_vertex_handle center_vh = it_tr->center_vertex();
if (tr.current_dimension() < m_intrinsic_dim)
continue;
diff --git a/include/gudhi/Witness_complex.h b/include/gudhi/Witness_complex.h
index 63f03687..bcfe8484 100644
--- a/include/gudhi/Witness_complex.h
+++ b/include/gudhi/Witness_complex.h
@@ -106,10 +106,6 @@ class Witness_complex {
std::cerr << "Witness complex cannot create complex - squared relaxation parameter must be non-negative.\n";
return false;
}
- if (limit_dimension < 0) {
- std::cerr << "Witness complex cannot create complex - limit dimension must be non-negative.\n";
- return false;
- }
ActiveWitnessList active_witnesses;
Landmark_id k = 0; /* current dimension in iterative construction */
for (auto w : nearest_landmark_table_)
diff --git a/include/gudhi/pick_n_random_points.h b/include/gudhi/pick_n_random_points.h
index f0e3f1f1..8c90b6bf 100644
--- a/include/gudhi/pick_n_random_points.h
+++ b/include/gudhi/pick_n_random_points.h
@@ -52,7 +52,7 @@ typename OutputIterator>
void pick_n_random_points(Point_container const &points,
std::size_t final_size,
OutputIterator output_it) {
-#ifdef GUDHI_SUBS_PROFILING
+#ifdef GUDHI_SUBSAMPLING_PROFILING
Gudhi::Clock t;
#endif
@@ -72,7 +72,7 @@ void pick_n_random_points(Point_container const &points,
for (int l : landmarks)
*output_it++ = points[l];
-#ifdef GUDHI_SUBS_PROFILING
+#ifdef GUDHI_SUBSAMPLING_PROFILING
t.end();
std::cerr << "Random landmark choice took " << t.num_seconds()
<< " seconds." << std::endl;
diff --git a/include/gudhi/random_point_generators.h b/include/gudhi/random_point_generators.h
index 2ec465ef..9df77760 100644
--- a/include/gudhi/random_point_generators.h
+++ b/include/gudhi/random_point_generators.h
@@ -282,6 +282,38 @@ std::vector<typename Kernel::Point_d> generate_points_on_sphere_d(std::size_t nu
}
template <typename Kernel>
+std::vector<typename Kernel::Point_d> generate_points_in_ball_d(std::size_t num_points, int dim, double radius) {
+ typedef typename Kernel::Point_d Point;
+ Kernel k;
+ CGAL::Random rng;
+ CGAL::Random_points_in_ball_d<Point> generator(dim, radius);
+ std::vector<Point> points;
+ points.reserve(num_points);
+ for (std::size_t i = 0; i < num_points;) {
+ Point p = *generator++;
+ points.push_back(p);
+ ++i;
+ }
+ return points;
+}
+
+template <typename Kernel>
+std::vector<typename Kernel::Point_d> generate_points_in_cube_d(std::size_t num_points, int dim, double radius) {
+ typedef typename Kernel::Point_d Point;
+ Kernel k;
+ CGAL::Random rng;
+ CGAL::Random_points_in_cube_d<Point> generator(dim, radius);
+ std::vector<Point> points;
+ points.reserve(num_points);
+ for (std::size_t i = 0; i < num_points;) {
+ Point p = *generator++;
+ points.push_back(p);
+ ++i;
+ }
+ return points;
+}
+
+template <typename Kernel>
std::vector<typename Kernel::Point_d> generate_points_on_two_spheres_d(std::size_t num_points, int dim, double radius,
double distance_between_centers,
double radius_noise_percentage = 0.) {
diff --git a/include/gudhi/reader_utils.h b/include/gudhi/reader_utils.h
index 97a87edd..90be4fc7 100644
--- a/include/gudhi/reader_utils.h
+++ b/include/gudhi/reader_utils.h
@@ -1,8 +1,8 @@
-/* This file is part of the Gudhi Library. The Gudhi library
- * (Geometric Understanding in Higher Dimensions) is a generic C++
+/* 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(s): Clement Maria, Pawel Dlotko
+ * Author(s): Clement Maria, Pawel Dlotko, Clement Jamin
*
* Copyright (C) 2014 INRIA
*
@@ -24,7 +24,9 @@
#define READER_UTILS_H_
#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/Debug_utils.h>
+#include <boost/function_output_iterator.hpp>
#include <boost/graph/adjacency_list.hpp>
#include <iostream>
@@ -34,6 +36,9 @@
#include <string>
#include <vector>
#include <utility> // for pair
+#include <tuple> // for std::make_tuple
+
+namespace Gudhi {
// Keep this file tag for Doxygen to parse the code, otherwise, functions are not documented.
// It is required for global functions and variables.
@@ -50,7 +55,7 @@
* X21 X22 ... X2d<br>
* etc<br>
*/
-inline void read_points(std::string file_name, std::vector< std::vector< double > > & points) {
+inline void read_points(std::string file_name, std::vector<std::vector<double>>& points) {
std::ifstream in_file(file_name.c_str(), std::ios::in);
if (!in_file.is_open()) {
std::cerr << "Unable to open file " << file_name << std::endl;
@@ -60,14 +65,13 @@ inline void read_points(std::string file_name, std::vector< std::vector< double
std::string line;
double x;
while (getline(in_file, line)) {
- std::vector< double > point;
+ std::vector<double> point;
std::istringstream iss(line);
while (iss >> x) {
point.push_back(x);
}
// Check for empty lines
- if (!point.empty())
- points.push_back(point);
+ if (!point.empty()) points.push_back(point);
}
in_file.close();
}
@@ -88,17 +92,20 @@ inline void read_points(std::string file_name, std::vector< std::vector< double
* Every simplex must appear exactly once.
* Simplices of dimension more than 1 are ignored.
*/
-template< typename Graph_t, typename Filtration_value, typename Vertex_handle >
+template <typename Graph_t, typename Filtration_value, typename Vertex_handle>
Graph_t read_graph(std::string file_name) {
std::ifstream in_(file_name.c_str(), std::ios::in);
if (!in_.is_open()) {
- std::cerr << "Unable to open file " << file_name << std::endl;
+ std::string error_str("read_graph - Unable to open file ");
+ error_str.append(file_name);
+ std::cerr << error_str << std::endl;
+ throw std::invalid_argument(error_str);
}
- typedef std::pair< Vertex_handle, Vertex_handle > Edge_t;
- std::vector< Edge_t > edges;
- std::vector< Filtration_value > edges_fil;
- std::map< Vertex_handle, Filtration_value > vertices;
+ typedef std::pair<Vertex_handle, Vertex_handle> Edge_t;
+ std::vector<Edge_t> edges;
+ std::vector<Filtration_value> edges_fil;
+ std::map<Vertex_handle, Filtration_value> vertices;
std::string line;
int dim;
@@ -108,8 +115,7 @@ Graph_t read_graph(std::string file_name) {
std::istringstream iss(line);
while (iss >> dim) {
switch (dim) {
- case 0:
- {
+ case 0: {
iss >> u;
iss >> fil;
vertices[u] = fil;
@@ -118,8 +124,7 @@ Graph_t read_graph(std::string file_name) {
}
break;
}
- case 1:
- {
+ case 1: {
iss >> u;
iss >> v;
iss >> fil;
@@ -127,16 +132,13 @@ Graph_t read_graph(std::string file_name) {
edges_fil.push_back(fil);
break;
}
- default:
- {
- break;
- }
+ default: { break; }
}
}
}
in_.close();
- if ((size_t) (max_h + 1) != vertices.size()) {
+ if ((size_t)(max_h + 1) != vertices.size()) {
std::cerr << "Error: vertices must be labeled from 0 to n-1 \n";
}
@@ -164,8 +166,8 @@ Graph_t read_graph(std::string file_name) {
* Every simplex must appear exactly once.
* Simplices of dimension more than 1 are ignored.
*/
-template< typename Vertex_handle, typename Filtration_value >
-bool read_simplex(std::istream & in_, std::vector< Vertex_handle > & simplex, Filtration_value & fil) {
+template <typename Vertex_handle, typename Filtration_value>
+bool read_simplex(std::istream& in_, std::vector<Vertex_handle>& simplex, Filtration_value& fil) {
int dim = 0;
if (!(in_ >> dim)) return false;
Vertex_handle v;
@@ -189,8 +191,8 @@ bool read_simplex(std::istream & in_, std::vector< Vertex_handle > & simplex, Fi
* The key of a simplex is its position in the filtration order and also the number of its row in the file.
* Dimi ki1 ki2 ... kiDimi Fili means that the ith simplex in the filtration has dimension Dimi, filtration value
* fil1 and simplices with key ki1 ... kiDimi in its boundary.*/
-template< typename Simplex_key, typename Filtration_value >
-bool read_hasse_simplex(std::istream & in_, std::vector< Simplex_key > & boundary, Filtration_value & fil) {
+template <typename Simplex_key, typename Filtration_value>
+bool read_hasse_simplex(std::istream& in_, std::vector<Simplex_key>& boundary, Filtration_value& fil) {
int dim;
if (!(in_ >> dim)) return false;
if (dim == 0) {
@@ -209,7 +211,7 @@ bool read_hasse_simplex(std::istream & in_, std::vector< Simplex_key > & boundar
/**
* @brief Read a lower triangular distance matrix from a csv file. We assume that the .csv store the whole
* (square) matrix.
- *
+ *
* @author Pawel Dlotko
*
* Square matrix file format:<br>
@@ -226,13 +228,13 @@ bool read_hasse_simplex(std::istream & in_, std::vector< Simplex_key > & boundar
* Dj1;Dj2;...;Dj(j-1);<br>
*
**/
-template< typename Filtration_value >
-std::vector< std::vector< Filtration_value > > read_lower_triangular_matrix_from_csv_file(const std::string& filename,
- const char separator = ';') {
+template <typename Filtration_value>
+std::vector<std::vector<Filtration_value>> read_lower_triangular_matrix_from_csv_file(const std::string& filename,
+ const char separator = ';') {
#ifdef DEBUG_TRACES
std::cout << "Using procedure read_lower_triangular_matrix_from_csv_file \n";
#endif // DEBUG_TRACES
- std::vector< std::vector< Filtration_value > > result;
+ std::vector<std::vector<Filtration_value>> result;
std::ifstream in;
in.open(filename.c_str());
if (!in.is_open()) {
@@ -243,7 +245,7 @@ std::vector< std::vector< Filtration_value > > read_lower_triangular_matrix_from
// the first line is emtpy, so we ignore it:
std::getline(in, line);
- std::vector< Filtration_value > values_in_this_line;
+ std::vector<Filtration_value> values_in_this_line;
result.push_back(values_in_this_line);
int number_of_line = 0;
@@ -251,11 +253,10 @@ std::vector< std::vector< Filtration_value > > read_lower_triangular_matrix_from
// first, read the file line by line to a string:
while (std::getline(in, line)) {
// if line is empty, break
- if (line.size() == 0)
- break;
+ if (line.size() == 0) break;
// if the last element of a string is comma:
- if (line[ line.size() - 1 ] == separator) {
+ if (line[line.size() - 1] == separator) {
// then shrink the string by one
line.pop_back();
}
@@ -268,7 +269,7 @@ std::vector< std::vector< Filtration_value > > read_lower_triangular_matrix_from
// and now read the doubles.
int number_of_entry = 0;
- std::vector< Filtration_value > values_in_this_line;
+ std::vector<Filtration_value> values_in_this_line;
while (iss.good()) {
double entry;
iss >> entry;
@@ -277,7 +278,7 @@ std::vector< std::vector< Filtration_value > > read_lower_triangular_matrix_from
}
++number_of_entry;
}
- if (!values_in_this_line.empty())result.push_back(values_in_this_line);
+ if (!values_in_this_line.empty()) result.push_back(values_in_this_line);
++number_of_line;
}
in.close();
@@ -295,4 +296,74 @@ std::vector< std::vector< Filtration_value > > read_lower_triangular_matrix_from
return result;
} // read_lower_triangular_matrix_from_csv_file
+/**
+Reads a file containing persistence intervals.
+Each line might contain 2, 3 or 4 values: [[field] dimension] birth death
+The output iterator `out` is used this way: `*out++ = std::make_tuple(dim, birth, death);`
+where `dim` is an `int`, `birth` a `double`, and `death` a `double`.
+Note: the function does not check that birth <= death.
+**/
+template <typename OutputIterator>
+void read_persistence_intervals_and_dimension(std::string const& filename, OutputIterator out) {
+ std::ifstream in(filename);
+ if (!in.is_open()) {
+ std::string error_str("read_persistence_intervals_and_dimension - Unable to open file ");
+ error_str.append(filename);
+ std::cerr << error_str << std::endl;
+ throw std::invalid_argument(error_str);
+ }
+
+ while (!in.eof()) {
+ std::string line;
+ getline(in, line);
+ if (line.length() != 0 && line[0] != '#') {
+ double numbers[4];
+ int n = sscanf(line.c_str(), "%lf %lf %lf %lf", &numbers[0], &numbers[1], &numbers[2], &numbers[3]);
+ if (n >= 2) {
+ int dim = (n >= 3 ? static_cast<int>(numbers[n - 3]) : -1);
+ *out++ = std::make_tuple(dim, numbers[n - 2], numbers[n - 1]);
+ }
+ }
+ }
+}
+
+/**
+Reads a file containing persistence intervals.
+Each line might contain 2, 3 or 4 values: [[field] dimension] birth death
+The return value is an `std::map<dim, std::vector<std::pair<birth, death>>>`
+where `dim` is an `int`, `birth` a `double`, and `death` a `double`.
+Note: the function does not check that birth <= death.
+**/
+inline std::map<int, std::vector<std::pair<double, double>>> read_persistence_intervals_grouped_by_dimension(
+ std::string const& filename) {
+ std::map<int, std::vector<std::pair<double, double>>> ret;
+ read_persistence_intervals_and_dimension(
+ filename, boost::make_function_output_iterator([&ret](std::tuple<int, double, double> t) {
+ ret[get<0>(t)].push_back(std::make_pair(get<1>(t), get<2>(t)));
+ }));
+ return ret;
+}
+
+/**
+Reads a file containing persistence intervals.
+Each line might contain 2, 3 or 4 values: [[field] dimension] birth death
+If `only_this_dim` = -1, dimension is ignored and all lines are returned.
+If `only_this_dim` is >= 0, only the lines where dimension = `only_this_dim`
+(or where dimension is not specified) are returned.
+The return value is an `std::vector<std::pair<birth, death>>`
+where `dim` is an `int`, `birth` a `double`, and `death` a `double`.
+Note: the function does not check that birth <= death.
+**/
+inline std::vector<std::pair<double, double>> read_persistence_intervals_in_dimension(std::string const& filename,
+ int only_this_dim = -1) {
+ std::vector<std::pair<double, double>> ret;
+ read_persistence_intervals_and_dimension(
+ filename, boost::make_function_output_iterator([only_this_dim, &ret](std::tuple<int, double, double> t) {
+ if (only_this_dim == get<0>(t) || only_this_dim == -1) ret.emplace_back(get<1>(t), get<2>(t));
+ }));
+ return ret;
+}
+
+} // namespace Gudhi
+
#endif // READER_UTILS_H_
diff --git a/include/gudhi/sparsify_point_set.h b/include/gudhi/sparsify_point_set.h
index 507f8c79..7d3b97fb 100644
--- a/include/gudhi/sparsify_point_set.h
+++ b/include/gudhi/sparsify_point_set.h
@@ -83,7 +83,7 @@ sparsify_point_set(
*output_it++ = *it_pt;
- auto ins_range = points_ds.query_incremental_nearest_neighbors(*it_pt);
+ auto ins_range = points_ds.incremental_nearest_neighbors(*it_pt);
// If another point Q is closer that min_squared_dist, mark Q to be dropped
for (auto const& neighbor : ins_range) {