From 88a36ffad6c11279990c1c96df32b95c1f6f526c Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Fri, 3 Jul 2020 13:57:49 +0200 Subject: A fix proposal for boudaries of a simplex python version --- src/python/include/Simplex_tree_interface.h | 11 +++++++++++ 1 file changed, 11 insertions(+) (limited to 'src/python/include') diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 56d7c41d..c4f18eeb 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -36,6 +36,7 @@ class Simplex_tree_interface : public Simplex_tree { using Skeleton_simplex_iterator = typename Base::Skeleton_simplex_iterator; using Complex_simplex_iterator = typename Base::Complex_simplex_iterator; using Extended_filtration_data = typename Base::Extended_filtration_data; + using Boundary_simplex_iterator = typename Base::Boundary_simplex_iterator; public: @@ -188,6 +189,16 @@ class Simplex_tree_interface : public Simplex_tree { // this specific case works because the range is just a pair of iterators - won't work if range was a vector return Base::skeleton_simplex_range(dimension).end(); } + + Boundary_simplex_iterator get_boundary_iterator_begin(const Simplex& simplex) { + // this specific case works because the range is just a pair of iterators - won't work if range was a vector + return Base::boundary_simplex_range(Base::find(simplex)).begin(); + } + + Boundary_simplex_iterator get_boundary_iterator_end(const Simplex& simplex) { + // this specific case works because the range is just a pair of iterators - won't work if range was a vector + return Base::boundary_simplex_range(Base::find(simplex)).end(); + } }; } // namespace Gudhi -- cgit v1.2.3 From 458bc2dcf5044e1d5fde5326b2be35e526abd457 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Wed, 12 Aug 2020 13:06:03 +0200 Subject: code review: boundaries uses only once find and return a pair of iterator. Exception is raised when not found. tested --- src/python/gudhi/simplex_tree.pxd | 3 +-- src/python/gudhi/simplex_tree.pyx | 14 ++++++++------ src/python/include/Simplex_tree_interface.h | 12 +++++------- src/python/test/test_simplex_tree.py | 9 +++++++++ 4 files changed, 23 insertions(+), 15 deletions(-) (limited to 'src/python/include') diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index a64599e7..e1b03391 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -71,8 +71,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": vector[Simplex_tree_simplex_handle].const_iterator get_filtration_iterator_end() nogil Simplex_tree_skeleton_iterator get_skeleton_iterator_begin(int dimension) nogil Simplex_tree_skeleton_iterator get_skeleton_iterator_end(int dimension) nogil - Simplex_tree_boundary_iterator get_boundary_iterator_begin(vector[int] simplex) nogil - Simplex_tree_boundary_iterator get_boundary_iterator_end(vector[int] simplex) nogil + pair[Simplex_tree_boundary_iterator, Simplex_tree_boundary_iterator] get_boundary_iterators(vector[int] simplex) nogil except + cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface>": diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index 3ebae923..bc5b43f4 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -293,12 +293,14 @@ cdef class SimplexTree: :returns: The (simplices of the) boundary of a simplex :rtype: generator with tuples(simplex, filtration) """ - cdef Simplex_tree_boundary_iterator it = self.get_ptr().get_boundary_iterator_begin(simplex) - cdef Simplex_tree_boundary_iterator end = self.get_ptr().get_boundary_iterator_end(simplex) - - while it != end: - yield self.get_ptr().get_simplex_and_filtration(dereference(it)) - preincrement(it) + cdef pair[Simplex_tree_boundary_iterator, Simplex_tree_boundary_iterator] it = self.get_ptr().get_boundary_iterators(simplex) + + try: + while it.first != it.second: + yield self.get_ptr().get_simplex_and_filtration(dereference(it.first)) + preincrement(it.first) + except RuntimeError: + print("simplex not found - cannot find boundaries") def remove_maximal_simplex(self, simplex): diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index c4f18eeb..2c695d1b 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -190,14 +190,12 @@ class Simplex_tree_interface : public Simplex_tree { return Base::skeleton_simplex_range(dimension).end(); } - Boundary_simplex_iterator get_boundary_iterator_begin(const Simplex& simplex) { + std::pair get_boundary_iterators(const Simplex& simplex) { + auto bd_sh = Base::find(simplex); + if (bd_sh == Base::null_simplex()) + throw std::runtime_error("simplex not found - cannot find boundaries"); // this specific case works because the range is just a pair of iterators - won't work if range was a vector - return Base::boundary_simplex_range(Base::find(simplex)).begin(); - } - - Boundary_simplex_iterator get_boundary_iterator_end(const Simplex& simplex) { - // this specific case works because the range is just a pair of iterators - won't work if range was a vector - return Base::boundary_simplex_range(Base::find(simplex)).end(); + return std::make_pair(Base::boundary_simplex_range(bd_sh).begin(), Base::boundary_simplex_range(bd_sh).end()); } }; diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 036e88df..828400fb 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -350,3 +350,12 @@ def test_boundaries_iterator(): assert list(st.get_boundaries([1, 2, 3])) == [([1, 2], 1.0), ([1, 3], 1.0), ([2, 3], 1.0)] assert list(st.get_boundaries([2, 3, 4])) == [([2, 3], 1.0), ([2, 4], 2.0), ([3, 4], 2.0)] assert list(st.get_boundaries([2])) == [] + + with pytest.raises(RuntimeError): + list(st.get_boundaries([])) + + with pytest.raises(RuntimeError): + list(st.get_boundaries([0, 4])) # (0, 4) does not exist + + with pytest.raises(RuntimeError): + list(st.get_boundaries([6])) # (6) does not exist -- cgit v1.2.3 From d9ae2cb822631d68e4ef1cec7eed0124080c0020 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Mon, 24 Aug 2020 09:17:23 +0200 Subject: code review: call boundary_simplex_range only once --- src/python/include/Simplex_tree_interface.h | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'src/python/include') diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index ab3f13f1..baff3850 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -226,7 +226,8 @@ class Simplex_tree_interface : public Simplex_tree { if (bd_sh == Base::null_simplex()) throw std::runtime_error("simplex not found - cannot find boundaries"); // this specific case works because the range is just a pair of iterators - won't work if range was a vector - return std::make_pair(Base::boundary_simplex_range(bd_sh).begin(), Base::boundary_simplex_range(bd_sh).end()); + auto boundary_srange = Base::boundary_simplex_range(bd_sh); + return std::make_pair(boundary_srange.begin(), boundary_srange.end()); } }; -- cgit v1.2.3 From 2bbba93e7f0837b42def9bed13a6fa790c0eabda Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sat, 31 Oct 2020 23:39:01 +0100 Subject: s/kernel/distance/ for choose_n_farthest_points argument --- .../include/gudhi/Sparse_rips_complex.h | 14 +---- src/Subsampling/example/CMakeLists.txt | 4 +- .../example/example_choose_n_farthest_points.cpp | 2 +- .../example/example_custom_distance.cpp | 44 +++++++++++++++ src/Subsampling/example/example_custom_kernel.cpp | 63 ---------------------- .../include/gudhi/choose_n_farthest_points.h | 26 +++++---- .../test/test_choose_n_farthest_points.cpp | 20 +++---- src/Witness_complex/doc/Witness_complex_doc.h | 6 +-- .../example/example_strong_witness_complex_off.cpp | 3 +- .../example/example_witness_complex_off.cpp | 3 +- .../example/example_witness_complex_sphere.cpp | 2 +- .../utilities/strong_witness_persistence.cpp | 3 +- .../utilities/weak_witness_persistence.cpp | 3 +- src/common/doc/examples.h | 2 +- src/common/doc/installation.h | 8 +-- src/python/include/Subsampling_interface.h | 10 ++-- 16 files changed, 93 insertions(+), 120 deletions(-) create mode 100644 src/Subsampling/example/example_custom_distance.cpp delete mode 100644 src/Subsampling/example/example_custom_kernel.cpp (limited to 'src/python/include') diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index 1b250818..a5501004 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -67,8 +67,7 @@ class Sparse_rips_complex { : epsilon_(epsilon) { GUDHI_CHECK(epsilon > 0, "epsilon must be positive"); auto dist_fun = [&](Vertex_handle i, Vertex_handle j) { return distance(points[i], points[j]); }; - Ker kernel(dist_fun); - subsampling::choose_n_farthest_points(kernel, boost::irange(0, boost::size(points)), -1, -1, + subsampling::choose_n_farthest_points(dist_fun, boost::irange(0, boost::size(points)), -1, -1, std::back_inserter(sorted_points), std::back_inserter(params)); compute_sparse_graph(dist_fun, epsilon, mini, maxi); } @@ -128,17 +127,6 @@ class Sparse_rips_complex { } private: - // choose_n_farthest_points wants the distance function in this form... - template - struct Ker { - typedef std::size_t Point_d; // index into point range - Ker(Distance& d) : dist(d) {} - // Despite the name, this is not squared... - typedef Distance Squared_distance_d; - Squared_distance_d& squared_distance_d_object() const { return dist; } - Distance& dist; - }; - // PointRange must be random access. template void compute_sparse_graph(Distance& dist, double epsilon, Filtration_value mini, Filtration_value maxi) { diff --git a/src/Subsampling/example/CMakeLists.txt b/src/Subsampling/example/CMakeLists.txt index 28aab103..fb6875e1 100644 --- a/src/Subsampling/example/CMakeLists.txt +++ b/src/Subsampling/example/CMakeLists.txt @@ -3,7 +3,7 @@ project(Subsampling_examples) if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) add_executable(Subsampling_example_pick_n_random_points example_pick_n_random_points.cpp) add_executable(Subsampling_example_choose_n_farthest_points example_choose_n_farthest_points.cpp) - add_executable(Subsampling_example_custom_kernel example_custom_kernel.cpp) + add_executable(Subsampling_example_custom_distance example_custom_distance.cpp) add_executable(Subsampling_example_sparsify_point_set example_sparsify_point_set.cpp) target_link_libraries(Subsampling_example_sparsify_point_set ${CGAL_LIBRARY}) @@ -16,7 +16,7 @@ if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) install(TARGETS Subsampling_example_pick_n_random_points DESTINATION bin) install(TARGETS Subsampling_example_choose_n_farthest_points DESTINATION bin) - install(TARGETS Subsampling_example_custom_kernel DESTINATION bin) + install(TARGETS Subsampling_example_custom_distance DESTINATION bin) install(TARGETS Subsampling_example_sparsify_point_set DESTINATION bin) endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) diff --git a/src/Subsampling/example/example_choose_n_farthest_points.cpp b/src/Subsampling/example/example_choose_n_farthest_points.cpp index 27cf5d4e..e8b3ce2d 100644 --- a/src/Subsampling/example/example_choose_n_farthest_points.cpp +++ b/src/Subsampling/example/example_choose_n_farthest_points.cpp @@ -20,7 +20,7 @@ int main(void) { K k; std::vector results; - Gudhi::subsampling::choose_n_farthest_points(k, points, 100, + Gudhi::subsampling::choose_n_farthest_points(k.squared_distance_d_object(), points, 100, Gudhi::subsampling::random_starting_point, std::back_inserter(results)); std::clog << "Before sparsification: " << points.size() << " points.\n"; diff --git a/src/Subsampling/example/example_custom_distance.cpp b/src/Subsampling/example/example_custom_distance.cpp new file mode 100644 index 00000000..3325b12d --- /dev/null +++ b/src/Subsampling/example/example_custom_distance.cpp @@ -0,0 +1,44 @@ +#include + +#include +#include +#include + + +typedef unsigned Point; + +/* The class Distance contains a distance function defined on the set of points {0, 1, 2, 3} + * and computes a distance according to the matrix: + * 0 1 2 4 + * 1 0 4 2 + * 2 4 0 1 + * 4 2 1 0 + */ +class Distance { + private: + std::vector> matrix_; + + public: + Distance() { + matrix_.push_back({0, 1, 2, 4}); + matrix_.push_back({1, 0, 4, 2}); + matrix_.push_back({2, 4, 0, 1}); + matrix_.push_back({4, 2, 1, 0}); + } + + double operator()(Point p1, Point p2) const { + return matrix_[p1][p2]; + } +}; + +int main(void) { + std::vector points = {0, 1, 2, 3}; + std::vector results; + + Gudhi::subsampling::choose_n_farthest_points(Distance(), points, 2, + Gudhi::subsampling::random_starting_point, + std::back_inserter(results)); + std::clog << "Before sparsification: " << points.size() << " points.\n"; + std::clog << "After sparsification: " << results.size() << " points.\n"; + std::clog << "Result table: {" << results[0] << "," << results[1] << "}\n"; +} diff --git a/src/Subsampling/example/example_custom_kernel.cpp b/src/Subsampling/example/example_custom_kernel.cpp deleted file mode 100644 index 535bf42a..00000000 --- a/src/Subsampling/example/example_custom_kernel.cpp +++ /dev/null @@ -1,63 +0,0 @@ -#include - -#include -#include -#include - - -/* The class Kernel contains a distance function defined on the set of points {0, 1, 2, 3} - * and computes a distance according to the matrix: - * 0 1 2 4 - * 1 0 4 2 - * 2 4 0 1 - * 4 2 1 0 - */ -class Kernel { - public: - typedef double FT; - typedef unsigned Point_d; - - // Class Squared_distance_d - class Squared_distance_d { - private: - std::vector> matrix_; - - public: - Squared_distance_d() { - matrix_.push_back(std::vector({0, 1, 2, 4})); - matrix_.push_back(std::vector({1, 0, 4, 2})); - matrix_.push_back(std::vector({2, 4, 0, 1})); - matrix_.push_back(std::vector({4, 2, 1, 0})); - } - - FT operator()(Point_d p1, Point_d p2) { - return matrix_[p1][p2]; - } - }; - - // Constructor - Kernel() {} - - // Object of type Squared_distance_d - Squared_distance_d squared_distance_d_object() const { - return Squared_distance_d(); - } -}; - -int main(void) { - typedef Kernel K; - typedef typename K::Point_d Point_d; - - K k; - std::vector points = {0, 1, 2, 3}; - std::vector results; - - Gudhi::subsampling::choose_n_farthest_points(k, points, 2, - Gudhi::subsampling::random_starting_point, - std::back_inserter(results)); - std::clog << "Before sparsification: " << points.size() << " points.\n"; - std::clog << "After sparsification: " << results.size() << " points.\n"; - std::clog << "Result table: {" << results[0] << "," << results[1] << "}\n"; - - return 0; -} diff --git a/src/Subsampling/include/gudhi/choose_n_farthest_points.h b/src/Subsampling/include/gudhi/choose_n_farthest_points.h index b70af8a0..561dcf3e 100644 --- a/src/Subsampling/include/gudhi/choose_n_farthest_points.h +++ b/src/Subsampling/include/gudhi/choose_n_farthest_points.h @@ -38,33 +38,33 @@ enum : std::size_t { * \ingroup subsampling * \brief Subsample by a greedy strategy of iteratively adding the farthest point from the * current chosen point set to the subsampling. - * The iteration starts with the landmark `starting point` or, if `starting point==random_starting_point`, with a random landmark. - * \tparam Kernel must provide a type Kernel::Squared_distance_d which is a model of the - * concept Kernel_d::Squared_distance_d (despite the name, taken from CGAL, this can be any kind of metric or proximity measure). - * It must also contain a public member `squared_distance_d_object()` that returns an object of this type. - * \tparam Point_range Range whose value type is Kernel::Point_d. It must provide random-access - * via `operator[]` and the points should be stored contiguously in memory. - * \tparam PointOutputIterator Output iterator whose value type is Kernel::Point_d. + * The iteration starts with the landmark `starting point` or, if `starting point==random_starting_point`, + * with a random landmark. + * \tparam Distance must provide an operator() that takes 2 points (value type of the range) + * and returns their distance (or some more general proximity measure). + * \tparam Point_range Random access range of points. + * \tparam PointOutputIterator Output iterator whose value type is the point type. * \tparam DistanceOutputIterator Output iterator for distances. * \details It chooses `final_size` points from a random access range * `input_pts` (or the number of distinct points if `final_size` is larger) * and outputs them in the output iterator `output_it`. It also * outputs the distance from each of those points to the set of previous * points in `dist_it`. - * @param[in] k A kernel object. + * @param[in] dist A distance function. * @param[in] input_pts Const reference to the input points. * @param[in] final_size The size of the subsample to compute. * @param[in] starting_point The seed in the farthest point algorithm. * @param[out] output_it The output iterator for points. * @param[out] dist_it The optional output iterator for distances. + * + * \warning Older versions of this function took a CGAL kernel as argument. Users need to replace `k` with `k.squared_distance_d_object()` in the first argument of every call to `choose_n_farthest_points`. * */ -template < typename Kernel, +template < typename Distance, typename Point_range, typename PointOutputIterator, typename DistanceOutputIterator = Null_output_iterator> -void choose_n_farthest_points(Kernel const &k, +void choose_n_farthest_points(Distance dist, Point_range const &input_pts, std::size_t final_size, std::size_t starting_point, @@ -86,8 +86,6 @@ void choose_n_farthest_points(Kernel const &k, starting_point = dis(gen); } - typename Kernel::Squared_distance_d sqdist = k.squared_distance_d_object(); - std::size_t current_number_of_landmarks = 0; // counter for landmarks const double infty = std::numeric_limits::infinity(); // infinity (see next entry) std::vector< double > dist_to_L(nb_points, infty); // vector of current distances to L from input_pts @@ -100,7 +98,7 @@ void choose_n_farthest_points(Kernel const &k, *dist_it++ = dist_to_L[curr_max_w]; std::size_t i = 0; for (auto&& p : input_pts) { - double curr_dist = sqdist(p, input_pts[curr_max_w]); + double curr_dist = dist(p, input_pts[curr_max_w]); if (curr_dist < dist_to_L[i]) dist_to_L[i] = curr_dist; ++i; diff --git a/src/Subsampling/test/test_choose_n_farthest_points.cpp b/src/Subsampling/test/test_choose_n_farthest_points.cpp index b318d58e..94793295 100644 --- a/src/Subsampling/test/test_choose_n_farthest_points.cpp +++ b/src/Subsampling/test/test_choose_n_farthest_points.cpp @@ -44,7 +44,8 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_choose_farthest_point, Kernel, list_of_tested landmarks.clear(); Kernel k; - Gudhi::subsampling::choose_n_farthest_points(k, points, 100, Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks)); + auto d = k.squared_distance_d_object(); + Gudhi::subsampling::choose_n_farthest_points(d, points, 100, Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks)); BOOST_CHECK(landmarks.size() == 100); for (auto landmark : landmarks) @@ -61,32 +62,33 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_choose_farthest_point_limits, Kernel, list_of std::vector< FT > distances; landmarks.clear(); Kernel k; + auto d = k.squared_distance_d_object(); // Choose -1 farthest points in an empty point cloud - Gudhi::subsampling::choose_n_farthest_points(k, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); + Gudhi::subsampling::choose_n_farthest_points(d, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); BOOST_CHECK(landmarks.size() == 0); landmarks.clear(); distances.clear(); // Choose 0 farthest points in an empty point cloud - Gudhi::subsampling::choose_n_farthest_points(k, points, 0, -1, std::back_inserter(landmarks), std::back_inserter(distances)); + Gudhi::subsampling::choose_n_farthest_points(d, points, 0, -1, std::back_inserter(landmarks), std::back_inserter(distances)); BOOST_CHECK(landmarks.size() == 0); landmarks.clear(); distances.clear(); // Choose 1 farthest points in an empty point cloud - Gudhi::subsampling::choose_n_farthest_points(k, points, 1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); + Gudhi::subsampling::choose_n_farthest_points(d, points, 1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); BOOST_CHECK(landmarks.size() == 0); landmarks.clear(); distances.clear(); std::vector point({0.0, 0.0, 0.0, 0.0}); points.emplace_back(point.begin(), point.end()); // Choose -1 farthest points in a one point cloud - Gudhi::subsampling::choose_n_farthest_points(k, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); + Gudhi::subsampling::choose_n_farthest_points(d, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); BOOST_CHECK(landmarks.size() == 1 && distances.size() == 1); BOOST_CHECK(distances[0] == std::numeric_limits::infinity()); landmarks.clear(); distances.clear(); // Choose 0 farthest points in a one point cloud - Gudhi::subsampling::choose_n_farthest_points(k, points, 0, -1, std::back_inserter(landmarks), std::back_inserter(distances)); + Gudhi::subsampling::choose_n_farthest_points(d, points, 0, -1, std::back_inserter(landmarks), std::back_inserter(distances)); BOOST_CHECK(landmarks.size() == 0 && distances.size() == 0); landmarks.clear(); distances.clear(); // Choose 1 farthest points in a one point cloud - Gudhi::subsampling::choose_n_farthest_points(k, points, 1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); + Gudhi::subsampling::choose_n_farthest_points(d, points, 1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); BOOST_CHECK(landmarks.size() == 1 && distances.size() == 1); BOOST_CHECK(distances[0] == std::numeric_limits::infinity()); landmarks.clear(); distances.clear(); @@ -94,7 +96,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_choose_farthest_point_limits, Kernel, list_of std::vector point2({1.0, 0.0, 0.0, 0.0}); points.emplace_back(point2.begin(), point2.end()); // Choose all farthest points among 2 points - Gudhi::subsampling::choose_n_farthest_points(k, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); + Gudhi::subsampling::choose_n_farthest_points(d, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); BOOST_CHECK(landmarks.size() == 2 && distances.size() == 2); BOOST_CHECK(distances[0] == std::numeric_limits::infinity()); BOOST_CHECK(distances[1] == 1); @@ -102,7 +104,7 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(test_choose_farthest_point_limits, Kernel, list_of // Ignore duplicated points points.emplace_back(point.begin(), point.end()); - Gudhi::subsampling::choose_n_farthest_points(k, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); + Gudhi::subsampling::choose_n_farthest_points(d, points, -1, -1, std::back_inserter(landmarks), std::back_inserter(distances)); BOOST_CHECK(landmarks.size() == 2 && distances.size() == 2); BOOST_CHECK(distances[0] == std::numeric_limits::infinity()); BOOST_CHECK(distances[1] == 1); diff --git a/src/Witness_complex/doc/Witness_complex_doc.h b/src/Witness_complex/doc/Witness_complex_doc.h index 62203054..202f4539 100644 --- a/src/Witness_complex/doc/Witness_complex_doc.h +++ b/src/Witness_complex/doc/Witness_complex_doc.h @@ -92,11 +92,11 @@ int main(int argc, char * const argv[]) { // Choose landmarks (one can choose either of the two methods below) // Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks)); - Gudhi::subsampling::choose_n_farthest_points(K(), point_vector, nbL, Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks)); + Gudhi::subsampling::choose_n_farthest_points(K().squared_distance_d_object(), point_vector, nbL, + Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks)); // Compute witness complex - Witness_complex witness_complex(landmarks, - point_vector); + Witness_complex witness_complex(landmarks, point_vector); witness_complex.create_complex(simplex_tree, alpha2, lim_dim); } diff --git a/src/Witness_complex/example/example_strong_witness_complex_off.cpp b/src/Witness_complex/example/example_strong_witness_complex_off.cpp index 583a04ab..2bb135bf 100644 --- a/src/Witness_complex/example/example_strong_witness_complex_off.cpp +++ b/src/Witness_complex/example/example_strong_witness_complex_off.cpp @@ -43,7 +43,8 @@ int main(int argc, char* const argv[]) { // Choose landmarks (decomment one of the following two lines) // Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks)); - Gudhi::subsampling::choose_n_farthest_points(K(), point_vector, nbL, Gudhi::subsampling::random_starting_point, + Gudhi::subsampling::choose_n_farthest_points(K().squared_distance_d_object(), point_vector, + nbL, Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks)); // Compute witness complex diff --git a/src/Witness_complex/example/example_witness_complex_off.cpp b/src/Witness_complex/example/example_witness_complex_off.cpp index 3635da78..e1384c73 100644 --- a/src/Witness_complex/example/example_witness_complex_off.cpp +++ b/src/Witness_complex/example/example_witness_complex_off.cpp @@ -47,7 +47,8 @@ int main(int argc, char * const argv[]) { // Choose landmarks (decomment one of the following two lines) // Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks)); - Gudhi::subsampling::choose_n_farthest_points(K(), point_vector, nbL, Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks)); + Gudhi::subsampling::choose_n_farthest_points(K().squared_distance_d_object(), point_vector, nbL, + Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks)); // Compute witness complex start = clock(); diff --git a/src/Witness_complex/example/example_witness_complex_sphere.cpp b/src/Witness_complex/example/example_witness_complex_sphere.cpp index 78d5db4f..12a56de4 100644 --- a/src/Witness_complex/example/example_witness_complex_sphere.cpp +++ b/src/Witness_complex/example/example_witness_complex_sphere.cpp @@ -53,7 +53,7 @@ int main(int argc, char* const argv[]) { // Choose landmarks start = clock(); // Gudhi::subsampling::pick_n_random_points(point_vector, number_of_landmarks, std::back_inserter(landmarks)); - Gudhi::subsampling::choose_n_farthest_points(K(), point_vector, number_of_landmarks, + Gudhi::subsampling::choose_n_farthest_points(K().squared_distance_d_object(), point_vector, number_of_landmarks, Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks)); diff --git a/src/Witness_complex/utilities/strong_witness_persistence.cpp b/src/Witness_complex/utilities/strong_witness_persistence.cpp index 1f61c77c..614de0d4 100644 --- a/src/Witness_complex/utilities/strong_witness_persistence.cpp +++ b/src/Witness_complex/utilities/strong_witness_persistence.cpp @@ -61,7 +61,8 @@ int main(int argc, char* argv[]) { // Choose landmarks (decomment one of the following two lines) // Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks)); - Gudhi::subsampling::choose_n_farthest_points(K(), witnesses, nbL, Gudhi::subsampling::random_starting_point, + Gudhi::subsampling::choose_n_farthest_points(K().squared_distance_d_object(), witnesses, nbL, + Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks)); // Compute witness complex diff --git a/src/Witness_complex/utilities/weak_witness_persistence.cpp b/src/Witness_complex/utilities/weak_witness_persistence.cpp index 93050af5..5ea31d6b 100644 --- a/src/Witness_complex/utilities/weak_witness_persistence.cpp +++ b/src/Witness_complex/utilities/weak_witness_persistence.cpp @@ -61,7 +61,8 @@ int main(int argc, char* argv[]) { // Choose landmarks (decomment one of the following two lines) // Gudhi::subsampling::pick_n_random_points(point_vector, nbL, std::back_inserter(landmarks)); - Gudhi::subsampling::choose_n_farthest_points(K(), witnesses, nbL, Gudhi::subsampling::random_starting_point, + Gudhi::subsampling::choose_n_farthest_points(K().squared_distance_d_object(), witnesses, nbL, + Gudhi::subsampling::random_starting_point, std::back_inserter(landmarks)); // Compute witness complex diff --git a/src/common/doc/examples.h b/src/common/doc/examples.h index c19b3444..474f8699 100644 --- a/src/common/doc/examples.h +++ b/src/common/doc/examples.h @@ -42,7 +42,7 @@ * @example Persistence_representations/persistence_landscape.cpp * @example Tangential_complex/example_basic.cpp * @example Tangential_complex/example_with_perturb.cpp - * @example Subsampling/example_custom_kernel.cpp + * @example Subsampling/example_custom_distance.cpp * @example Subsampling/example_choose_n_farthest_points.cpp * @example Subsampling/example_sparsify_point_set.cpp * @example Subsampling/example_pick_n_random_points.cpp diff --git a/src/common/doc/installation.h b/src/common/doc/installation.h index 9df1c2f0..a6b9292b 100644 --- a/src/common/doc/installation.h +++ b/src/common/doc/installation.h @@ -113,8 +113,8 @@ make doxygen * Spatial_searching/example_spatial_searching.cpp * \li * Subsampling/example_choose_n_farthest_points.cpp - * \li - * Subsampling/example_custom_kernel.cpp + * \li + * Subsampling/example_custom_distance.cpp * \li * Subsampling/example_pick_n_random_points.cpp * \li @@ -153,8 +153,8 @@ make doxygen * Spatial_searching/example_spatial_searching.cpp * \li * Subsampling/example_choose_n_farthest_points.cpp - * \li - * Subsampling/example_custom_kernel.cpp + * \li + * Subsampling/example_custom_distance.cpp * \li * Subsampling/example_pick_n_random_points.cpp * \li diff --git a/src/python/include/Subsampling_interface.h b/src/python/include/Subsampling_interface.h index cdda851f..6aee7231 100644 --- a/src/python/include/Subsampling_interface.h +++ b/src/python/include/Subsampling_interface.h @@ -11,6 +11,7 @@ #ifndef INCLUDE_SUBSAMPLING_INTERFACE_H_ #define INCLUDE_SUBSAMPLING_INTERFACE_H_ +#include #include #include #include @@ -27,14 +28,13 @@ namespace subsampling { using Subsampling_dynamic_kernel = CGAL::Epick_d< CGAL::Dynamic_dimension_tag >; using Subsampling_point_d = Subsampling_dynamic_kernel::Point_d; -using Subsampling_ft = Subsampling_dynamic_kernel::FT; // ------ choose_n_farthest_points ------ std::vector> subsampling_n_farthest_points(const std::vector>& points, unsigned nb_points) { std::vector> landmarks; - Subsampling_dynamic_kernel k; - choose_n_farthest_points(k, points, nb_points, random_starting_point, std::back_inserter(landmarks)); + choose_n_farthest_points(Euclidean_distance(), points, nb_points, + random_starting_point, std::back_inserter(landmarks)); return landmarks; } @@ -42,8 +42,8 @@ std::vector> subsampling_n_farthest_points(const std::vector std::vector> subsampling_n_farthest_points(const std::vector>& points, unsigned nb_points, unsigned starting_point) { std::vector> landmarks; - Subsampling_dynamic_kernel k; - choose_n_farthest_points(k, points, nb_points, starting_point, std::back_inserter(landmarks)); + choose_n_farthest_points(Euclidean_distance(), points, nb_points, + starting_point, std::back_inserter(landmarks)); return landmarks; } -- cgit v1.2.3 From 6279ac91dd7e2d3206e8a380d38cb2e5d503e9dc Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Tue, 8 Dec 2020 17:54:46 +0100 Subject: Fix #388 for the tests to pass --- src/cmake/modules/GUDHI_third_party_libraries.cmake | 1 + src/python/gudhi/simplex_tree.pyx | 3 +++ src/python/include/Simplex_tree_interface.h | 6 ++++++ 3 files changed, 10 insertions(+) (limited to 'src/python/include') diff --git a/src/cmake/modules/GUDHI_third_party_libraries.cmake b/src/cmake/modules/GUDHI_third_party_libraries.cmake index e1566877..e2684aa4 100644 --- a/src/cmake/modules/GUDHI_third_party_libraries.cmake +++ b/src/cmake/modules/GUDHI_third_party_libraries.cmake @@ -58,6 +58,7 @@ endif(WITH_GUDHI_USE_TBB) set(CGAL_WITH_EIGEN3_VERSION 0.0.0) find_package(Eigen3 3.1.0) if (EIGEN3_FOUND) + add_definitions(-DGUDHI_USE_EIGEN3) include( ${EIGEN3_USE_FILE} ) set(CGAL_WITH_EIGEN3_VERSION ${CGAL_VERSION}) endif (EIGEN3_FOUND) diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index cdd2e87b..7d6ab89a 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -627,6 +627,9 @@ cdef class SimplexTree: :param nb_iterations: The number of edge collapse iterations to perform. Default is 1. :type nb_iterations: int + + :note: collapse_edges function requires `Eigen `_ >= 3.1.0, otherwise no action is + performed. """ # Backup old pointer cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr() diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index baff3850..2bd704b4 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -15,7 +15,9 @@ #include #include #include +#ifdef GUDHI_USE_EIGEN3 #include +#endif #include #include @@ -162,6 +164,7 @@ class Simplex_tree_interface : public Simplex_tree { } Simplex_tree_interface* collapse_edges(int nb_collapse_iteration) { +#ifdef GUDHI_USE_EIGEN3 using Filtered_edge = std::tuple; std::vector edges; for (Simplex_handle sh : Base::skeleton_simplex_range(1)) { @@ -187,6 +190,9 @@ class Simplex_tree_interface : public Simplex_tree { collapsed_stree_ptr->insert({std::get<0>(remaining_edge), std::get<1>(remaining_edge)}, std::get<2>(remaining_edge)); } return collapsed_stree_ptr; +#else + return this; +#endif } // Iterator over the simplex tree -- cgit v1.2.3 From fda0084941ece5d41a258d19ca4eb0b3d87384a4 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Wed, 9 Dec 2020 09:41:13 +0100 Subject: Fix #388 --- src/cmake/modules/GUDHI_third_party_libraries.cmake | 1 - src/python/CMakeLists.txt | 1 + src/python/gudhi/simplex_tree.pxd | 2 ++ src/python/gudhi/simplex_tree.pyx | 3 +++ src/python/include/Simplex_tree_interface.h | 10 +++++++++- src/python/test/test_simplex_tree.py | 14 +++++++++----- 6 files changed, 24 insertions(+), 7 deletions(-) (limited to 'src/python/include') diff --git a/src/cmake/modules/GUDHI_third_party_libraries.cmake b/src/cmake/modules/GUDHI_third_party_libraries.cmake index e2684aa4..e1566877 100644 --- a/src/cmake/modules/GUDHI_third_party_libraries.cmake +++ b/src/cmake/modules/GUDHI_third_party_libraries.cmake @@ -58,7 +58,6 @@ endif(WITH_GUDHI_USE_TBB) set(CGAL_WITH_EIGEN3_VERSION 0.0.0) find_package(Eigen3 3.1.0) if (EIGEN3_FOUND) - add_definitions(-DGUDHI_USE_EIGEN3) include( ${EIGEN3_USE_FILE} ) set(CGAL_WITH_EIGEN3_VERSION ${CGAL_VERSION}) endif (EIGEN3_FOUND) diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 2d5f793a..45c89609 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -133,6 +133,7 @@ if(PYTHONINTERP_FOUND) add_gudhi_debug_info("Eigen3 version ${EIGEN3_VERSION}") # No problem, even if no CGAL found set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_EIGEN3_ENABLED', ") + set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DGUDHI_USE_EIGEN3', ") endif (EIGEN3_FOUND) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'off_reader', ") diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 3c4cbed3..283830ff 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -74,6 +74,8 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": Simplex_tree_skeleton_iterator get_skeleton_iterator_begin(int dimension) nogil Simplex_tree_skeleton_iterator get_skeleton_iterator_end(int dimension) nogil pair[Simplex_tree_boundary_iterator, Simplex_tree_boundary_iterator] get_boundary_iterators(vector[int] simplex) nogil except + + + cdef int _GUDHI_USE_EIGEN3 cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface>": diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index 7d6ab89a..eb44c0fc 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -17,6 +17,9 @@ __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" __license__ = "MIT" +# For unitary tests purpose +__GUDHI_USE_EIGEN3 = _GUDHI_USE_EIGEN3 + # SimplexTree python interface cdef class SimplexTree: """The simplex tree is an efficient and flexible data structure for diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 2bd704b4..50592e25 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -27,6 +27,12 @@ namespace Gudhi { +#ifdef GUDHI_USE_EIGEN3 +const int _GUDHI_USE_EIGEN3 = 1; +#else +const int _GUDHI_USE_EIGEN3 = 0; +#endif + template class Simplex_tree_interface : public Simplex_tree { public: @@ -191,7 +197,9 @@ class Simplex_tree_interface : public Simplex_tree { } return collapsed_stree_ptr; #else - return this; + // If no Eigen3, return a copy, as it will be deleted in pyx + Simplex_tree_interface* collapsed_stree_ptr = new Simplex_tree_interface(*this); + return collapsed_stree_ptr; #endif } diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 3b23fa0b..15b472ee 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -8,7 +8,7 @@ - YYYY/MM Author: Description of the modification """ -from gudhi import SimplexTree +from gudhi import SimplexTree, simplex_tree import pytest __author__ = "Vincent Rouvreau" @@ -353,11 +353,15 @@ def test_collapse_edges(): assert st.num_simplices() == 10 + # If no Eigen3, collapse_edges just return a copy, no action. Maybe it would require some user warning st.collapse_edges() - assert st.num_simplices() == 9 - assert st.find([1, 3]) == False - for simplex in st.get_skeleton(0): - assert simplex[1] == 1. + if simplex_tree.__GUDHI_USE_EIGEN3: + assert st.num_simplices() == 9 + assert st.find([1, 3]) == False + for simplex in st.get_skeleton(0): + assert simplex[1] == 1. + else: + assert st.num_simplices() == 10 def test_reset_filtration(): st = SimplexTree() -- cgit v1.2.3 From 315b0e1b4301ffab98bbc36d064344fba72076c5 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Thu, 10 Dec 2020 09:24:05 +0100 Subject: Fix memory leak --- src/python/include/Alpha_complex_factory.h | 9 ++++++--- 1 file changed, 6 insertions(+), 3 deletions(-) (limited to 'src/python/include') diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h index d699ad9b..3405fdd6 100644 --- a/src/python/include/Alpha_complex_factory.h +++ b/src/python/include/Alpha_complex_factory.h @@ -48,11 +48,14 @@ static CgalPointType pt_cython_to_cgal(std::vector const& vec) { class Abstract_alpha_complex { public: virtual std::vector get_point(int vh) = 0; + virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) = 0; + + virtual ~Abstract_alpha_complex() = default; }; -class Exact_Alphacomplex_dD : public Abstract_alpha_complex { +class Exact_Alphacomplex_dD final : public Abstract_alpha_complex { private: using Kernel = CGAL::Epeck_d; using Point = typename Kernel::Point_d; @@ -78,7 +81,7 @@ class Exact_Alphacomplex_dD : public Abstract_alpha_complex { Alpha_complex alpha_complex_; }; -class Inexact_Alphacomplex_dD : public Abstract_alpha_complex { +class Inexact_Alphacomplex_dD final : public Abstract_alpha_complex { private: using Kernel = CGAL::Epick_d; using Point = typename Kernel::Point_d; @@ -104,7 +107,7 @@ class Inexact_Alphacomplex_dD : public Abstract_alpha_complex { }; template -class Alphacomplex_3D : public Abstract_alpha_complex { +class Alphacomplex_3D final : public Abstract_alpha_complex { private: using Point = typename Alpha_complex_3d::Bare_point_3; -- cgit v1.2.3 From 957da77f9484972ce34d0415502887f92080878e Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Fri, 11 Dec 2020 09:38:12 +0100 Subject: code review: GUDHI_USE_EIGEN3 generated by CMake in __init__.py as suggested and roll back the other solution --- src/python/CMakeLists.txt | 2 ++ src/python/gudhi/__init__.py.in | 4 ++++ src/python/gudhi/simplex_tree.pxd | 2 -- src/python/gudhi/simplex_tree.pyx | 3 --- src/python/include/Simplex_tree_interface.h | 6 ------ src/python/test/test_simplex_tree.py | 4 ++-- 6 files changed, 8 insertions(+), 13 deletions(-) (limited to 'src/python/include') diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 45c89609..5a245aac 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -129,11 +129,13 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DDEBUG_TRACES', ") endif() + set(GUDHI_USE_EIGEN3 "False") if (EIGEN3_FOUND) add_gudhi_debug_info("Eigen3 version ${EIGEN3_VERSION}") # No problem, even if no CGAL found set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DCGAL_EIGEN3_ENABLED', ") set(GUDHI_PYTHON_EXTRA_COMPILE_ARGS "${GUDHI_PYTHON_EXTRA_COMPILE_ARGS}'-DGUDHI_USE_EIGEN3', ") + set(GUDHI_USE_EIGEN3 "True") endif (EIGEN3_FOUND) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'off_reader', ") diff --git a/src/python/gudhi/__init__.py.in b/src/python/gudhi/__init__.py.in index 79e12fbc..3043201a 100644 --- a/src/python/gudhi/__init__.py.in +++ b/src/python/gudhi/__init__.py.in @@ -23,6 +23,10 @@ __all__ = [@GUDHI_PYTHON_MODULES@ @GUDHI_PYTHON_MODULES_EXTRA@] __available_modules = '' __missing_modules = '' +# For unitary tests purpose +# could use "if 'collapse_edges' in gudhi.__all__" when collapse edges will have a python module +__GUDHI_USE_EIGEN3 = @GUDHI_USE_EIGEN3@ + # Try to import * from gudhi.__module_name for default modules. # Extra modules require an explicit import by the user (mostly because of # unusual dependencies, but also to avoid cluttering namespace gudhi and diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 283830ff..3c4cbed3 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -74,8 +74,6 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": Simplex_tree_skeleton_iterator get_skeleton_iterator_begin(int dimension) nogil Simplex_tree_skeleton_iterator get_skeleton_iterator_end(int dimension) nogil pair[Simplex_tree_boundary_iterator, Simplex_tree_boundary_iterator] get_boundary_iterators(vector[int] simplex) nogil except + - - cdef int _GUDHI_USE_EIGEN3 cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi": cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface>": diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index eb44c0fc..7d6ab89a 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -17,9 +17,6 @@ __author__ = "Vincent Rouvreau" __copyright__ = "Copyright (C) 2016 Inria" __license__ = "MIT" -# For unitary tests purpose -__GUDHI_USE_EIGEN3 = _GUDHI_USE_EIGEN3 - # SimplexTree python interface cdef class SimplexTree: """The simplex tree is an efficient and flexible data structure for diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 50592e25..82444609 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -27,12 +27,6 @@ namespace Gudhi { -#ifdef GUDHI_USE_EIGEN3 -const int _GUDHI_USE_EIGEN3 = 1; -#else -const int _GUDHI_USE_EIGEN3 = 0; -#endif - template class Simplex_tree_interface : public Simplex_tree { public: diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 15b472ee..0c072baa 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -8,7 +8,7 @@ - YYYY/MM Author: Description of the modification """ -from gudhi import SimplexTree, simplex_tree +from gudhi import SimplexTree, __GUDHI_USE_EIGEN3 import pytest __author__ = "Vincent Rouvreau" @@ -355,7 +355,7 @@ def test_collapse_edges(): # If no Eigen3, collapse_edges just return a copy, no action. Maybe it would require some user warning st.collapse_edges() - if simplex_tree.__GUDHI_USE_EIGEN3: + if __GUDHI_USE_EIGEN3: assert st.num_simplices() == 9 assert st.find([1, 3]) == False for simplex in st.get_skeleton(0): -- cgit v1.2.3 From 40e0976e8dae27f0e30f3c9df1fd7de1a7343948 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Fri, 11 Dec 2020 10:31:09 +0100 Subject: code review: throw an eception if collapse_edges when no Eigen3 --- src/python/gudhi/simplex_tree.pxd | 2 +- src/python/gudhi/simplex_tree.pyx | 4 ++-- src/python/include/Simplex_tree_interface.h | 4 +--- src/python/test/test_simplex_tree.py | 7 ++++--- 4 files changed, 8 insertions(+), 9 deletions(-) (limited to 'src/python/include') diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd index 3c4cbed3..000323af 100644 --- a/src/python/gudhi/simplex_tree.pxd +++ b/src/python/gudhi/simplex_tree.pxd @@ -63,7 +63,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi": bool make_filtration_non_decreasing() nogil void compute_extended_filtration() nogil vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(vector[pair[int, pair[double, double]]] dgm, double min_persistence) nogil - Simplex_tree_interface_full_featured* collapse_edges(int nb_collapse_iteration) nogil + Simplex_tree_interface_full_featured* collapse_edges(int nb_collapse_iteration) nogil except + void reset_filtration(double filtration, int dimension) nogil # Iterators over Simplex tree pair[vector[int], double] get_simplex_and_filtration(Simplex_tree_simplex_handle f_simplex) nogil diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx index 7d6ab89a..665d41e6 100644 --- a/src/python/gudhi/simplex_tree.pyx +++ b/src/python/gudhi/simplex_tree.pyx @@ -628,8 +628,8 @@ cdef class SimplexTree: :param nb_iterations: The number of edge collapse iterations to perform. Default is 1. :type nb_iterations: int - :note: collapse_edges function requires `Eigen `_ >= 3.1.0, otherwise no action is - performed. + :note: collapse_edges method requires `Eigen `_ >= 3.1.0 and an exception is thrown + if this method is not available. """ # Backup old pointer cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr() diff --git a/src/python/include/Simplex_tree_interface.h b/src/python/include/Simplex_tree_interface.h index 82444609..629f6083 100644 --- a/src/python/include/Simplex_tree_interface.h +++ b/src/python/include/Simplex_tree_interface.h @@ -191,9 +191,7 @@ class Simplex_tree_interface : public Simplex_tree { } return collapsed_stree_ptr; #else - // If no Eigen3, return a copy, as it will be deleted in pyx - Simplex_tree_interface* collapsed_stree_ptr = new Simplex_tree_interface(*this); - return collapsed_stree_ptr; + throw std::runtime_error("Unable to collapse edges as it requires Eigen3 >= 3.1.0."); #endif } diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py index 0c072baa..a3eacaa9 100755 --- a/src/python/test/test_simplex_tree.py +++ b/src/python/test/test_simplex_tree.py @@ -353,15 +353,16 @@ def test_collapse_edges(): assert st.num_simplices() == 10 - # If no Eigen3, collapse_edges just return a copy, no action. Maybe it would require some user warning - st.collapse_edges() if __GUDHI_USE_EIGEN3: + st.collapse_edges() assert st.num_simplices() == 9 assert st.find([1, 3]) == False for simplex in st.get_skeleton(0): assert simplex[1] == 1. else: - assert st.num_simplices() == 10 + # If no Eigen3, collapse_edges throws an exception + with pytest.raises(RuntimeError): + st.collapse_edges() def test_reset_filtration(): st = SimplexTree() -- cgit v1.2.3 From af98f16e12ec9d1af7d925ecdc53b4daefea6ebe Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Thu, 18 Mar 2021 16:16:48 +0100 Subject: Add weight support --- src/python/doc/alpha_complex_user.rst | 104 ++++++++++++++++++++++----- src/python/gudhi/alpha_complex.pyx | 48 ++++++++----- src/python/include/Alpha_complex_factory.h | 65 +++++++++++++++-- src/python/include/Alpha_complex_interface.h | 28 +++++--- 4 files changed, 194 insertions(+), 51 deletions(-) (limited to 'src/python/include') diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst index b116944c..6f35cc15 100644 --- a/src/python/doc/alpha_complex_user.rst +++ b/src/python/doc/alpha_complex_user.rst @@ -44,16 +44,15 @@ This example builds the alpha-complex from the given points: .. testcode:: - import gudhi - alpha_complex = gudhi.AlphaComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]]) - - simplex_tree = alpha_complex.create_simplex_tree() - result_str = 'Alpha complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \ - repr(simplex_tree.num_simplices()) + ' simplices - ' + \ - repr(simplex_tree.num_vertices()) + ' vertices.' - print(result_str) + from gudhi import AlphaComplex + ac = AlphaComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]]) + + stree = ac.create_simplex_tree() + print('Alpha complex is of dimension ', stree.dimension(), ' - ', + stree.num_simplices(), ' simplices - ', stree.num_vertices(), ' vertices.') + fmt = '%s -> %.2f' - for filtered_value in simplex_tree.get_filtration(): + for filtered_value in stree.get_filtration(): print(fmt % tuple(filtered_value)) The output is: @@ -174,6 +173,78 @@ of speed-up, since we always first build the full filtered complex, so it is rec :paramref:`~gudhi.AlphaComplex.create_simplex_tree.max_alpha_square`. In the following example, a threshold of :math:`\alpha^2 = 32.0` is used. +Weighted specific version +^^^^^^^^^^^^^^^^^^^^^^^^^ + +:Requires: `Eigen `_ :math:`\geq` 3.1.0 and `CGAL `_ :math:`\geq` 5.1.0. + +A weighted version for Alpha complex is available. It is like a usual Alpha complex, but based on a +`CGAL regular triangulation `_ +of Delaunay. + +This example builds the CGAL weighted alpha shapes from a small molecule, and initializes the alpha complex with +it. This example is taken from +`CGAL 3d weighted alpha shapes `_. + +Then, it is asked to display information about the alpha complex. + +.. testcode:: + + from gudhi import AlphaComplex + wgt_ac = AlphaComplex(weighted_points=[[[ 1., -1., -1.], 4.], + [[-1., 1., -1.], 4.], + [[-1., -1., 1.], 4.], + [[ 1., 1., 1.], 4.], + [[ 2., 2., 2.], 1.]]) + # equivalent to: + # wgt_ac = AlphaComplex(points=[[ 1., -1., -1.], + # [-1., 1., -1.], + # [-1., -1., 1.], + # [ 1., 1., 1.], + # [ 2., 2., 2.]], + # weights = [4., 4., 4., 4., 1.]) + + stree = wgt_ac.create_simplex_tree() + print('Weighted alpha complex is of dimension ', stree.dimension(), ' - ', + stree.num_simplices(), ' simplices - ', stree.num_vertices(), ' vertices.') + fmt = '%s -> %.2f' + for filtered_value in stree.get_filtration(): + print(fmt % tuple(filtered_value)) + +The output is: + +.. testoutput:: + + Weighted alpha complex is of dimension 3 - 29 simplices - 5 vertices. + [ 0 ] -> [-4] + [ 1 ] -> [-4] + [ 2 ] -> [-4] + [ 3 ] -> [-4] + [ 1, 0 ] -> [-2] + [ 2, 0 ] -> [-2] + [ 2, 1 ] -> [-2] + [ 3, 0 ] -> [-2] + [ 3, 1 ] -> [-2] + [ 3, 2 ] -> [-2] + [ 2, 1, 0 ] -> [-1.33333] + [ 3, 1, 0 ] -> [-1.33333] + [ 3, 2, 0 ] -> [-1.33333] + [ 3, 2, 1 ] -> [-1.33333] + [ 3, 2, 1, 0 ] -> [-1] + [ 4 ] -> [-1] + [ 4, 2 ] -> [-1] + [ 4, 0 ] -> [23] + [ 4, 1 ] -> [23] + [ 4, 2, 0 ] -> [23] + [ 4, 2, 1 ] -> [23] + [ 4, 3 ] -> [23] + [ 4, 3, 2 ] -> [23] + [ 4, 1, 0 ] -> [95] + [ 4, 2, 1, 0 ] -> [95] + [ 4, 3, 0 ] -> [95] + [ 4, 3, 1 ] -> [95] + [ 4, 3, 2, 0 ] -> [95] + [ 4, 3, 2, 1 ] -> [95] Example from OFF file ^^^^^^^^^^^^^^^^^^^^^ @@ -186,14 +257,9 @@ Then, it computes the persistence diagram and displays it: :include-source: import matplotlib.pyplot as plt - import gudhi - alpha_complex = gudhi.AlphaComplex(off_file=gudhi.__root_source_dir__ + \ - '/data/points/tore3D_300.off') - simplex_tree = alpha_complex.create_simplex_tree() - result_str = 'Alpha complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \ - repr(simplex_tree.num_simplices()) + ' simplices - ' + \ - repr(simplex_tree.num_vertices()) + ' vertices.' - print(result_str) - diag = simplex_tree.persistence() - gudhi.plot_persistence_diagram(diag) + import gudhi as gd + ac = gd.AlphaComplex(off_file=gd.__root_source_dir__ + '/data/points/tore3D_300.off') + stree = ac.create_simplex_tree() + dgm = stree.persistence() + gd.plot_persistence_diagram(dgm, legend = True) plt.show() diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx index b7c20f74..681faebe 100644 --- a/src/python/gudhi/alpha_complex.pyx +++ b/src/python/gudhi/alpha_complex.pyx @@ -29,7 +29,7 @@ __license__ = "GPL v3" cdef extern from "Alpha_complex_interface.h" namespace "Gudhi": cdef cppclass Alpha_complex_interface "Gudhi::alpha_complex::Alpha_complex_interface": - Alpha_complex_interface(vector[vector[double]] points, bool fast_version, bool exact_version) nogil except + + Alpha_complex_interface(vector[vector[double]] points, vector[double] weights, bool fast_version, bool exact_version) nogil except + vector[double] get_point(int vertex) nogil except + void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool default_filtration_value) nogil except + @@ -56,26 +56,35 @@ cdef class AlphaComplex: cdef Alpha_complex_interface * this_ptr # Fake constructor that does nothing but documenting the constructor - def __init__(self, points=[], off_file='', weights=[], weight_file='', precision='safe'): + def __init__(self, points=[], off_file='', weights=[], weight_file='', weighted_points=[], + precision='safe'): """AlphaComplex constructor. :param points: A list of points in d-Dimension. - :type points: list of list of double + :type points: Iterable[Iterable[float]] - :param off_file: An `OFF file style `_ name. `points` are - read and overwritten by the points in the `off_file`. + :param off_file: An `OFF file style `_ name. + If an `off_file` is given with `points` or `weighted_points`, only points from the + file are taken into account. :type off_file: string :param weights: A list of weights. If set, the number of weights must correspond to the number of points. - :type weights: list of double + :type weights: Iterable[float] :param weight_file: A file containing a list of weights (one per line). - `weights` are read and overwritten by the weights in the `weight_file`. - If set, the number of weights must correspond to the number of points. + If a `weight_file` is given with `weights` or `weighted_points`, only weights from the + file are taken into account. + :type weight_file: string - :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. + :param weighted_points: A list of points in d-Dimension and its weight. + If `weighted_points` are given with `weights` or `points`, these last ones will + not be taken into account. + :type weighted_points: Iterable[Iterable[float], float] + + :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is + 'safe'. :type precision: string :raises FileNotFoundError: If `off_file` and/or `weight_file` is set but not found. @@ -83,11 +92,16 @@ cdef class AlphaComplex: """ # The real cython constructor - def __cinit__(self, points = [], off_file = '', weights=[], weight_file='', precision = 'safe'): + def __cinit__(self, points = [], off_file = '', weights=[], weight_file='', weighted_points=[], + precision = 'safe'): assert precision in ['fast', 'safe', 'exact'], "Alpha complex precision can only be 'fast', 'safe' or 'exact'" cdef bool fast = precision == 'fast' cdef bool exact = precision == 'exact' + if len(weighted_points) > 0: + points = [wpt[0] for wpt in weighted_points] + weights = [wpt[1] for wpt in weighted_points] + if off_file: if os.path.isfile(off_file): points = read_points_from_off_file(off_file = off_file) @@ -100,20 +114,18 @@ cdef class AlphaComplex: else: raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), weight_file) + # weights are set but is inconsistent with the number of points + if len(weights) != 0 and len(weights) != len(points): + raise ValueError("Inconsistency between the number of points and weights") + # need to copy the points to use them without the gil cdef vector[vector[double]] pts cdef vector[double] wgts pts = points + wgts = weights if len(weights) == 0: with nogil: - self.this_ptr = new Alpha_complex_interface(pts, fast, exact) - else: - if len(weights) == len(points): - wgts = weights - with nogil: - self.this_ptr = new Alpha_complex_interface(pts, fast, exact) - else: - raise ValueError("Inconsistency between the number of points and weights") + self.this_ptr = new Alpha_complex_interface(pts, wgts, fast, exact) def __dealloc__(self): if self.this_ptr != NULL: diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h index 3405fdd6..36e98615 100644 --- a/src/python/include/Alpha_complex_factory.h +++ b/src/python/include/Alpha_complex_factory.h @@ -55,13 +55,13 @@ class Abstract_alpha_complex { virtual ~Abstract_alpha_complex() = default; }; -class Exact_Alphacomplex_dD final : public Abstract_alpha_complex { +class Exact_alpha_complex_dD final : public Abstract_alpha_complex { private: using Kernel = CGAL::Epeck_d; using Point = typename Kernel::Point_d; public: - Exact_Alphacomplex_dD(const std::vector>& points, bool exact_version) + Exact_alpha_complex_dD(const std::vector>& points, bool exact_version) : exact_version_(exact_version), alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal)) { } @@ -81,13 +81,13 @@ class Exact_Alphacomplex_dD final : public Abstract_alpha_complex { Alpha_complex alpha_complex_; }; -class Inexact_Alphacomplex_dD final : public Abstract_alpha_complex { +class Inexact_alpha_complex_dD final : public Abstract_alpha_complex { private: using Kernel = CGAL::Epick_d; using Point = typename Kernel::Point_d; public: - Inexact_Alphacomplex_dD(const std::vector>& points, bool exact_version) + Inexact_alpha_complex_dD(const std::vector>& points, bool exact_version) : exact_version_(exact_version), alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal)) { } @@ -106,8 +106,61 @@ class Inexact_Alphacomplex_dD final : public Abstract_alpha_complex { Alpha_complex alpha_complex_; }; +class Exact_weighted_alpha_complex_dD final : public Abstract_alpha_complex { + private: + using Kernel = CGAL::Epeck_d; + using Point = typename Kernel::Point_d; + + public: + Exact_weighted_alpha_complex_dD(const std::vector>& points, + const std::vector& weights, bool exact_version) + : exact_version_(exact_version), + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal), weights) { + } + + virtual std::vector get_point(int vh) override { + Point const& point = alpha_complex_.get_point(vh).point(); + return pt_cgal_to_cython(point); + } + + virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, + bool default_filtration_value) override { + return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value); + } + + private: + bool exact_version_; + Alpha_complex alpha_complex_; +}; + +class Inexact_weighted_alpha_complex_dD final : public Abstract_alpha_complex { + private: + using Kernel = CGAL::Epick_d; + using Point = typename Kernel::Point_d; + + public: + Inexact_weighted_alpha_complex_dD(const std::vector>& points, + const std::vector& weights, bool exact_version) + : exact_version_(exact_version), + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal), weights) { + } + + virtual std::vector get_point(int vh) override { + Point const& point = alpha_complex_.get_point(vh).point(); + return pt_cgal_to_cython(point); + } + virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, + bool default_filtration_value) override { + return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value); + } + + private: + bool exact_version_; + Alpha_complex alpha_complex_; +}; + template -class Alphacomplex_3D final : public Abstract_alpha_complex { +class Alpha_complex_3D final : public Abstract_alpha_complex { private: using Point = typename Alpha_complex_3d::Bare_point_3; @@ -116,7 +169,7 @@ class Alphacomplex_3D final : public Abstract_alpha_complex { } public: - Alphacomplex_3D(const std::vector>& points) + Alpha_complex_3D(const std::vector>& points) : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3)) { } diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index 23be194d..43c96b2f 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -27,8 +27,11 @@ namespace alpha_complex { class Alpha_complex_interface { public: - Alpha_complex_interface(const std::vector>& points, bool fast_version, bool exact_version) + Alpha_complex_interface(const std::vector>& points, + const std::vector& weights, + bool fast_version, bool exact_version) : points_(points), + weights_(weights), fast_version_(fast_version), exact_version_(exact_version) { } @@ -41,13 +44,13 @@ class Alpha_complex_interface { bool default_filtration_value) { if (points_.size() > 0) { std::size_t dimension = points_[0].size(); - if (dimension == 3 && !default_filtration_value) { + if (dimension == 3 && weights_.size() == 0 && !default_filtration_value) { if (fast_version_) - alpha_ptr_ = std::make_unique>(points_); + alpha_ptr_ = std::make_unique>(points_); else if (exact_version_) - alpha_ptr_ = std::make_unique>(points_); + alpha_ptr_ = std::make_unique>(points_); else - alpha_ptr_ = std::make_unique>(points_); + alpha_ptr_ = std::make_unique>(points_); if (!alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value)) { // create_simplex_tree will fail if all points are on a plane - Retry with dD by setting dimension to 2 dimension--; @@ -55,11 +58,19 @@ class Alpha_complex_interface { } } // Not ** else ** because we have to take into account if 3d fails - if (dimension != 3 || default_filtration_value) { + if (dimension != 3 || weights_.size() != 0 || default_filtration_value) { if (fast_version_) { - alpha_ptr_ = std::make_unique(points_, exact_version_); + if (weights_.size() == 0) { + alpha_ptr_ = std::make_unique(points_, exact_version_); + } else { + alpha_ptr_ = std::make_unique(points_, weights_, exact_version_); + } } else { - alpha_ptr_ = std::make_unique(points_, exact_version_); + if (weights_.size() == 0) { + alpha_ptr_ = std::make_unique(points_, exact_version_); + } else { + alpha_ptr_ = std::make_unique(points_, weights_, exact_version_); + } } alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value); } @@ -69,6 +80,7 @@ class Alpha_complex_interface { private: std::unique_ptr alpha_ptr_; std::vector> points_; + std::vector weights_; bool fast_version_; bool exact_version_; }; -- cgit v1.2.3 From e4381a3e2ad79d3150cd03704bef3fc006e7c54b Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Sat, 3 Apr 2021 10:27:08 +0200 Subject: Python alpha complex specific 3d with weighted version and functor to get points --- src/python/CMakeLists.txt | 2 + src/python/gudhi/alpha_complex_3d.pyx | 129 ++++++++++++++++++++++++ src/python/include/Alpha_complex_factory.h | 48 +++++++-- src/python/include/Alpha_complex_interface.h | 59 ++++------- src/python/include/Alpha_complex_interface_3d.h | 71 +++++++++++++ 5 files changed, 261 insertions(+), 48 deletions(-) create mode 100644 src/python/gudhi/alpha_complex_3d.pyx create mode 100644 src/python/include/Alpha_complex_interface_3d.h (limited to 'src/python/include') diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 73303a24..307181b7 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -61,6 +61,7 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'subsampling', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'tangential_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'alpha_complex', ") + set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'alpha_complex_3d', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'euclidean_witness_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'euclidean_strong_witness_complex', ") # Modules that should not be auto-imported in __init__.py @@ -156,6 +157,7 @@ if(PYTHONINTERP_FOUND) endif () if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ") + set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex_3d', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'subsampling', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'tangential_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'euclidean_witness_complex', ") diff --git a/src/python/gudhi/alpha_complex_3d.pyx b/src/python/gudhi/alpha_complex_3d.pyx new file mode 100644 index 00000000..3959004a --- /dev/null +++ b/src/python/gudhi/alpha_complex_3d.pyx @@ -0,0 +1,129 @@ +# 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): Vincent Rouvreau +# +# Copyright (C) 2021 Inria +# +# Modification(s): +# - YYYY/MM Author: Description of the modification + +from __future__ import print_function +from cython cimport numeric +from libcpp.vector cimport vector +from libcpp.utility cimport pair +from libcpp.string cimport string +from libcpp cimport bool +from libc.stdint cimport intptr_t +import errno +import os +import warnings + +from gudhi.simplex_tree cimport * +from gudhi.simplex_tree import SimplexTree +from gudhi import read_points_from_off_file + +__author__ = "Vincent Rouvreau" +__copyright__ = "Copyright (C) 2021 Inria" +__license__ = "GPL v3" + +cdef extern from "Alpha_complex_interface_3d.h" namespace "Gudhi": + cdef cppclass Alpha_complex_interface_3d "Gudhi::alpha_complex::Alpha_complex_interface_3d": + Alpha_complex_interface_3d(vector[vector[double]] points, vector[double] weights, bool fast_version, bool exact_version) nogil except + + vector[double] get_point(int vertex) nogil except + + void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square) nogil except + + +# AlphaComplex3D python interface +cdef class AlphaComplex3D: + """AlphaComplex3D is a simplicial complex constructed from the finite cells + of a Delaunay Triangulation. + + The filtration value of each simplex is computed as the square of the + circumradius of the simplex if the circumsphere is empty (the simplex is + then said to be Gabriel), and as the minimum of the filtration values of + the codimension 1 cofaces that make it not Gabriel otherwise. + + All simplices that have a filtration value strictly greater than a given + alpha squared value are not inserted into the complex. + + .. note:: + + When AlphaComplex3D is constructed with an infinite value of alpha, the + complex is a Delaunay complex. + + """ + + cdef Alpha_complex_interface_3d * this_ptr + + # Fake constructor that does nothing but documenting the constructor + def __init__(self, points=[], weights=[], precision='safe'): + """AlphaComplex3D constructor. + + :param points: A list of points in d-Dimension. + :type points: Iterable[Iterable[float]] + + :param weights: A list of weights. If set, the number of weights must correspond to the + number of points. + :type weights: Iterable[float] + + :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is + 'safe'. + :type precision: string + + :raises ValueError: In case of inconsistency between the number of points and weights. + """ + + # The real cython constructor + def __cinit__(self, points = [], weights=[], precision = 'safe'): + assert precision in ['fast', 'safe', 'exact'], \ + "Alpha complex precision can only be 'fast', 'safe' or 'exact'" + cdef bool fast = precision == 'fast' + cdef bool exact = precision == 'exact' + + # weights are set but is inconsistent with the number of points + if len(weights) != 0 and len(weights) != len(points): + raise ValueError("Inconsistency between the number of points and weights") + + # need to copy the points to use them without the gil + cdef vector[vector[double]] pts + cdef vector[double] wgts + pts = points + wgts = weights + with nogil: + self.this_ptr = new Alpha_complex_interface_3d(pts, wgts, fast, exact) + + def __dealloc__(self): + if self.this_ptr != NULL: + del self.this_ptr + + def __is_defined(self): + """Returns true if AlphaComplex3D pointer is not NULL. + """ + return self.this_ptr != NULL + + def get_point(self, vertex): + """This function returns the point corresponding to a given vertex from the :class:`~gudhi.SimplexTree`. + + :param vertex: The vertex. + :type vertex: int + :rtype: list of float + :returns: the point. + """ + return self.this_ptr.get_point(vertex) + + def create_simplex_tree(self, max_alpha_square = float('inf')): + """ + :param max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to + infinity, and there is very little point using anything else since it does not save time. + :type max_alpha_square: float + :returns: A simplex tree created from the Delaunay Triangulation. + :rtype: SimplexTree + """ + stree = SimplexTree() + cdef double mas = max_alpha_square + cdef intptr_t stree_int_ptr=stree.thisptr + with nogil: + self.this_ptr.create_simplex_tree(stree_int_ptr, + mas) + return stree diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h index 36e98615..5d3bfb65 100644 --- a/src/python/include/Alpha_complex_factory.h +++ b/src/python/include/Alpha_complex_factory.h @@ -31,6 +31,34 @@ namespace Gudhi { namespace alpha_complex { +template +struct Point_cgal_to_cython; + +template +struct Point_cgal_to_cython { + std::vector operator()(CgalPointType const& point) const + { + std::vector vd; + vd.reserve(point.dimension()); + for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++) + vd.push_back(CGAL::to_double(*coord)); + return vd; + } +}; + +template +struct Point_cgal_to_cython { + std::vector operator()(CgalPointType const& weighted_point) const + { + auto point = weighted_point.point(); + std::vector vd; + vd.reserve(point.dimension()); + for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++) + vd.push_back(CGAL::to_double(*coord)); + return vd; + } +}; + template std::vector pt_cgal_to_cython(CgalPointType const& point) { std::vector vd; @@ -159,13 +187,14 @@ class Inexact_weighted_alpha_complex_dD final : public Abstract_alpha_complex { Alpha_complex alpha_complex_; }; -template +template class Alpha_complex_3D final : public Abstract_alpha_complex { private: - using Point = typename Alpha_complex_3d::Bare_point_3; + using Bare_point = typename Alpha_complex_3d::Bare_point_3; + using Point = typename Alpha_complex_3d::Point_3; - static Point pt_cython_to_cgal_3(std::vector const& vec) { - return Point(vec[0], vec[1], vec[2]); + static Bare_point pt_cython_to_cgal_3(std::vector const& vec) { + return Bare_point(vec[0], vec[1], vec[2]); } public: @@ -173,18 +202,23 @@ class Alpha_complex_3D final : public Abstract_alpha_complex { : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3)) { } + Alpha_complex_3D(const std::vector>& points, const std::vector& weights) + : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3), weights) { + } + virtual std::vector get_point(int vh) override { Point const& point = alpha_complex_.get_point(vh); - return pt_cgal_to_cython(point); + return Point_cgal_to_cython()(point); } virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) override { - return alpha_complex_.create_complex(*simplex_tree, max_alpha_square); + alpha_complex_.create_complex(*simplex_tree, max_alpha_square); + return true; } private: - Alpha_complex_3d alpha_complex_; + Alpha_complex_3d alpha_complex_; }; diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index 43c96b2f..31a8147b 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -30,10 +30,20 @@ class Alpha_complex_interface { Alpha_complex_interface(const std::vector>& points, const std::vector& weights, bool fast_version, bool exact_version) - : points_(points), - weights_(weights), - fast_version_(fast_version), - exact_version_(exact_version) { + : empty_point_set_(points.size() == 0) { + if (fast_version) { + if (weights.size() == 0) { + alpha_ptr_ = std::make_unique(points, exact_version); + } else { + alpha_ptr_ = std::make_unique(points, weights, exact_version); + } + } else { + if (weights.size() == 0) { + alpha_ptr_ = std::make_unique(points, exact_version); + } else { + alpha_ptr_ = std::make_unique(points, weights, exact_version); + } + } } std::vector get_point(int vh) { @@ -42,47 +52,14 @@ class Alpha_complex_interface { void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) { - if (points_.size() > 0) { - std::size_t dimension = points_[0].size(); - if (dimension == 3 && weights_.size() == 0 && !default_filtration_value) { - if (fast_version_) - alpha_ptr_ = std::make_unique>(points_); - else if (exact_version_) - alpha_ptr_ = std::make_unique>(points_); - else - alpha_ptr_ = std::make_unique>(points_); - if (!alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value)) { - // create_simplex_tree will fail if all points are on a plane - Retry with dD by setting dimension to 2 - dimension--; - alpha_ptr_.reset(); - } - } - // Not ** else ** because we have to take into account if 3d fails - if (dimension != 3 || weights_.size() != 0 || default_filtration_value) { - if (fast_version_) { - if (weights_.size() == 0) { - alpha_ptr_ = std::make_unique(points_, exact_version_); - } else { - alpha_ptr_ = std::make_unique(points_, weights_, exact_version_); - } - } else { - if (weights_.size() == 0) { - alpha_ptr_ = std::make_unique(points_, exact_version_); - } else { - alpha_ptr_ = std::make_unique(points_, weights_, exact_version_); - } - } - alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value); - } - } + // Nothing to be done in case of an empty point set + if (!empty_point_set_) + alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value); } private: std::unique_ptr alpha_ptr_; - std::vector> points_; - std::vector weights_; - bool fast_version_; - bool exact_version_; + bool empty_point_set_; }; } // namespace alpha_complex diff --git a/src/python/include/Alpha_complex_interface_3d.h b/src/python/include/Alpha_complex_interface_3d.h new file mode 100644 index 00000000..bb66b8e1 --- /dev/null +++ b/src/python/include/Alpha_complex_interface_3d.h @@ -0,0 +1,71 @@ +/* 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): Vincent Rouvreau + * + * Copyright (C) 2021 Inria + * + * Modification(s): + * - YYYY/MM Author: Description of the modification + */ + +#ifndef INCLUDE_ALPHA_COMPLEX_INTERFACE_3D_H_ +#define INCLUDE_ALPHA_COMPLEX_INTERFACE_3D_H_ + +#include "Alpha_complex_factory.h" +#include + +#include "Simplex_tree_interface.h" + +#include +#include +#include +#include // for std::unique_ptr + +namespace Gudhi { + +namespace alpha_complex { + +class Alpha_complex_interface_3d { + public: + Alpha_complex_interface_3d(const std::vector>& points, + const std::vector& weights, + bool fast_version, bool exact_version) + : empty_point_set_(points.size() == 0) { + const bool weighted = (weights.size() > 0); + if (fast_version) + if (weighted) + alpha_ptr_ = std::make_unique>(points, weights); + else + alpha_ptr_ = std::make_unique>(points); + else if (exact_version) + if (weighted) + alpha_ptr_ = std::make_unique>(points, weights); + else + alpha_ptr_ = std::make_unique>(points); + else + if (weighted) + alpha_ptr_ = std::make_unique>(points, weights); + else + alpha_ptr_ = std::make_unique>(points); + } + + std::vector get_point(int vh) { + return alpha_ptr_->get_point(vh); + } + + void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square) { + // Nothing to be done in case of an empty point set + if (!empty_point_set_) + alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, false); + } + + private: + std::unique_ptr alpha_ptr_; + bool empty_point_set_; +}; + +} // namespace alpha_complex + +} // namespace Gudhi + +#endif // INCLUDE_ALPHA_COMPLEX_INTERFACE_3D_H_ -- cgit v1.2.3 From d79ef2f79c33d8e8b8fd88e16f30e16433ab6644 Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Wed, 7 Apr 2021 16:47:34 +0200 Subject: alpha complex python module requires cgal 5.1 for weighted version. Alpha complex factory does not require a weighted specific version for alpha complex dD --- src/python/CMakeLists.txt | 12 ++-- src/python/include/Alpha_complex_factory.h | 102 ++++++++------------------- src/python/include/Alpha_complex_interface.h | 13 ++-- 3 files changed, 46 insertions(+), 81 deletions(-) (limited to 'src/python/include') diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 307181b7..0739d7b5 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -156,13 +156,15 @@ if(PYTHONINTERP_FOUND) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'nerve_gic', ") endif () if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) - set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex_3d', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'subsampling', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'tangential_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'euclidean_witness_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'euclidean_strong_witness_complex', ") endif () + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) + set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex', ") + endif () if(CGAL_FOUND) # Add CGAL compilation args @@ -345,13 +347,15 @@ if(PYTHONINTERP_FOUND) # Test examples - if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) # Bottleneck and Alpha add_test(NAME alpha_rips_persistence_bottleneck_distance_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}" ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/alpha_rips_persistence_bottleneck_distance.py" -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -t 0.15 -d 3) + endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) # Tangential add_test(NAME tangential_complex_plain_homology_from_off_file_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} @@ -421,7 +425,7 @@ if(PYTHONINTERP_FOUND) add_gudhi_py_test(test_cover_complex) endif (NOT CGAL_VERSION VERSION_LESS 4.11.0) - if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) # Alpha add_test(NAME alpha_complex_from_points_example_py_test WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR} @@ -433,7 +437,7 @@ if(PYTHONINTERP_FOUND) ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/alpha_complex_diagram_persistence_from_off_file_example.py" --no-diagram -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off) add_gudhi_py_test(test_alpha_complex) - endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) + endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0) if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) # Euclidean witness diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h index 5d3bfb65..7d45af7c 100644 --- a/src/python/include/Alpha_complex_factory.h +++ b/src/python/include/Alpha_complex_factory.h @@ -31,9 +31,11 @@ namespace Gudhi { namespace alpha_complex { +// template Functor that transforms a CGAL point to a vector of double as expected by cython template struct Point_cgal_to_cython; +// Specialized Unweighted Functor template struct Point_cgal_to_cython { std::vector operator()(CgalPointType const& point) const @@ -46,6 +48,7 @@ struct Point_cgal_to_cython { } }; +// Specialized Weighted Functor template struct Point_cgal_to_cython { std::vector operator()(CgalPointType const& weighted_point) const @@ -59,15 +62,7 @@ struct Point_cgal_to_cython { } }; -template -std::vector pt_cgal_to_cython(CgalPointType const& point) { - std::vector vd; - vd.reserve(point.dimension()); - for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++) - vd.push_back(CGAL::to_double(*coord)); - return vd; -} - +// Function that transforms a cython point (aka. a vector of double) to a CGAL point template static CgalPointType pt_cython_to_cgal(std::vector const& vec) { return CgalPointType(vec.size(), vec.begin(), vec.end()); @@ -83,20 +78,29 @@ class Abstract_alpha_complex { virtual ~Abstract_alpha_complex() = default; }; +template class Exact_alpha_complex_dD final : public Abstract_alpha_complex { private: using Kernel = CGAL::Epeck_d; - using Point = typename Kernel::Point_d; + using Bare_point = typename Kernel::Point_d; + using Point = std::conditional_t; public: Exact_alpha_complex_dD(const std::vector>& points, bool exact_version) : exact_version_(exact_version), - alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal)) { + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal)) { + } + + Exact_alpha_complex_dD(const std::vector>& points, + const std::vector& weights, bool exact_version) + : exact_version_(exact_version), + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal), weights) { } virtual std::vector get_point(int vh) override { - Point const& point = alpha_complex_.get_point(vh); - return pt_cgal_to_cython(point); + // Can be a Weighted or a Bare point in function of Weighted + return Point_cgal_to_cython()(alpha_complex_.get_point(vh)); } virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, @@ -106,76 +110,32 @@ class Exact_alpha_complex_dD final : public Abstract_alpha_complex { private: bool exact_version_; - Alpha_complex alpha_complex_; + Alpha_complex alpha_complex_; }; +template class Inexact_alpha_complex_dD final : public Abstract_alpha_complex { private: using Kernel = CGAL::Epick_d; - using Point = typename Kernel::Point_d; + using Bare_point = typename Kernel::Point_d; + using Point = std::conditional_t; public: Inexact_alpha_complex_dD(const std::vector>& points, bool exact_version) : exact_version_(exact_version), - alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal)) { - } - - virtual std::vector get_point(int vh) override { - Point const& point = alpha_complex_.get_point(vh); - return pt_cgal_to_cython(point); - } - virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, - bool default_filtration_value) override { - return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value); - } - - private: - bool exact_version_; - Alpha_complex alpha_complex_; -}; - -class Exact_weighted_alpha_complex_dD final : public Abstract_alpha_complex { - private: - using Kernel = CGAL::Epeck_d; - using Point = typename Kernel::Point_d; - - public: - Exact_weighted_alpha_complex_dD(const std::vector>& points, - const std::vector& weights, bool exact_version) - : exact_version_(exact_version), - alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal), weights) { - } - - virtual std::vector get_point(int vh) override { - Point const& point = alpha_complex_.get_point(vh).point(); - return pt_cgal_to_cython(point); - } - - virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, - bool default_filtration_value) override { - return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value); + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal)) { } - private: - bool exact_version_; - Alpha_complex alpha_complex_; -}; - -class Inexact_weighted_alpha_complex_dD final : public Abstract_alpha_complex { - private: - using Kernel = CGAL::Epick_d; - using Point = typename Kernel::Point_d; - - public: - Inexact_weighted_alpha_complex_dD(const std::vector>& points, - const std::vector& weights, bool exact_version) + Inexact_alpha_complex_dD(const std::vector>& points, + const std::vector& weights, bool exact_version) : exact_version_(exact_version), - alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal), weights) { + alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal), weights) { } virtual std::vector get_point(int vh) override { - Point const& point = alpha_complex_.get_point(vh).point(); - return pt_cgal_to_cython(point); + // Can be a Weighted or a Bare point in function of Weighted + return Point_cgal_to_cython()(alpha_complex_.get_point(vh)); } virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) override { @@ -184,7 +144,7 @@ class Inexact_weighted_alpha_complex_dD final : public Abstract_alpha_complex { private: bool exact_version_; - Alpha_complex alpha_complex_; + Alpha_complex alpha_complex_; }; template @@ -207,8 +167,8 @@ class Alpha_complex_3D final : public Abstract_alpha_complex { } virtual std::vector get_point(int vh) override { - Point const& point = alpha_complex_.get_point(vh); - return Point_cgal_to_cython()(point); + // Can be a Weighted or a Bare point in function of Weighted + return Point_cgal_to_cython()(alpha_complex_.get_point(vh)); } virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index 31a8147b..ed243f19 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -31,17 +31,18 @@ class Alpha_complex_interface { const std::vector& weights, bool fast_version, bool exact_version) : empty_point_set_(points.size() == 0) { + const bool weighted = (weights.size() > 0); if (fast_version) { - if (weights.size() == 0) { - alpha_ptr_ = std::make_unique(points, exact_version); + if (weighted) { + alpha_ptr_ = std::make_unique>(points, weights, exact_version); } else { - alpha_ptr_ = std::make_unique(points, weights, exact_version); + alpha_ptr_ = std::make_unique>(points, exact_version); } } else { - if (weights.size() == 0) { - alpha_ptr_ = std::make_unique(points, exact_version); + if (weighted) { + alpha_ptr_ = std::make_unique>(points, weights, exact_version); } else { - alpha_ptr_ = std::make_unique(points, weights, exact_version); + alpha_ptr_ = std::make_unique>(points, exact_version); } } } -- cgit v1.2.3 From 2f507728035c76f16c732b911c52a9118a9b52dd Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Wed, 9 Jun 2021 09:57:26 +0200 Subject: code review: constify auto --- src/python/include/Alpha_complex_factory.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/python/include') diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h index 7d45af7c..fbbf8896 100644 --- a/src/python/include/Alpha_complex_factory.h +++ b/src/python/include/Alpha_complex_factory.h @@ -53,7 +53,7 @@ template struct Point_cgal_to_cython { std::vector operator()(CgalPointType const& weighted_point) const { - auto point = weighted_point.point(); + const auto& point = weighted_point.point(); std::vector vd; vd.reserve(point.dimension()); for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++) -- cgit v1.2.3 From f70e6f26f329428184fc5cf935ad4dfc20648bfb Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Thu, 1 Jul 2021 17:30:38 +0200 Subject: Revert AlphaComplex3D. To be done with periodic --- .../example/Alpha_complex_3d_from_points.cpp | 29 ++-- .../Weighted_alpha_complex_3d_from_points.cpp | 29 ++-- src/Alpha_complex/include/gudhi/Alpha_complex_3d.h | 19 +-- .../test/Alpha_complex_3d_unit_test.cpp | 30 ---- .../test/Weighted_alpha_complex_unit_test.cpp | 2 +- src/python/CMakeLists.txt | 2 - src/python/doc/alpha_complex_ref.rst | 6 - src/python/doc/alpha_complex_user.rst | 14 -- src/python/gudhi/alpha_complex_3d.pyx | 135 ----------------- src/python/include/Alpha_complex_factory.h | 35 ----- src/python/test/test_alpha_complex_3d.py | 159 --------------------- 11 files changed, 42 insertions(+), 418 deletions(-) delete mode 100644 src/python/gudhi/alpha_complex_3d.pyx delete mode 100755 src/python/test/test_alpha_complex_3d.py (limited to 'src/python/include') diff --git a/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp b/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp index dd3c0225..a2c85138 100644 --- a/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp +++ b/src/Alpha_complex/example/Alpha_complex_3d_from_points.cpp @@ -34,22 +34,23 @@ int main(int argc, char **argv) { Alpha_complex_3d alpha_complex_from_points(points); Gudhi::Simplex_tree<> simplex; - alpha_complex_from_points.create_complex(simplex); - // ---------------------------------------------------------------------------- - // Display information about the alpha complex - // ---------------------------------------------------------------------------- - std::clog << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() - << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; + if (alpha_complex_from_points.create_complex(simplex)) { + // ---------------------------------------------------------------------------- + // Display information about the alpha complex + // ---------------------------------------------------------------------------- + std::clog << "Alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() + << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; - std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; - for (auto f_simplex : simplex.filtration_simplex_range()) { - std::clog << " ( "; - for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { - std::clog << vertex << " "; + std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + for (auto f_simplex : simplex.filtration_simplex_range()) { + std::clog << " ( "; + for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { + std::clog << vertex << " "; + } + std::clog << ") -> " + << "[" << simplex.filtration(f_simplex) << "] "; + std::clog << std::endl; } - std::clog << ") -> " - << "[" << simplex.filtration(f_simplex) << "] "; - std::clog << std::endl; } return 0; } diff --git a/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp b/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp index 507d6413..ee12d418 100644 --- a/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp +++ b/src/Alpha_complex/example/Weighted_alpha_complex_3d_from_points.cpp @@ -30,22 +30,23 @@ int main(int argc, char **argv) { Weighted_alpha_complex_3d alpha_complex_from_points(weighted_points); Gudhi::Simplex_tree<> simplex; - alpha_complex_from_points.create_complex(simplex); - // ---------------------------------------------------------------------------- - // Display information about the alpha complex - // ---------------------------------------------------------------------------- - std::clog << "Weighted alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() - << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; + if (alpha_complex_from_points.create_complex(simplex)) { + // ---------------------------------------------------------------------------- + // Display information about the alpha complex + // ---------------------------------------------------------------------------- + std::clog << "Weighted alpha complex is of dimension " << simplex.dimension() << " - " << simplex.num_simplices() + << " simplices - " << simplex.num_vertices() << " vertices." << std::endl; - std::clog << "Iterator on weighted alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; - for (auto f_simplex : simplex.filtration_simplex_range()) { - std::clog << " ( "; - for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { - std::clog << vertex << " "; + std::clog << "Iterator on weighted alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; + for (auto f_simplex : simplex.filtration_simplex_range()) { + std::clog << " ( "; + for (auto vertex : simplex.simplex_vertex_range(f_simplex)) { + std::clog << vertex << " "; + } + std::clog << ") -> " + << "[" << simplex.filtration(f_simplex) << "] "; + std::clog << std::endl; } - std::clog << ") -> " - << "[" << simplex.filtration(f_simplex) << "] "; - std::clog << std::endl; } return 0; } diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h index 73f9dd41..4e5fc933 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex_3d.h @@ -48,7 +48,6 @@ #include // for std::unique_ptr #include // for std::conditional and std::enable_if #include // for numeric_limits<> -#include // for domain_error and invalid_argument // Make compilation fail - required for external projects - https://github.com/GUDHI/gudhi-devel/issues/10 #if CGAL_VERSION_NR < 1041101000 @@ -429,18 +428,19 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ * @param[in] max_alpha_square maximum for alpha square value. Default value is +\f$\infty\f$, and there is very * little point using anything else since it does not save time. * - * @exception invalid_argument In debug mode, if `complex` given as argument is not empty. - * @exception domain_error If `points` given in the constructor are on a 2d plane. + * @return true if creation succeeds, false otherwise. * * @pre The simplicial complex must be empty (no vertices). * */ template - void create_complex(SimplicialComplexForAlpha3d& complex, + bool create_complex(SimplicialComplexForAlpha3d& complex, Filtration_value max_alpha_square = std::numeric_limits::infinity()) { - GUDHI_CHECK(complex.num_vertices() == 0, - std::invalid_argument("Alpha_complex_3d create_complex - The complex given as argument is not empty")); + if (complex.num_vertices() > 0) { + std::cerr << "Alpha_complex_3d create_complex - complex is not empty\n"; + return false; // ----- >> + } using Complex_vertex_handle = typename SimplicialComplexForAlpha3d::Vertex_handle; using Simplex_tree_vector_vertex = std::vector; @@ -461,8 +461,10 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ #ifdef DEBUG_TRACES std::clog << "filtration_with_alpha_values returns : " << objects.size() << " objects" << std::endl; #endif // DEBUG_TRACES - if (objects.size() == 0) - throw std::domain_error("Alpha_complex_3d create_complex - no triangulation as points are on a 2d plane"); + if (objects.size() == 0) { + std::cerr << "Alpha_complex_3d create_complex - no triangulation as points are on a 2d plane\n"; + return false; // ----- >> + } using Alpha_value_iterator = typename std::vector::const_iterator; Alpha_value_iterator alpha_value_iterator = alpha_values.begin(); @@ -557,6 +559,7 @@ Weighted_alpha_complex_3d::Weighted_point_3 wp0(Weighted_alpha_complex_3d::Bare_ // Remove all simplices that have a filtration value greater than max_alpha_square complex.prune_above_filtration(max_alpha_square); // -------------------------------------------------------------------------------------------- + return true; } /** \brief get_point returns the point corresponding to the vertex given as parameter. diff --git a/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp index 94021954..a4ecb6ad 100644 --- a/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp +++ b/src/Alpha_complex/test/Alpha_complex_3d_unit_test.cpp @@ -11,7 +11,6 @@ #define BOOST_TEST_DYN_LINK #define BOOST_TEST_MODULE "alpha_complex_3d" #include -#include #include // float comparison #include @@ -37,7 +36,6 @@ using Safe_alpha_complex_3d = using Exact_alpha_complex_3d = Gudhi::alpha_complex::Alpha_complex_3d; - template std::vector get_points() { std::vector points; @@ -199,31 +197,3 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_3d_from_points) { ++safe_sh; } } - -typedef boost::mpl::list list_of_alpha_variants; - -BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_3d_exceptions_points_on_plane, Alpha, list_of_alpha_variants) { - std::vector points_on_plane; - points_on_plane.emplace_back(1.0, 1.0 , 0.0); - points_on_plane.emplace_back(7.0, 0.0 , 0.0); - points_on_plane.emplace_back(4.0, 6.0 , 0.0); - points_on_plane.emplace_back(9.0, 6.0 , 0.0); - points_on_plane.emplace_back(0.0, 14.0, 0.0); - points_on_plane.emplace_back(2.0, 19.0, 0.0); - points_on_plane.emplace_back(9.0, 17.0, 0.0); - - Alpha alpha_complex(points_on_plane); - Gudhi::Simplex_tree<> stree; - - BOOST_CHECK_THROW(alpha_complex.create_complex(stree), std::domain_error); -} - -BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_3d_exceptions_non_empty_simplex_tree, Alpha, list_of_alpha_variants) { - Alpha alpha_complex(get_points()); - Gudhi::Simplex_tree<> stree; - stree.insert_simplex_and_subfaces({2,1,0}, 3.0); - -#ifdef GUDHI_DEBUG - BOOST_CHECK_THROW(alpha_complex.create_complex(stree), std::invalid_argument); -#endif // GUDHI_DEBUG -} \ No newline at end of file diff --git a/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp index 4e1a38df..875704ee 100644 --- a/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp +++ b/src/Alpha_complex/test/Weighted_alpha_complex_unit_test.cpp @@ -83,7 +83,7 @@ BOOST_AUTO_TEST_CASE(Weighted_alpha_complex_3d_comparison) { // Weighted alpha complex for 3D version Exact_weighted_alpha_complex_3d alpha_complex_3D_from_weighted_points(w_points_3); Gudhi::Simplex_tree<> w_simplex_3; - alpha_complex_3D_from_weighted_points.create_complex(w_simplex_3); + BOOST_CHECK(alpha_complex_3D_from_weighted_points.create_complex(w_simplex_3)); std::clog << "Iterator on weighted alpha complex 3D simplices in the filtration order, with [filtration value]:" << std::endl; diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt index 669239b8..bfa78131 100644 --- a/src/python/CMakeLists.txt +++ b/src/python/CMakeLists.txt @@ -62,7 +62,6 @@ if(PYTHONINTERP_FOUND) set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'subsampling', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'tangential_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'alpha_complex', ") - set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'alpha_complex_3d', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'euclidean_witness_complex', ") set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'euclidean_strong_witness_complex', ") # Modules that should not be auto-imported in __init__.py @@ -158,7 +157,6 @@ if(PYTHONINTERP_FOUND) set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'nerve_gic', ") endif () if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0) - set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'alpha_complex_3d', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'subsampling', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'tangential_complex', ") set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'euclidean_witness_complex', ") diff --git a/src/python/doc/alpha_complex_ref.rst b/src/python/doc/alpha_complex_ref.rst index 49321368..eaa72551 100644 --- a/src/python/doc/alpha_complex_ref.rst +++ b/src/python/doc/alpha_complex_ref.rst @@ -11,9 +11,3 @@ Alpha complex reference manual :undoc-members: .. automethod:: gudhi.AlphaComplex.__init__ - -.. autoclass:: gudhi.AlphaComplex3D - :members: - :undoc-members: - - .. automethod:: gudhi.AlphaComplex3D.__init__ diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst index d7b09246..db0ccdc9 100644 --- a/src/python/doc/alpha_complex_user.rst +++ b/src/python/doc/alpha_complex_user.rst @@ -254,17 +254,3 @@ Then, it computes the persistence diagram and displays it: dgm = stree.persistence() gd.plot_persistence_diagram(dgm, legend = True) plt.show() - -3d specific version -^^^^^^^^^^^^^^^^^^^ - -:Requires: `Eigen `_ :math:`\geq` 3.1.0 and `CGAL `_ :math:`\geq` 4.11.0. - -A specific module for Alpha complex is available in 3d (cf. :class:`~gudhi.AlphaComplex3D`) and -allows to construct standard and weighted versions of alpha complexes. - -Remark -"""""" - -* Contrary to the dD version, with the 3d version, the vertices in the output simplex tree are not guaranteed to match - the order of the input points. One can use :func:`~gudhi.AlphaComplex3D.get_point` to get the initial point back. diff --git a/src/python/gudhi/alpha_complex_3d.pyx b/src/python/gudhi/alpha_complex_3d.pyx deleted file mode 100644 index 578011a7..00000000 --- a/src/python/gudhi/alpha_complex_3d.pyx +++ /dev/null @@ -1,135 +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): Vincent Rouvreau -# -# Copyright (C) 2021 Inria -# -# Modification(s): -# - YYYY/MM Author: Description of the modification - -from __future__ import print_function -from cython cimport numeric -from libcpp.vector cimport vector -from libcpp.utility cimport pair -from libcpp.string cimport string -from libcpp cimport bool -from libc.stdint cimport intptr_t -import errno -import os -import warnings - -from gudhi.simplex_tree cimport * -from gudhi.simplex_tree import SimplexTree -from gudhi import read_points_from_off_file - -__author__ = "Vincent Rouvreau" -__copyright__ = "Copyright (C) 2021 Inria" -__license__ = "GPL v3" - -cdef extern from "Alpha_complex_interface_3d.h" namespace "Gudhi": - cdef cppclass Alpha_complex_interface_3d "Gudhi::alpha_complex::Alpha_complex_interface_3d": - Alpha_complex_interface_3d(vector[vector[double]] points, vector[double] weights, bool fast_version, bool exact_version) nogil except + - vector[double] get_point(int vertex) nogil except + - void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square) nogil except + - -# AlphaComplex3D python interface -cdef class AlphaComplex3D: - """AlphaComplex3D is a simplicial complex constructed from the finite cells of a Delaunay Triangulation. - - The filtration value of each simplex is computed as the square of the circumradius of the simplex if the - circumsphere is empty (the simplex is then said to be Gabriel), and as the minimum of the filtration values of the - codimension 1 cofaces that make it not Gabriel otherwise. - - All simplices that have a filtration value strictly greater than a given alpha squared value are not inserted into - the complex. - - .. note:: - - When AlphaComplex3D is constructed with an infinite value of alpha, the complex is a Delaunay complex. - - .. warning:: - - Contrary to the dD version, with the 3d version, the vertices in the output simplex tree are not guaranteed to - match the order of the input points. One can use :func:`~gudhi.AlphaComplex3D.get_point` to get the initial - point back. - """ - - cdef Alpha_complex_interface_3d * this_ptr - - # Fake constructor that does nothing but documenting the constructor - def __init__(self, points=[], weights=[], precision='safe'): - """AlphaComplex3D constructor. - - :param points: A list of points in 3d. - :type points: Iterable[Iterable[float]] - - :param weights: A list of weights. If set, the number of weights must correspond to the number of points. - :type weights: Iterable[float] - - :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'. - :type precision: string - - :raises ValueError: If the given points are not in 3d. - :raises ValueError: In case of inconsistency between the number of points and weights. - """ - - # The real cython constructor - def __cinit__(self, points = [], weights=[], precision = 'safe'): - assert precision in ['fast', 'safe', 'exact'], "Alpha complex precision can only be 'fast', 'safe' or 'exact'" - cdef bool fast = precision == 'fast' - cdef bool exact = precision == 'exact' - - if len(points) > 0: - if len(points[0]) != 3: - raise ValueError("AlphaComplex3D only accepts 3d points as an input") - - # weights are set but is inconsistent with the number of points - if len(weights) != 0 and len(weights) != len(points): - raise ValueError("Inconsistency between the number of points and weights") - - # need to copy the points to use them without the gil - cdef vector[vector[double]] pts - cdef vector[double] wgts - pts = points - wgts = weights - with nogil: - self.this_ptr = new Alpha_complex_interface_3d(pts, wgts, fast, exact) - - def __dealloc__(self): - if self.this_ptr != NULL: - del self.this_ptr - - def __is_defined(self): - """Returns true if AlphaComplex3D pointer is not NULL. - """ - return self.this_ptr != NULL - - def get_point(self, vertex): - """This function returns the point corresponding to a given vertex from the :class:`~gudhi.SimplexTree`. - - :param vertex: The vertex. - :type vertex: int - :rtype: list of float - :returns: the point. - """ - return self.this_ptr.get_point(vertex) - - def create_simplex_tree(self, max_alpha_square = float('inf')): - """ - :param max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to - infinity, and there is very little point using anything else since it does not save time. - :type max_alpha_square: float - :returns: A simplex tree created from the Delaunay Triangulation. - :rtype: SimplexTree - - :raises ValueError: If the points given at construction time are on a plane. - """ - stree = SimplexTree() - cdef double mas = max_alpha_square - cdef intptr_t stree_int_ptr=stree.thisptr - with nogil: - self.this_ptr.create_simplex_tree(stree_int_ptr, - mas) - return stree diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h index fbbf8896..298469fe 100644 --- a/src/python/include/Alpha_complex_factory.h +++ b/src/python/include/Alpha_complex_factory.h @@ -147,41 +147,6 @@ class Inexact_alpha_complex_dD final : public Abstract_alpha_complex { Alpha_complex alpha_complex_; }; -template -class Alpha_complex_3D final : public Abstract_alpha_complex { - private: - using Bare_point = typename Alpha_complex_3d::Bare_point_3; - using Point = typename Alpha_complex_3d::Point_3; - - static Bare_point pt_cython_to_cgal_3(std::vector const& vec) { - return Bare_point(vec[0], vec[1], vec[2]); - } - - public: - Alpha_complex_3D(const std::vector>& points) - : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3)) { - } - - Alpha_complex_3D(const std::vector>& points, const std::vector& weights) - : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal_3), weights) { - } - - virtual std::vector get_point(int vh) override { - // Can be a Weighted or a Bare point in function of Weighted - return Point_cgal_to_cython()(alpha_complex_.get_point(vh)); - } - - virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, - bool default_filtration_value) override { - alpha_complex_.create_complex(*simplex_tree, max_alpha_square); - return true; - } - - private: - Alpha_complex_3d alpha_complex_; -}; - - } // namespace alpha_complex } // namespace Gudhi diff --git a/src/python/test/test_alpha_complex_3d.py b/src/python/test/test_alpha_complex_3d.py deleted file mode 100755 index a5d9373b..00000000 --- a/src/python/test/test_alpha_complex_3d.py +++ /dev/null @@ -1,159 +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): Vincent Rouvreau - - Copyright (C) 2021 Inria - - Modification(s): - - YYYY/MM Author: Description of the modification -""" - -from gudhi import AlphaComplex3D -import pytest -import numpy as np - -try: - # python3 - from itertools import zip_longest -except ImportError: - # python2 - from itertools import izip_longest as zip_longest - - - -def _empty_alpha(precision): - alpha_complex = AlphaComplex3D(precision = precision) - assert alpha_complex.__is_defined() == True - -def _one_3d_point_alpha(precision): - alpha_complex = AlphaComplex3D(points=[[0, 0, 0]], precision = precision) - assert alpha_complex.__is_defined() == True - -def test_empty_alpha(): - for precision in ['fast', 'safe', 'exact']: - _empty_alpha(precision) - _one_3d_point_alpha(precision) - -def _infinite_alpha(precision): - point_list = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0], [0, 0, 1], [1, 0, 1], [0, 1, 1], [1, 1, 1]] - alpha_complex = AlphaComplex3D(points=point_list, precision = precision) - assert alpha_complex.__is_defined() == True - - stree = alpha_complex.create_simplex_tree() - assert stree.__is_persistence_defined() == False - - assert stree.num_simplices() == 51 - assert stree.num_vertices() == len(point_list) - - for filtration in stree.get_filtration(): - if len(filtration[0]) == 1: - assert filtration[1] == 0. - if len(filtration[0]) == 4: - assert filtration[1] == 0.75 - - for idx in range(len(point_list)): - pt_idx = point_list.index(alpha_complex.get_point(idx)) - assert pt_idx >= 0 - assert pt_idx < len(point_list) - - with pytest.raises(IndexError): - alpha_complex.get_point(len(point_list)) - -def test_infinite_alpha(): - for precision in ['fast', 'safe', 'exact']: - _infinite_alpha(precision) - -def _filtered_alpha(precision): - point_list = [[0, 0, 0], [1, 0, 0], [0, 1, 0], [1, 1, 0], [0, 0, 1], [1, 0, 1], [0, 1, 1], [1, 1, 1]] - filtered_alpha = AlphaComplex3D(points=point_list, precision = precision) - - stree = filtered_alpha.create_simplex_tree(max_alpha_square=0.25) - - assert stree.num_simplices() == 20 - assert stree.num_vertices() == len(point_list) - - for filtration in stree.get_filtration(): - if len(filtration[0]) == 1: - assert filtration[1] == 0. - elif len(filtration[0]) == 2: - assert filtration[1] == 0.25 - else: - assert False - - for idx in range(len(point_list)): - pt_idx = point_list.index(filtered_alpha.get_point(idx)) - assert pt_idx >= 0 - assert pt_idx < len(point_list) - - with pytest.raises(IndexError): - filtered_alpha.get_point(len(point_list)) - -def test_filtered_alpha(): - for precision in ['fast', 'safe', 'exact']: - _filtered_alpha(precision) - -def _3d_points_on_a_plane(precision): - alpha = AlphaComplex3D(points = [[1.0, 1.0 , 0.0], - [7.0, 0.0 , 0.0], - [4.0, 6.0 , 0.0], - [9.0, 6.0 , 0.0], - [0.0, 14.0, 0.0], - [2.0, 19.0, 0.0], - [9.0, 17.0, 0.0]], precision = precision) - - with pytest.raises(ValueError): - stree = alpha.create_simplex_tree() - -def test_3d_points_on_a_plane(): - for precision in ['fast', 'safe', 'exact']: - _3d_points_on_a_plane(precision) - -def test_inconsistency_points_and_weights(): - points = [[1.0, 1.0 , 1.0], - [7.0, 0.0 , 2.0], - [4.0, 6.0 , 0.0], - [9.0, 6.0 , 1.0], - [0.0, 14.0, 2.0], - [2.0, 19.0, 0.0], - [9.0, 17.0, 1.0]] - with pytest.raises(ValueError): - # 7 points, 8 weights, on purpose - alpha = AlphaComplex3D(points = points, - weights = [1., 2., 3., 4., 5., 6., 7., 8.]) - - with pytest.raises(ValueError): - # 7 points, 6 weights, on purpose - alpha = AlphaComplex3D(points = points, - weights = [1., 2., 3., 4., 5., 6.]) - -def _weighted_doc_example(precision): - pts = [[ 1., -1., -1.], - [-1., 1., -1.], - [-1., -1., 1.], - [ 1., 1., 1.], - [ 2., 2., 2.]] - wgts = [4., 4., 4., 4., 1.] - alpha = AlphaComplex3D(points = pts, weights = wgts, precision = precision) - stree = alpha.create_simplex_tree() - - # Needs to retrieve points as points are shuffled - get_idx = lambda idx: pts.index(alpha.get_point(idx)) - indices = [get_idx(x) for x in range(len(pts))] - - assert stree.filtration([indices[x] for x in [0, 1, 2, 3]]) == pytest.approx(-1.) - assert stree.filtration([indices[x] for x in [0, 1, 3, 4]]) == pytest.approx(95.) - assert stree.filtration([indices[x] for x in [0, 2, 3, 4]]) == pytest.approx(95.) - assert stree.filtration([indices[x] for x in [1, 2, 3, 4]]) == pytest.approx(95.) - -def test_weighted_doc_example(): - for precision in ['fast', 'safe', 'exact']: - _weighted_doc_example(precision) - -def test_points_not_in_3d(): - with pytest.raises(ValueError): - alpha = AlphaComplex3D(points = np.random.rand(6,2)) - with pytest.raises(ValueError): - alpha = AlphaComplex3D(points = np.random.rand(6,4)) - - alpha = AlphaComplex3D(points = np.random.rand(6,3)) - assert alpha.__is_defined() == True \ No newline at end of file -- cgit v1.2.3 From 1617cac6fe6333b2e63f9931399a27643dad45ac Mon Sep 17 00:00:00 2001 From: VincentRouvreau Date: Mon, 6 Sep 2021 15:51:20 +0200 Subject: Only dD version of Alpha complex for python --- src/python/include/Alpha_complex_interface_3d.h | 71 ------------------------- 1 file changed, 71 deletions(-) delete mode 100644 src/python/include/Alpha_complex_interface_3d.h (limited to 'src/python/include') diff --git a/src/python/include/Alpha_complex_interface_3d.h b/src/python/include/Alpha_complex_interface_3d.h deleted file mode 100644 index bb66b8e1..00000000 --- a/src/python/include/Alpha_complex_interface_3d.h +++ /dev/null @@ -1,71 +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): Vincent Rouvreau - * - * Copyright (C) 2021 Inria - * - * Modification(s): - * - YYYY/MM Author: Description of the modification - */ - -#ifndef INCLUDE_ALPHA_COMPLEX_INTERFACE_3D_H_ -#define INCLUDE_ALPHA_COMPLEX_INTERFACE_3D_H_ - -#include "Alpha_complex_factory.h" -#include - -#include "Simplex_tree_interface.h" - -#include -#include -#include -#include // for std::unique_ptr - -namespace Gudhi { - -namespace alpha_complex { - -class Alpha_complex_interface_3d { - public: - Alpha_complex_interface_3d(const std::vector>& points, - const std::vector& weights, - bool fast_version, bool exact_version) - : empty_point_set_(points.size() == 0) { - const bool weighted = (weights.size() > 0); - if (fast_version) - if (weighted) - alpha_ptr_ = std::make_unique>(points, weights); - else - alpha_ptr_ = std::make_unique>(points); - else if (exact_version) - if (weighted) - alpha_ptr_ = std::make_unique>(points, weights); - else - alpha_ptr_ = std::make_unique>(points); - else - if (weighted) - alpha_ptr_ = std::make_unique>(points, weights); - else - alpha_ptr_ = std::make_unique>(points); - } - - std::vector get_point(int vh) { - return alpha_ptr_->get_point(vh); - } - - void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square) { - // Nothing to be done in case of an empty point set - if (!empty_point_set_) - alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, false); - } - - private: - std::unique_ptr alpha_ptr_; - bool empty_point_set_; -}; - -} // namespace alpha_complex - -} // namespace Gudhi - -#endif // INCLUDE_ALPHA_COMPLEX_INTERFACE_3D_H_ -- cgit v1.2.3 From 11b8a288946c953472165ffbb5dbfbeb5ca64d26 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Mon, 8 Nov 2021 10:03:24 +0100 Subject: code review: factorization - Point_cgal_to_cython weighted version calls unweighted version --- src/python/include/Alpha_complex_factory.h | 6 +----- 1 file changed, 1 insertion(+), 5 deletions(-) (limited to 'src/python/include') diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h index 298469fe..8f820fc9 100644 --- a/src/python/include/Alpha_complex_factory.h +++ b/src/python/include/Alpha_complex_factory.h @@ -54,11 +54,7 @@ struct Point_cgal_to_cython { std::vector operator()(CgalPointType const& weighted_point) const { const auto& point = weighted_point.point(); - std::vector vd; - vd.reserve(point.dimension()); - for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++) - vd.push_back(CGAL::to_double(*coord)); - return vd; + return Point_cgal_to_cython()(point); } }; -- cgit v1.2.3 From 72a56067c5a6f105ff60996e7257063505030e73 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Wed, 19 Jan 2022 09:50:13 +0100 Subject: Code review: exact_version was useless in Inexact_alpha_complex_dD --- src/python/include/Alpha_complex_factory.h | 14 +++++--------- src/python/include/Alpha_complex_interface.h | 4 ++-- 2 files changed, 7 insertions(+), 11 deletions(-) (limited to 'src/python/include') diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h index 8f820fc9..14ec1f86 100644 --- a/src/python/include/Alpha_complex_factory.h +++ b/src/python/include/Alpha_complex_factory.h @@ -118,15 +118,12 @@ class Inexact_alpha_complex_dD final : public Abstract_alpha_complex { typename Kernel::Point_d>; public: - Inexact_alpha_complex_dD(const std::vector>& points, bool exact_version) - : exact_version_(exact_version), - alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal)) { + Inexact_alpha_complex_dD(const std::vector>& points) + : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal)) { } - Inexact_alpha_complex_dD(const std::vector>& points, - const std::vector& weights, bool exact_version) - : exact_version_(exact_version), - alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal), weights) { + Inexact_alpha_complex_dD(const std::vector>& points, const std::vector& weights) + : alpha_complex_(boost::adaptors::transform(points, pt_cython_to_cgal), weights) { } virtual std::vector get_point(int vh) override { @@ -135,11 +132,10 @@ class Inexact_alpha_complex_dD final : public Abstract_alpha_complex { } virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) override { - return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value); + return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, false, default_filtration_value); } private: - bool exact_version_; Alpha_complex alpha_complex_; }; diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index ed243f19..188036ed 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -34,9 +34,9 @@ class Alpha_complex_interface { const bool weighted = (weights.size() > 0); if (fast_version) { if (weighted) { - alpha_ptr_ = std::make_unique>(points, weights, exact_version); + alpha_ptr_ = std::make_unique>(points, weights); } else { - alpha_ptr_ = std::make_unique>(points, exact_version); + alpha_ptr_ = std::make_unique>(points); } } else { if (weighted) { -- cgit v1.2.3 From e88833431fbdd2b58b00fe3d3cf84973700477b3 Mon Sep 17 00:00:00 2001 From: Vincent Rouvreau Date: Wed, 26 Jan 2022 09:49:20 +0100 Subject: Code review: Remove empty_point_set_ that can be deduced from num_vertices (new method in Alpha_complex -> Abstract_alpha_complex -> In/Exact_alpha_complex_dD -> Alpha_complex_interface --- src/Alpha_complex/include/gudhi/Alpha_complex.h | 12 +++++++++++- src/Alpha_complex/test/Alpha_complex_unit_test.cpp | 15 +++++++++++++++ src/python/include/Alpha_complex_factory.h | 10 ++++++++++ src/python/include/Alpha_complex_interface.h | 6 ++---- 4 files changed, 38 insertions(+), 5 deletions(-) (limited to 'src/python/include') diff --git a/src/Alpha_complex/include/gudhi/Alpha_complex.h b/src/Alpha_complex/include/gudhi/Alpha_complex.h index e03bb161..028ec9bb 100644 --- a/src/Alpha_complex/include/gudhi/Alpha_complex.h +++ b/src/Alpha_complex/include/gudhi/Alpha_complex.h @@ -20,6 +20,7 @@ #include #include // isnan, fmax #include // for std::unique_ptr +#include // for std::size_t #include #include // aka. Weighted Delaunay triangulation @@ -213,6 +214,15 @@ class Alpha_complex { Alpha_complex (Alpha_complex&& other) = delete; Alpha_complex& operator= (Alpha_complex&& other) = delete; + /** \brief Returns the number of finite vertices in the triangulation. + */ + std::size_t num_vertices() const { + if (triangulation_ == nullptr) + return 0; + else + return triangulation_->number_of_vertices(); + } + /** \brief get_point returns the point corresponding to the vertex given as parameter. * * @param[in] vertex Vertex handle of the point to retrieve. @@ -373,7 +383,7 @@ class Alpha_complex { // -------------------------------------------------------------------------------------------- // Simplex_tree construction from loop on triangulation finite full cells list - if (triangulation_->number_of_vertices() > 0) { + if (num_vertices() > 0) { for (auto cit = triangulation_->finite_full_cells_begin(); cit != triangulation_->finite_full_cells_end(); ++cit) { diff --git a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp index 4b37e4bd..f74ad217 100644 --- a/src/Alpha_complex/test/Alpha_complex_unit_test.cpp +++ b/src/Alpha_complex/test/Alpha_complex_unit_test.cpp @@ -56,6 +56,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of Gudhi::Simplex_tree<> simplex_tree_60; BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_60, max_alpha_square_value)); + std::clog << "alpha_complex_from_file.num_vertices()=" << alpha_complex_from_file.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_file.num_vertices() == 7); + std::clog << "simplex_tree_60.dimension()=" << simplex_tree_60.dimension() << std::endl; BOOST_CHECK(simplex_tree_60.dimension() == 2); @@ -72,6 +75,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_OFF_file, TestedKernel, list_of Gudhi::Simplex_tree<> simplex_tree_59; BOOST_CHECK(alpha_complex_from_file.create_complex(simplex_tree_59, max_alpha_square_value)); + std::clog << "alpha_complex_from_file.num_vertices()=" << alpha_complex_from_file.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_file.num_vertices() == 7); + std::clog << "simplex_tree_59.dimension()=" << simplex_tree_59.dimension() << std::endl; BOOST_CHECK(simplex_tree_59.dimension() == 2); @@ -120,6 +126,9 @@ BOOST_AUTO_TEST_CASE(Alpha_complex_from_points) { Gudhi::Simplex_tree<> simplex_tree; BOOST_CHECK(alpha_complex_from_points.create_complex(simplex_tree)); + std::clog << "alpha_complex_from_points.num_vertices()=" << alpha_complex_from_points.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_points.num_vertices() == points.size()); + // Another way to check num_simplices std::clog << "Iterator on alpha complex simplices in the filtration order, with [filtration value]:" << std::endl; int num_simplices = 0; @@ -240,6 +249,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_from_empty_points, TestedKernel, lis // ---------------------------------------------------------------------------- Gudhi::alpha_complex::Alpha_complex alpha_complex_from_points(points); + std::clog << "alpha_complex_from_points.num_vertices()=" << alpha_complex_from_points.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_points.num_vertices() == points.size()); + // Test to the limit BOOST_CHECK_THROW (alpha_complex_from_points.get_point(0), std::out_of_range); @@ -291,6 +303,9 @@ BOOST_AUTO_TEST_CASE_TEMPLATE(Alpha_complex_with_duplicated_points, TestedKernel std::clog << "create_complex" << std::endl; BOOST_CHECK(alpha_complex_from_points.create_complex(simplex_tree)); + std::clog << "alpha_complex_from_points.num_vertices()=" << alpha_complex_from_points.num_vertices() << std::endl; + BOOST_CHECK(alpha_complex_from_points.num_vertices() < points.size()); + std::clog << "simplex_tree.num_vertices()=" << simplex_tree.num_vertices() << std::endl; BOOST_CHECK(simplex_tree.num_vertices() < points.size()); diff --git a/src/python/include/Alpha_complex_factory.h b/src/python/include/Alpha_complex_factory.h index 14ec1f86..3d20aa8f 100644 --- a/src/python/include/Alpha_complex_factory.h +++ b/src/python/include/Alpha_complex_factory.h @@ -70,6 +70,8 @@ class Abstract_alpha_complex { virtual bool create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) = 0; + + virtual std::size_t num_vertices() const = 0; virtual ~Abstract_alpha_complex() = default; }; @@ -104,6 +106,10 @@ class Exact_alpha_complex_dD final : public Abstract_alpha_complex { return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, exact_version_, default_filtration_value); } + virtual std::size_t num_vertices() const { + return alpha_complex_.num_vertices(); + } + private: bool exact_version_; Alpha_complex alpha_complex_; @@ -135,6 +141,10 @@ class Inexact_alpha_complex_dD final : public Abstract_alpha_complex { return alpha_complex_.create_complex(*simplex_tree, max_alpha_square, false, default_filtration_value); } + virtual std::size_t num_vertices() const { + return alpha_complex_.num_vertices(); + } + private: Alpha_complex alpha_complex_; }; diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h index 188036ed..671af4a4 100644 --- a/src/python/include/Alpha_complex_interface.h +++ b/src/python/include/Alpha_complex_interface.h @@ -29,8 +29,7 @@ class Alpha_complex_interface { public: Alpha_complex_interface(const std::vector>& points, const std::vector& weights, - bool fast_version, bool exact_version) - : empty_point_set_(points.size() == 0) { + bool fast_version, bool exact_version) { const bool weighted = (weights.size() > 0); if (fast_version) { if (weighted) { @@ -54,13 +53,12 @@ class Alpha_complex_interface { void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool default_filtration_value) { // Nothing to be done in case of an empty point set - if (!empty_point_set_) + if (alpha_ptr_->num_vertices() > 0) alpha_ptr_->create_simplex_tree(simplex_tree, max_alpha_square, default_filtration_value); } private: std::unique_ptr alpha_ptr_; - bool empty_point_set_; }; } // namespace alpha_complex -- cgit v1.2.3