From f5c97bc9fd1c247045d35ddf261f9afe4d024406 Mon Sep 17 00:00:00 2001 From: glisse Date: Fri, 20 Jul 2018 15:03:46 +0000 Subject: First try at interfacing the sparse rips in python. Needs at least documentation and tests. git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/branches/sparserips-python@3697 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 75bef59e90355853ee24807ca7453c4bb0a38f43 --- src/cython/include/Rips_complex_interface.h | 39 +++++++++++++++++++---------- 1 file changed, 26 insertions(+), 13 deletions(-) (limited to 'src/cython/include') diff --git a/src/cython/include/Rips_complex_interface.h b/src/cython/include/Rips_complex_interface.h index 8b6c9c35..f6fc79a1 100644 --- a/src/cython/include/Rips_complex_interface.h +++ b/src/cython/include/Rips_complex_interface.h @@ -25,8 +25,11 @@ #include #include +#include #include +#include + #include "Simplex_tree_interface.h" #include @@ -43,28 +46,38 @@ class Rips_complex_interface { using Distance_matrix = std::vector::Filtration_value>>; public: - Rips_complex_interface(const std::vector>& values, double threshold, bool euclidean) { - if (euclidean) { - // Rips construction where values is a vector of points - rips_complex_ = new Rips_complex::Filtration_value>(values, threshold, - Gudhi::Euclidean_distance()); - } else { - // Rips construction where values is a distance matrix - rips_complex_ = new Rips_complex::Filtration_value>(values, threshold); - } + void init_points(const std::vector>& points, double threshold) { + rips_complex_.emplace(points, threshold, Gudhi::Euclidean_distance()); + } + void init_matrix(const std::vector>& matrix, double threshold) { + rips_complex_.emplace(matrix, threshold); } - ~Rips_complex_interface() { - delete rips_complex_; + void init_points_sparse(const std::vector>& points, double threshold, double epsilon) { + sparse_rips_complex_.emplace(points, Gudhi::Euclidean_distance(), epsilon); + threshold_ = threshold; + } + void init_matrix_sparse(const std::vector>& matrix, double threshold, double epsilon) { + sparse_rips_complex_.emplace(matrix, epsilon); + threshold_ = threshold; } void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, int dim_max) { - rips_complex_->create_complex(*simplex_tree, dim_max); + if (rips_complex_) + rips_complex_->create_complex(*simplex_tree, dim_max); + else { + sparse_rips_complex_->create_complex(*simplex_tree, dim_max); + // This pruning should be done much earlier! It isn't that useful for sparse Rips, but it would be inconsistent not to do it. + simplex_tree->prune_above_filtration(threshold_); + } simplex_tree->initialize_filtration(); } private: - Rips_complex::Filtration_value>* rips_complex_; + // std::variant would work, but we don't require C++17 yet, and boost::variant is not super convenient. Anyway, storing a graph would make more sense. Or changing the interface completely so there is no such storage. + boost::optional::Filtration_value>> rips_complex_; + boost::optional::Filtration_value>> sparse_rips_complex_; + double threshold_; }; } // namespace rips_complex -- cgit v1.2.3 From cae07e49a1700613bf709f7d0bbf3d1c640d4a2c Mon Sep 17 00:00:00 2001 From: glisse Date: Tue, 2 Oct 2018 11:09:47 +0000 Subject: Catch exception by const& git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/trunk@3922 636b058d-ea47-450e-bf9e-a15bfbe3eedb Former-commit-id: 1b0cea99516d6b42e5ab3b0fdb21468b321e6d6c --- src/cython/include/Alpha_complex_interface.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src/cython/include') diff --git a/src/cython/include/Alpha_complex_interface.h b/src/cython/include/Alpha_complex_interface.h index 8cf527fc..faa059d1 100644 --- a/src/cython/include/Alpha_complex_interface.h +++ b/src/cython/include/Alpha_complex_interface.h @@ -60,7 +60,7 @@ class Alpha_complex_interface { Point_d ph = alpha_complex_->get_point(vh); for (auto coord = ph.cartesian_begin(); coord < ph.cartesian_end(); coord++) vd.push_back(*coord); - } catch (std::out_of_range outofrange) { + } catch (std::out_of_range const&) { // std::out_of_range is thrown in case not found. Other exceptions must be re-thrown } return vd; -- cgit v1.2.3 From 6112bd15d73b00a07edd54eeaf9c94cdff93aa9e Mon Sep 17 00:00:00 2001 From: ROUVREAU Vincent Date: Tue, 12 Mar 2019 23:12:51 +0100 Subject: Add some documentation for Fix #31 Cythonization of Fix #31 Add some tests --- .../include/gudhi/Tangential_complex.h | 13 +++++++++- .../test/test_tangential_complex.cpp | 30 ++++++++++++++++++++++ src/cython/cython/tangential_complex.pyx | 26 +++++++++++++++++++ src/cython/doc/tangential_complex_user.rst | 2 ++ ...complex_plain_homology_from_off_file_example.py | 1 + src/cython/include/Tangential_complex_interface.h | 15 +++++++---- src/cython/test/test_tangential_complex.py | 9 +++++++ 7 files changed, 90 insertions(+), 6 deletions(-) (limited to 'src/cython/include') diff --git a/src/Tangential_complex/include/gudhi/Tangential_complex.h b/src/Tangential_complex/include/gudhi/Tangential_complex.h index 37cdf1b4..4a78127c 100644 --- a/src/Tangential_complex/include/gudhi/Tangential_complex.h +++ b/src/Tangential_complex/include/gudhi/Tangential_complex.h @@ -322,7 +322,11 @@ class Tangential_complex { for (std::size_t i = 0; i < m_points.size(); ++i) m_are_tangent_spaces_computed[i] = true; } - /// Computes the tangential complex. + /** \brief Computes the tangential complex. + * \exception std::invalid_argument In debug mode, if the computed star dimension is too low. Try to set a bigger + * maximal edge length value with `Tangential_complex::set_max_squared_edge_length` if + * this happens. + */ void compute_tangential_complex() { #ifdef GUDHI_TC_PERFORM_EXTRA_CHECKS std::cerr << red << "WARNING: GUDHI_TC_PERFORM_EXTRA_CHECKS is defined. " @@ -1984,6 +1988,13 @@ class Tangential_complex { return os; } + /** \brief Sets the maximal possible squared edge length for the edges in the triangulations. + * + * @param[in] max_squared_edge_length Maximal possible squared edge length. + * + * If the maximal edge length value is too low `Tangential_complex::compute_tangential_complex` will throw an + * exception in debug mode. + */ void set_max_squared_edge_length(FT max_squared_edge_length) { m_max_squared_edge_length = max_squared_edge_length; } private: diff --git a/src/Tangential_complex/test/test_tangential_complex.cpp b/src/Tangential_complex/test/test_tangential_complex.cpp index 4e2d4f65..103b8b30 100644 --- a/src/Tangential_complex/test/test_tangential_complex.cpp +++ b/src/Tangential_complex/test/test_tangential_complex.cpp @@ -126,3 +126,33 @@ BOOST_AUTO_TEST_CASE(test_mini_tangential) { BOOST_CHECK(stree.num_vertices() == 4); BOOST_CHECK(stree.num_simplices() == 6); } + +#ifdef GUDHI_DEBUG +BOOST_AUTO_TEST_CASE(test_basic_example_throw) { + typedef CGAL::Epick_d Kernel; + typedef Kernel::FT FT; + typedef Kernel::Point_d Point; + typedef Kernel::Vector_d Vector; + typedef tc::Tangential_complex TC; + + const int INTRINSIC_DIM = 2; + const int AMBIENT_DIM = 3; + const int NUM_POINTS = 1000; + + Kernel k; + + // Generate points on a 2-sphere + CGAL::Random_points_on_sphere_d generator(AMBIENT_DIM, 3.); + std::vector points; + points.reserve(NUM_POINTS); + for (int i = 0; i < NUM_POINTS; ++i) + points.push_back(*generator++); + + // Compute the TC + TC tc(points, INTRINSIC_DIM, k); + tc.set_max_squared_edge_length(0.01); + std::cout << "test_basic_example_throw - set_max_squared_edge_length(0.01) to make GUDHI_CHECK fail" << std::endl; + BOOST_CHECK_THROW(tc.compute_tangential_complex(), std::invalid_argument); + +} +#endif diff --git a/src/cython/cython/tangential_complex.pyx b/src/cython/cython/tangential_complex.pyx index 4bb07076..293ef8cb 100644 --- a/src/cython/cython/tangential_complex.pyx +++ b/src/cython/cython/tangential_complex.pyx @@ -36,6 +36,7 @@ cdef extern from "Tangential_complex_interface.h" namespace "Gudhi": Tangential_complex_interface(int intrisic_dim, vector[vector[double]] points) # bool from_file is a workaround for cython to find the correct signature Tangential_complex_interface(int intrisic_dim, string off_file, bool from_file) + void compute_tangential_complex() except + vector[double] get_point(unsigned vertex) unsigned number_of_vertices() unsigned number_of_simplices() @@ -43,6 +44,7 @@ cdef extern from "Tangential_complex_interface.h" namespace "Gudhi": unsigned number_of_inconsistent_stars() void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree) void fix_inconsistencies_using_perturbation(double max_perturb, double time_limit) + void set_max_squared_edge_length(double max_squared_edge_length) # TangentialComplex python interface cdef class TangentialComplex: @@ -92,6 +94,17 @@ cdef class TangentialComplex: """ return self.thisptr != NULL + def compute_tangential_complex(self): + """This function computes the tangential complex. + + Raises: + ValueError: In debug mode, if the computed star dimension is too + low. Try to set a bigger maximal edge length value with + :func:`~gudhi.Tangential_complex.set_max_squared_edge_length` + if this happens. + """ + self.thisptr.compute_tangential_complex() + def get_point(self, vertex): """This function returns the point corresponding to a given vertex. @@ -152,3 +165,16 @@ cdef class TangentialComplex: """ self.thisptr.fix_inconsistencies_using_perturbation(max_perturb, time_limit) + + def set_max_squared_edge_length(self, max_squared_edge_length): + """Sets the maximal possible squared edge length for the edges in the + triangulations. + + :param max_squared_edge_length: Maximal possible squared edge length. + :type max_squared_edge_length: double + + If the maximal edge length value is too low + :func:`~gudhi.Tangential_complex.compute_tangential_complex` + will throw an exception in debug mode. + """ + self.thisptr.set_max_squared_edge_length(max_squared_edge_length) diff --git a/src/cython/doc/tangential_complex_user.rst b/src/cython/doc/tangential_complex_user.rst index 5ce69e86..97471baf 100644 --- a/src/cython/doc/tangential_complex_user.rst +++ b/src/cython/doc/tangential_complex_user.rst @@ -128,6 +128,7 @@ This example builds the Tangential complex of point set read in an OFF file. import gudhi tc = gudhi.TangentialComplex(intrisic_dim = 1, off_file=gudhi.__root_source_dir__ + '/data/points/alphacomplexdoc.off') + tc.compute_tangential_complex() result_str = 'Tangential contains ' + repr(tc.num_simplices()) + \ ' simplices - ' + repr(tc.num_vertices()) + ' vertices.' print(result_str) @@ -175,6 +176,7 @@ simplices. import gudhi tc = gudhi.TangentialComplex(intrisic_dim = 1, points=[[0.0, 0.0], [1.0, 0.0], [0.0, 1.0], [1.0, 1.0]]) + tc.compute_tangential_complex() result_str = 'Tangential contains ' + repr(tc.num_vertices()) + ' vertices.' print(result_str) diff --git a/src/cython/example/tangential_complex_plain_homology_from_off_file_example.py b/src/cython/example/tangential_complex_plain_homology_from_off_file_example.py index 0f8f5e80..536517d1 100755 --- a/src/cython/example/tangential_complex_plain_homology_from_off_file_example.py +++ b/src/cython/example/tangential_complex_plain_homology_from_off_file_example.py @@ -50,6 +50,7 @@ with open(args.file, 'r') as f: print("TangentialComplex creation from points read in a OFF file") tc = gudhi.TangentialComplex(intrisic_dim = args.intrisic_dim, off_file=args.file) + tc.compute_tangential_complex() st = tc.create_simplex_tree() message = "Number of simplices=" + repr(st.num_simplices()) diff --git a/src/cython/include/Tangential_complex_interface.h b/src/cython/include/Tangential_complex_interface.h index 71418886..c4ddbdbe 100644 --- a/src/cython/include/Tangential_complex_interface.h +++ b/src/cython/include/Tangential_complex_interface.h @@ -49,8 +49,6 @@ class Tangential_complex_interface { Dynamic_kernel k; tangential_complex_ = new TC(points, intrisic_dim, k); - tangential_complex_->compute_tangential_complex(); - num_inconsistencies_ = tangential_complex_->number_of_inconsistent_simplices(); } Tangential_complex_interface(int intrisic_dim, const std::string& off_file_name, bool from_file = true) { @@ -60,14 +58,17 @@ class Tangential_complex_interface { std::vector points = off_reader.get_point_cloud(); tangential_complex_ = new TC(points, intrisic_dim, k); - tangential_complex_->compute_tangential_complex(); - num_inconsistencies_ = tangential_complex_->number_of_inconsistent_simplices(); } ~Tangential_complex_interface() { delete tangential_complex_; } + void compute_tangential_complex() { + tangential_complex_->compute_tangential_complex(); + num_inconsistencies_ = tangential_complex_->number_of_inconsistent_simplices(); + } + std::vector get_point(unsigned vh) { std::vector vd; if (vh < tangential_complex_->number_of_vertices()) { @@ -104,7 +105,11 @@ class Tangential_complex_interface { simplex_tree->initialize_filtration(); } - private: + void set_max_squared_edge_length(double max_squared_edge_length) { + tangential_complex_->set_max_squared_edge_length(max_squared_edge_length); + } + +private: TC* tangential_complex_; TC::Num_inconsistencies num_inconsistencies_; }; diff --git a/src/cython/test/test_tangential_complex.py b/src/cython/test/test_tangential_complex.py index 5385a0d3..5c62f278 100755 --- a/src/cython/test/test_tangential_complex.py +++ b/src/cython/test/test_tangential_complex.py @@ -32,6 +32,15 @@ def test_tangential(): tc = TangentialComplex(intrisic_dim = 1, points=point_list) assert tc.__is_defined() == True assert tc.num_vertices() == 4 + assert tc.num_simplices() == 0 + assert tc.num_inconsistent_simplices() == 0 + assert tc.num_inconsistent_stars() == 0 + + tc.compute_tangential_complex() + assert tc.num_vertices() == 4 + assert tc.num_simplices() == 4 + assert tc.num_inconsistent_simplices() == 0 + assert tc.num_inconsistent_stars() == 0 st = tc.create_simplex_tree() assert st.__is_defined() == True -- cgit v1.2.3 From 6525c78704489b0c8cb62b2e3f882ce6113c0f0d Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Tue, 16 Apr 2019 23:19:19 +0200 Subject: Add min and max filtration values for the sparse Rips --- src/Rips_complex/doc/Intro_rips_complex.h | 11 +++++---- .../include/gudhi/Sparse_rips_complex.h | 27 ++++++++++++++-------- src/cython/doc/rips_complex_user.rst | 4 ++-- src/cython/include/Rips_complex_interface.h | 13 +++-------- 4 files changed, 30 insertions(+), 25 deletions(-) (limited to 'src/cython/include') diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index 1aac804b..1d8cd2ba 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -99,10 +99,13 @@ namespace rips_complex { * \cite cavanna15geometric, and in a video \cite cavanna15visualizing. * * The interface of `Sparse_rips_complex` is similar to the one for the usual - * `Rips_complex`, except that one has to specify the approximation factor, and - * there is no option to limit the maximum filtration value (the way the - * approximation is done means that larger filtration values are much cheaper - * to handle than low filtration values, so the gain would be too small). + * `Rips_complex`, except that one has to specify the approximation factor. + * There is an option to limit the minimum and maximum filtration values, but + * they are not recommended: the way the approximation is done means that + * larger filtration values are much cheaper to handle than low filtration + * values, so the gain in ignoring the large ones is small, and + * `Gudhi::subsampling::sparsify_point_set()` is a more efficient way of + * ignoring small filtration values. * * Theoretical guarantees are only available for \f$\epsilon<1\f$. The * construction accepts larger values of ε, and the size of the complex diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index 1e8fb6a3..a249cd8e 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -70,17 +70,19 @@ class Sparse_rips_complex { * @param[in] points Range of points. * @param[in] distance Distance function that returns a `Filtration_value` from 2 given points. * @param[in] epsilon Approximation parameter. epsilon must be positive. + * @param[in] mini Minimal filtration value. Ignore anything below this scale. This is a less efficient version of `Gudhi::subsampling::sparsify_point_set()`. + * @param[in] maxi Maximal filtration value. Ignore anything above this scale. * */ template - Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon) + Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon, Filtration_value mini=-std::numeric_limits::infinity(), Filtration_value maxi=std::numeric_limits::infinity()) : 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, std::back_inserter(sorted_points), std::back_inserter(params)); - compute_sparse_graph(dist_fun, epsilon); + compute_sparse_graph(dist_fun, epsilon, mini, maxi); } /** \brief Sparse_rips_complex constructor from a distance matrix. @@ -90,11 +92,14 @@ class Sparse_rips_complex { * \f$j\f$ as long as \f$ 0 \leqslant j < i \leqslant * distance\_matrix.size().\f$ * @param[in] epsilon Approximation parameter. epsilon must be positive. + * @param[in] mini Minimal filtration value. Ignore anything below this scale. This is a less efficient version of `Gudhi::subsampling::sparsify_point_set()`. + * @param[in] maxi Maximal filtration value. Ignore anything above this scale. */ template - Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon) + Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon, Filtration_value mini=-std::numeric_limits::infinity(), Filtration_value maxi=std::numeric_limits::infinity()) : Sparse_rips_complex(boost::irange(0, boost::size(distance_matrix)), - [&](Vertex_handle i, Vertex_handle j) { return (i==j) ? 0 : (i - void compute_sparse_graph(Distance& dist, double epsilon) { + void compute_sparse_graph(Distance& dist, double epsilon, Filtration_value mini, Filtration_value maxi) { const auto& points = sorted_points; // convenience alias const int n = boost::size(points); graph_.~Graph(); @@ -164,13 +169,15 @@ class Sparse_rips_complex { // TODO(MG): // - make it parallel // - only test near-enough neighbors - for (int i = 0; i < n; ++i) + for (int i = 0; i < n; ++i) { + auto&& pi = points[i]; + auto li = params[i]; + if (li < mini) break; for (int j = i + 1; j < n; ++j) { - auto&& pi = points[i]; auto&& pj = points[j]; auto d = dist(pi, pj); - auto li = params[i]; auto lj = params[j]; + if (lj < mini) break; GUDHI_CHECK(lj <= li, "Bad furthest point sorting"); Filtration_value alpha; @@ -182,8 +189,10 @@ class Sparse_rips_complex { else continue; - add_edge(pi, pj, alpha, graph_); + if (alpha <= maxi) + add_edge(pi, pj, alpha, graph_); } + } } Graph graph_; diff --git a/src/cython/doc/rips_complex_user.rst b/src/cython/doc/rips_complex_user.rst index e814b4c3..1d340dbe 100644 --- a/src/cython/doc/rips_complex_user.rst +++ b/src/cython/doc/rips_complex_user.rst @@ -50,8 +50,8 @@ by more than the length used to define "too close". A more general technique is to use a sparse approximation of the Rips introduced by Don Sheehy :cite:`sheehy13linear`. We are using the version described in :cite:`buchet16efficient` (except that we multiply all filtration -values by 2, to match the usual Rips complex), which proves a -:math:`\frac{1+\varepsilon}{1-\varepsilon}`-interleaving, although in practice the +values by 2, to match the usual Rips complex). :cite:`cavanna15geometric` proves +a :math:`\frac{1}{1-\varepsilon}`-interleaving, although in practice the error is usually smaller. A more intuitive presentation of the idea is available in :cite:`cavanna15geometric`, and in a video :cite:`cavanna15visualizing`. Passing an extra argument `sparse=0.3` at the diff --git a/src/cython/include/Rips_complex_interface.h b/src/cython/include/Rips_complex_interface.h index 1a6e2477..40aff299 100644 --- a/src/cython/include/Rips_complex_interface.h +++ b/src/cython/include/Rips_complex_interface.h @@ -54,23 +54,17 @@ class Rips_complex_interface { } void init_points_sparse(const std::vector>& points, double threshold, double epsilon) { - sparse_rips_complex_.emplace(points, Gudhi::Euclidean_distance(), epsilon); - threshold_ = threshold; + sparse_rips_complex_.emplace(points, Gudhi::Euclidean_distance(), epsilon, -std::numeric_limits::infinity(), threshold); } void init_matrix_sparse(const std::vector>& matrix, double threshold, double epsilon) { - sparse_rips_complex_.emplace(matrix, epsilon); - threshold_ = threshold; + sparse_rips_complex_.emplace(matrix, epsilon, -std::numeric_limits::infinity(), threshold); } void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, int dim_max) { if (rips_complex_) rips_complex_->create_complex(*simplex_tree, dim_max); - else { + else sparse_rips_complex_->create_complex(*simplex_tree, dim_max); - // This pruning should be done much earlier! It isn't that useful for sparse Rips, - // but it would be inconsistent not to do it. - simplex_tree->prune_above_filtration(threshold_); - } simplex_tree->initialize_filtration(); } @@ -79,7 +73,6 @@ class Rips_complex_interface { // Anyway, storing a graph would make more sense. Or changing the interface completely so there is no such storage. boost::optional::Filtration_value>> rips_complex_; boost::optional::Filtration_value>> sparse_rips_complex_; - double threshold_; }; } // namespace rips_complex -- cgit v1.2.3