Start fixing the sparse rips to match the true definition (not a clique complex)

3 files changed, 47 insertions, 5 deletions

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<class OneSkeletonGraph>
   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 <CODE>bool block_simplex(Simplex_handle sh)</CODE>
+   *
+   * 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 &epsilon;, 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 &epsilon; 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 <typename RandomAccessPointRange, typename Distance>
-  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<Vertex_handle> sorted_points;
     std::vector<Filtration_value> 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