From baf00e2d21884bd3cc711e281ae77fe31e794b32 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Tue, 19 Mar 2019 12:05:05 +0100 Subject: Start fixing the sparse rips to match the true definition (not a clique complex) --- .../concept/SimplicialComplexForRips.h | 21 +++++++++++++++++++ src/Rips_complex/doc/Intro_rips_complex.h | 7 +++++-- .../include/gudhi/Sparse_rips_complex.h | 24 +++++++++++++++++++--- 3 files changed, 47 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/Rips_complex/concept/SimplicialComplexForRips.h b/src/Rips_complex/concept/SimplicialComplexForRips.h index 3c5acecf..36ab1b0c 100644 --- a/src/Rips_complex/concept/SimplicialComplexForRips.h +++ b/src/Rips_complex/concept/SimplicialComplexForRips.h @@ -34,6 +34,9 @@ struct SimplicialComplexForRips { /** \brief Type used to store the filtration values of the simplicial complex. */ typedef unspecified Filtration_value; + /** \brief Handle type to a simplex contained in the simplicial complex. */ + typedef unspecified Simplex_handle; + /** \brief Inserts a given `Gudhi::rips_complex::Rips_complex::OneSkeletonGraph` in the simplicial complex. */ template void insert_graph(const OneSkeletonGraph& skel_graph); @@ -42,6 +45,24 @@ struct SimplicialComplexForRips { * explained in \ref ripsdefinition. */ void expansion(int max_dim); + /** \brief Expands a simplicial complex containing only a graph. Simplices corresponding to cliques in the graph are added + * incrementally, faces before cofaces, unless the simplex has dimension larger than `max_dim` or `block_simplex` + * returns true for this simplex. + * + * @param[in] max_dim Expansion maximal dimension value. + * @param[in] block_simplex Blocker oracle. Its concept is bool block_simplex(Simplex_handle sh) + * + * The function identifies a candidate simplex whose faces are all already in the complex, inserts + * it with a filtration value corresponding to the maximum of the filtration values of the faces, then calls + * `block_simplex` on a `Simplex_handle` for this new simplex. If `block_simplex` returns true, the simplex is + * removed, otherwise it is kept. + */ + template< typename Blocker > + void expansion_with_blockers(int max_dim, Blocker block_simplex); + + /** \brief Returns a range over the vertices of a simplex. */ + unspecified simplex_vertex_range(Simplex_handle sh); + /** \brief Returns the number of vertices in the simplicial complex. */ std::size_t num_vertices(); diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index a2537036..1aac804b 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -92,8 +92,8 @@ namespace rips_complex { * The sparse Rips filtration was 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 - * \f$\frac{1+\epsilon}{1-\epsilon}\f$-interleaving, although in practice the + * by 2, to match the usual Rips complex), for which \cite cavanna15geometric proves a + * \f$(1,\frac{1}{1-\epsilon})\f$-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. @@ -107,6 +107,9 @@ namespace rips_complex { * Theoretical guarantees are only available for \f$\epsilon<1\f$. The * construction accepts larger values of ε, and the size of the complex * keeps decreasing, but there is no guarantee on the quality of the result. + * Note that while the number of edges decreases when ε increases, the + * number of higher-dimensional simplices may not be monotonous when + * \f$\frac12\leq\epsilon\leq 1\f$. * * \section ripspointsdistance Point cloud and distance function * diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index 00da148f..87d267d0 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -47,7 +47,9 @@ namespace rips_complex { * * \details * This class is used to construct a sparse \f$(1+O(\epsilon))\f$-approximation of `Rips_complex`, i.e. a filtered - * simplicial complex that is multiplicatively \f$(1+O(\epsilon))\f$-interleaved with the Rips filtration. + * simplicial complex that is multiplicatively + * \f$(1+O(\epsilon))\f$-interleaved with the Rips filtration. More precisely, + * this is a \f$(1,\frac{1}{1-\epsilon}\f$-interleaving. * * \tparam Filtration_value is the type used to store the filtration values of the simplicial complex. */ @@ -71,7 +73,8 @@ class Sparse_rips_complex { * */ template - Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon) { + Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon) + : epsilon_(epsilon) { GUDHI_CHECK(epsilon > 0, "epsilon must be positive"); std::vector sorted_points; std::vector params; @@ -111,7 +114,21 @@ class Sparse_rips_complex { std::invalid_argument("Sparse_rips_complex::create_complex - simplicial complex is not empty")); complex.insert_graph(graph_); - complex.expansion(dim_max); + if(epsilon_ >= 1) { + complex.expansion(dim_max); + return; + } + double cst = epsilon_ * (1 - epsilon_); + auto block = [=cst,&complex](typename SimplicialComplexForRips::Simplex_handle sh){ + auto filt = complex.filtration(sh); + auto mini = file * cst; + for(auto v : complex.simplex_vertex_range(sh)){ + if(lambda[v] < mini) // FIXME: store lambda/params somewhere!!! + return true; // v died before this simplex could be born + } + return false; + }; + complex.expansion_with_blockers(dim_max, block); } private: @@ -166,6 +183,7 @@ class Sparse_rips_complex { } Graph graph_; + double epsilon_; }; } // namespace rips_complex -- cgit v1.2.3 From ea4d23b73f74873723aade040a041da639eea09f Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Tue, 19 Mar 2019 12:28:31 +0100 Subject: Filter during expansion. Completely untested. --- src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 14 ++++++++++---- 1 file changed, 10 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index 87d267d0..7686d666 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -76,13 +76,11 @@ class Sparse_rips_complex { Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon) : epsilon_(epsilon) { GUDHI_CHECK(epsilon > 0, "epsilon must be positive"); - std::vector sorted_points; - std::vector params; 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(sorted_points, params, dist_fun, epsilon); + compute_sparse_graph(dist_fun, epsilon); } /** \brief Sparse_rips_complex constructor from a distance matrix. @@ -118,6 +116,10 @@ class Sparse_rips_complex { complex.expansion(dim_max); return; } + const int n = boost::size(params); + std::vector lambda(n); + for(int i=0;i - void compute_sparse_graph(const PointRange& points, const ParamRange& params, Distance& dist, double epsilon) { + void compute_sparse_graph(Distance& dist, double epsilon) { + const auto& points = sorted_points; // convenience alias const int n = boost::size(points); graph_.~Graph(); new (&graph_) Graph(n); @@ -184,6 +187,9 @@ class Sparse_rips_complex { Graph graph_; double epsilon_; + // Because of the arbitrary split between constructor and create_complex + std::vector sorted_points; + std::vector params; }; } // namespace rips_complex -- cgit v1.2.3 From 79c5d2acf5a9af5b11dac611fe7b81de75e53b89 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Tue, 19 Mar 2019 16:54:25 +0100 Subject: Fixes --- src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index 7686d666..11d0b4c4 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -121,9 +121,9 @@ class Sparse_rips_complex { for(int i=0;i + template void compute_sparse_graph(Distance& dist, double epsilon) { const auto& points = sorted_points; // convenience alias const int n = boost::size(points); -- cgit v1.2.3 From 9ceaa84d52d939a117d78b49fd19c8900387dadc Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Thu, 11 Apr 2019 22:58:00 +0200 Subject: comments --- src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 5 ++++- 1 file changed, 4 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index 11d0b4c4..015c94d9 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -118,6 +118,7 @@ class Sparse_rips_complex { } const int n = boost::size(params); std::vector lambda(n); + // lambda[original_order]=params[sorted_order] for(int i=0;i sorted_points; + // params[sorted_order]=distance to previous points std::vector params; }; -- cgit v1.2.3 From 0191e039a08ccd68b9916fafb2da98260e5b4efb Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Sat, 13 Apr 2019 00:08:23 +0200 Subject: Handle triangular matrix --- src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index 015c94d9..3aafc5ff 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -94,7 +94,7 @@ class Sparse_rips_complex { template Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon) : Sparse_rips_complex(boost::irange(0, boost::size(distance_matrix)), - [&](Vertex_handle i, Vertex_handle j) { return distance_matrix[j][i]; }, epsilon) {} + [&](Vertex_handle i, Vertex_handle j) { return (i==j) ? 0 : (i Date: Sat, 13 Apr 2019 00:11:01 +0200 Subject: Fix the documentation of what relation i and j satisfy in the lower triangular part of a matrix. --- src/Rips_complex/include/gudhi/Rips_complex.h | 2 +- src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/Rips_complex/include/gudhi/Rips_complex.h b/src/Rips_complex/include/gudhi/Rips_complex.h index e902e52c..ee100867 100644 --- a/src/Rips_complex/include/gudhi/Rips_complex.h +++ b/src/Rips_complex/include/gudhi/Rips_complex.h @@ -90,7 +90,7 @@ class Rips_complex { * @param[in] threshold Rips value. * * \tparam DistanceMatrix must have a `size()` method and on which `distance_matrix[i][j]` returns - * the distance between points \f$i\f$ and \f$j\f$ as long as \f$ 0 \leqslant i < j \leqslant + * the distance between points \f$i\f$ and \f$j\f$ as long as \f$ 0 \leqslant j < i \leqslant * distance\_matrix.size().\f$ */ template diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index 3aafc5ff..1e8fb6a3 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -87,7 +87,7 @@ class Sparse_rips_complex { * * @param[in] distance_matrix Range of range of distances. * `distance_matrix[i][j]` returns the distance between points \f$i\f$ and - * \f$j\f$ as long as \f$ 0 \leqslant i < j \leqslant + * \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. */ -- 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') 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 From a3f4f17948eb073623cfb4e50c6f6d36bb4c47ee Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Thu, 18 Apr 2019 11:32:26 +0200 Subject: Doc reformat. --- src/Rips_complex/doc/Intro_rips_complex.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h index 1d8cd2ba..97d66fbd 100644 --- a/src/Rips_complex/doc/Intro_rips_complex.h +++ b/src/Rips_complex/doc/Intro_rips_complex.h @@ -125,8 +125,8 @@ namespace rips_complex { * * \include Rips_complex/example_one_skeleton_rips_from_points.cpp * - * When launching (Rips maximal distance between 2 points is 12.0, is expanded until dimension 1 - one skeleton graph - * in other words): + * When launching (Rips maximal distance between 2 points is 12.0, is expanded + * until dimension 1 - one skeleton graph in other words): * * \code $> ./Rips_complex_example_one_skeleton_from_points * \endcode -- cgit v1.2.3 From d697f28afbcc6be4d5346b6a390b2549d7bf481f Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Fri, 10 May 2019 16:35:59 +0200 Subject: max_dim=INT_MAX for sparse Rips by default. --- src/Rips_complex/utilities/ripscomplex.md | 4 ++-- src/Rips_complex/utilities/sparse_rips_persistence.cpp | 4 ++-- 2 files changed, 4 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/Rips_complex/utilities/ripscomplex.md b/src/Rips_complex/utilities/ripscomplex.md index 108cdd50..03838085 100644 --- a/src/Rips_complex/utilities/ripscomplex.md +++ b/src/Rips_complex/utilities/ripscomplex.md @@ -85,7 +85,7 @@ properly, this is a known issue. ## sparse_rips_persistence ## This program computes the persistent homology with coefficient field *Z/pZ* -of a sparse (1+epsilon)/(1-epsilon)-approximation of the Rips complex defined on a set of input Euclidean points. The output diagram contains one bar per line, written with the convention: +of a sparse 1/(1-epsilon)-approximation of the Rips complex defined on a set of input Euclidean points. The output diagram contains one bar per line, written with the convention: `p dim birth death` @@ -100,7 +100,7 @@ where `dim` is the dimension of the homological feature, `birth` and `death` are * `-h [ --help ]` Produce help message * `-o [ --output-file ]` Name of file in which the persistence diagram is written. Default print in standard output. * `-e [ --approximation ]` (default = .5) Epsilon, where the sparse Rips complex is a (1+epsilon)/(1-epsilon)-approximation of the Rips complex. -* `-d [ --cpx-dimension ]` (default = 1) Maximal dimension of the Rips complex we want to compute. +* `-d [ --cpx-dimension ]` (default = INT_MAX) Maximal dimension of the Rips complex we want to compute. * `-p [ --field-charac ]` (default = 11) Characteristic p of the coefficient field Z/pZ for computing homology. * `-m [ --min-persistence ]` (default = 0) Minimal lifetime of homology feature to be recorded. Enter a negative value to see zero length intervals. diff --git a/src/Rips_complex/utilities/sparse_rips_persistence.cpp b/src/Rips_complex/utilities/sparse_rips_persistence.cpp index 6d4d86fd..3840d9f7 100644 --- a/src/Rips_complex/utilities/sparse_rips_persistence.cpp +++ b/src/Rips_complex/utilities/sparse_rips_persistence.cpp @@ -98,7 +98,7 @@ void program_options(int argc, char* argv[], std::string& off_file_points, std:: "Name of file in which the persistence diagram is written. Default print in std::cout")( "approximation,e", po::value(&epsilon)->default_value(.5), "Epsilon, where the sparse Rips complex is a (1+epsilon)-approximation of the Rips complex.")( - "cpx-dimension,d", po::value(&dim_max)->default_value(1), + "cpx-dimension,d", po::value(&dim_max)->default_value(std::numeric_limits::max()), "Maximal dimension of the Rips complex we want to compute.")( "field-charac,p", po::value(&p)->default_value(11), "Characteristic p of the coefficient field Z/pZ for computing homology.")( @@ -119,7 +119,7 @@ void program_options(int argc, char* argv[], std::string& off_file_points, std:: if (vm.count("help") || !vm.count("input-file")) { std::cout << std::endl; std::cout << "Compute the persistent homology with coefficient field Z/pZ \n"; - std::cout << "of a sparse (1+epsilon)-approximation of the Rips complex \ndefined on a set of input points.\n \n"; + std::cout << "of a sparse 1/(1-epsilon)-approximation of the Rips complex \ndefined on a set of input points.\n \n"; std::cout << "The output diagram contains one bar per line, written with the convention: \n"; std::cout << " p dim b d \n"; std::cout << "where dim is the dimension of the homological feature,\n"; -- cgit v1.2.3 From fe2f26629481faf316c74dafa7eae892bbc7a556 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Fri, 10 May 2019 22:00:54 +0200 Subject: Fix blocker formula I had copied it straight from the paper, without adding the factor 2 we have everywhere. --- src/Rips_complex/include/gudhi/Sparse_rips_complex.h | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) (limited to 'src') diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h index a249cd8e..8df6e387 100644 --- a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h +++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h @@ -126,7 +126,7 @@ class Sparse_rips_complex { // lambda[original_order]=params[sorted_order] for(int i=0;i= 1 || d * epsilon <= lj * (1 + 1 / (1 - epsilon)))) - alpha = (d - lj / epsilon) * 2; - else + else if (d * epsilon > li + lj) continue; + else { + alpha = (d - lj / epsilon) * 2; + // Keep the test exactly the same as in block to avoid inconsistencies + if (epsilon < 1 && alpha * cst > lj) + continue; + } if (alpha <= maxi) add_edge(pi, pj, alpha, graph_); -- cgit v1.2.3