Merge sparse rips branch in trunk

git-svn-id: svn+ssh://scm.gforge.inria.fr/svnroot/gudhi/trunk@3336 636b058d-ea47-450e-bf9e-a15bfbe3eedb

Former-commit-id: 01dcab006197637ccba9c29f8990375fe05ad7b4

12 files changed, 576 insertions, 42 deletions

diff --git a/biblio/bibliography.bib b/biblio/bibliography.bib
index b101cb76..16f43d6f 100644
--- a/biblio/bibliography.bib
+++ b/biblio/bibliography.bib
@@ -1,12 +1,12 @@
 @inproceedings{gudhilibrary_ICMS14,
-  author    = {Cl\'ement Maria and Jean-Daniel Boissonnat and 
+  author    = {Cl\'ement Maria and Jean-Daniel Boissonnat and
               Marc Glisse and Mariette Yvinec},
   title     = {The {G}udhi Library: Simplicial Complexes and Persistent Homology},
   booktitle = {ICMS},
   year      = {2014},
 }
 
-@article{Carriere17c,	
+@article{Carriere17c,
   author    = {Carri\`ere, Mathieu and Michel, Bertrand and Oudot, Steve},
   title     = {{Statistical Analysis and Parameter Selection for Mapper}},
   journal   = {CoRR},
@@ -20,7 +20,7 @@
  booktitle = {Proceedings of the Twenty-ninth Annual Symposium on Computational Geometry},
  year = {2013},
  pages = {107--116},
-} 
+}
 
 @article{Carriere16,
   title={{Structure and Stability of the 1-Dimensional Mapper}},
@@ -70,7 +70,7 @@ language={English}
  publisher = {ACM Press/Addison-Wesley Publishing Co.},
  address = {New York, NY, USA},
  keywords = {level of detail, mutiresolution modeling, non-manifold, pair contraction, surface simplification},
-} 
+}
 
 
 @inproceedings{Lindstrom,
@@ -85,7 +85,7 @@ language={English}
  numpages = {8},
  publisher = {IEEE Computer Society Press},
  address = {Los Alamitos, CA, USA},
-} 
+}
 
 
 @article{boissonnatmariacamalgorithmica,
@@ -201,13 +201,13 @@ language={English},
   bibsource = {DBLP, http://dblp.uni-trier.de}
 }
 @article{RS62,
-	author={J. B. Rosser and L. Schoenfeld},
-	title={Approximate Formulas for some Functions of Prime Numbers},
-	journal= {ijm},
-	volume= 6,
-	year= 1962,
-	pages={64-94},
-	mrnumber={25 \#1139}
+        author={J. B. Rosser and L. Schoenfeld},
+        title={Approximate Formulas for some Functions of Prime Numbers},
+        journal= {ijm},
+        volume= 6,
+        year= 1962,
+        pages={64-94},
+        mrnumber={25 \#1139}
 }
 @article{DBLP:journals/siamcomp/Furer09,
   author    = {Martin F{\"u}rer},
@@ -410,7 +410,7 @@ url="http://dx.doi.org/10.1007/s00454-013-9557-2"
  edition = {2},
  publisher = {Cambridge University Press},
  address = {New York, NY, USA},
-} 
+}
 @book{DBLP:books/aw/AhoHU74,
   author    = {Alfred V. Aho and
                John E. Hopcroft and
@@ -508,7 +508,7 @@ proceedings{DBLP:conf/issac/1996,
  pages = {501--509},
  numpages = {9},
  keywords = {clique complexes, data structure, edge contraction, flag complexes, high dimensions, homotopy equivalence, shape reconstruction, shape simplification, simplicial complexes, vietoris-rips complexes},
-} 
+}
 
 
 
@@ -571,7 +571,7 @@ note = "http://gmplib.org/",
  publisher = {ACM},
  address = {New York, NY, USA},
  keywords = {discrete and computational geometry, embedding theorem, persistent homology},
-} 
+}
 @article{Fredman:1984:SST:828.1884,
  author = {Fredman, Michael L. and Koml\'{o}s, J\'{a}nos and Szemer{\'e}di, Endre},
  title = {Storing a Sparse Table with {O}(1) Worst Case Access Time},
@@ -589,7 +589,7 @@ note = "http://gmplib.org/",
  acmid = {1884},
  publisher = {ACM},
  address = {New York, NY, USA},
-} 
+}
 %------------------------------------------------------------------
 
 
@@ -739,7 +739,7 @@ inproceedings{DBLP:conf/soda/1997,
  publisher = {Kluwer Academic Publishers},
  address = {Hingham, MA, USA},
  keywords = {Filtration, Klein bottle, Manifold, Natural images, Persistent homology, Topology},
-} 
+}
 
 article{DBLP:journals/ijcv/LeePM03,
   author    = {Ann B. Lee and
@@ -767,7 +767,7 @@ pages={234115},
 year={2010},
 abstract={Understanding energy landscapes is a major challenge in chemistry and biology. Although a wide variety of methods have been invented and applied to this problem, very little is understood about the actual mathematical structures underlying such landscapes. Perhaps the most general assumption is the idea that energy landscapes are low-dimensional manifolds embedded in high-dimensional Euclidean space. While this is a very mild assumption, we have discovered an example of an energy landscape which is nonmanifold, demonstrating previously unknown mathematical complexity. The example occurs in the energy landscape of cyclo-octane, which was found to have the structure of a reducible algebraic variety, composed of the union of a sphere and a Klein bottle, intersecting in two rings.}
 }
- 
+
 
 book{hatcher2002algebraic,
   title={Algebraic Topology},
@@ -803,7 +803,7 @@ book{hatcher2002algebraic,
  isbn = {0070131511},
  edition = {2nd},
  publisher = {McGraw-Hill Higher Education},
-} 
+}
 
 @article{DBLP:journals/focm/CarlssonS10,
   author    = {Gunnar E. Carlsson and
@@ -1017,7 +1017,7 @@ language={English}
     journal = {Journal of Symbolic Computation.},
     year = {2016}
 }
-                         
+
 @ARTICLE{Fasy_Kim_Lecci_Maria_tda,
     author = {B. Fasy and J. Kim and F. Lecci and C. Maria},
     title = {Introduction to the R package TDA.},
@@ -1064,11 +1064,42 @@ language={English}
 }
 
 @ARTICLE{Carriere_Oudot_Ovsjanikov_top_signatures_3d,
-    author = {M. Carrière and S. Oudot and M. Ovsjanikov},
+    author = {M. Carri\`ere and S. Oudot and M. Ovsjanikov},
     title = {Stable Topological Signatures for Points on 3D Shapes.},
     journal = {Proc. Sympos. on Geometry Processing},
     year = {2015}
 }
-           
 
+@article{buchet16efficient,
+  title = {Efficient and Robust Persistent Homology for Measures},
+  author = {Micka\"{e}l Buchet and Fr\'{e}d\'{e}ric Chazal and Steve Y. Oudot and Donald Sheehy},
+  booktitle = {Computational Geometry: Theory and Applications},
+  volume = {58},
+  pages = {70--96},
+  year = {2016}
+}
+
+@inproceedings{cavanna15geometric,
+  author = {Nicholas J. Cavanna and Mahmoodreza Jahanseir and Donald R. Sheehy},
+  booktitle = {Proceedings of the Canadian Conference on Computational Geometry},
+  title = {A Geometric Perspective on Sparse Filtrations},
+  year = {2015}
+}
+
+@inproceedings{cavanna15visualizing,
+  author = {Nicholas J. Cavanna and Mahmoodreza Jahanseir and Donald R. Sheehy},
+  booktitle = {Proceedings of the 31st International Symposium on Computational Geometry},
+  title = {Visualizing Sparse Filtrations},
+  year = {2015},
+  doi = {10.4230/LIPIcs.SOCG.2015.23}
+}
 
+@article{sheehy13linear,
+  title = {Linear-Size Approximations to the {V}ietoris-{R}ips Filtration},
+  author = {Donald R. Sheehy},
+  journal = {Discrete \& Computational Geometry},
+  volume = {49},
+  number = {4},
+  pages = {778--796},
+  year = {2013}
+}
diff --git a/biblio/how_to_cite_gudhi.bib b/biblio/how_to_cite_gudhi.bib
index 5994124a..942f8d7e 100644
--- a/biblio/how_to_cite_gudhi.bib
+++ b/biblio/how_to_cite_gudhi.bib
@@ -97,7 +97,7 @@
 }
 
 @incollection{gudhi:RipsComplex
-, author =  "Cl\'ement Maria, Pawel Dlotko, Vincent Rouvreau"
+, author =  "Cl\'ement Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse"
 , title =   "Rips complex"
 , publisher =  "{GUDHI Editorial Board}"
 , booktitle =   "{GUDHI} User and Reference Manual"
diff --git a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
index e0a147b3..a8c9afa3 100644
--- a/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
+++ b/src/Persistent_cohomology/include/gudhi/Persistent_cohomology.h
@@ -285,7 +285,7 @@ class Persistent_cohomology {
         }
       }
       cpx_->assign_key(sigma, cpx_->null_key());
-    } else {  // If ku == kv, same connected component: create a 1-cocycle class.
+    } else if (dim_max_ > 1) {  // If ku == kv, same connected component: create a 1-cocycle class.
       create_cocycle(sigma, coeff_field_.multiplicative_identity(), coeff_field_.characteristic());
     }
   }
diff --git a/src/Rips_complex/doc/Intro_rips_complex.h b/src/Rips_complex/doc/Intro_rips_complex.h
index 496c4218..5a551e60 100644
--- a/src/Rips_complex/doc/Intro_rips_complex.h
+++ b/src/Rips_complex/doc/Intro_rips_complex.h
@@ -29,28 +29,41 @@ namespace rips_complex {
 
 /**  \defgroup rips_complex Rips complex
  * 
- * \author    Clément Maria, Pawel Dlotko, Vincent Rouvreau
+ * \author    Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse
  * 
  * @{
  * 
  * \section ripsdefinition Rips complex definition
  * 
- * Rips_complex
- * <a target="_blank" href="https://en.wikipedia.org/wiki/Vietoris%E2%80%93Rips_complex">(Wikipedia)</a> is a
- * one skeleton graph that allows to construct a
- * <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a>
- * from it.
- * The input can be a point cloud with a given distance function, or a distance matrix.
- * 
- * The filtration value of each edge is computed from a user-given distance function, or directly from the distance
- * matrix.
+ * The Vietoris-Rips complex
+ * <a target="_blank" href="https://en.wikipedia.org/wiki/Vietoris%E2%80%93Rips_complex">(Wikipedia)</a>
+ * is an abstract simplicial complex
+ * defined on a finite metric space, where each simplex corresponds to a subset
+ * of point whose diameter is smaller that some given threshold.
+ * Varying the threshold, we can also see the Rips complex as a filtration of
+ * the \f$(n-1)-\f$dimensional simplex, where the filtration value of each
+ * simplex is the diameter of the corresponding subset of points.
+ *
+ * This filtered complex is most often used as an approximation of the
+ * Čech complex. After rescaling (Rips using the length of the edges and Čech
+ * the half-length), they share the same 1-skeleton and are multiplicatively
+ * 2-interleaved or better. While it is slightly bigger, it is also much
+ * easier to compute.
+ *
+ * The number of simplices in the full Rips complex is exponential in the
+ * number of vertices, it is thus usually restricted, by excluding all the
+ * simplices with filtration value larger than some threshold, and keeping only
+ * the dim_max-skeleton.
+ *
+ * In order to build this complex, the algorithm first builds the graph.
+ * The filtration value of each edge is computed from a user-given distance
+ * function, or directly read from the distance matrix.
+ * In a second step, this graph is inserted in a simplicial complex, which then
+ * gets expanded to a flag complex.
  * 
- * All edges that have a filtration value strictly greater than a given threshold value are not inserted into
- * the complex.
+ * The input can be given as a range of points and a distance function, or as a
+ * distance matrix.
  * 
- * When creating a simplicial complex from this one skeleton graph, Rips inserts the one skeleton graph into the data
- * structure, and then expands the simplicial complex when required.
- *
  * Vertex name correspond to the index of the point in the given range (aka. the point cloud).
  * 
  * \image html "rips_complex_representation.png" "Rips-complex one skeleton graph representation"
@@ -61,7 +74,36 @@ namespace rips_complex {
  * 
  * If the Rips_complex interfaces are not detailed enough for your need, please refer to
  * <a href="_persistent_cohomology_2rips_persistence_step_by_step_8cpp-example.html">
- * rips_persistence_step_by_step.cpp</a> example, where the graph construction over the Simplex_tree is more detailed.
+ * rips_persistence_step_by_step.cpp</a> example, where the constructions of the graph and
+ * the Simplex_tree are more detailed.
+ *
+ * \section sparserips Sparse Rips complex
+ *
+ * Even truncated in filtration value and dimension, the Rips complex remains
+ * quite large. However, it is possible to approximate it by a much smaller
+ * filtered simplicial complex (linear size, with constants that depend on
+ * &epsilon; and the doubling dimension of the space) that is
+ * \f$(1+O(\epsilon))-\f$interleaved with it (in particular, their persistence
+ * diagrams are at log-bottleneck distance at most \f$O(\epsilon)\f$).
+ *
+ * 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
+ * error is usually smaller.
+ * A more intuitive presentation of the idea is available in
+ * \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).
+ *
+ * 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.
  *
  * \section ripspointsdistance Point cloud and distance function
  * 
@@ -104,6 +146,24 @@ namespace rips_complex {
  * 
  * \include Rips_complex/full_skeleton_rips_for_doc.txt
  * 
+ *
+ * \subsection sparseripspointscloudexample Example of a sparse Rips from a point cloud
+ * 
+ * This example builds the full sparse Rips of a set of 2D Euclidean points, then prints some minimal
+ * information about the complex.
+ * 
+ * \include Rips_complex/example_sparse_rips.cpp
+ * 
+ * When launching:
+ * 
+ * \code $> ./Rips_complex_example_sparse
+ * \endcode
+ *
+ * the program output may be (the exact output varies from one run to the next):
+ *
+ * \code Sparse Rips complex is of dimension 2 - 19 simplices - 7 vertices.
+ * \endcode
+ * 
  * 
  * 
  * \section ripsdistancematrix Distance matrix
diff --git a/src/Rips_complex/example/CMakeLists.txt b/src/Rips_complex/example/CMakeLists.txt
index fcb1eaee..af86636b 100644
--- a/src/Rips_complex/example/CMakeLists.txt
+++ b/src/Rips_complex/example/CMakeLists.txt
@@ -11,6 +11,9 @@ add_executable ( Rips_complex_example_one_skeleton_from_distance_matrix example_
 
 add_executable ( Rips_complex_example_from_csv_distance_matrix example_rips_complex_from_csv_distance_matrix_file.cpp )
 
+# Sparse rips from points
+add_executable ( Rips_complex_example_sparse example_sparse_rips.cpp )
+
 # Correlation matrix
 add_executable (  Rips_complex_example_one_skeleton_rips_from_correlation_matrix example_one_skeleton_rips_from_correlation_matrix.cpp )
 
@@ -19,6 +22,7 @@ if (TBB_FOUND)
   target_link_libraries(Rips_complex_example_one_skeleton_from_points ${TBB_LIBRARIES})
   target_link_libraries(Rips_complex_example_one_skeleton_from_distance_matrix ${TBB_LIBRARIES})
   target_link_libraries(Rips_complex_example_from_csv_distance_matrix ${TBB_LIBRARIES})
+  target_link_libraries(Rips_complex_example_sparse ${TBB_LIBRARIES})
   target_link_libraries(Rips_complex_example_one_skeleton_rips_from_correlation_matrix ${TBB_LIBRARIES})
 endif()
 
@@ -26,6 +30,8 @@ add_test(NAME Rips_complex_example_one_skeleton_from_points
     COMMAND $<TARGET_FILE:Rips_complex_example_one_skeleton_from_points>)
 add_test(NAME Rips_complex_example_one_skeleton_from_distance_matrix
     COMMAND $<TARGET_FILE:Rips_complex_example_one_skeleton_from_distance_matrix>)
+add_test(NAME Rips_complex_example_sparse
+    COMMAND $<TARGET_FILE:Rips_complex_example_sparse>)
 add_test(NAME Rips_complex_example_one_skeleton_rips_from_correlation_matrix
     COMMAND $<TARGET_FILE:Rips_complex_example_one_skeleton_rips_from_correlation_matrix>)
 
diff --git a/src/Rips_complex/example/example_sparse_rips.cpp b/src/Rips_complex/example/example_sparse_rips.cpp
new file mode 100644
index 00000000..1c95b48c
--- /dev/null
+++ b/src/Rips_complex/example/example_sparse_rips.cpp
@@ -0,0 +1,30 @@
+#include <gudhi/Sparse_rips_complex.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/distance_functions.h>
+
+#include <iostream>
+#include <vector>
+
+int main() {
+  using Point = std::vector<double>;
+  using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
+  using Filtration_value = Simplex_tree::Filtration_value;
+  using Sparse_rips = Gudhi::rips_complex::Sparse_rips_complex<Filtration_value>;
+
+  Point points[] = {{1.0, 1.0}, {7.0, 0.0}, {4.0, 6.0}, {9.0, 6.0}, {0.0, 14.0}, {2.0, 19.0}, {9.0, 17.0}};
+
+  // ----------------------------------------------------------------------------
+  // Init from Euclidean points
+  // ----------------------------------------------------------------------------
+  double epsilon = 2;  // very rough, no guarantees
+  Sparse_rips sparse_rips(points, Gudhi::Euclidean_distance(), epsilon);
+
+  Simplex_tree stree;
+  sparse_rips.create_complex(stree, 10);
+
+  // ----------------------------------------------------------------------------
+  // Display information about the complex
+  // ----------------------------------------------------------------------------
+  std::cout << "Sparse Rips complex is of dimension " << stree.dimension() << " - " << stree.num_simplices()
+            << " simplices - " << stree.num_vertices() << " vertices." << std::endl;
+}
diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h
new file mode 100644
index 00000000..1a9d6ebb
--- /dev/null
+++ b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h
@@ -0,0 +1,175 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s): Marc Glisse
+ *
+ *    Copyright (C) 2018  INRIA
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#ifndef SPARSE_RIPS_COMPLEX_H_
+#define SPARSE_RIPS_COMPLEX_H_
+
+#include <gudhi/Debug_utils.h>
+#include <gudhi/graph_simplicial_complex.h>
+#include <gudhi/choose_n_farthest_points.h>
+
+#include <boost/graph/adjacency_list.hpp>
+#include <boost/range/metafunctions.hpp>
+
+#include <vector>
+
+namespace Gudhi {
+
+namespace rips_complex {
+
+// The whole interface is copied on Rips_complex. A redesign should be discussed with all complex creation classes in
+// mind.
+
+/**
+ * \class Sparse_rips_complex
+ * \brief Sparse Rips complex data structure.
+ *
+ * \ingroup 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.
+ *
+ * \tparam Filtration_value is the type used to store the filtration values of the simplicial complex.
+ */
+template <typename Filtration_value>
+class Sparse_rips_complex {
+ private:
+  // TODO: use a different graph where we know we can safely insert in parallel.
+  typedef typename boost::adjacency_list<boost::vecS, boost::vecS, boost::undirectedS,
+                                         boost::property<vertex_filtration_t, Filtration_value>,
+                                         boost::property<edge_filtration_t, Filtration_value>>
+      Graph;
+
+  typedef int Vertex_handle;
+
+ public:
+  /** \brief Sparse_rips_complex constructor from a list of points.
+   *
+   * @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.
+   *
+   */
+  template <typename RandomAccessPointRange, typename Distance>
+  Sparse_rips_complex(const RandomAccessPointRange& points, Distance distance, double epsilon) {
+    GUDHI_CHECK(epsilon > 0, "epsilon must be positive");
+    std::vector<Vertex_handle> sorted_points;
+    std::vector<Filtration_value> params;
+    auto dist_fun = [&](Vertex_handle i, Vertex_handle j) { return distance(points[i], points[j]); };
+    Ker<decltype(dist_fun)> kernel(dist_fun);
+    subsampling::choose_n_farthest_points(kernel, boost::irange<Vertex_handle>(0, boost::size(points)), -1, -1,
+                                          std::back_inserter(sorted_points), std::back_inserter(params));
+    compute_sparse_graph(sorted_points, params, dist_fun, epsilon);
+  }
+
+  /** \brief Sparse_rips_complex constructor from a distance matrix.
+   *
+   * @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
+   * distance\_matrix.size().\f$
+   * @param[in] epsilon Approximation parameter. epsilon must be positive.
+   */
+  template <typename DistanceMatrix>
+  Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon)
+      : Sparse_rips_complex(boost::irange<Vertex_handle>(0, boost::size(distance_matrix)),
+                            [&](Vertex_handle i, Vertex_handle j) { return distance_matrix[j][i]; }, epsilon) {}
+
+  /** \brief Fills the simplicial complex with the sparse Rips graph and
+   * expands it with all the cliques, stopping at a given maximal dimension.
+   *
+   * \tparam SimplicialComplexForRips must meet `SimplicialComplexForRips` concept.
+   *
+   * @param[in] complex the complex to fill
+   * @param[in] dim_max maximal dimension of the simplicial complex.
+   * @exception std::invalid_argument In debug mode, if `complex.num_vertices()` does not return 0.
+   *
+   */
+  template <typename SimplicialComplexForRips>
+  void create_complex(SimplicialComplexForRips& complex, int dim_max) {
+    GUDHI_CHECK(complex.num_vertices() == 0,
+                std::invalid_argument("Sparse_rips_complex::create_complex - simplicial complex is not empty"));
+
+    complex.insert_graph(graph_);
+    complex.expansion(dim_max);
+  }
+
+ private:
+  // choose_n_farthest_points wants the distance function in this form...
+  template <class Distance>
+  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 <typename PointRange, typename ParamRange, typename Distance>
+  void compute_sparse_graph(const PointRange& points, const ParamRange& params, Distance& dist, double epsilon) {
+    const int n = boost::size(points);
+    graph_.~Graph();
+    new (&graph_) Graph(n);
+    // for(auto v : vertices(g)) // doesn't work :-(
+    typename boost::graph_traits<Graph>::vertex_iterator v_i, v_e;
+    for (std::tie(v_i, v_e) = vertices(graph_); v_i != v_e; ++v_i) {
+      auto v = *v_i;
+      // This whole loop might not be necessary, leave it until someone investigates if it is safe to remove.
+      put(vertex_filtration_t(), graph_, v, 0);
+    }
+
+    // TODO:
+    // - make it parallel
+    // - only test near-enough neighbors
+    for (int i = 0; i < n; ++i)
+      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];
+        GUDHI_CHECK(lj <= li, "Bad furthest point sorting");
+        Filtration_value alpha;
+
+        // The paper has d/2 and d-lj/e to match the Cech, but we use doubles to match the Rips
+        if (d * epsilon <= 2 * lj)
+          alpha = d;
+        else if (d * epsilon <= li + lj && (epsilon >= 1 || d * epsilon <= lj * (1 + 1 / (1 - epsilon))))
+          alpha = (d - lj / epsilon) * 2;
+        else
+          continue;
+
+        add_edge(pi, pj, alpha, graph_);
+      }
+  }
+
+  Graph graph_;
+};
+
+}  // namespace rips_complex
+
+}  // namespace Gudhi
+
+#endif  // SPARSE_RIPS_COMPLEX_H_
diff --git a/src/Rips_complex/test/test_rips_complex.cpp b/src/Rips_complex/test/test_rips_complex.cpp
index 89afbc25..4e7b79d2 100644
--- a/src/Rips_complex/test/test_rips_complex.cpp
+++ b/src/Rips_complex/test/test_rips_complex.cpp
@@ -31,6 +31,7 @@
 #include <algorithm>    // std::max
 
 #include <gudhi/Rips_complex.h>
+#include <gudhi/Sparse_rips_complex.h>
 // to construct Rips_complex from a OFF file of points
 #include <gudhi/Points_off_io.h>
 #include <gudhi/Simplex_tree.h>
@@ -43,6 +44,7 @@ using Point = std::vector<double>;
 using Simplex_tree = Gudhi::Simplex_tree<>;
 using Filtration_value = Simplex_tree::Filtration_value;
 using Rips_complex = Gudhi::rips_complex::Rips_complex<Simplex_tree::Filtration_value>;
+using Sparse_rips_complex = Gudhi::rips_complex::Sparse_rips_complex<Simplex_tree::Filtration_value>;
 using Distance_matrix = std::vector<std::vector<Filtration_value>>;
 
 BOOST_AUTO_TEST_CASE(RIPS_DOC_OFF_file) {
@@ -230,6 +232,71 @@ BOOST_AUTO_TEST_CASE(Rips_complex_from_points) {
   }
 }
 
+BOOST_AUTO_TEST_CASE(Sparse_rips_complex_from_points) {
+  // This is a clone of the test above
+  // ----------------------------------------------------------------------------
+  // Init of a list of points
+  // ----------------------------------------------------------------------------
+  Vector_of_points points;
+  std::vector<double> coords = { 0.0, 0.0, 0.0, 1.0 };
+  points.push_back(Point(coords.begin(), coords.end()));
+  coords = { 0.0, 0.0, 1.0, 0.0 };
+  points.push_back(Point(coords.begin(), coords.end()));
+  coords = { 0.0, 1.0, 0.0, 0.0 };
+  points.push_back(Point(coords.begin(), coords.end()));
+  coords = { 1.0, 0.0, 0.0, 0.0 };
+  points.push_back(Point(coords.begin(), coords.end()));
+
+  // ----------------------------------------------------------------------------
+  // Init of a Rips complex from the list of points
+  // ----------------------------------------------------------------------------
+  // .001 is small enough that we get a deterministic result matching the exact Rips
+  Sparse_rips_complex sparse_rips(points, Custom_square_euclidean_distance(), .001);
+
+  std::cout << "========== Sparse_rips_complex_from_points ==========" << std::endl;
+  Simplex_tree st;
+  const int DIMENSION = 3;
+  sparse_rips.create_complex(st, DIMENSION);
+
+  // Another way to check num_simplices
+  std::cout << "Iterator on Rips complex simplices in the filtration order, with [filtration value]:" << std::endl;
+  int num_simplices = 0;
+  for (auto f_simplex : st.filtration_simplex_range()) {
+    num_simplices++;
+    std::cout << "   ( ";
+    for (auto vertex : st.simplex_vertex_range(f_simplex)) {
+      std::cout << vertex << " ";
+    }
+    std::cout << ") -> " << "[" << st.filtration(f_simplex) << "] ";
+    std::cout << std::endl;
+  }
+  BOOST_CHECK(num_simplices == 15);
+  std::cout << "st.num_simplices()=" << st.num_simplices() << std::endl;
+  BOOST_CHECK(st.num_simplices() == 15);
+
+  std::cout << "st.dimension()=" << st.dimension() << std::endl;
+  BOOST_CHECK(st.dimension() == DIMENSION);
+  std::cout << "st.num_vertices()=" << st.num_vertices() << std::endl;
+  BOOST_CHECK(st.num_vertices() == 4);
+
+  for (auto f_simplex : st.filtration_simplex_range()) {
+    std::cout << "dimension(" << st.dimension(f_simplex) << ") - f = " << st.filtration(f_simplex) << std::endl;
+    switch (st.dimension(f_simplex)) {
+      case 0:
+        GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 0.0);
+        break;
+      case 1:
+      case 2:
+      case 3:
+        GUDHI_TEST_FLOAT_EQUALITY_CHECK(st.filtration(f_simplex), 2.0);
+        break;
+      default:
+        BOOST_CHECK(false);  // Shall not happen
+        break;
+    }
+  }
+}
+
 BOOST_AUTO_TEST_CASE(Rips_doc_csv_file) {
   // ----------------------------------------------------------------------------
   //
diff --git a/src/Rips_complex/utilities/CMakeLists.txt b/src/Rips_complex/utilities/CMakeLists.txt
index b99fc808..deb73ff0 100644
--- a/src/Rips_complex/utilities/CMakeLists.txt
+++ b/src/Rips_complex/utilities/CMakeLists.txt
@@ -10,20 +10,26 @@ target_link_libraries(rips_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY})
 add_executable(rips_correlation_matrix_persistence rips_correlation_matrix_persistence.cpp)
 target_link_libraries(rips_correlation_matrix_persistence ${Boost_SYSTEM_LIBRARY} ${Boost_PROGRAM_OPTIONS_LIBRARY})
 
+add_executable(sparse_rips_persistence sparse_rips_persistence.cpp)
+target_link_libraries(sparse_rips_persistence ${Boost_PROGRAM_OPTIONS_LIBRARY})
+
 if (TBB_FOUND)
   target_link_libraries(rips_distance_matrix_persistence ${TBB_LIBRARIES})
   target_link_libraries(rips_persistence ${TBB_LIBRARIES})
   target_link_libraries(rips_correlation_matrix_persistence ${TBB_LIBRARIES})
+  target_link_libraries(sparse_rips_persistence ${TBB_LIBRARIES})
 endif()
 
 add_test(NAME Rips_complex_utility_from_rips_distance_matrix COMMAND $<TARGET_FILE:rips_distance_matrix_persistence>
     "${CMAKE_SOURCE_DIR}/data/distance_matrix/full_square_distance_matrix.csv" "-r" "1.0" "-d" "3" "-p" "3" "-m" "0")
 add_test(NAME Rips_complex_utility_from_rips_on_tore_3D COMMAND $<TARGET_FILE:rips_persistence>
     "${CMAKE_SOURCE_DIR}/data/points/tore3D_1307.off" "-r" "0.25" "-m" "0.5" "-d" "3" "-p" "3")
-
 add_test(NAME Rips_complex_utility_from_rips_correlation_matrix COMMAND $<TARGET_FILE:rips_correlation_matrix_persistence>
     "${CMAKE_SOURCE_DIR}/data/correlation_matrix/lower_triangular_correlation_matrix.csv" "-c" "0.3" "-d" "3" "-p" "3" "-m" "0")
+add_test(NAME Sparse_rips_complex_utility_on_tore_3D COMMAND $<TARGET_FILE:sparse_rips_persistence>
+    "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off" "-e" "0.5" "-m" "0.2" "-d" "3" "-p" "2")
 
 install(TARGETS rips_distance_matrix_persistence DESTINATION bin)
 install(TARGETS rips_persistence DESTINATION bin)
 install(TARGETS rips_correlation_matrix_persistence DESTINATION bin)
+install(TARGETS sparse_rips_persistence DESTINATION bin)
diff --git a/src/Rips_complex/utilities/ripscomplex.md b/src/Rips_complex/utilities/ripscomplex.md
index d4dbc36e..84f40cc1 100644
--- a/src/Rips_complex/utilities/ripscomplex.md
+++ b/src/Rips_complex/utilities/ripscomplex.md
@@ -83,3 +83,29 @@ Please refer to data/correlation_matrix/lower_triangular_correlation_matrix.csv
 
 As persistence diagrams points will be under the diagonal, bottleneck distance and persistence graphical tool will not work
 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)-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`
+
+where `dim` is the dimension of the homological feature, `birth` and `death` are respectively the birth and death of the feature, and `p` is the characteristic of the field *Z/pZ* used for homology coefficients (`p` must be a prime number).
+
+**Usage**
+
+`sparse_rips_persistence [options] <OFF input file>`
+
+**Allowed options**
+
+* `-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)-approximation of the Rips complex.
+* `-d [ --cpx-dimension ]` (default = 1) 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.
+
+**Example with Z/2Z coefficients**
+
+`sparse_rips_persistence ../../data/points/tore3D_1307.off -e .5 -m .2 -d 3 -p 2`
diff --git a/src/Rips_complex/utilities/sparse_rips_persistence.cpp b/src/Rips_complex/utilities/sparse_rips_persistence.cpp
new file mode 100644
index 00000000..d4bae3ba
--- /dev/null
+++ b/src/Rips_complex/utilities/sparse_rips_persistence.cpp
@@ -0,0 +1,133 @@
+/*    This file is part of the Gudhi Library. The Gudhi library
+ *    (Geometric Understanding in Higher Dimensions) is a generic C++
+ *    library for computational topology.
+ *
+ *    Author(s): Marc Glisse, Clément Maria
+ *
+ *    Copyright (C) 2018  INRIA
+ *
+ *    This program is free software: you can redistribute it and/or modify
+ *    it under the terms of the GNU General Public License as published by
+ *    the Free Software Foundation, either version 3 of the License, or
+ *    (at your option) any later version.
+ *
+ *    This program is distributed in the hope that it will be useful,
+ *    but WITHOUT ANY WARRANTY; without even the implied warranty of
+ *    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+ *    GNU General Public License for more details.
+ *
+ *    You should have received a copy of the GNU General Public License
+ *    along with this program.  If not, see <http://www.gnu.org/licenses/>.
+ */
+
+#include <gudhi/Sparse_rips_complex.h>
+#include <gudhi/distance_functions.h>
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Persistent_cohomology.h>
+#include <gudhi/Points_off_io.h>
+
+#include <boost/program_options.hpp>
+
+#include <string>
+#include <vector>
+
+// Types definition
+using Simplex_tree = Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_fast_persistence>;
+using Filtration_value = Simplex_tree::Filtration_value;
+using Sparse_rips = Gudhi::rips_complex::Sparse_rips_complex<Filtration_value>;
+using Field_Zp = Gudhi::persistent_cohomology::Field_Zp;
+using Persistent_cohomology = Gudhi::persistent_cohomology::Persistent_cohomology<Simplex_tree, Field_Zp>;
+using Point = std::vector<double>;
+using Points_off_reader = Gudhi::Points_off_reader<Point>;
+
+void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, double& epsilon,
+                     int& dim_max, int& p, Filtration_value& min_persistence);
+
+int main(int argc, char* argv[]) {
+  std::string off_file_points;
+  std::string filediag;
+  double epsilon;
+  int dim_max;
+  int p;
+  Filtration_value min_persistence;
+
+  program_options(argc, argv, off_file_points, filediag, epsilon, dim_max, p, min_persistence);
+
+  Points_off_reader off_reader(off_file_points);
+  Sparse_rips sparse_rips(off_reader.get_point_cloud(), Gudhi::Euclidean_distance(), epsilon);
+
+  // Construct the Rips complex in a Simplex Tree
+  Simplex_tree simplex_tree;
+
+  sparse_rips.create_complex(simplex_tree, dim_max);
+  std::cout << "The complex contains " << simplex_tree.num_simplices() << " simplices \n";
+  std::cout << "   and has dimension " << simplex_tree.dimension() << " \n";
+
+  // Sort the simplices in the order of the filtration
+  simplex_tree.initialize_filtration();
+
+  // Compute the persistence diagram of the complex
+  Persistent_cohomology pcoh(simplex_tree);
+  // initializes the coefficient field for homology
+  pcoh.init_coefficients(p);
+
+  pcoh.compute_persistent_cohomology(min_persistence);
+
+  // Output the diagram in filediag
+  if (filediag.empty()) {
+    pcoh.output_diagram();
+  } else {
+    std::ofstream out(filediag);
+    pcoh.output_diagram(out);
+    out.close();
+  }
+
+  return 0;
+}
+
+void program_options(int argc, char* argv[], std::string& off_file_points, std::string& filediag, double& epsilon,
+                     int& dim_max, int& p, Filtration_value& min_persistence) {
+  namespace po = boost::program_options;
+  po::options_description hidden("Hidden options");
+  hidden.add_options()("input-file", po::value<std::string>(&off_file_points),
+                       "Name of an OFF file containing a point set.\n");
+
+  po::options_description visible("Allowed options", 100);
+  visible.add_options()("help,h", "produce help message")(
+      "output-file,o", po::value<std::string>(&filediag)->default_value(std::string()),
+      "Name of file in which the persistence diagram is written. Default print in std::cout")(
+      "approximation,e", po::value<double>(&epsilon)->default_value(.5),
+      "Epsilon, where the sparse Rips complex is a (1+epsilon)-approximation of the Rips complex.")(
+      "cpx-dimension,d", po::value<int>(&dim_max)->default_value(1),
+      "Maximal dimension of the Rips complex we want to compute.")(
+      "field-charac,p", po::value<int>(&p)->default_value(11),
+      "Characteristic p of the coefficient field Z/pZ for computing homology.")(
+      "min-persistence,m", po::value<Filtration_value>(&min_persistence),
+      "Minimal lifetime of homology feature to be recorded. Default is 0. Enter a negative value to see zero length "
+      "intervals");
+
+  po::positional_options_description pos;
+  pos.add("input-file", 1);
+
+  po::options_description all;
+  all.add(visible).add(hidden);
+
+  po::variables_map vm;
+  po::store(po::command_line_parser(argc, argv).options(all).positional(pos).run(), vm);
+  po::notify(vm);
+
+  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 << "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";
+    std::cout << "b and d are respectively the birth and death of the feature and \n";
+    std::cout << "p is the characteristic of the field Z/pZ used for homology coefficients." << std::endl << std::endl;
+
+    std::cout << "Usage: " << argv[0] << " [options] input-file" << std::endl << std::endl;
+    std::cout << visible << std::endl;
+    std::abort();
+  }
+}
diff --git a/src/cython/doc/cubical_complex_user.rst b/src/cython/doc/cubical_complex_user.rst
index 34598f02..dd82ad93 100644
--- a/src/cython/doc/cubical_complex_user.rst
+++ b/src/cython/doc/cubical_complex_user.rst
@@ -152,6 +152,6 @@ End user programs are available in cython/example/ folder.
 Bibliography
 ============
 
-.. bibliography:: ../../bibliography.bib
+.. bibliography:: ../../biblio/bibliography.bib
    :filter: docnames
    :style: unsrt