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/Rips_complex') 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