From 91b0ff839b8058d3f5767e6ed80b93c23be2c98a Mon Sep 17 00:00:00 2001 From: Hind-M Date: Mon, 9 Aug 2021 16:21:24 +0200 Subject: First version of cech enhancement --- src/Cech_complex/benchmark/cech_complex_benchmark.cpp | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) (limited to 'src/Cech_complex/benchmark/cech_complex_benchmark.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index e489e8a4..c332c656 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -17,6 +17,8 @@ #include #include +#include // For EXACT or SAFE version + #include "boost/filesystem.hpp" // includes all needed Boost.Filesystem declarations #include @@ -30,7 +32,8 @@ using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; using Proximity_graph = Gudhi::Proximity_graph; using Rips_complex = Gudhi::rips_complex::Rips_complex; -using Cech_complex = Gudhi::cech_complex::Cech_complex; +using Kernel = CGAL::Epeck_d; +using Cech_complex = Gudhi::cech_complex::Cech_complex; class Minimal_enclosing_ball_radius { public: -- cgit v1.2.3 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 --- .../benchmark/cech_complex_benchmark.cpp | 2 +- .../example/cech_complex_example_from_points.cpp | 42 +------ .../example/cech_complex_step_by_step.cpp | 8 -- src/Cech_complex/include/gudhi/Cech_complex.h | 38 +++---- .../include/gudhi/Cech_complex_blocker.h | 121 +++------------------ src/Cech_complex/test/test_cech_complex.cpp | 37 ++----- src/Cech_complex/utilities/cech_persistence.cpp | 4 +- src/Simplex_tree/include/gudhi/Simplex_tree.h | 4 - src/common/include/gudhi/distance_functions.h | 20 +--- .../include/gudhi/graph_simplicial_complex.h | 3 - 10 files changed, 54 insertions(+), 225 deletions(-) (limited to 'src/Cech_complex/benchmark/cech_complex_benchmark.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index c332c656..4a1aa06e 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -33,7 +33,7 @@ using Points_off_reader = Gudhi::Points_off_reader; using Proximity_graph = Gudhi::Proximity_graph; using Rips_complex = Gudhi::rips_complex::Rips_complex; using Kernel = CGAL::Epeck_d; -using Cech_complex = Gudhi::cech_complex::Cech_complex; +using Cech_complex = Gudhi::cech_complex::Cech_complex; class Minimal_enclosing_ball_radius { public: diff --git a/src/Cech_complex/example/cech_complex_example_from_points.cpp b/src/Cech_complex/example/cech_complex_example_from_points.cpp index e78ad51d..78861951 100644 --- a/src/Cech_complex/example/cech_complex_example_from_points.cpp +++ b/src/Cech_complex/example/cech_complex_example_from_points.cpp @@ -10,25 +10,15 @@ int main() { // Type definitions -// using Point_cloud = std::vector>; using Simplex_tree = Gudhi::Simplex_tree; using Filtration_value = Simplex_tree::Filtration_value; using Kernel = CGAL::Epeck_d; using FT = typename Kernel::FT; using Point = typename Kernel::Point_d; using Point_cloud = std::vector; - using Cech_complex = Gudhi::cech_complex::Cech_complex; + using Cech_complex = Gudhi::cech_complex::Cech_complex; Point_cloud points; -// points.push_back({1., 0.}); // 0 -// points.push_back({0., 1.}); // 1 -// points.push_back({2., 1.}); // 2 -// points.push_back({3., 2.}); // 3 -// points.push_back({0., 3.}); // 4 -// points.push_back({3. + std::sqrt(3.), 3.}); // 5 - -// std::vector point({0.0, 0.0, 0.0, 0.0}); -// points.emplace_back(point.begin(), point.end()); std::vector point0({1., 0.}); points.emplace_back(point0.begin(), point0.end()); @@ -42,8 +32,6 @@ int main() { points.emplace_back(point4.begin(), point4.end()); std::vector point5({3. + std::sqrt(3.), 3.}); points.emplace_back(point5.begin(), point5.end()); - -/* std::vector point6({1., 4.}); points.emplace_back(point6.begin(), point6.end()); std::vector point7({3., 4.}); @@ -54,39 +42,15 @@ int main() { points.emplace_back(point9.begin(), point9.end()); std::vector point10({-0.5, 2.}); points.emplace_back(point10.begin(), point10.end()); -*/ -// points.emplace_back(Point(std::vector({1., 0.}))); -// points.emplace_back(Point(std::vector({0., 1.}))); -// points.emplace_back(Point(std::vector({2., 1.}))); -// points.emplace_back(Point(std::vector({3., 2.}))); -// points.emplace_back(Point(std::vector({0., 3.}))); -// points.emplace_back(Point(std::vector({3. + std::sqrt(3.), 3.}))); - - -// points.push_back(Point(1.0, 0.0)); -// points.push_back(Point(0.0, 1.0)); -// points.push_back(Point(2.0, 1.0)); -// points.push_back(Point(3.0, 2.0)); -// points.push_back(Point(0.0, 3.0)); -// points.push_back(Point(3.0 + std::sqrt(3.0), 3.0)); - - -// points.push_back({1., 4.}); // 6 -// points.push_back({3., 4.}); // 7 -// points.push_back({2., 4. + std::sqrt(3.)}); // 8 -// points.push_back({0., 4.}); // 9 -// points.push_back({-0.5, 2.}); // 10 // ---------------------------------------------------------------------------- // Init of a Cech complex from points // ---------------------------------------------------------------------------- - Filtration_value max_radius = 10.; //100.; - std::clog << "Hind: Just before the Cech constructor" << std::endl; + Filtration_value max_radius = 100.; //100.; Cech_complex cech_complex_from_points(points, max_radius); - std::clog << "Hind: Just after the Cech constructor" << std::endl; Simplex_tree stree; - cech_complex_from_points.create_complex(stree, 3); //6 + cech_complex_from_points.create_complex(stree, 6); //6 // ---------------------------------------------------------------------------- // Display information about the one skeleton Cech complex // ---------------------------------------------------------------------------- diff --git a/src/Cech_complex/example/cech_complex_step_by_step.cpp b/src/Cech_complex/example/cech_complex_step_by_step.cpp index ac08e6cc..44e7f945 100644 --- a/src/Cech_complex/example/cech_complex_step_by_step.cpp +++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp @@ -13,8 +13,6 @@ #include #include -#include // TODO to remove ? - #include #include @@ -36,7 +34,6 @@ using Simplex_tree = Gudhi::Simplex_tree<>; using Simplex_handle = Simplex_tree::Simplex_handle; using Filtration_value = Simplex_tree::Filtration_value; -// using Point = std::vector; using Kernel = CGAL::Epeck_d; using Point = typename Kernel::Point_d; using Points_off_reader = Gudhi::Points_off_reader; @@ -45,9 +42,6 @@ using Proximity_graph = Gudhi::Proximity_graph; class Cech_blocker { private: using Point_cloud = std::vector; -// using Point_iterator = Point_cloud::const_iterator; -// using Coordinate_iterator = Point::const_iterator; -// using Min_sphere = Gudhi::Miniball::Miniball>; public: bool operator()(Simplex_handle sh) { @@ -67,14 +61,12 @@ class Cech_blocker { } Cech_blocker(Simplex_tree& simplex_tree, Filtration_value max_radius, const std::vector& point_cloud) : simplex_tree_(simplex_tree), max_radius_(max_radius), point_cloud_(point_cloud) { -// dimension_ = point_cloud_[0].size(); } private: Simplex_tree simplex_tree_; Filtration_value max_radius_; std::vector point_cloud_; -// int dimension_; }; void program_options(int argc, char* argv[], std::string& off_file_points, Filtration_value& max_radius, int& dim_max); 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 diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 81efd6ae..51b466da 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -24,26 +24,19 @@ #include #include #include -#include #include // For EXACT or SAFE version - // Type definitions using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; -//using Point = std::vector; using Kernel = CGAL::Epeck_d; using FT = typename Kernel::FT; using Point = typename Kernel::Point_d; using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; -using Cech_complex = Gudhi::cech_complex::Cech_complex; - -// using Point_iterator = Point_cloud::const_iterator; -// using Coordinate_iterator = Point::const_iterator; -// using Min_sphere = Gudhi::Miniball::Miniball>; +using Cech_complex = Gudhi::cech_complex::Cech_complex; BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { // ---------------------------------------------------------------------------- @@ -52,18 +45,6 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { // // ---------------------------------------------------------------------------- Point_cloud points; -// points.push_back({1., 0.}); // 0 -// points.push_back({0., 1.}); // 1 -// points.push_back({2., 1.}); // 2 -// points.push_back({3., 2.}); // 3 -// points.push_back({0., 3.}); // 4 -// points.push_back({3. + std::sqrt(3.), 3.}); // 5 -// points.push_back({1., 4.}); // 6 -// points.push_back({3., 4.}); // 7 -// points.push_back({2., 4. + std::sqrt(3.)}); // 8 -// points.push_back({0., 4.}); // 9 -// points.push_back({-0.5, 2.}); // 10 - std::vector point0({1., 0.}); points.emplace_back(point0.begin(), point0.end()); @@ -156,36 +137,32 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { for (std::size_t vertex = 0; vertex <= 2; vertex++) { points012.push_back(cech_complex_for_doc.get_point(vertex)); } -// std::size_t dimension = points[0].end() - points[0].begin(); -// Min_sphere ms012(dimension, points012.begin(), points012.end()); Kernel kern; Simplex_tree::Filtration_value f012 = st2.filtration(st2.find({0, 1, 2})); - CGAL::NT_converter cast_to_double; - std::clog << "f012= " << f012 << " | points012_radius= " << std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points012.begin(), points012.end()))) << std::endl; - + std::clog << "f012= " << f012 << std::endl; + CGAL::NT_converter cast_to_double; GUDHI_TEST_FLOAT_EQUALITY_CHECK(f012, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points012.begin(), points012.end())))); Point_cloud points1410; points1410.push_back(cech_complex_for_doc.get_point(1)); points1410.push_back(cech_complex_for_doc.get_point(4)); points1410.push_back(cech_complex_for_doc.get_point(10)); -// Min_sphere ms1410(dimension, points1410.begin(), points1410.end()); Simplex_tree::Filtration_value f1410 = st2.filtration(st2.find({1, 4, 10})); - std::clog << "f1410= " << f1410 << " | points1410_radius= " << std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points1410.begin(), points1410.end()))) << std::endl; + std::clog << "f1410= " << f1410 << std::endl; - GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points1410.begin(), points1410.end())))); + // In this case, the computed sphere using CGAL kernel does not match the minimal enclosing ball; the filtration value check is therefore done against a hardcoded value + GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, 1.); Point_cloud points469; points469.push_back(cech_complex_for_doc.get_point(4)); points469.push_back(cech_complex_for_doc.get_point(6)); points469.push_back(cech_complex_for_doc.get_point(9)); -// Min_sphere ms469(dimension, points469.begin(), points469.end()); Simplex_tree::Filtration_value f469 = st2.filtration(st2.find({4, 6, 9})); - std::clog << "f469= " << f469 << " | points469_radius= " << std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points469.begin(), points469.end()))) << std::endl; + std::clog << "f469= " << f469 << std::endl; GUDHI_TEST_FLOAT_EQUALITY_CHECK(f469, std::sqrt(cast_to_double(kern.compute_squared_radius_d_object()(points469.begin(), points469.end())))); diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp index ccd7d453..0c945cad 100644 --- a/src/Cech_complex/utilities/cech_persistence.cpp +++ b/src/Cech_complex/utilities/cech_persistence.cpp @@ -25,14 +25,12 @@ // Types definition using Simplex_tree = Gudhi::Simplex_tree; using Filtration_value = Simplex_tree::Filtration_value; -// using Point = std::vector; -// using Point_cloud = std::vector; using Kernel = CGAL::Epeck_d; using Point = typename Kernel::Point_d; using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; -using Cech_complex = Gudhi::cech_complex::Cech_complex; +using Cech_complex = Gudhi::cech_complex::Cech_complex; using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; diff --git a/src/Simplex_tree/include/gudhi/Simplex_tree.h b/src/Simplex_tree/include/gudhi/Simplex_tree.h index f69ed6ec..85790baf 100644 --- a/src/Simplex_tree/include/gudhi/Simplex_tree.h +++ b/src/Simplex_tree/include/gudhi/Simplex_tree.h @@ -1278,10 +1278,8 @@ class Simplex_tree { intersection.emplace_back(next->first, Node(nullptr, filt)); } } - std::clog << "Hind: after intersection insertion" << std::endl; if (intersection.size() != 0) { // Reverse the order to insert - std::clog << "Hind: declare new siblings" << std::endl; Siblings * new_sib = new Siblings(siblings, // oncles simplex->first, // parent boost::adaptors::reverse(intersection)); // boost::container::ordered_unique_range_t @@ -1290,12 +1288,10 @@ class Simplex_tree { for (auto new_sib_member = new_sib->members().begin(); new_sib_member != new_sib->members().end(); new_sib_member++) { - std::clog << "Hind: check the blocker result" << std::endl; bool blocker_result = block_simplex(new_sib_member); // new_sib member has been blocked by the blocker function // add it to the list to be removed - do not perform it while looping on it if (blocker_result) { - std::clog << "Hind: add to list of blocked sib to be removed" << std::endl; blocked_new_sib_vertex_list.push_back(new_sib_member->first); } } diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index ae5168aa..a8ee4a75 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -65,19 +65,17 @@ class Euclidean_distance { * The points are assumed to have the same dimension. */ class Minimal_enclosing_ball_radius { public: - /** \brief TODO + /** \brief Enclosing ball radius from two points using CGAL. * * @param[in] point_1 * @param[in] point_2 - * @return - * \tparam Point + * @return Enclosing ball radius for the two points. + * \tparam Point must be a Kernel::Point_d from CGAL. * */ - //typename FT = typename Kernel::FT, template< typename Kernel = CGAL::Epeck_d, typename Point= typename Kernel::Point_d> double operator()(const Point& point_1, const Point& point_2) const { - std::clog << "Added template: distance betw points 1 and 2" << std::endl; Kernel kernel_; return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; } @@ -94,25 +92,21 @@ class Minimal_enclosing_ball_radius { template< typename Point > typename std::iterator_traits::type>::value_type operator()(const Point& point_1, const Point& point_2) const { - std::clog << "Hind: Minimal_enclosing_ball_radius point1 et 2; Euclidean" << std::endl; std::clog << "#" << *point_1.begin() << "##" << *point_2.begin() << std::endl; return Euclidean_distance()(point_1, point_2) / 2.; } - - /** \brief TODO + /** \brief Enclosing ball radius from a point cloud using CGAL. * * @param[in] point_cloud The points. - * @return - * \tparam Point_cloud + * @return Enclosing ball radius for the points. + * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. * */ - //typename FT = typename Kernel::FT, template< typename Kernel = CGAL::Epeck_d, typename Point= typename Kernel::Point_d, typename Point_cloud = std::vector> double operator()(const Point_cloud& point_cloud) const { - std::clog << "Added template: distance in point cloud" << std::endl; Kernel kernel_; return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); } @@ -133,8 +127,6 @@ class Minimal_enclosing_ball_radius { typename Coordinate = typename std::iterator_traits::value_type> Coordinate operator()(const Point_cloud& point_cloud) const { - std::clog << "Hind: Minimal_enclosing_ball_radius point cloud; Miniball" << std::endl; - using Min_sphere = Miniball::Miniball>; Min_sphere ms(boost::size(*point_cloud.begin()), point_cloud.begin(), point_cloud.end()); diff --git a/src/common/include/gudhi/graph_simplicial_complex.h b/src/common/include/gudhi/graph_simplicial_complex.h index 9190182c..da9dee7d 100644 --- a/src/common/include/gudhi/graph_simplicial_complex.h +++ b/src/common/include/gudhi/graph_simplicial_complex.h @@ -18,8 +18,6 @@ #include #include // for std::tie -#include - namespace Gudhi { /** @file * @brief Graph simplicial complex methods @@ -78,7 +76,6 @@ Proximity_graph compute_proximity_graph( for (auto it_u = points.begin(); it_u != points.end(); ++it_u) { idx_v = idx_u + 1; for (auto it_v = it_u + 1; it_v != points.end(); ++it_v, ++idx_v) { - std::clog << "#idx_u" << idx_u << "#idx_v " << idx_v << std::endl; fil = distance(*it_u, *it_v); if (fil <= threshold) { edges.emplace_back(idx_u, idx_v); -- cgit v1.2.3 From d3970dbbc16993d348092899eb8fcd1ea1aceb8d Mon Sep 17 00:00:00 2001 From: Hind-M Date: Fri, 1 Oct 2021 15:47:58 +0200 Subject: Separate Minimal_enclosing_ball_radius class which uses CGAL from common distance_functions.h file --- .../benchmark/cech_complex_benchmark.cpp | 2 +- .../example/cech_complex_step_by_step.cpp | 2 +- src/Cech_complex/include/gudhi/Cech_complex.h | 2 +- .../include/gudhi/Cech_complex/Cech_kernel.h | 149 +++++++++++++++++++++ src/Cech_complex/test/test_cech_complex.cpp | 2 +- src/Cech_complex/utilities/cech_persistence.cpp | 2 +- src/common/include/gudhi/distance_functions.h | 83 ------------ 7 files changed, 154 insertions(+), 88 deletions(-) create mode 100644 src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h (limited to 'src/Cech_complex/benchmark/cech_complex_benchmark.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 4a1aa06e..cfeb0725 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -9,7 +9,7 @@ */ #include -#include +#include #include #include #include diff --git a/src/Cech_complex/example/cech_complex_step_by_step.cpp b/src/Cech_complex/example/cech_complex_step_by_step.cpp index 44e7f945..2d8321b1 100644 --- a/src/Cech_complex/example/cech_complex_step_by_step.cpp +++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp @@ -9,7 +9,7 @@ */ #include -#include +#include #include #include diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 32a78aec..7bbf97d1 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -11,7 +11,7 @@ #ifndef CECH_COMPLEX_H_ #define CECH_COMPLEX_H_ -#include // for Gudhi::Minimal_enclosing_ball_radius +#include // for Gudhi::Minimal_enclosing_ball_radius #include // for Gudhi::Proximity_graph #include // for GUDHI_CHECK #include // for Gudhi::cech_complex::Cech_blocker diff --git a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h new file mode 100644 index 00000000..93af90d2 --- /dev/null +++ b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h @@ -0,0 +1,149 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Hind Montassif + * + * Copyright (C) 2021 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef CECH_KERNEL_H_ +#define CECH_KERNEL_H_ + +#include + +#include // for std::sqrt + +namespace Gudhi { + +// namespace cech_complex { + +/** @brief Compute the radius of the minimal enclosing ball between Points given by a range of coordinates. + * The points are assumed to have the same dimension. */ +class Minimal_enclosing_ball_radius { + public: + /** \brief Enclosing ball radius from two points using CGAL. + * + * @param[in] point_1 + * @param[in] point_2 + * @return Enclosing ball radius for the two points. + * \tparam Point must be a Kernel::Point_d from CGAL. + * + */ + template< typename Kernel = CGAL::Epeck_d, + typename Point= typename Kernel::Point_d> + double operator()(const Point& point_1, const Point& point_2) const { + Kernel kernel_; + return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; + } + + /** \brief Minimal_enclosing_ball_radius from two points. + * + * @param[in] point_1 First point. + * @param[in] point_2 second point. + * @return The minimal enclosing ball radius for the two points (aka. Euclidean distance / 2.). + * + * \tparam Point must be a range of Cartesian coordinates. + * + */ +// template< typename Point > +// typename std::iterator_traits::type>::value_type +// operator()(const Point& point_1, const Point& point_2) const { +// std::clog << "#" << *point_1.begin() << "##" << *point_2.begin() << std::endl; +// return Euclidean_distance()(point_1, point_2) / 2.; +// } + + /** \brief Enclosing ball radius from a point cloud using CGAL. + * + * @param[in] point_cloud The points. + * @return Enclosing ball radius for the points. + * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. + * + */ + template< typename Kernel = CGAL::Epeck_d, + typename Point= typename Kernel::Point_d, + typename Point_cloud = std::vector> + double operator()(const Point_cloud& point_cloud) const { + Kernel kernel_; + return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); + } + + /** \brief Minimal_enclosing_ball_radius from a point cloud. + * + * @param[in] point_cloud The points. + * @return The minimal enclosing ball radius for the points. + * + * \tparam Point_cloud must be a range of points with Cartesian coordinates. + * Point_cloud is a range over a range of Coordinate. + * + */ +// template< typename Point_cloud, +// typename Point_iterator = typename boost::range_const_iterator::type, +// typename Point = typename std::iterator_traits::value_type, +// typename Coordinate_iterator = typename boost::range_const_iterator::type, +// typename Coordinate = typename std::iterator_traits::value_type> +// Coordinate +// operator()(const Point_cloud& point_cloud) const { +// using Min_sphere = Miniball::Miniball>; +// +// Min_sphere ms(boost::size(*point_cloud.begin()), point_cloud.begin(), point_cloud.end()); +// #ifdef DEBUG_TRACES +// std::clog << "Minimal_enclosing_ball_radius = " << std::sqrt(ms.squared_radius()) << " | nb points = " +// << boost::size(point_cloud) << " | dimension = " +// << boost::size(*point_cloud.begin()) << std::endl; +// #endif // DEBUG_TRACES +// +// return std::sqrt(ms.squared_radius()); +// } +}; + +/** + * \class Cech_kernel + * \brief Cech complex kernel container. + * + * \details + * The Cech complex kernel container stores CGAL Kernel and dispatch basic computations. + */ + +// template < typename Kernel > +// class Cech_kernel { +// private: +// // Kernel for functions access. +// Kernel kernel_; +// public: +// using Point_d = typename Kernel::Point_d; +// // Numeric type of coordinates in the kernel +// using FT = typename Kernel::FT; +// // Sphere is a pair of point and squared radius. +// using Sphere = typename std::pair; +// +// int get_dimension(const Point_d& p0) const { +// return kernel_.point_dimension_d_object()(p0); +// } +// +// 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)); +// } +// +// template +// FT get_squared_radius(PointIterator begin, PointIterator end) const { +// return kernel_.compute_squared_radius_d_object()(begin, end); +// } +// +// FT get_squared_radius(const Sphere& sph) const { +// return sph.second; +// } +// }; + + +//} // namespace cech_complex + +// namespace cechcomplex = cech_complex; + +} // namespace Gudhi + +#endif // CECH_KERNEL_H_ diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 51b466da..7d8c3c22 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -22,7 +22,7 @@ // to construct Cech_complex from a OFF file of points #include #include -#include +#include #include #include // For EXACT or SAFE version diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp index 0c945cad..ccf63e3e 100644 --- a/src/Cech_complex/utilities/cech_persistence.cpp +++ b/src/Cech_complex/utilities/cech_persistence.cpp @@ -9,7 +9,7 @@ */ #include -#include +#include #include #include #include diff --git a/src/common/include/gudhi/distance_functions.h b/src/common/include/gudhi/distance_functions.h index a8ee4a75..5e5a1e31 100644 --- a/src/common/include/gudhi/distance_functions.h +++ b/src/common/include/gudhi/distance_functions.h @@ -13,13 +13,9 @@ #include -#include - #include #include -#include - #include // for std::sqrt #include // for std::decay #include // for std::begin, std::end @@ -61,85 +57,6 @@ class Euclidean_distance { } }; -/** @brief Compute the radius of the minimal enclosing ball between Points given by a range of coordinates. - * The points are assumed to have the same dimension. */ -class Minimal_enclosing_ball_radius { - public: - /** \brief Enclosing ball radius from two points using CGAL. - * - * @param[in] point_1 - * @param[in] point_2 - * @return Enclosing ball radius for the two points. - * \tparam Point must be a Kernel::Point_d from CGAL. - * - */ - template< typename Kernel = CGAL::Epeck_d, - typename Point= typename Kernel::Point_d> - double operator()(const Point& point_1, const Point& point_2) const { - Kernel kernel_; - return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; - } - - /** \brief Minimal_enclosing_ball_radius from two points. - * - * @param[in] point_1 First point. - * @param[in] point_2 second point. - * @return The minimal enclosing ball radius for the two points (aka. Euclidean distance / 2.). - * - * \tparam Point must be a range of Cartesian coordinates. - * - */ - template< typename Point > - typename std::iterator_traits::type>::value_type - operator()(const Point& point_1, const Point& point_2) const { - std::clog << "#" << *point_1.begin() << "##" << *point_2.begin() << std::endl; - return Euclidean_distance()(point_1, point_2) / 2.; - } - - /** \brief Enclosing ball radius from a point cloud using CGAL. - * - * @param[in] point_cloud The points. - * @return Enclosing ball radius for the points. - * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. - * - */ - template< typename Kernel = CGAL::Epeck_d, - typename Point= typename Kernel::Point_d, - typename Point_cloud = std::vector> - double operator()(const Point_cloud& point_cloud) const { - Kernel kernel_; - return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); - } - - /** \brief Minimal_enclosing_ball_radius from a point cloud. - * - * @param[in] point_cloud The points. - * @return The minimal enclosing ball radius for the points. - * - * \tparam Point_cloud must be a range of points with Cartesian coordinates. - * Point_cloud is a range over a range of Coordinate. - * - */ - template< typename Point_cloud, - typename Point_iterator = typename boost::range_const_iterator::type, - typename Point = typename std::iterator_traits::value_type, - typename Coordinate_iterator = typename boost::range_const_iterator::type, - typename Coordinate = typename std::iterator_traits::value_type> - Coordinate - operator()(const Point_cloud& point_cloud) const { - using Min_sphere = Miniball::Miniball>; - - Min_sphere ms(boost::size(*point_cloud.begin()), point_cloud.begin(), point_cloud.end()); -#ifdef DEBUG_TRACES - std::clog << "Minimal_enclosing_ball_radius = " << std::sqrt(ms.squared_radius()) << " | nb points = " - << boost::size(point_cloud) << " | dimension = " - << boost::size(*point_cloud.begin()) << std::endl; -#endif // DEBUG_TRACES - - return std::sqrt(ms.squared_radius()); - } -}; - } // namespace Gudhi #endif // DISTANCE_FUNCTIONS_H_ -- cgit v1.2.3 From 7d3d9e57c7c72b0762db910c2638b08e596199df Mon Sep 17 00:00:00 2001 From: Hind-M Date: Tue, 19 Oct 2021 16:16:21 +0200 Subject: Make Cech benchmark work --- .../benchmark/cech_complex_benchmark.cpp | 24 ++++++++++++++-------- 1 file changed, 15 insertions(+), 9 deletions(-) (limited to 'src/Cech_complex/benchmark/cech_complex_benchmark.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index cfeb0725..06d90757 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -9,7 +9,7 @@ */ #include -#include +#include #include #include #include @@ -33,7 +33,10 @@ using Points_off_reader = Gudhi::Points_off_reader; using Proximity_graph = Gudhi::Proximity_graph; using Rips_complex = Gudhi::rips_complex::Rips_complex; using Kernel = CGAL::Epeck_d; -using Cech_complex = Gudhi::cech_complex::Cech_complex; +using Point_cgal = typename Kernel::Point_d; +using Point_cloud_cgal = std::vector; +using Points_off_reader_cgal = Gudhi::Points_off_reader; +using Cech_complex = Gudhi::cech_complex::Cech_complex; class Minimal_enclosing_ball_radius { public: @@ -65,6 +68,7 @@ int main(int argc, char* argv[]) { // Extract the points from the file filepoints Points_off_reader off_reader(off_file_points); + Points_off_reader_cgal off_reader_cgal(off_file_points); Gudhi::Clock euclidean_clock("Gudhi::Euclidean_distance"); // Compute the proximity graph of the points @@ -79,16 +83,16 @@ int main(int argc, char* argv[]) { off_reader.get_point_cloud(), threshold, Minimal_enclosing_ball_radius()); std::clog << miniball_clock << std::endl; - Gudhi::Clock common_miniball_clock("Gudhi::Minimal_enclosing_ball_radius()"); + Gudhi::Clock cgal_miniball_clock("Gudhi::Minimal_enclosing_ball_radius_cgal()"); // Compute the proximity graph of the points - Proximity_graph common_miniball_prox_graph = Gudhi::compute_proximity_graph( - off_reader.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius()); - std::clog << common_miniball_clock << std::endl; + Proximity_graph cgal_miniball_prox_graph = Gudhi::compute_proximity_graph( + off_reader_cgal.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius()); + std::clog << cgal_miniball_clock << std::endl; boost::filesystem::path full_path(boost::filesystem::current_path()); std::clog << "Current path is : " << full_path << std::endl; - std::clog << "File name;Radius;Rips time;Cech time; Ratio Rips/Cech time;Rips nb simplices;Cech nb simplices;" + std::clog << "File name; Radius; Rips time; Cech time; Ratio Rips/Cech time; Rips nb simplices; Cech nb simplices;" << std::endl; boost::filesystem::directory_iterator end_itr; // default construction yields past-the-end for (boost::filesystem::directory_iterator itr(boost::filesystem::current_path()); itr != end_itr; ++itr) { @@ -96,13 +100,15 @@ int main(int argc, char* argv[]) { if (itr->path().extension() == ".off") // see below { Points_off_reader off_reader(itr->path().string()); + Points_off_reader_cgal off_reader_cgal(itr->path().string()); + Point p0 = off_reader.get_point_cloud()[0]; for (Filtration_value radius = 0.1; radius < 0.4; radius += 0.1) { std::clog << itr->path().stem() << ";"; std::clog << radius << ";"; Gudhi::Clock rips_clock("Rips computation"); - Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), radius, + Rips_complex rips_complex_from_points(off_reader_cgal.get_point_cloud(), radius, Gudhi::Minimal_enclosing_ball_radius()); Simplex_tree rips_stree; rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); @@ -113,7 +119,7 @@ int main(int argc, char* argv[]) { std::clog << rips_sec << ";"; Gudhi::Clock cech_clock("Cech computation"); - Cech_complex cech_complex_from_points(off_reader.get_point_cloud(), radius); + Cech_complex cech_complex_from_points(off_reader_cgal.get_point_cloud(), radius); Simplex_tree cech_stree; cech_complex_from_points.create_complex(cech_stree, p0.size() - 1); // ------------------------------------------ -- cgit v1.2.3 From b1f40dd2c4397c1975533c54a54538160c727d55 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Thu, 6 Jan 2022 11:39:05 +0100 Subject: Make kernel a parameter of Minimal_enclosing_ball_radius class Use Epick in cech benchmark instead of Epeck --- src/Cech_complex/benchmark/cech_complex_benchmark.cpp | 8 ++++---- src/Cech_complex/include/gudhi/Cech_complex.h | 3 +-- .../include/gudhi/Cech_complex/Cech_kernel.h | 16 ++++++++-------- src/Cech_complex/test/test_cech_complex.cpp | 4 ++-- 4 files changed, 15 insertions(+), 16 deletions(-) (limited to 'src/Cech_complex/benchmark/cech_complex_benchmark.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 06d90757..e715b513 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -17,7 +17,7 @@ #include #include -#include // For EXACT or SAFE version +#include #include "boost/filesystem.hpp" // includes all needed Boost.Filesystem declarations @@ -32,7 +32,7 @@ using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; using Proximity_graph = Gudhi::Proximity_graph; using Rips_complex = Gudhi::rips_complex::Rips_complex; -using Kernel = CGAL::Epeck_d; +using Kernel = CGAL::Epick_d>; using Point_cgal = typename Kernel::Point_d; using Point_cloud_cgal = std::vector; using Points_off_reader_cgal = Gudhi::Points_off_reader; @@ -86,7 +86,7 @@ int main(int argc, char* argv[]) { Gudhi::Clock cgal_miniball_clock("Gudhi::Minimal_enclosing_ball_radius_cgal()"); // Compute the proximity graph of the points Proximity_graph cgal_miniball_prox_graph = Gudhi::compute_proximity_graph( - off_reader_cgal.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius()); + off_reader_cgal.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius()); std::clog << cgal_miniball_clock << std::endl; boost::filesystem::path full_path(boost::filesystem::current_path()); @@ -109,7 +109,7 @@ int main(int argc, char* argv[]) { std::clog << radius << ";"; Gudhi::Clock rips_clock("Rips computation"); Rips_complex rips_complex_from_points(off_reader_cgal.get_point_cloud(), radius, - Gudhi::Minimal_enclosing_ball_radius()); + Gudhi::Minimal_enclosing_ball_radius()); Simplex_tree rips_stree; rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); // ------------------------------------------ diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 7bbf97d1..0031d861 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -18,7 +18,6 @@ #include #include // for exception management -#include namespace Gudhi { @@ -78,7 +77,7 @@ class Cech_complex { point_cloud_.assign(points.begin(), points.end()); cech_skeleton_graph_ = Gudhi::compute_proximity_graph( - point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); + point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); } /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal diff --git a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h index 348bb57d..89012206 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h +++ b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h @@ -11,9 +11,10 @@ #ifndef CECH_KERNEL_H_ #define CECH_KERNEL_H_ -#include +#include // for #include #include // for std::sqrt +#include namespace Gudhi { @@ -21,8 +22,14 @@ namespace Gudhi { /** @brief Compute the radius of the minimal enclosing ball between Points given by a range of coordinates. * The points are assumed to have the same dimension. */ +template class Minimal_enclosing_ball_radius { + private: + Kernel kernel_; public: + using Point = typename Kernel::Point_d; + using Point_cloud = typename std::vector; + /** \brief Enclosing ball radius from two points using CGAL. * * @param[in] point_1 @@ -31,10 +38,7 @@ class Minimal_enclosing_ball_radius { * \tparam Point must be a Kernel::Point_d from CGAL. * */ - template< typename Kernel = CGAL::Epeck_d, - typename Point= typename Kernel::Point_d> double operator()(const Point& point_1, const Point& point_2) const { - Kernel kernel_; return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; } @@ -46,11 +50,7 @@ class Minimal_enclosing_ball_radius { * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. * */ - template< typename Kernel = CGAL::Epeck_d, - typename Point= typename Kernel::Point_d, - typename Point_cloud = std::vector> double operator()(const Point_cloud& point_cloud) const { - Kernel kernel_; return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); } diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index 7d8c3c22..ca7a9778 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -108,11 +108,11 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { std::clog << vertex << ","; vp.push_back(points.at(vertex)); } - std::clog << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1)) + std::clog << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1)) << " - filtration =" << st.filtration(f_simplex) << std::endl; BOOST_CHECK(vp.size() == 2); GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), - Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1))); + Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1))); } } -- cgit v1.2.3 From 2d1fb6b63f0ca0c7e027cc298fc16198a6283df1 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Fri, 11 Feb 2022 15:00:27 +0100 Subject: Do some code clean up/renaming --- .../benchmark/cech_complex_benchmark.cpp | 13 ++- .../example/cech_complex_example_from_points.cpp | 2 +- .../example/cech_complex_step_by_step.cpp | 6 +- src/Cech_complex/include/gudhi/Cech_complex.h | 40 +++----- .../include/gudhi/Cech_complex/Cech_kernel.h | 107 --------------------- .../include/gudhi/Cech_complex_blocker.h | 22 ++--- .../include/gudhi/sphere_circumradius.h | 62 ++++++++++++ src/Cech_complex/test/test_cech_complex.cpp | 10 +- src/Cech_complex/utilities/cech_persistence.cpp | 5 +- 9 files changed, 106 insertions(+), 161 deletions(-) delete mode 100644 src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h create mode 100644 src/Cech_complex/include/gudhi/sphere_circumradius.h (limited to 'src/Cech_complex/benchmark/cech_complex_benchmark.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 94c5fa4f..b283e1a8 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -34,9 +34,8 @@ using Proximity_graph = Gudhi::Proximity_graph; using Rips_complex = Gudhi::rips_complex::Rips_complex; using Kernel = CGAL::Epick_d>; using Point_cgal = typename Kernel::Point_d; -using Point_cloud_cgal = std::vector; using Points_off_reader_cgal = Gudhi::Points_off_reader; -using Cech_complex = Gudhi::cech_complex::Cech_complex; +using Cech_complex = Gudhi::cech_complex::Cech_complex; class Minimal_enclosing_ball_radius { public: @@ -83,11 +82,11 @@ int main(int argc, char* argv[]) { off_reader.get_point_cloud(), threshold, Minimal_enclosing_ball_radius()); std::clog << miniball_clock << std::endl; - Gudhi::Clock cgal_miniball_clock("Gudhi::Minimal_enclosing_ball_radius_cgal()"); + Gudhi::Clock cgal_circumsphere_clock("Gudhi::cech_complex::Sphere_circumradius_cgal()"); // Compute the proximity graph of the points - Proximity_graph cgal_miniball_prox_graph = Gudhi::compute_proximity_graph( - off_reader_cgal.get_point_cloud(), threshold, Gudhi::Minimal_enclosing_ball_radius()); - std::clog << cgal_miniball_clock << std::endl; + Proximity_graph cgal_circumsphere_prox_graph = Gudhi::compute_proximity_graph( + off_reader_cgal.get_point_cloud(), threshold, Gudhi::cech_complex::Sphere_circumradius()); + std::clog << cgal_circumsphere_clock << std::endl; boost::filesystem::path full_path(boost::filesystem::current_path()); std::clog << "Current path is : " << full_path << std::endl; @@ -109,7 +108,7 @@ int main(int argc, char* argv[]) { std::clog << radius << ";"; Gudhi::Clock rips_clock("Rips computation"); Rips_complex rips_complex_from_points(off_reader_cgal.get_point_cloud(), radius, - Gudhi::Minimal_enclosing_ball_radius()); + Gudhi::cech_complex::Sphere_circumradius()); Simplex_tree rips_stree; rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); // ------------------------------------------ diff --git a/src/Cech_complex/example/cech_complex_example_from_points.cpp b/src/Cech_complex/example/cech_complex_example_from_points.cpp index 78861951..38021e4a 100644 --- a/src/Cech_complex/example/cech_complex_example_from_points.cpp +++ b/src/Cech_complex/example/cech_complex_example_from_points.cpp @@ -16,7 +16,7 @@ int main() { using FT = typename Kernel::FT; using Point = typename Kernel::Point_d; using Point_cloud = std::vector; - using Cech_complex = Gudhi::cech_complex::Cech_complex; + using Cech_complex = Gudhi::cech_complex::Cech_complex; Point_cloud points; diff --git a/src/Cech_complex/example/cech_complex_step_by_step.cpp b/src/Cech_complex/example/cech_complex_step_by_step.cpp index c8dd1585..4401f6af 100644 --- a/src/Cech_complex/example/cech_complex_step_by_step.cpp +++ b/src/Cech_complex/example/cech_complex_step_by_step.cpp @@ -9,7 +9,7 @@ */ #include -#include +#include #include #include @@ -52,7 +52,7 @@ class Cech_blocker { std::clog << "#(" << vertex << ")#"; #endif // DEBUG_TRACES } - Filtration_value radius = Gudhi::Minimal_enclosing_ball_radius()(points); + Filtration_value radius = Gudhi::cech_complex::Sphere_circumradius()(points); #ifdef DEBUG_TRACES std::clog << "radius = " << radius << " - " << (radius > max_radius_) << std::endl; #endif // DEBUG_TRACES @@ -83,7 +83,7 @@ int main(int argc, char* argv[]) { // Compute the proximity graph of the points Proximity_graph prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), max_radius, - Gudhi::Minimal_enclosing_ball_radius()); + Gudhi::cech_complex::Sphere_circumradius()); // Construct the Cech complex in a Simplex Tree Simplex_tree st; diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 0031d861..375be1d2 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -11,7 +11,7 @@ #ifndef CECH_COMPLEX_H_ #define CECH_COMPLEX_H_ -#include // for Gudhi::Minimal_enclosing_ball_radius +#include // for Gudhi::cech_complex::Sphere_circumradius #include // for Gudhi::Proximity_graph #include // for GUDHI_CHECK #include // for Gudhi::cech_complex::Cech_blocker @@ -31,23 +31,21 @@ namespace cech_complex { * * \details * The data structure is a proximity graph, containing edges when the edge length is less or equal - * to a given max_radius. Edge length is computed from `Gudhi::Minimal_enclosing_ball_radius` distance function. + * to a given max_radius. Edge length is computed from `Gudhi::cech_complex::Sphere_circumradius` distance function. * - * \tparam SimplicialComplexForProximityGraph furnishes `Vertex_handle` and `Filtration_value` type definition required - * by `Gudhi::Proximity_graph`. + * \tparam Kernel CGAL kernel. + * + * \tparam SimplicialComplexForCechComplex furnishes `Vertex_handle` and `Filtration_value` type definition required + * by `Gudhi::Proximity_graph` and Cech blocker. * - * \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 - using Vertex_handle = typename SimplicialComplexForProximityGraph::Vertex_handle; - using Filtration_value = typename SimplicialComplexForProximityGraph::Filtration_value; - using Proximity_graph = Gudhi::Proximity_graph; - - public: + using Vertex_handle = typename SimplicialComplexForCechComplex::Vertex_handle; + using Filtration_value = typename SimplicialComplexForCechComplex::Filtration_value; + using Proximity_graph = Gudhi::Proximity_graph; using cech_blocker = Cech_blocker; @@ -57,27 +55,21 @@ class Cech_complex { // 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; + using Sphere = typename cech_blocker::Sphere; - public: + public: /** \brief Cech_complex constructor from a list of points. * - * @param[in] points Range of points. + * @param[in] points Vector of points where each point is defined as `kernel::Point_d`. * @param[in] max_radius Maximal radius value. * - * \tparam ForwardPointRange must be a range of Point. Point must be a range of copyable Cartesian coordinates. - * */ - Cech_complex(const ForwardPointRange& points, Filtration_value max_radius) : max_radius_(max_radius) { - // Point cloud deep copy - -// point_cloud_.reserve(boost::size(points)); -// for (auto&& point : points) point_cloud_.emplace_back(point.cartesian_begin(), point.cartesian_end()); + Cech_complex(const Point_cloud & points, Filtration_value max_radius) : max_radius_(max_radius) { point_cloud_.assign(points.begin(), points.end()); - cech_skeleton_graph_ = Gudhi::compute_proximity_graph( - point_cloud_, max_radius_, Gudhi::Minimal_enclosing_ball_radius()); + cech_skeleton_graph_ = Gudhi::compute_proximity_graph( + point_cloud_, max_radius_, Sphere_circumradius()); } /** \brief Initializes the simplicial complex from the proximity graph and expands it until a given maximal diff --git a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h b/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h deleted file mode 100644 index 89012206..00000000 --- a/src/Cech_complex/include/gudhi/Cech_complex/Cech_kernel.h +++ /dev/null @@ -1,107 +0,0 @@ -/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. - * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. - * Author(s): Hind Montassif - * - * Copyright (C) 2021 Inria - * - * Modification(s): - * - YYYY/MM Author: Description of the modification - */ - -#ifndef CECH_KERNEL_H_ -#define CECH_KERNEL_H_ - -#include // for #include - -#include // for std::sqrt -#include - -namespace Gudhi { - -// namespace cech_complex { - -/** @brief Compute the radius of the minimal enclosing ball between Points given by a range of coordinates. - * The points are assumed to have the same dimension. */ -template -class Minimal_enclosing_ball_radius { - private: - Kernel kernel_; - public: - using Point = typename Kernel::Point_d; - using Point_cloud = typename std::vector; - - /** \brief Enclosing ball radius from two points using CGAL. - * - * @param[in] point_1 - * @param[in] point_2 - * @return Enclosing ball radius for the two points. - * \tparam Point must be a Kernel::Point_d from CGAL. - * - */ - double operator()(const Point& point_1, const Point& point_2) const { - return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; - } - - - /** \brief Enclosing ball radius from a point cloud using CGAL. - * - * @param[in] point_cloud The points. - * @return Enclosing ball radius for the points. - * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. - * - */ - double operator()(const Point_cloud& point_cloud) const { - return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); - } - -}; - -/** - * \class Cech_kernel - * \brief Cech complex kernel container. - * - * \details - * The Cech complex kernel container stores CGAL Kernel and dispatch basic computations. - */ - -// template < typename Kernel > -// class Cech_kernel { -// private: -// // Kernel for functions access. -// Kernel kernel_; -// public: -// using Point_d = typename Kernel::Point_d; -// // Numeric type of coordinates in the kernel -// using FT = typename Kernel::FT; -// // Sphere is a pair of point and squared radius. -// using Sphere = typename std::pair; -// -// int get_dimension(const Point_d& p0) const { -// return kernel_.point_dimension_d_object()(p0); -// } -// -// 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)); -// } -// -// template -// FT get_squared_radius(PointIterator begin, PointIterator end) const { -// return kernel_.compute_squared_radius_d_object()(begin, end); -// } -// -// FT get_squared_radius(const Sphere& sph) const { -// return sph.second; -// } -// }; - - -//} // namespace cech_complex - -// namespace cechcomplex = cech_complex; - -} // namespace Gudhi - -#endif // CECH_KERNEL_H_ diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index f7f86534..1a696422 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -31,17 +31,15 @@ namespace cech_complex { * \details * Čech blocker is an oracle constructed from a Cech_complex and a simplicial complex. * - * \tparam SimplicialComplexForProximityGraph furnishes `Simplex_handle` and `Filtration_value` type definition, + * \tparam SimplicialComplexForCech furnishes `Simplex_handle` and `Filtration_value` type definition, * `simplex_vertex_range(Simplex_handle sh)`and `assign_filtration(Simplex_handle sh, Filtration_value filt)` methods. * - * \tparam Chech_complex is required by the blocker. + * \tparam Cech_complex is required by the blocker. + * + * \tparam Kernel CGAL kernel. */ template class Cech_blocker { - private: - - using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; - using Filtration_value = typename SimplicialComplexForCech::Filtration_value; public: @@ -51,10 +49,10 @@ class Cech_blocker { // Sphere is a pair of point and squared radius. using Sphere = typename std::pair; - template - FT get_squared_radius(PointIterator begin, PointIterator end) const { - return kernel_.compute_squared_radius_d_object()(begin, end); - } + private: + + using Simplex_handle = typename SimplicialComplexForCech::Simplex_handle; + using Filtration_value = typename SimplicialComplexForCech::Filtration_value; template Sphere get_sphere(PointIterator begin, PointIterator end) const { @@ -63,6 +61,7 @@ class Cech_blocker { return std::make_pair(std::move(c), std::move(r)); } + public: /** \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. @@ -108,10 +107,11 @@ class Cech_blocker { 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; + std::clog << "center: " << sph.first << ", radius: " << radius << std::endl; #endif // DEBUG_TRACES if (cast_to_double(sph.second) < cast_to_double(min_enclos_ball.second)) min_enclos_ball = sph; + break; } } // Get the minimal radius of all faces enclosing balls if exists diff --git a/src/Cech_complex/include/gudhi/sphere_circumradius.h b/src/Cech_complex/include/gudhi/sphere_circumradius.h new file mode 100644 index 00000000..a6dec3dc --- /dev/null +++ b/src/Cech_complex/include/gudhi/sphere_circumradius.h @@ -0,0 +1,62 @@ +/* This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT. + * See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details. + * Author(s): Hind Montassif + * + * Copyright (C) 2021 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef SPHERE_CIRCUMRADIUS_H_ +#define SPHERE_CIRCUMRADIUS_H_ + +#include // for #include + +#include // for std::sqrt +#include + +namespace Gudhi { + +namespace cech_complex { + +/** @brief Compute the circumradius of the sphere passing through points given by a range of coordinates. + * The points are assumed to have the same dimension. */ +template +class Sphere_circumradius { + private: + Kernel kernel_; + public: + using Point = typename Kernel::Point_d; + using Point_cloud = typename std::vector; + + /** \brief Circumradius of sphere passing through two points using CGAL. + * + * @param[in] point_1 + * @param[in] point_2 + * @return Sphere circumradius passing through two points. + * \tparam Point must be a Kernel::Point_d from CGAL. + * + */ + double operator()(const Point& point_1, const Point& point_2) const { + return std::sqrt(CGAL::to_double(kernel_.squared_distance_d_object()(point_1, point_2))) / 2.; + } + + /** \brief Circumradius of sphere passing through point cloud using CGAL. + * + * @param[in] point_cloud The points. + * @return Sphere circumradius passing through the points. + * \tparam Point_cloud must be a range of Kernel::Point_d points from CGAL. + * + */ + double operator()(const Point_cloud& point_cloud) const { + return std::sqrt(CGAL::to_double(kernel_.compute_squared_radius_d_object()(point_cloud.begin(), point_cloud.end()))); + } + +}; + +} // namespace cech_complex + +} // namespace Gudhi + +#endif // SPHERE_CIRCUMRADIUS_H_ diff --git a/src/Cech_complex/test/test_cech_complex.cpp b/src/Cech_complex/test/test_cech_complex.cpp index ca7a9778..4cf8b68f 100644 --- a/src/Cech_complex/test/test_cech_complex.cpp +++ b/src/Cech_complex/test/test_cech_complex.cpp @@ -22,7 +22,7 @@ // to construct Cech_complex from a OFF file of points #include #include -#include +#include #include #include // For EXACT or SAFE version @@ -36,7 +36,7 @@ using Point = typename Kernel::Point_d; using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; -using Cech_complex = Gudhi::cech_complex::Cech_complex; +using Cech_complex = Gudhi::cech_complex::Cech_complex; BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { // ---------------------------------------------------------------------------- @@ -108,11 +108,11 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { std::clog << vertex << ","; vp.push_back(points.at(vertex)); } - std::clog << ") - distance =" << Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1)) + std::clog << ") - distance =" << Gudhi::cech_complex::Sphere_circumradius()(vp.at(0), vp.at(1)) << " - filtration =" << st.filtration(f_simplex) << std::endl; BOOST_CHECK(vp.size() == 2); GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), - Gudhi::Minimal_enclosing_ball_radius()(vp.at(0), vp.at(1))); + Gudhi::cech_complex::Sphere_circumradius()(vp.at(0), vp.at(1))); } } @@ -153,7 +153,7 @@ BOOST_AUTO_TEST_CASE(Cech_complex_for_documentation) { Simplex_tree::Filtration_value f1410 = st2.filtration(st2.find({1, 4, 10})); std::clog << "f1410= " << f1410 << std::endl; - // In this case, the computed sphere using CGAL kernel does not match the minimal enclosing ball; the filtration value check is therefore done against a hardcoded value + // In this case, the computed circumsphere using CGAL kernel does not match the minimal enclosing ball; the filtration value check is therefore done against a hardcoded value GUDHI_TEST_FLOAT_EQUALITY_CHECK(f1410, 1.); Point_cloud points469; diff --git a/src/Cech_complex/utilities/cech_persistence.cpp b/src/Cech_complex/utilities/cech_persistence.cpp index ccf63e3e..82992f2d 100644 --- a/src/Cech_complex/utilities/cech_persistence.cpp +++ b/src/Cech_complex/utilities/cech_persistence.cpp @@ -9,7 +9,7 @@ */ #include -#include +#include #include #include #include @@ -28,9 +28,8 @@ using Filtration_value = Simplex_tree::Filtration_value; using Kernel = CGAL::Epeck_d; using Point = typename Kernel::Point_d; -using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; -using Cech_complex = Gudhi::cech_complex::Cech_complex; +using Cech_complex = Gudhi::cech_complex::Cech_complex; using Field_Zp = Gudhi::persistent_cohomology::Field_Zp; using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology; -- cgit v1.2.3 From af74316148165cb01c9f28b8b05e1b9764e4579a Mon Sep 17 00:00:00 2001 From: Hind-M Date: Fri, 4 Mar 2022 14:05:15 +0100 Subject: Use Euclidean_distance instead of CGAL dependant Sphere_circumradius for Rips in Cech benchmark --- src/Cech_complex/benchmark/cech_complex_benchmark.cpp | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) (limited to 'src/Cech_complex/benchmark/cech_complex_benchmark.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index b283e1a8..bf013a81 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -107,8 +107,7 @@ int main(int argc, char* argv[]) { std::clog << itr->path().stem() << ";"; std::clog << radius << ";"; Gudhi::Clock rips_clock("Rips computation"); - Rips_complex rips_complex_from_points(off_reader_cgal.get_point_cloud(), radius, - Gudhi::cech_complex::Sphere_circumradius()); + Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), radius, Gudhi::Euclidean_distance()); Simplex_tree rips_stree; rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); // ------------------------------------------ -- cgit v1.2.3 From d5e1760353e7d6ba66975a90f9a2768a48f0abf8 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Tue, 8 Mar 2022 15:27:19 +0100 Subject: Modify cech benchmark to include both Epick and Epeck cases --- .../benchmark/cech_complex_benchmark.cpp | 177 +++++++++++++-------- 1 file changed, 109 insertions(+), 68 deletions(-) (limited to 'src/Cech_complex/benchmark/cech_complex_benchmark.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index bf013a81..9cf24542 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -18,6 +18,7 @@ #include #include +#include #include "boost/filesystem.hpp" // includes all needed Boost.Filesystem declarations @@ -32,10 +33,6 @@ using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; using Proximity_graph = Gudhi::Proximity_graph; using Rips_complex = Gudhi::rips_complex::Rips_complex; -using Kernel = CGAL::Epick_d>; -using Point_cgal = typename Kernel::Point_d; -using Points_off_reader_cgal = Gudhi::Points_off_reader; -using Cech_complex = Gudhi::cech_complex::Cech_complex; class Minimal_enclosing_ball_radius { public: @@ -61,79 +58,123 @@ class Minimal_enclosing_ball_radius { } }; -int main(int argc, char* argv[]) { - std::string off_file_points = "tore3D_1307.off"; - Filtration_value threshold = 1e20; +enum distance_type { Euclidean_dist, Minimal_enclosing_ball_dist, CGAL_dist }; - // Extract the points from the file filepoints - Points_off_reader off_reader(off_file_points); - Points_off_reader_cgal off_reader_cgal(off_file_points); +template> +void benchmark_prox_graph(const std::string& off_file_points, const Filtration_value& threshold, const std::string& msg, distance_type dist = CGAL_dist) { + if (dist != CGAL_dist) { + std::cerr << "Error: when CGAL is used, the distance should be CGAL_dist" << std::endl; + exit(-1); + } + if (!use_cgal) { + std::cerr << "Warning: if kernel is given, CGAL will be used" << std::endl; + } + using Point_cgal = typename Kernel::Point_d; + using Points_off_reader_cgal = Gudhi::Points_off_reader; - Gudhi::Clock euclidean_clock("Gudhi::Euclidean_distance"); - // Compute the proximity graph of the points - Proximity_graph euclidean_prox_graph = Gudhi::compute_proximity_graph( - off_reader.get_point_cloud(), threshold, Gudhi::Euclidean_distance()); + // Extract the points from the file filepoints + Points_off_reader_cgal off_reader_cgal(off_file_points); - std::clog << euclidean_clock << std::endl; + Gudhi::Clock cgal_circumsphere_clock("Gudhi::cech_complex::Sphere_circumradius_cgal()"); + // Compute the proximity graph of the points + Proximity_graph cgal_circumsphere_prox_graph = Gudhi::compute_proximity_graph(off_reader_cgal.get_point_cloud(), threshold, + Gudhi::cech_complex::Sphere_circumradius()); + std::clog << msg << " - " << cgal_circumsphere_clock << std::endl; +} - Gudhi::Clock miniball_clock("Minimal_enclosing_ball_radius"); - // Compute the proximity graph of the points - Proximity_graph miniball_prox_graph = Gudhi::compute_proximity_graph( - off_reader.get_point_cloud(), threshold, Minimal_enclosing_ball_radius()); - std::clog << miniball_clock << std::endl; +template +void benchmark_prox_graph(const std::string& off_file_points, const Filtration_value& threshold, const std::string& msg, distance_type dist) { + // Extract the points from the file filepoints + Points_off_reader off_reader(off_file_points); + + if (dist == Euclidean_dist) { + Gudhi::Clock euclidean_clock("Gudhi::Euclidean_distance"); + // Compute the proximity graph of the points + Proximity_graph euclidean_prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), threshold, + Gudhi::Euclidean_distance()); + std::clog << msg << " - " << euclidean_clock << std::endl; + } + else if (dist == Minimal_enclosing_ball_dist) { + Gudhi::Clock miniball_clock("Minimal_enclosing_ball_radius"); + // Compute the proximity graph of the points + Proximity_graph miniball_prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), threshold, + Minimal_enclosing_ball_radius()); + std::clog << msg << " - " << miniball_clock << std::endl; + } + else { + std::cerr << "Error: when CGAL is not used, the distance should be either Euclidean_dist or Minimal_enclosing_ball_dist" << std::endl; + exit(-1); + } +} + +template +void benchmark_cech(const std::string& off_file_points, const Filtration_value& radius, const int& dim_max) { + using Point_cgal = typename Kernel::Point_d; + using Points_off_reader_cgal = Gudhi::Points_off_reader; + using Cech_complex = Gudhi::cech_complex::Cech_complex; + + // Extract the points from the file filepoints + Points_off_reader_cgal off_reader_cgal(off_file_points); + + Gudhi::Clock cech_clock("Cech computation"); + Cech_complex cech_complex_from_points(off_reader_cgal.get_point_cloud(), radius); + Simplex_tree cech_stree; + cech_complex_from_points.create_complex(cech_stree, dim_max); + + // ------------------------------------------ + // Display information about the Cech complex + // ------------------------------------------ + double cech_sec = cech_clock.num_seconds(); + std::clog << cech_sec << " ; "; + std::clog << cech_stree.num_simplices() << " ; "; +} - Gudhi::Clock cgal_circumsphere_clock("Gudhi::cech_complex::Sphere_circumradius_cgal()"); - // Compute the proximity graph of the points - Proximity_graph cgal_circumsphere_prox_graph = Gudhi::compute_proximity_graph( - off_reader_cgal.get_point_cloud(), threshold, Gudhi::cech_complex::Sphere_circumradius()); - std::clog << cgal_circumsphere_clock << std::endl; +int main(int argc, char* argv[]) { + std::string off_file_points = "tore3D_1307.off"; + Filtration_value threshold = 1e20; + + benchmark_prox_graph(off_file_points, threshold, "Euclidean distance", Euclidean_dist); + benchmark_prox_graph(off_file_points, threshold, "Minimal_enclosing_ball", Minimal_enclosing_ball_dist); + benchmark_prox_graph>>(off_file_points, threshold, "Epick"); + benchmark_prox_graph>>(off_file_points, threshold, "Epeck"); - boost::filesystem::path full_path(boost::filesystem::current_path()); - std::clog << "Current path is : " << full_path << std::endl; + boost::filesystem::path full_path(boost::filesystem::current_path()); + std::clog << "Current path is : " << full_path << std::endl; - std::clog << "File name; Radius; Rips time; Cech time; Ratio Rips/Cech time; Rips nb simplices; Cech nb simplices;" + std::clog << "File name ; Radius ; Rips time ; Epick Cech time ; Epick Cech nb simplices ; Epeck Cech time ; Epeck Cech nb simplices ; Rips nb simplices;" << std::endl; - boost::filesystem::directory_iterator end_itr; // default construction yields past-the-end - for (boost::filesystem::directory_iterator itr(boost::filesystem::current_path()); itr != end_itr; ++itr) { - if (!boost::filesystem::is_directory(itr->status())) { - if (itr->path().extension() == ".off") // see below - { - Points_off_reader off_reader(itr->path().string()); - Points_off_reader_cgal off_reader_cgal(itr->path().string()); - - Point p0 = off_reader.get_point_cloud()[0]; - - for (Filtration_value radius = 0.1; radius < 0.4; radius += 0.1) { - std::clog << itr->path().stem() << ";"; - std::clog << radius << ";"; - Gudhi::Clock rips_clock("Rips computation"); - Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), radius, Gudhi::Euclidean_distance()); - Simplex_tree rips_stree; - rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); - // ------------------------------------------ - // Display information about the Rips complex - // ------------------------------------------ - double rips_sec = rips_clock.num_seconds(); - std::clog << rips_sec << ";"; - - Gudhi::Clock cech_clock("Cech computation"); - Cech_complex cech_complex_from_points(off_reader_cgal.get_point_cloud(), radius); - Simplex_tree cech_stree; - cech_complex_from_points.create_complex(cech_stree, p0.size() - 1); - // ------------------------------------------ - // Display information about the Cech complex - // ------------------------------------------ - double cech_sec = cech_clock.num_seconds(); - std::clog << cech_sec << ";"; - std::clog << cech_sec / rips_sec << ";"; - - assert(rips_stree.num_simplices() >= cech_stree.num_simplices()); - std::clog << rips_stree.num_simplices() << ";"; - std::clog << cech_stree.num_simplices() << ";" << std::endl; + boost::filesystem::directory_iterator end_itr; // default construction yields past-the-end + for (boost::filesystem::directory_iterator itr(boost::filesystem::current_path()); itr != end_itr; ++itr) { + if (!boost::filesystem::is_directory(itr->status())) { + if (itr->path().extension() == ".off") { + Points_off_reader off_reader(itr->path().string()); + Point p0 = off_reader.get_point_cloud()[0]; + + for (Filtration_value radius = 0.1; radius < 0.4; radius += 0.1) { + std::clog << itr->path().stem() << " ; "; + std::clog << radius << " ; "; + + Gudhi::Clock rips_clock("Rips computation"); + Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), radius, Gudhi::Euclidean_distance()); + Simplex_tree rips_stree; + rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); + // ------------------------------------------ + // Display information about the Rips complex + // ------------------------------------------ + double rips_sec = rips_clock.num_seconds(); + std::clog << rips_sec << " ; "; + + // -------------- + // Cech complex + // -------------- + benchmark_cech>>(itr->path().string(), radius, p0.size() - 1); + benchmark_cech>>(itr->path().string(), radius, p0.size() - 1); + + std::clog << rips_stree.num_simplices() << ";" << std::endl; + } + } } - } } - } - return 0; + return 0; } -- cgit v1.2.3 From dbb65c3f3eb82d080e47b40b52deb03814d8da31 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Mon, 25 Apr 2022 13:18:11 +0200 Subject: Remove proximity_graph computation benchmark Add Dynamic_dimension_tag case --- .../benchmark/cech_complex_benchmark.cpp | 97 ++-------------------- 1 file changed, 8 insertions(+), 89 deletions(-) (limited to 'src/Cech_complex/benchmark/cech_complex_benchmark.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 9cf24542..d2a71879 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -10,12 +10,10 @@ #include #include -#include #include #include #include #include -#include #include #include @@ -29,86 +27,11 @@ using Simplex_tree = Gudhi::Simplex_tree<>; using Filtration_value = Simplex_tree::Filtration_value; using Point = std::vector; -using Point_cloud = std::vector; using Points_off_reader = Gudhi::Points_off_reader; -using Proximity_graph = Gudhi::Proximity_graph; using Rips_complex = Gudhi::rips_complex::Rips_complex; -class Minimal_enclosing_ball_radius { - public: - // boost::range_value is not SFINAE-friendly so we cannot use it in the return type - template - typename std::iterator_traits::type>::value_type operator()( - const Point& p1, const Point& p2) const { - // Type def - using Point_cloud = std::vector; - using Point_iterator = typename Point_cloud::const_iterator; - using Coordinate_iterator = typename Point::const_iterator; - using Min_sphere = - typename Gudhi::Miniball::Miniball>; - - Point_cloud point_cloud; - point_cloud.push_back(p1); - point_cloud.push_back(p2); - - GUDHI_CHECK((p1.end() - p1.begin()) == (p2.end() - p2.begin()), "inconsistent point dimensions"); - Min_sphere min_sphere(p1.end() - p1.begin(), point_cloud.begin(), point_cloud.end()); - - return std::sqrt(min_sphere.squared_radius()); - } -}; - -enum distance_type { Euclidean_dist, Minimal_enclosing_ball_dist, CGAL_dist }; - -template> -void benchmark_prox_graph(const std::string& off_file_points, const Filtration_value& threshold, const std::string& msg, distance_type dist = CGAL_dist) { - if (dist != CGAL_dist) { - std::cerr << "Error: when CGAL is used, the distance should be CGAL_dist" << std::endl; - exit(-1); - } - if (!use_cgal) { - std::cerr << "Warning: if kernel is given, CGAL will be used" << std::endl; - } - using Point_cgal = typename Kernel::Point_d; - using Points_off_reader_cgal = Gudhi::Points_off_reader; - - // Extract the points from the file filepoints - Points_off_reader_cgal off_reader_cgal(off_file_points); - - Gudhi::Clock cgal_circumsphere_clock("Gudhi::cech_complex::Sphere_circumradius_cgal()"); - // Compute the proximity graph of the points - Proximity_graph cgal_circumsphere_prox_graph = Gudhi::compute_proximity_graph(off_reader_cgal.get_point_cloud(), threshold, - Gudhi::cech_complex::Sphere_circumradius()); - std::clog << msg << " - " << cgal_circumsphere_clock << std::endl; -} - -template -void benchmark_prox_graph(const std::string& off_file_points, const Filtration_value& threshold, const std::string& msg, distance_type dist) { - // Extract the points from the file filepoints - Points_off_reader off_reader(off_file_points); - - if (dist == Euclidean_dist) { - Gudhi::Clock euclidean_clock("Gudhi::Euclidean_distance"); - // Compute the proximity graph of the points - Proximity_graph euclidean_prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), threshold, - Gudhi::Euclidean_distance()); - std::clog << msg << " - " << euclidean_clock << std::endl; - } - else if (dist == Minimal_enclosing_ball_dist) { - Gudhi::Clock miniball_clock("Minimal_enclosing_ball_radius"); - // Compute the proximity graph of the points - Proximity_graph miniball_prox_graph = Gudhi::compute_proximity_graph(off_reader.get_point_cloud(), threshold, - Minimal_enclosing_ball_radius()); - std::clog << msg << " - " << miniball_clock << std::endl; - } - else { - std::cerr << "Error: when CGAL is not used, the distance should be either Euclidean_dist or Minimal_enclosing_ball_dist" << std::endl; - exit(-1); - } -} - template -void benchmark_cech(const std::string& off_file_points, const Filtration_value& radius, const int& dim_max) { +Simplex_tree benchmark_cech(const std::string& off_file_points, const Filtration_value& radius, const int& dim_max) { using Point_cgal = typename Kernel::Point_d; using Points_off_reader_cgal = Gudhi::Points_off_reader; using Cech_complex = Gudhi::cech_complex::Cech_complex; @@ -126,23 +49,16 @@ void benchmark_cech(const std::string& off_file_points, const Filtration_value& // ------------------------------------------ double cech_sec = cech_clock.num_seconds(); std::clog << cech_sec << " ; "; - std::clog << cech_stree.num_simplices() << " ; "; + return cech_stree; } int main(int argc, char* argv[]) { - std::string off_file_points = "tore3D_1307.off"; - Filtration_value threshold = 1e20; - - benchmark_prox_graph(off_file_points, threshold, "Euclidean distance", Euclidean_dist); - benchmark_prox_graph(off_file_points, threshold, "Minimal_enclosing_ball", Minimal_enclosing_ball_dist); - benchmark_prox_graph>>(off_file_points, threshold, "Epick"); - benchmark_prox_graph>>(off_file_points, threshold, "Epeck"); - boost::filesystem::path full_path(boost::filesystem::current_path()); std::clog << "Current path is : " << full_path << std::endl; - std::clog << "File name ; Radius ; Rips time ; Epick Cech time ; Epick Cech nb simplices ; Epeck Cech time ; Epeck Cech nb simplices ; Rips nb simplices;" - << std::endl; + std::clog << "File name ; Radius ; Rips time ; Dim-3 Epick Cech time ; Dynamic_dim Epick Cech time ; " + "Dim-3 Epeck Cech time ; Dynamic_dim Epeck Cech time ; Cech nb simplices ; Rips nb simplices;" + << std::endl; boost::filesystem::directory_iterator end_itr; // default construction yields past-the-end for (boost::filesystem::directory_iterator itr(boost::filesystem::current_path()); itr != end_itr; ++itr) { if (!boost::filesystem::is_directory(itr->status())) { @@ -168,8 +84,11 @@ int main(int argc, char* argv[]) { // Cech complex // -------------- benchmark_cech>>(itr->path().string(), radius, p0.size() - 1); + benchmark_cech>(itr->path().string(), radius, p0.size() - 1); benchmark_cech>>(itr->path().string(), radius, p0.size() - 1); + auto cech_stree = benchmark_cech>(itr->path().string(), radius, p0.size() - 1); + std::clog << cech_stree.num_simplices() << " ; "; std::clog << rips_stree.num_simplices() << ";" << std::endl; } } -- cgit v1.2.3 From 868369dd61fb6ef475ffa3af724907927121b6bb Mon Sep 17 00:00:00 2001 From: Hind-M Date: Thu, 16 Jun 2022 15:54:21 +0200 Subject: Add exact option for exact cech variant --- .../benchmark/cech_complex_benchmark.cpp | 22 ++++++++++++++-------- src/Cech_complex/include/gudhi/Cech_complex.h | 6 ++++-- .../include/gudhi/Cech_complex_blocker.h | 21 ++++++++++++++------- 3 files changed, 32 insertions(+), 17 deletions(-) (limited to 'src/Cech_complex/benchmark/cech_complex_benchmark.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index d2a71879..19142780 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -31,7 +31,7 @@ using Points_off_reader = Gudhi::Points_off_reader; using Rips_complex = Gudhi::rips_complex::Rips_complex; template -Simplex_tree benchmark_cech(const std::string& off_file_points, const Filtration_value& radius, const int& dim_max) { +Simplex_tree benchmark_cech(const std::string& off_file_points, const Filtration_value& radius, const int& dim_max, const bool exact) { using Point_cgal = typename Kernel::Point_d; using Points_off_reader_cgal = Gudhi::Points_off_reader; using Cech_complex = Gudhi::cech_complex::Cech_complex; @@ -42,7 +42,7 @@ Simplex_tree benchmark_cech(const std::string& off_file_points, const Filtration Gudhi::Clock cech_clock("Cech computation"); Cech_complex cech_complex_from_points(off_reader_cgal.get_point_cloud(), radius); Simplex_tree cech_stree; - cech_complex_from_points.create_complex(cech_stree, dim_max); + cech_complex_from_points.create_complex(cech_stree, dim_max, exact); // ------------------------------------------ // Display information about the Cech complex @@ -56,8 +56,9 @@ int main(int argc, char* argv[]) { boost::filesystem::path full_path(boost::filesystem::current_path()); std::clog << "Current path is : " << full_path << std::endl; - std::clog << "File name ; Radius ; Rips time ; Dim-3 Epick Cech time ; Dynamic_dim Epick Cech time ; " - "Dim-3 Epeck Cech time ; Dynamic_dim Epeck Cech time ; Cech nb simplices ; Rips nb simplices;" + std::clog << "File name ; Radius ; Rips time ; Dim-3 Fast Cech time ; Dynamic_dim Fast Cech time ; " + "Dim-3 Safe Cech time ; Dynamic_dim Safe Cech time ; Dim-3 Exact Cech time ; Dynamic_dim Exact Cech time ; " + "Cech nb simplices ; Rips nb simplices;" << std::endl; boost::filesystem::directory_iterator end_itr; // default construction yields past-the-end for (boost::filesystem::directory_iterator itr(boost::filesystem::current_path()); itr != end_itr; ++itr) { @@ -83,10 +84,15 @@ int main(int argc, char* argv[]) { // -------------- // Cech complex // -------------- - benchmark_cech>>(itr->path().string(), radius, p0.size() - 1); - benchmark_cech>(itr->path().string(), radius, p0.size() - 1); - benchmark_cech>>(itr->path().string(), radius, p0.size() - 1); - auto cech_stree = benchmark_cech>(itr->path().string(), radius, p0.size() - 1); + // Fast + benchmark_cech>>(itr->path().string(), radius, p0.size() - 1, false); + benchmark_cech>(itr->path().string(), radius, p0.size() - 1, false); + // Safe + benchmark_cech>>(itr->path().string(), radius, p0.size() - 1, false); + benchmark_cech>(itr->path().string(), radius, p0.size() - 1, false); + // Exact + benchmark_cech>>(itr->path().string(), radius, p0.size() - 1, true); + auto cech_stree = benchmark_cech>(itr->path().string(), radius, p0.size() - 1, true); std::clog << cech_stree.num_simplices() << " ; "; std::clog << rips_stree.num_simplices() << ";" << std::endl; diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index fc39f75b..2c6d3df5 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -78,17 +78,19 @@ class Cech_complex { * * @param[in] complex SimplicialComplexForCech to be created. * @param[in] dim_max graph expansion until this given maximal dimension. + * @param[in] exact Exact filtration values computation. Not exact if `Kernel` is not CGAL::Epeck_d. * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0. * */ - void create_complex(SimplicialComplexForCechComplex& complex, int dim_max) { + void create_complex(SimplicialComplexForCechComplex& complex, int dim_max, const bool exact = false) { 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 complex.insert_graph(cech_skeleton_graph_); // expand the graph until dimension dim_max - complex.expansion_with_blockers(dim_max, cech_blocker(&complex, this)); + complex.expansion_with_blockers(dim_max, cech_blocker(&complex, this, exact)); } /** @return max_radius value given at construction. */ diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 3141d27a..087390b6 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -94,9 +94,9 @@ class Cech_blocker { Point_cloud face_points; for (auto vertex : sc_ptr_->simplex_vertex_range(face)) { face_points.push_back(cc_ptr_->get_point(vertex)); - #ifdef DEBUG_TRACES - std::clog << "#(" << vertex << ")#"; - #endif // DEBUG_TRACES +#ifdef DEBUG_TRACES + std::clog << "#(" << vertex << ")#"; +#endif // DEBUG_TRACES } sph = get_sphere(face_points.cbegin(), face_points.cend()); // Put edge sphere in cache @@ -107,10 +107,13 @@ class Cech_blocker { } // Check if the minimal enclosing ball of current face contains the extra point if (kernel_.squared_distance_d_object()(sph.first, cc_ptr_->get_point(extra)) <= sph.second) { - #ifdef DEBUG_TRACES - std::clog << "center: " << sph.first << ", radius: " << radius << std::endl; - #endif // DEBUG_TRACES +#ifdef DEBUG_TRACES + std::clog << "center: " << sph.first << ", radius: " << radius << std::endl; +#endif // DEBUG_TRACES is_min_enclos_ball = true; +#if CGAL_VERSION_NR >= 1050000000 + if(exact_) CGAL::exact(sph.second); +#endif radius = std::sqrt(cast_to_fv(sph.second)); sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); cc_ptr_->get_cache().push_back(sph); @@ -124,6 +127,9 @@ class Cech_blocker { points.push_back(cc_ptr_->get_point(vertex)); } Sphere sph = get_sphere(points.cbegin(), points.cend()); +#if CGAL_VERSION_NR >= 1050000000 + if(exact_) CGAL::exact(sph.second); +#endif radius = std::sqrt(cast_to_fv(sph.second)); sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); @@ -138,12 +144,13 @@ class Cech_blocker { } /** \internal \brief Čech complex blocker constructor. */ - Cech_blocker(SimplicialComplexForCech* sc_ptr, Cech_complex* cc_ptr) : sc_ptr_(sc_ptr), cc_ptr_(cc_ptr) {} + Cech_blocker(SimplicialComplexForCech* sc_ptr, Cech_complex* cc_ptr, const bool exact) : sc_ptr_(sc_ptr), cc_ptr_(cc_ptr), exact_(exact) {} private: SimplicialComplexForCech* sc_ptr_; Cech_complex* cc_ptr_; Kernel kernel_; + const bool exact_; }; } // namespace cech_complex -- cgit v1.2.3 From 3fa972970514333d4db22ec7628c5c1a4de3c6e8 Mon Sep 17 00:00:00 2001 From: Hind-M Date: Tue, 21 Jun 2022 15:04:27 +0200 Subject: -Add/modify some comments -Some other minor changes -Change license to LGPL --- .../benchmark/cech_complex_benchmark.cpp | 20 +++++++++++--------- src/Cech_complex/include/gudhi/Cech_complex.h | 2 +- .../include/gudhi/Cech_complex_blocker.h | 2 +- src/common/doc/main_page.md | 2 +- 4 files changed, 14 insertions(+), 12 deletions(-) (limited to 'src/Cech_complex/benchmark/cech_complex_benchmark.cpp') diff --git a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp index 19142780..a9dc5d0d 100644 --- a/src/Cech_complex/benchmark/cech_complex_benchmark.cpp +++ b/src/Cech_complex/benchmark/cech_complex_benchmark.cpp @@ -61,20 +61,22 @@ int main(int argc, char* argv[]) { "Cech nb simplices ; Rips nb simplices;" << std::endl; boost::filesystem::directory_iterator end_itr; // default construction yields past-the-end + // For every ".off" file in the current directory, and for 3 predefined thresholds, compare Rips and various Cech constructions for (boost::filesystem::directory_iterator itr(boost::filesystem::current_path()); itr != end_itr; ++itr) { if (!boost::filesystem::is_directory(itr->status())) { if (itr->path().extension() == ".off") { Points_off_reader off_reader(itr->path().string()); Point p0 = off_reader.get_point_cloud()[0]; - - for (Filtration_value radius = 0.1; radius < 0.4; radius += 0.1) { + // Loop over the different thresholds + for (Filtration_value radius = 0.1; radius < 0.35; radius += 0.1) { std::clog << itr->path().stem() << " ; "; std::clog << radius << " ; "; Gudhi::Clock rips_clock("Rips computation"); Rips_complex rips_complex_from_points(off_reader.get_point_cloud(), radius, Gudhi::Euclidean_distance()); Simplex_tree rips_stree; - rips_complex_from_points.create_complex(rips_stree, p0.size() - 1); + int dim_max = p0.size() - 1; + rips_complex_from_points.create_complex(rips_stree, dim_max); // ------------------------------------------ // Display information about the Rips complex // ------------------------------------------ @@ -85,14 +87,14 @@ int main(int argc, char* argv[]) { // Cech complex // -------------- // Fast - benchmark_cech>>(itr->path().string(), radius, p0.size() - 1, false); - benchmark_cech>(itr->path().string(), radius, p0.size() - 1, false); + benchmark_cech>>(itr->path().string(), radius, dim_max, false); + benchmark_cech>(itr->path().string(), radius, dim_max, false); // Safe - benchmark_cech>>(itr->path().string(), radius, p0.size() - 1, false); - benchmark_cech>(itr->path().string(), radius, p0.size() - 1, false); + benchmark_cech>>(itr->path().string(), radius, dim_max, false); + benchmark_cech>(itr->path().string(), radius, dim_max, false); // Exact - benchmark_cech>>(itr->path().string(), radius, p0.size() - 1, true); - auto cech_stree = benchmark_cech>(itr->path().string(), radius, p0.size() - 1, true); + benchmark_cech>>(itr->path().string(), radius, dim_max, true); + auto cech_stree = benchmark_cech>(itr->path().string(), radius, dim_max, true); std::clog << cech_stree.num_simplices() << " ; "; std::clog << rips_stree.num_simplices() << ";" << std::endl; diff --git a/src/Cech_complex/include/gudhi/Cech_complex.h b/src/Cech_complex/include/gudhi/Cech_complex.h index 2c6d3df5..bae21d28 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex.h +++ b/src/Cech_complex/include/gudhi/Cech_complex.h @@ -30,7 +30,7 @@ namespace cech_complex { * \ingroup cech_complex * * \details - * Cech complex is a simplicial complex constructed from a proximity graph, where the set of all simplices is filtered + * Cech complex is a simplicial complex where the set of all simplices is filtered * by the radius of their minimal enclosing ball and bounded by the given max_radius. * * \tparam Kernel CGAL kernel: either Epick_d or Epeck_d. diff --git a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h index 087390b6..9cd49a52 100644 --- a/src/Cech_complex/include/gudhi/Cech_complex_blocker.h +++ b/src/Cech_complex/include/gudhi/Cech_complex_blocker.h @@ -133,7 +133,7 @@ class Cech_blocker { radius = std::sqrt(cast_to_fv(sph.second)); sc_ptr_->assign_key(sh, cc_ptr_->get_cache().size()); - cc_ptr_->get_cache().push_back(sph); + cc_ptr_->get_cache().push_back(std::move(sph)); } #ifdef DEBUG_TRACES diff --git a/src/common/doc/main_page.md b/src/common/doc/main_page.md index 2cb02e3f..ce903405 100644 --- a/src/common/doc/main_page.md +++ b/src/common/doc/main_page.md @@ -180,7 +180,7 @@ Author: Vincent Rouvreau
Introduced in: GUDHI 2.2.0
- Copyright: MIT [(GPL v3)](../../licensing/)
+ Copyright: MIT [(LGPL v3)](../../licensing/)
Requires: \ref cgal -- cgit v1.2.3