From 866f6ce614e9c09c97fed12c8c0c2c9fb84fad3f Mon Sep 17 00:00:00 2001 From: Gard Spreemann Date: Sun, 8 Oct 2017 11:15:17 +0200 Subject: GUDHI 2.0.1 as released by upstream in a tarball. --- include/gudhi/Bitmap_cubical_complex.h | 27 +++-- include/gudhi/Bottleneck.h | 4 +- include/gudhi/Clock.h | 20 ++-- include/gudhi/Euclidean_strong_witness_complex.h | 2 +- include/gudhi/Euclidean_witness_complex.h | 2 +- include/gudhi/Hasse_complex.h | 11 +- include/gudhi/Kd_tree_search.h | 81 ++++++++----- include/gudhi/Neighbors_finder.h | 8 +- include/gudhi/Persistence_graph.h | 8 +- include/gudhi/Persistent_cohomology.h | 13 +- include/gudhi/Simplex_tree.h | 26 +--- include/gudhi/Strong_witness_complex.h | 4 - include/gudhi/Tangential_complex.h | 21 ++-- include/gudhi/Witness_complex.h | 4 - include/gudhi/pick_n_random_points.h | 4 +- include/gudhi/random_point_generators.h | 32 +++++ include/gudhi/reader_utils.h | 145 +++++++++++++++++------ include/gudhi/sparsify_point_set.h | 2 +- 18 files changed, 253 insertions(+), 161 deletions(-) (limited to 'include') 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& dimensions, - const std::vector& top_dimensional_cells) : + const std::vector& 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& dimensions, - const std::vector& top_dimensional_cells, + const std::vector& 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::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::max()) return this->data[sh]; - return std::numeric_limits::max(); + if (sh != null_simplex()) return this->data[sh]; + return std::numeric_limits::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::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::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::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::Simplex_handle& sh1, const typename Bitmap_cubical_complex::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 double bottleneck_distance(const Persistence_diagram1 &diag1, const Persistence_diagram2 &diag2, - double e = std::numeric_limits::min()) { + double e = (std::numeric_limits::min)()) { Persistence_graph g(diag1, diag2, e); if (g.bottleneck_alive() == std::numeric_limits::infinity()) return std::numeric_limits::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 - +#include #include +#include 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(end-startTime).count() / 1000.; } else { - return (endTime - startTime).total_milliseconds() / 1000.; + return std::chrono::duration_cast(endTime-startTime).count() / 1000.; } } private: - mutable boost::posix_time::ptime startTime, endTime; + mutable std::chrono::time_point 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 #include // for pair #include +#include // for infinity value #ifdef GUDHI_USE_TBB #include @@ -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::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 > complex_; std::vector 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 #include #include +#include #include #include @@ -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 * CGAL dD spatial searching data structures. - * 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 k-nearest or k-farthest neighbor query, where k is fixed and the k nearest - * or farthest points are computed right away, - * and the incremental nearest or farthest neighbor query, where no number of neighbors is provided during the call, as the + * There are two types of queries: the k-nearest or k-furthest neighbor query, where k is fixed and the k nearest + * or furthest points are computed right away, + * and the incremental nearest or furthest neighbor query, 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 STraits; + typedef CGAL::Distance_adapter< + std::ptrdiff_t, + Point_property_map, + CGAL::Euclidean_distance > Orthogonal_distance; typedef CGAL::Orthogonal_k_neighbor_search 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` 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, 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` 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 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::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::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::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::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 + 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 const&r) const { + bool inner_range_intersects(CGAL::Kd_tree_rectangle 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 const&r) const { + bool outer_range_contains(CGAL::Kd_tree_rectangle 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::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::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::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::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 & st) { typedef Simplex_tree ST; std::vector 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 & 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 @@ -281,6 +281,38 @@ std::vector generate_points_on_sphere_d(std::size_t nu return points; } +template +std::vector 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 generator(dim, radius); + std::vector 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 +std::vector 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 generator(dim, radius); + std::vector 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 std::vector generate_points_on_two_spheres_d(std::size_t num_points, int dim, double radius, double distance_between_centers, 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 +#include +#include #include #include @@ -34,6 +36,9 @@ #include #include #include // for pair +#include // 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
* etc
*/ -inline void read_points(std::string file_name, std::vector< std::vector< double > > & points) { +inline void read_points(std::string file_name, std::vector>& 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 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 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 Edge_t; + std::vector edges; + std::vector edges_fil; + std::map 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 +bool read_simplex(std::istream& in_, std::vector& 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 +bool read_hasse_simplex(std::istream& in_, std::vector& 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:
@@ -226,13 +228,13 @@ bool read_hasse_simplex(std::istream & in_, std::vector< Simplex_key > & boundar * Dj1;Dj2;...;Dj(j-1);
* **/ -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 +std::vector> 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> 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 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 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 +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(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>>` +where `dim` is an `int`, `birth` a `double`, and `death` a `double`. +Note: the function does not check that birth <= death. +**/ +inline std::map>> read_persistence_intervals_grouped_by_dimension( + std::string const& filename) { + std::map>> ret; + read_persistence_intervals_and_dimension( + filename, boost::make_function_output_iterator([&ret](std::tuple 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>` +where `dim` is an `int`, `birth` a `double`, and `death` a `double`. +Note: the function does not check that birth <= death. +**/ +inline std::vector> read_persistence_intervals_in_dimension(std::string const& filename, + int only_this_dim = -1) { + std::vector> ret; + read_persistence_intervals_and_dimension( + filename, boost::make_function_output_iterator([only_this_dim, &ret](std::tuple 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) { -- cgit v1.2.3