From cb838b2ea4a4db9c54f71103001bdafb90766306 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Tue, 24 Mar 2020 06:37:00 +0100 Subject: merge https://github.com/mglisse/gudhi-devel/tree/alpha-cache and fix conflicts --- src/Alpha_complex/include/gudhi/Alpha_complex.h | 89 ++++++++++--------------- 1 file changed, 37 insertions(+), 52 deletions(-) (limited to 'src/Alpha_complex') diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 1b5d6997..eb4ef427 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -132,6 +132,8 @@ class Alpha_complex { Delaunay_triangulation* triangulation_; /** \brief Kernel for triangulation_ functions access.*/ Kernel kernel_; + /** \brief Cache for geometric constructions: circumcenter and squared radius of a simplex.*/ + std::vector> cache_; public: /** \brief Alpha_complex constructor from an OFF file name. @@ -246,6 +248,24 @@ class Alpha_complex { } } + template + auto& get_cache(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) { + auto k = cplx.key(s); + if(k==cplx.null_key()){ + k = cache_.size(); + cplx.assign_key(s, k); + // Use a transform_range? Check the impact on perf. + thread_local std::vector v; + v.clear(); + for (auto vertex : cplx.simplex_vertex_range(s)) + v.push_back(get_point(vertex)); + Point_d c = kernel_.construct_circumcenter_d_object()(v.cbegin(), v.cend()); + typename Kernel::FT r = kernel_.squared_distance_d_object()(c, v[0]); + cache_.emplace_back(std::move(c), std::move(r)); + } + return cache_[k]; + } + public: /** \brief Inserts all Delaunay triangulation into the simplicial complex. * It also computes the filtration values accordingly to the \ref createcomplexalgorithm if default_filtration_value @@ -324,46 +344,28 @@ class Alpha_complex { if (!default_filtration_value) { // -------------------------------------------------------------------------------------------- - // Will be re-used many times - Vector_of_CGAL_points pointVector; // ### For i : d -> 0 for (int decr_dim = triangulation_->maximal_dimension(); decr_dim >= 0; decr_dim--) { // ### Foreach Sigma of dim i for (Simplex_handle f_simplex : complex.skeleton_simplex_range(decr_dim)) { int f_simplex_dim = complex.dimension(f_simplex); if (decr_dim == f_simplex_dim) { - pointVector.clear(); - #ifdef DEBUG_TRACES - std::clog << "Sigma of dim " << decr_dim << " is"; - #endif // DEBUG_TRACES - for (auto vertex : complex.simplex_vertex_range(f_simplex)) { - pointVector.push_back(get_point(vertex)); - #ifdef DEBUG_TRACES - std::clog << " " << vertex; - #endif // DEBUG_TRACES - } - #ifdef DEBUG_TRACES - std::clog << std::endl; - #endif // DEBUG_TRACES // ### If filt(Sigma) is NaN : filt(Sigma) = alpha(Sigma) if (std::isnan(complex.filtration(f_simplex))) { Filtration_value alpha_complex_filtration = 0.0; // No need to compute squared_radius on a single point - alpha is 0.0 if (f_simplex_dim > 0) { - // squared_radius function initialization - Squared_Radius squared_radius = kernel_.compute_squared_radius_d_object(); - - CGAL::NT_converter cv; - auto sqrad = squared_radius(pointVector.begin(), pointVector.end()); - #if CGAL_VERSION_NR >= 1050000000 + auto const& sqrad = get_cache(complex, f_simplex).second; +#if CGAL_VERSION_NR >= 1050000000 if(exact) CGAL::exact(sqrad); - #endif +#endif + CGAL::NT_converter cv; alpha_complex_filtration = cv(sqrad); } complex.assign_filtration(f_simplex, alpha_complex_filtration); - #ifdef DEBUG_TRACES +#ifdef DEBUG_TRACES std::clog << "filt(Sigma) is NaN : filt(Sigma) =" << complex.filtration(f_simplex) << std::endl; - #endif // DEBUG_TRACES +#endif // DEBUG_TRACES } // No need to propagate further, unweighted points all have value 0 if (decr_dim > 1) @@ -388,9 +390,7 @@ class Alpha_complex { void propagate_alpha_filtration(SimplicialComplexForAlpha& complex, Simplex_handle f_simplex) { // From SimplicialComplexForAlpha type required to assign filtration values. typedef typename SimplicialComplexForAlpha::Filtration_value Filtration_value; -#ifdef DEBUG_TRACES typedef typename SimplicialComplexForAlpha::Vertex_handle Vertex_handle; -#endif // DEBUG_TRACES // ### Foreach Tau face of Sigma for (auto f_boundary : complex.boundary_simplex_range(f_simplex)) { @@ -414,33 +414,18 @@ class Alpha_complex { #endif // DEBUG_TRACES // ### Else } else { - // insert the Tau points in a vector for is_gabriel function - Vector_of_CGAL_points pointVector; -#ifdef DEBUG_TRACES - Vertex_handle vertexForGabriel = Vertex_handle(); -#endif // DEBUG_TRACES - for (auto vertex : complex.simplex_vertex_range(f_boundary)) { - pointVector.push_back(get_point(vertex)); - } - // Retrieve the Sigma point that is not part of Tau - parameter for is_gabriel function - Point_d point_for_gabriel; - for (auto vertex : complex.simplex_vertex_range(f_simplex)) { - point_for_gabriel = get_point(vertex); - if (std::find(pointVector.begin(), pointVector.end(), point_for_gabriel) == pointVector.end()) { -#ifdef DEBUG_TRACES - // vertex is not found in Tau - vertexForGabriel = vertex; -#endif // DEBUG_TRACES - // No need to continue loop - break; - } - } - // is_gabriel function initialization - Is_Gabriel is_gabriel = kernel_.side_of_bounded_sphere_d_object(); - bool is_gab = is_gabriel(pointVector.begin(), pointVector.end(), point_for_gabriel) - != CGAL::ON_BOUNDED_SIDE; + // Find which vertex of f_simplex is missing in f_boundary. We could actually write a variant of boundary_simplex_range that gives pairs (f_boundary, vertex). We rely on the fact that simplex_vertex_range is sorted. + auto longlist = complex.simplex_vertex_range(f_simplex); + auto shortlist = complex.simplex_vertex_range(f_boundary); + auto longiter = std::begin(longlist); + auto shortiter = std::begin(shortlist); + auto enditer = std::end(shortlist); + while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; } + Vertex_handle extra = *longiter; + auto const& cache=get_cache(complex, f_boundary); + bool is_gab = kernel_.squared_distance_d_object()(cache.first, get_point(extra)) >= cache.second; #ifdef DEBUG_TRACES - std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << vertexForGabriel << std::endl; + std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << extra << std::endl; #endif // DEBUG_TRACES // ### If Tau is not Gabriel of Sigma if (false == is_gab) { -- cgit v1.2.3 From 9d89f57376e619515d99ad88c2cdeef35daaedd5 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Wed, 1 Apr 2020 09:04:18 +0200 Subject: code review: use operator[] instead of at() --- src/Alpha_complex/include/gudhi/Alpha_complex.h | 14 ++++++++++++-- 1 file changed, 12 insertions(+), 2 deletions(-) (limited to 'src/Alpha_complex') diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index eb4ef427..4369071c 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -248,6 +248,16 @@ class Alpha_complex { } } + /** \brief get_point_ returns the point corresponding to the vertex given as parameter. + * Only for internal use for faster access. + * + * @param[in] vertex Vertex handle of the point to retrieve. + * @return The point found. + */ + const Point_d& get_point_(std::size_t vertex) const { + return vertex_handle_to_iterator_[vertex]->point(); + } + template auto& get_cache(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) { auto k = cplx.key(s); @@ -258,7 +268,7 @@ class Alpha_complex { thread_local std::vector v; v.clear(); for (auto vertex : cplx.simplex_vertex_range(s)) - v.push_back(get_point(vertex)); + v.push_back(get_point_(vertex)); Point_d c = kernel_.construct_circumcenter_d_object()(v.cbegin(), v.cend()); typename Kernel::FT r = kernel_.squared_distance_d_object()(c, v[0]); cache_.emplace_back(std::move(c), std::move(r)); @@ -423,7 +433,7 @@ class Alpha_complex { while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; } Vertex_handle extra = *longiter; auto const& cache=get_cache(complex, f_boundary); - bool is_gab = kernel_.squared_distance_d_object()(cache.first, get_point(extra)) >= cache.second; + bool is_gab = kernel_.squared_distance_d_object()(cache.first, get_point_(extra)) >= cache.second; #ifdef DEBUG_TRACES std::clog << " | Tau is_gabriel(Sigma)=" << is_gab << " - vertexForGabriel=" << extra << std::endl; #endif // DEBUG_TRACES -- cgit v1.2.3 From 0b19e1f991fdebbdb622d3101135eaee65c4ed5d Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Wed, 1 Apr 2020 14:45:37 +0200 Subject: Split the cache per dimension Try to reduce slightly the memory use. --- src/Alpha_complex/include/gudhi/Alpha_complex.h | 33 +++++++++++++++++++------ 1 file changed, 25 insertions(+), 8 deletions(-) (limited to 'src/Alpha_complex') diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index 4369071c..ba91998d 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -112,9 +112,6 @@ class Alpha_complex { typedef typename Kernel::Side_of_bounded_sphere_d Is_Gabriel; typedef typename Kernel::Point_dimension_d Point_Dimension; - // Type required to compute squared radius, or side of bounded sphere on a vector of points. - typedef typename std::vector Vector_of_CGAL_points; - // Vertex_iterator type from CGAL. typedef typename Delaunay_triangulation::Vertex_iterator CGAL_vertex_iterator; @@ -124,6 +121,9 @@ class Alpha_complex { // Structure to switch from simplex tree vertex handle to CGAL vertex iterator. typedef typename std::vector< CGAL_vertex_iterator > Vector_vertex_iterator; + // Numeric type of coordinates in the kernel + typedef typename Kernel::FT FT; + private: /** \brief Vertex iterator vector to switch from simplex tree vertex handle to CGAL vertex iterator. * Vertex handles are inserted sequentially, starting at 0.*/ @@ -133,7 +133,7 @@ class Alpha_complex { /** \brief Kernel for triangulation_ functions access.*/ Kernel kernel_; /** \brief Cache for geometric constructions: circumcenter and squared radius of a simplex.*/ - std::vector> cache_; + std::vector> cache_, old_cache_; public: /** \brief Alpha_complex constructor from an OFF file name. @@ -258,24 +258,39 @@ class Alpha_complex { return vertex_handle_to_iterator_[vertex]->point(); } + /// Return a reference to the circumcenter and circumradius, writing them in the cache if necessary. template auto& get_cache(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) { auto k = cplx.key(s); if(k==cplx.null_key()){ k = cache_.size(); cplx.assign_key(s, k); - // Use a transform_range? Check the impact on perf. + // Using a transform_range is slower, currently. thread_local std::vector v; v.clear(); for (auto vertex : cplx.simplex_vertex_range(s)) v.push_back(get_point_(vertex)); Point_d c = kernel_.construct_circumcenter_d_object()(v.cbegin(), v.cend()); - typename Kernel::FT r = kernel_.squared_distance_d_object()(c, v[0]); + FT r = kernel_.squared_distance_d_object()(c, v[0]); cache_.emplace_back(std::move(c), std::move(r)); } return cache_[k]; } + /// Return the circumradius, either from the old cache or computed, without writing to the cache. + template + auto radius(SimplicialComplexForAlpha& cplx, typename SimplicialComplexForAlpha::Simplex_handle s) { + auto k = cplx.key(s); + if(k!=cplx.null_key()) + return old_cache_[k].second; + // Using a transform_range is slower, currently. + thread_local std::vector v; + v.clear(); + for (auto vertex : cplx.simplex_vertex_range(s)) + v.push_back(get_point_(vertex)); + return kernel_.compute_squared_radius_d_object()(v.cbegin(), v.cend()); + } + public: /** \brief Inserts all Delaunay triangulation into the simplicial complex. * It also computes the filtration values accordingly to the \ref createcomplexalgorithm if default_filtration_value @@ -365,11 +380,11 @@ class Alpha_complex { Filtration_value alpha_complex_filtration = 0.0; // No need to compute squared_radius on a single point - alpha is 0.0 if (f_simplex_dim > 0) { - auto const& sqrad = get_cache(complex, f_simplex).second; + auto const& sqrad = radius(complex, f_simplex); #if CGAL_VERSION_NR >= 1050000000 if(exact) CGAL::exact(sqrad); #endif - CGAL::NT_converter cv; + CGAL::NT_converter cv; alpha_complex_filtration = cv(sqrad); } complex.assign_filtration(f_simplex, alpha_complex_filtration); @@ -382,6 +397,8 @@ class Alpha_complex { propagate_alpha_filtration(complex, f_simplex); } } + old_cache_ = std::move(cache_); + cache_.clear(); } // -------------------------------------------------------------------------------------------- -- cgit v1.2.3 From d302e90dcf4b284e6dc8b3ab21e8a67fb9cf5179 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Thu, 16 Apr 2020 15:40:45 +0200 Subject: Update the concept of the simplicial complex We use the key now. It wouldn't be hard to use an unordered_map, but since we usually have an unused field key... --- src/Alpha_complex/concept/SimplicialComplexForAlpha.h | 18 ++++++++++++++++++ 1 file changed, 18 insertions(+) (limited to 'src/Alpha_complex') diff --git a/src/Alpha_complex/concept/SimplicialComplexForAlpha.h b/src/Alpha_complex/concept/SimplicialComplexForAlpha.h index 1c6c3b0c..c20c3201 100644 --- a/src/Alpha_complex/concept/SimplicialComplexForAlpha.h +++ b/src/Alpha_complex/concept/SimplicialComplexForAlpha.h @@ -72,6 +72,24 @@ struct SimplicialComplexForAlpha { /** \brief Return type of an insertion of a simplex */ typedef unspecified Insertion_result_type; + + /** \name Map interface + * Conceptually a `std::unordered_map`. + * @{ */ + /** \brief Data stored for each simplex. + * + * Must be an integer type. */ + typedef unspecified Simplex_key; + /** \brief Returns a constant dummy number that is either negative, + * or at least as large as the number of simplices. Suggested value: -1. */ + Simplex_key null_key (); + /** \brief Returns the number stored for a simplex by `assign_key()`. + * + * If `assign_key()` has not been called, it must return `null_key()`. */ + Simplex_key key ( Simplex_handle sh ); + /** \brief Store a number for a simplex, which can later be retrieved with `key()`. */ + void assign_key(Simplex_handle sh, Simplex_key n); + /** @} */ }; } // namespace alpha_complex -- cgit v1.2.3