From 767d9fca5da2d3dd9698a5c27e9bedc159271f67 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Thu, 30 Sep 2021 11:05:17 +0200 Subject: Move minimal enclosing balls cache to Cech_complex.h instead of the blocker Modify cech test to be more relevant regarding the current algorithm Do some cleaning --- src/Cech_complex/include/gudhi/Cech_complex.h | 38 +++---- .../include/gudhi/Cech_complex_blocker.h | 121 +++------------------ 2 files changed, 36 insertions(+), 123 deletions(-) (limited to 'src/Cech_complex/include/gudhi') diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 2c6f202a..32a78aec 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -40,7 +40,7 @@ namespace cech_complex { * \tparam ForwardPointRange must be a range for which `std::begin()` and `std::end()` methods return input * iterators on a point. `std::begin()` and `std::end()` methods are also required for a point. */ -template +template class Cech_complex { private: // Required by compute_proximity_graph @@ -48,17 +48,17 @@ class Cech_complex { using Filtration_value = typename SimplicialComplexForProximityGraph::Filtration_value; using Proximity_graph = Gudhi::Proximity_graph; - // Retrieve Coordinate type from ForwardPointRange -// using Point_from_range_iterator = typename boost::range_const_iterator::type; -// using Point_from_range = typename std::iterator_traits::value_type; -// using Coordinate_iterator = typename boost::range_const_iterator::type; -// using Coordinate = typename std::iterator_traits::value_type; - public: - // Point and Point_cloud type definition - //using Point = std::vector; - using Point = typename Kernel::Point_d; - using Point_cloud = std::vector; + + using cech_blocker = Cech_blocker; + + using Point_d = typename cech_blocker::Point_d; + using Point_cloud = std::vector; + + // Numeric type of coordinates in the kernel + using FT = typename cech_blocker::FT; + // Sphere is a pair of point and squared radius. + using Sphere = typename std::pair; public: /** \brief Cech_complex constructor from a list of points. @@ -77,7 +77,6 @@ class Cech_complex { point_cloud_.assign(points.begin(), points.end()); - std::clog << "Hind: Just before the graph compute" << std::endl; cech_skeleton_graph_ = Gudhi::compute_proximity_graph( point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); } @@ -90,19 +89,14 @@ class Cech_complex { * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0. * */ - template void create_complex(SimplicialComplexForCechComplex& complex, int dim_max) { - std::clog << "Hind: in create complex" << std::endl; GUDHI_CHECK(complex.num_vertices() == 0, std::invalid_argument("Cech_complex::create_complex - simplicial complex is not empty")); // insert the proximity graph in the simplicial complex - std::clog << "Hind: before insert_graph" << std::endl; complex.insert_graph(cech_skeleton_graph_); // expand the graph until dimension dim_max - std::clog << "Hind: before expansion_with_blockers" << std::endl; - complex.expansion_with_blockers(dim_max, - Cech_blocker(&complex, this)); + complex.expansion_with_blockers(dim_max, cech_blocker(&complex, this)); } /** @return max_radius value given at construction. */ @@ -111,12 +105,18 @@ class Cech_complex { /** @param[in] vertex Point position in the range. * @return The point. */ - const Point& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; } + const Point_d& get_point(Vertex_handle vertex) const { return point_cloud_[vertex]; } + + /** + * @return Vector of cached spheres. + */ + std::vector & get_cache() { return cache_; } private: Proximity_graph cech_skeleton_graph_; Filtration_value max_radius_; Point_cloud point_cloud_; + std::vector cache_; }; } // namespace cech_complex diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 3cac9ee2..5edd005d 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -11,20 +11,11 @@ #ifndef CECH_COMPLEX_BLOCKER_H_ #define CECH_COMPLEX_BLOCKER_H_ -// TODO to remove -#include // for Gudhi::Minimal_enclosing_ball_radius - -#include // for CGAL::to_double -#include // -// #include -#include // For EXACT or SAFE version -#include +#include // for casting from FT to double #include #include #include -#include -#include #include // for std::sqrt namespace Gudhi { @@ -48,7 +39,6 @@ namespace cech_complex { template class Cech_blocker { private: -// using Point_cloud = typename Cech_complex::Point_cloud; using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; using Filtration_value = typename SimplicialComplexForCech::Filtration_value; @@ -61,28 +51,18 @@ class Cech_blocker { // Sphere is a pair of point and squared radius. using Sphere = typename std::pair; - - /** \internal \brief TODO - * \param[in] - * \return */ template FT get_squared_radius(PointIterator begin, PointIterator end) const { return kernel_.compute_squared_radius_d_object()(begin, end); } - - /** \internal \brief TODO - * \param[in] - * \return */ + template Sphere get_sphere(PointIterator begin, PointIterator end) const { Point_d c = kernel_.construct_circumcenter_d_object()(begin, end); FT r = kernel_.squared_distance_d_object()(c, *begin); return std::make_pair(std::move(c), std::move(r)); } - - /** \internal \brief TODO - * \param[in] - * \return */ + template class CompareSpheresRadii { @@ -93,7 +73,7 @@ class Cech_blocker { return cast_to_double(firstSphere.second) < cast_to_double(secondSphere.second); } }; - + /** \internal \brief Čech complex blocker operator() - the oracle - assigns the filtration value from the simplex * radius and returns if the simplex expansion must be blocked. * \param[in] sh The Simplex_handle. @@ -102,104 +82,54 @@ class Cech_blocker { using Point_cloud = std::vector; CGAL::NT_converter cast_to_double; Filtration_value radius = 0.; -// std::string key_to_permute; -// std::vector faces_keys; - + // for each face of simplex sh, test outsider point is indeed inside enclosing ball, if yes, take it and exit loop, otherwise, new sphere is circumsphere of all vertices - // std::set enclosing_ball_radii; std::set > enclosing_ball_spheres; for (auto face : sc_ptr_->boundary_simplex_range(sh)) { - // Find which vertex of sh is missing in face. We rely on the fact that simplex_vertex_range is sorted. auto longlist = sc_ptr_->simplex_vertex_range(sh); auto shortlist = sc_ptr_->simplex_vertex_range(face); -// std::clog << "Hind debug: within FACE loop "<< std::endl; - // TODO to remove -// for (auto i = std::begin(longlist); i != std::end(longlist);++i) -// std::clog << "Hind debug: longlist: " << cc_ptr_->get_point(*i) << std::endl; -// for (auto i = std::begin(shortlist); i != std::end(shortlist);++i) -// std::clog << "Hind debug: shortlist: " << cc_ptr_->get_point(*i) << std::endl; - auto longiter = std::begin(longlist); auto shortiter = std::begin(shortlist); auto enditer = std::end(shortlist); while(shortiter != enditer && *longiter == *shortiter) { ++longiter; ++shortiter; } auto extra = *longiter; // Vertex_handle -// std::clog << "Hind debug: extra vertex: " << cc_ptr_->get_point(extra) << std::endl; - Point_cloud face_points; -// std::string key, key_extra; for (auto vertex : sc_ptr_->simplex_vertex_range(face)) { face_points.push_back(cc_ptr_->get_point(vertex)); -// key.append(std::to_string(vertex)); #ifdef DEBUG_TRACES std::clog << "#(" << vertex << ")#"; #endif // DEBUG_TRACES } -// key_extra = key; -// key_extra.append(std::to_string(extra)); -// faces_keys.push_back(key_extra); -// key_to_permute = key_extra; -// std::clog << "END OF VERTICES " << std::endl; -// std::clog << "KEY is: " << key << std::endl; -// std::clog << "KEY extra is: " << key_extra << std::endl; Sphere sph; - auto k = sc_ptr_->key(sh); - if(k != sc_ptr_->null_key()) - sph = cache_[k]; + auto face_sh = sc_ptr_->find(sc_ptr_->simplex_vertex_range(face)); + auto k = sc_ptr_->key(face_sh); + if(k != sc_ptr_->null_key()) { + sph = cc_ptr_->get_cache().at(k); + } else { sph = get_sphere(face_points.cbegin(), face_points.cend()); } -// auto it = cache_.find(key); -// if(it != cache_.end()) -// sph = it->second; -// else { -// sph = get_sphere(face_points.cbegin(), face_points.cend()); -// } if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) <= sph.second) { radius = std::sqrt(cast_to_double(sph.second)); #ifdef DEBUG_TRACES std::clog << "circumcenter: " << sph.first << ", radius: " << radius << std::endl; #endif // DEBUG_TRACES -// std::clog << "distance FYI: " << kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " < " << cast_to_double(sph.second) << std::endl; - // enclosing_ball_radii.insert(radius); enclosing_ball_spheres.insert(sph); -// cache_[key_extra] = sph; -// sc_ptr_->assign_key(sh, cache_.size()); -// cache_.push_back(sph); } -// else {// TODO to remove -// std::clog << "vertex not included BECAUSE DISTANCE: "<< kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) << " AND RAD SPHERE: " << sph.second << std::endl; -// } } // Get the minimal radius of all faces enclosing balls if exists if (!enclosing_ball_spheres.empty()) { - // radius = *enclosing_ball_radii.begin(); Sphere sph_min = *enclosing_ball_spheres.begin(); radius = std::sqrt(cast_to_double(sph_min.second)); - // std::clog << "CHECK that radius of min sphere is min radius: " << std::sqrt(cast_to_double(sph_min.second)) << "; and RADIUS min: " << radius << std::endl; - // Set all key_to_permute permutations to min sphere in cache -// do -// { -// if (cache_.find(key_to_permute) != cache_.end()) { -// if (cast_to_double(cache_[key_to_permute].second) > cast_to_double(sph_min.second)) -// cache_[key_to_permute] = sph_min; -// } -// else { -// cache_[key_to_permute] = sph_min; -// } -// } while(std::next_permutation(key_to_permute.begin(), key_to_permute.end())); -// for (auto k : faces_keys) { -// cache_[k] = sph_min; -// } - sc_ptr_->assign_key(sh, cache_.size()); - cache_.push_back(sph_min); + + sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); + cc_ptr_->get_cache().push_back(sph_min); } -// std::clog << "END OF FACES ; radius = " << radius << std::endl; - + if (radius == 0.) { // Spheres of each face don't contain the whole simplex Point_cloud points; for (auto vertex : sc_ptr_->simplex_vertex_range(sh)) { @@ -207,25 +137,10 @@ class Cech_blocker { } Sphere sph = get_sphere(points.cbegin(), points.cend()); radius = std::sqrt(cast_to_double(sph.second)); -// std::clog << "GLOBAL SPHERE radius = " << radius << std::endl; - // Set all key_to_permute permutations to sphere in cache -// do -// { -// // if (cache_.find(key_to_permute) != cache_.end()) { -// // if (cast_to_double(cache_[key_to_permute].second) > cast_to_double(sph.second)) -// // cache_[key_to_permute] = sph; -// // } -// // else { -// // cache_[key_to_permute] = sph; -// // } -// } while(std::next_permutation(key_to_permute.begin(), key_to_permute.end())); -// for (auto k : faces_keys) { -// cache_[k] = sph; -// } -// sc_ptr_->assign_key(sh, cache_.size()); -// cache_.push_back(sph); - } + sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); + cc_ptr_->get_cache().push_back(sph); + } #ifdef DEBUG_TRACES if (radius > cc_ptr_->max_radius()) std::clog << "radius > max_radius => expansion is blocked\n"; @@ -241,8 +156,6 @@ class Cech_blocker { SimplicialComplexForCech* sc_ptr_; Cech_complex* cc_ptr_; Kernel kernel_; -// std::map cache_; - std::vector cache_; }; } // namespace cech_complex -- cgit v1.2.3