Merge branch 'master' into simplex_tree_insert_duplicated_vertices_fix_vincent

38 files changed, 1605 insertions, 118 deletions

diff --git a/src/Alpha_complex/utilities/CMakeLists.txt b/src/Alpha_complex/utilities/CMakeLists.txt
index b12c9690..e76edc5f 100644
--- a/src/Alpha_complex/utilities/CMakeLists.txt
+++ b/src/Alpha_complex/utilities/CMakeLists.txt
@@ -23,7 +23,7 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
 
   add_test(NAME Alpha_complex_utilities_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
           "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off"
-          "-p" "2" "-m" "0.45" "-o" "alpha.pers")
+          "-p" "2" "-m" "0.45" "-o" "safe.pers")
 
   add_test(NAME Alpha_complex_utilities_exact_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
           "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off"
@@ -31,13 +31,13 @@ if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
 
   add_test(NAME Alpha_complex_utilities_safe_alpha_complex_3d COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
           "${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off"
-          "-p" "2" "-m" "0.45" "-o" "safe.pers" "-s")
+          "-p" "2" "-m" "0.45" "-o" "fast.pers" "-f")
 
   if (DIFF_PATH)
     add_test(Alpha_complex_utilities_diff_alpha_complex_3d ${DIFF_PATH}
-            "exact.pers" "alpha.pers")
+            "exact.pers" "safe.pers")
     add_test(Alpha_complex_utilities_diff_alpha_complex_3d ${DIFF_PATH}
-            "safe.pers" "alpha.pers")
+            "fast.pers" "safe.pers")
   endif()
 
   add_test(NAME Alpha_complex_utilities_periodic_alpha_complex_3d_persistence COMMAND $<TARGET_FILE:alpha_complex_3d_persistence>
diff --git a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
index 19e608ad..09c84eb3 100644
--- a/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
+++ b/src/Alpha_complex/utilities/alpha_complex_3d_persistence.cpp
@@ -94,11 +94,11 @@ int main(int argc, char **argv) {
   int coeff_field_characteristic = 0;
   Filtration_value min_persistence = 0.;
   bool exact_version = false;
-  bool safe_version = false;
+  bool fast_version = false;
   bool weighted_version = false;
   bool periodic_version = false;
 
-  program_options(argc, argv, off_file_points, exact_version, safe_version, weight_file, cuboid_file, output_file_diag,
+  program_options(argc, argv, off_file_points, exact_version, fast_version, weight_file, cuboid_file, output_file_diag,
                   alpha_square_max_value, coeff_field_characteristic, min_persistence);
 
   std::vector<double> weights;
@@ -120,16 +120,16 @@ int main(int argc, char **argv) {
     periodic_version = true;
   }
 
-  Gudhi::alpha_complex::complexity complexity = Gudhi::alpha_complex::complexity::FAST;
+  Gudhi::alpha_complex::complexity complexity = Gudhi::alpha_complex::complexity::SAFE;
   if (exact_version) {
-    if (safe_version) {
-      std::cerr << "You cannot set the exact and the safe version." << std::endl;
+    if (fast_version) {
+      std::cerr << "You cannot set the exact and the fast version." << std::endl;
       exit(-1);
     }
     complexity = Gudhi::alpha_complex::complexity::EXACT;
   }
-  if (safe_version) {
-    complexity = Gudhi::alpha_complex::complexity::SAFE;
+  if (fast_version) {
+    complexity = Gudhi::alpha_complex::complexity::FAST;
   }
 
   Simplex_tree simplex_tree;
@@ -258,7 +258,7 @@ int main(int argc, char **argv) {
   return 0;
 }
 
-void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &safe,
+void program_options(int argc, char *argv[], std::string &off_file_points, bool &exact, bool &fast,
                      std::string &weight_file, std::string &cuboid_file, std::string &output_file_diag,
                      Filtration_value &alpha_square_max_value, int &coeff_field_characteristic,
                      Filtration_value &min_persistence) {
@@ -270,9 +270,9 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool
   po::options_description visible("Allowed options", 100);
   visible.add_options()("help,h", "produce help message")(
       "exact,e", po::bool_switch(&exact),
-      "To activate exact version of Alpha complex 3d (default is false, not available if safe is set)")(
-      "safe,s", po::bool_switch(&safe),
-      "To activate safe version of Alpha complex 3d (default is false, not available if exact is set)")(
+      "To activate exact version of Alpha complex 3d (default is false, not available if fast is set)")(
+      "fast,f", po::bool_switch(&fast),
+      "To activate fast version of Alpha complex 3d (default is false, not available if exact is set)")(
       "weight-file,w", po::value<std::string>(&weight_file)->default_value(std::string()),
       "Name of file containing a point weights. Format is one weight per line:\n  W1\n  ...\n  Wn ")(
       "cuboid-file,c", po::value<std::string>(&cuboid_file),
@@ -303,7 +303,7 @@ void program_options(int argc, char *argv[], std::string &off_file_points, bool
     std::cout << std::endl;
     std::cout << "Compute the persistent homology with coefficient field Z/pZ \n";
     std::cout << "of a 3D Alpha complex defined on a set of input points.\n";
-    std::cout << "3D Alpha complex can be exact or safe, weighted and/or periodic\n\n";
+    std::cout << "3D Alpha complex can be safe (by default) exact or fast, weighted and/or periodic\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";
diff --git a/src/Alpha_complex/utilities/alphacomplex.md b/src/Alpha_complex/utilities/alphacomplex.md
index 98f56802..50a39d32 100644
--- a/src/Alpha_complex/utilities/alphacomplex.md
+++ b/src/Alpha_complex/utilities/alphacomplex.md
@@ -107,6 +107,7 @@ It must be in the format described
 points (one value per line).

 * `-e [ --exact ]` for the exact computation version (not compatible with

 weight and periodic version).

+* `-f [ --fast ]` for the fast computation version.

 

 **Example**

 

diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt
index b40d506a..40fdcf2b 100644
--- a/src/CMakeLists.txt
+++ b/src/CMakeLists.txt
@@ -35,6 +35,7 @@ add_gudhi_module(Skeleton_blocker)
 add_gudhi_module(Spatial_searching)
 add_gudhi_module(Subsampling)
 add_gudhi_module(Tangential_complex)
+add_gudhi_module(Toplex_map)
 add_gudhi_module(Witness_complex)
 add_gudhi_module(Nerve_GIC)
 
diff --git a/src/Doxyfile.in b/src/Doxyfile.in
index 858a9299..1c293d1c 100644
--- a/src/Doxyfile.in
+++ b/src/Doxyfile.in
@@ -785,7 +785,8 @@ EXCLUDE                = data/                   \
 	                     GudhUI/                 \
                          cmake/                  \
                          src/cython/             \
-                         include/gudhi_patches/
+                         include/gudhi_patches/  \
+                         README.md
 
 # The EXCLUDE_SYMLINKS tag can be used to select whether or not files or
 # directories that are symbolic links (a Unix file system feature) are excluded
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..97d66fbd 100644
--- a/src/Rips_complex/doc/Intro_rips_complex.h
+++ b/src/Rips_complex/doc/Intro_rips_complex.h
@@ -92,21 +92,27 @@ 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.
  *
  * 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 &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
  * 
@@ -119,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
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<typename DistanceMatrix>
diff --git a/src/Rips_complex/include/gudhi/Sparse_rips_complex.h b/src/Rips_complex/include/gudhi/Sparse_rips_complex.h
index 00da148f..8df6e387 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.
  */
@@ -68,32 +70,36 @@ 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 <typename RandomAccessPointRange, typename Distance>
-  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<Filtration_value>::infinity(), Filtration_value maxi=std::numeric_limits<Filtration_value>::infinity())
+      : epsilon_(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);
+    compute_sparse_graph(dist_fun, epsilon, mini, maxi);
   }
 
   /** \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
+   * \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 <typename DistanceMatrix>
-  Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon)
+  Sparse_rips_complex(const DistanceMatrix& distance_matrix, double epsilon, Filtration_value mini=-std::numeric_limits<Filtration_value>::infinity(), Filtration_value maxi=std::numeric_limits<Filtration_value>::infinity())
       : Sparse_rips_complex(boost::irange<Vertex_handle>(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<j) ? distance_matrix[j][i] : distance_matrix[i][j]; },
+                            epsilon, mini, maxi) {}
 
   /** \brief Fills the simplicial complex with the sparse Rips graph and
    * expands it with all the cliques, stopping at a given maximal dimension.
@@ -111,7 +117,26 @@ 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;
+    }
+    const int n = boost::size(params);
+    std::vector<Filtration_value> lambda(n);
+    // lambda[original_order]=params[sorted_order]
+    for(int i=0;i<n;++i)
+      lambda[sorted_points[i]] = params[i];
+    double cst = epsilon_ * (1 - epsilon_) / 2;
+    auto block = [cst,&complex,&lambda](typename SimplicialComplexForRips::Simplex_handle sh){
+      auto filt = complex.filtration(sh);
+      auto mini = filt * cst;
+      for(auto v : complex.simplex_vertex_range(sh)){
+        if(lambda[v] < mini)
+          return true; // v died before this simplex could be born
+      }
+      return false;
+    };
+    complex.expansion_with_blockers(dim_max, block);
   }
 
  private:
@@ -127,9 +152,11 @@ class Sparse_rips_complex {
   };
 
   // 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) {
+  template <typename Distance>
+  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);
+    double cst = epsilon * (1 - epsilon) / 2;
     graph_.~Graph();
     new (&graph_) Graph(n);
     // for(auto v : vertices(g)) // doesn't work :-(
@@ -143,29 +170,43 @@ 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;
 
         // 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
+        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;
+        }
 
-        add_edge(pi, pj, alpha, graph_);
+        if (alpha <= maxi)
+          add_edge(pi, pj, alpha, graph_);
       }
+    }
   }
 
   Graph graph_;
+  double epsilon_;
+  // Because of the arbitrary split between constructor and create_complex
+  // sorted_points[sorted_order]=original_order
+  std::vector<Vertex_handle> sorted_points;
+  // params[sorted_order]=distance to previous points
+  std::vector<Filtration_value> params;
 };
 
 }  // namespace rips_complex
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<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),
+      "cpx-dimension,d", po::value<int>(&dim_max)->default_value(std::numeric_limits<int>::max()),
       "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.")(
@@ -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";
diff --git a/src/Tangential_complex/doc/Intro_tangential_complex.h b/src/Tangential_complex/doc/Intro_tangential_complex.h
index f4fc8ac7..501f4a8b 100644
--- a/src/Tangential_complex/doc/Intro_tangential_complex.h
+++ b/src/Tangential_complex/doc/Intro_tangential_complex.h
@@ -35,9 +35,11 @@ namespace tangential_complex {
 
 \section tangentialdefinition Definition
 
-A Tangential Delaunay complex is a <a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a>
+A Tangential Delaunay complex is a
+<a target="_blank" href="https://en.wikipedia.org/wiki/Simplicial_complex">simplicial complex</a>
 designed to reconstruct a \f$k\f$-dimensional smooth manifold embedded in \f$d\f$-dimensional Euclidean space. 
-The input is a point sample coming from an unknown manifold, which means that the points lie close to a structure of "small" intrinsic dimension.
+The input is a point sample coming from an unknown manifold, which means that the points lie close to a structure of
+"small" intrinsic dimension.
 The running time depends only linearly on the extrinsic dimension \f$ d \f$
 and exponentially on the intrinsic dimension \f$ k \f$.
 
@@ -46,17 +48,19 @@ An extensive description of the Tangential complex can be found in \cite tangent
 \subsection whatisthetc What is a Tangential Complex?
 
 Let us start with the description of the Tangential complex of a simple example, with \f$ k=1 \f$ and \f$ d=2 \f$.
-The input data is 4 points \f$ P \f$ located on a curve embedded in 2D.
+The point set \f$ \mathscr P \f$ is located on a closed curve embedded in 2D.
+Only 4 points will be displayed (more are required for PCA) to simplify the figures.
 \image html "tc_example_01.png" "The input"
-For each point \f$ p \f$, estimate its tangent subspace \f$ T_p \f$ (e.g. using PCA).
+For each point \f$ P \f$, estimate its tangent subspace \f$ T_P \f$ using PCA.
 \image html "tc_example_02.png" "The estimated normals"
-Let us add the Voronoi diagram of the points in orange. For each point \f$ p \f$, construct its star in the Delaunay triangulation of \f$ P \f$ restricted to \f$ T_p \f$.
+Let us add the Voronoi diagram of the points in orange. For each point \f$ P \f$, construct its star in the Delaunay
+triangulation of \f$ \mathscr P \f$ restricted to \f$ T_P \f$.
 \image html "tc_example_03.png" "The Voronoi diagram"
 The Tangential Delaunay complex is the union of those stars.
 
 In practice, neither the ambient Voronoi diagram nor the ambient Delaunay triangulation is computed.
-Instead, local \f$ k \f$-dimensional regular triangulations are computed with a limited number of points as we only need the star of each point.
-More details can be found in \cite tangentialcomplex2014.
+Instead, local \f$ k \f$-dimensional regular triangulations are computed with a limited number of points as we only
+need the star of each point. More details can be found in \cite tangentialcomplex2014.
 
 \subsection inconsistencies Inconsistencies
 
@@ -65,7 +69,7 @@ An inconsistency occurs when a simplex is not in the star of all its vertices.
 
 Let us take the same example.
 \image html "tc_example_07_before.png" "Before"
-Let us slightly move the tangent subspace \f$ T_q \f$
+Let us slightly move the tangent subspace \f$ T_Q \f$
 \image html "tc_example_07_after.png" "After"
 Now, the star of \f$ Q \f$ contains \f$ QP \f$, but the star of \f$ P \f$ does not contain \f$ QP \f$. We have an inconsistency.
 \image html "tc_example_08.png" "After"
diff --git a/src/Tangential_complex/example/example_basic.cpp b/src/Tangential_complex/example/example_basic.cpp
index 4f2b859e..ab35edf0 100644
--- a/src/Tangential_complex/example/example_basic.cpp
+++ b/src/Tangential_complex/example/example_basic.cpp
@@ -1,5 +1,7 @@
 #include <gudhi/Tangential_complex.h>
 #include <gudhi/sparsify_point_set.h>
+//#include <gudhi/Fake_simplex_tree.h>
+
 
 #include <CGAL/Epick_d.h>
 #include <CGAL/Random.h>
@@ -20,7 +22,7 @@ CGAL::Parallel_tag> TC;
 int main(void) {
   const int INTRINSIC_DIM = 2;
   const int AMBIENT_DIM = 3;
-  const int NUM_POINTS = 1000;
+  const int NUM_POINTS = 100;
 
   Kernel k;
 
@@ -37,6 +39,7 @@ int main(void) {
 
   // Export the TC into a Simplex_tree
   Gudhi::Simplex_tree<> stree;
+  //Gudhi::Fake_simplex_tree stree;
   tc.create_complex(stree);
 
   // Display stats about inconsistencies
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<CGAL::Dynamic_dimension_tag> Kernel;
+  typedef Kernel::FT FT;
+  typedef Kernel::Point_d Point;
+  typedef Kernel::Vector_d Vector;
+  typedef tc::Tangential_complex<Kernel, CGAL::Dynamic_dimension_tag,CGAL::Parallel_tag> 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<Point> generator(AMBIENT_DIM, 3.);
+  std::vector<Point> 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/Toplex_map/benchmark/CMakeLists.txt b/src/Toplex_map/benchmark/CMakeLists.txt
new file mode 100644
index 00000000..2d58a156
--- /dev/null
+++ b/src/Toplex_map/benchmark/CMakeLists.txt
@@ -0,0 +1,3 @@
+project(Toplex_map_benchmark)
+
+add_executable(Toplex_map_benchmark benchmark_tm.cpp)
diff --git a/src/Toplex_map/benchmark/benchmark_tm.cpp b/src/Toplex_map/benchmark/benchmark_tm.cpp
new file mode 100644
index 00000000..eedb442b
--- /dev/null
+++ b/src/Toplex_map/benchmark/benchmark_tm.cpp
@@ -0,0 +1,137 @@
+/*    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:       François Godi, Vincent Rouvreau
+ *
+ *    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 <iostream>
+#include <random>
+#include <chrono>
+
+#include <gudhi/Simplex_tree.h>
+#include <gudhi/Toplex_map.h>
+#include <gudhi/Lazy_toplex_map.h>
+
+using namespace Gudhi;
+
+typedef Toplex_map::Simplex Simplex;
+typedef Toplex_map::Vertex Vertex;
+typedef std::pair<Simplex_tree<>::Simplex_handle, bool> typePairSimplexBool;
+
+class ST_wrapper {
+ public:
+  void insert_simplex(const Simplex& tau) {
+    /*std::cout << "insert_simplex - " << simplexTree.num_simplices() << " - ";
+    for (auto v : tau)
+      std::cout << v << ", ";
+    std::cout << std::endl;
+    */
+    simplexTree.insert_simplex_and_subfaces(tau);
+  }
+
+  bool membership(const Simplex& tau) { return simplexTree.find(tau) != simplexTree.null_simplex(); }
+
+  Vertex contraction(const Vertex x, const Vertex y) {
+    // TODO (VR): edge contraction is not yet available for Simplex_tree
+    return y;
+  }
+
+  std::size_t num_maximal_simplices() { return simplexTree.num_simplices(); }
+
+ private:
+  Simplex_tree<> simplexTree;
+  void erase_max(const Simplex& sigma) {
+    if (membership(sigma)) simplexTree.remove_maximal_simplex(simplexTree.find(sigma));
+  }
+};
+
+int n = 300;
+
+int nb_insert_simplex1 = 3000;
+int nb_membership1 = 4000;
+int nb_contraction = 300;
+int nb_insert_simplex2 = 3000;
+int nb_membership2 = 400000;
+
+Simplex random_simplex(int n, std::size_t d) {
+  std::random_device rd;
+  std::mt19937 gen(rd());
+  std::uniform_int_distribution<std::size_t> dis(1, n);
+  Simplex s;
+  while (s.size() < d) s.insert(dis(gen));
+  return s;
+}
+
+std::vector<Simplex> r_vector_simplices(int n, int max_d, int m) {
+  std::random_device rd;
+  std::mt19937 gen(rd());
+  std::uniform_int_distribution<std::size_t> dis(1, max_d);
+  std::vector<Simplex> v;
+  for (int i = 0; i < m; i++) v.push_back(random_simplex(n, dis(gen)));
+  return v;
+}
+
+template <typename complex_type>
+void chrono(int n, int d) {
+  complex_type K;
+  std::vector<Simplex> simplices_insert_simplex1 = r_vector_simplices(n, d, nb_insert_simplex1);
+  std::vector<Simplex> simplices_membership1 = r_vector_simplices(n, d, nb_membership1);
+  std::vector<Simplex> simplices_insert_simplex2 = r_vector_simplices(n - 2 * nb_contraction, d, nb_insert_simplex2);
+  std::vector<Simplex> simplices_membership2 = r_vector_simplices(n - 2 * nb_contraction, d, nb_membership2);
+  std::chrono::time_point<std::chrono::system_clock> start, end;
+
+  for (const Simplex& s : simplices_insert_simplex1) K.insert_simplex(s);
+
+  for (const Simplex& s : simplices_membership1) K.membership(s);
+
+  start = std::chrono::system_clock::now();
+  for (int i = 1; i <= nb_contraction; i++) K.contraction(n - 2 * i, n - 2 * i - 1);
+  end = std::chrono::system_clock::now();
+  auto c3 = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+
+  start = std::chrono::system_clock::now();
+  for (const Simplex& s : simplices_insert_simplex2) K.insert_simplex(s);
+  end = std::chrono::system_clock::now();
+  auto c1 = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+
+  start = std::chrono::system_clock::now();
+  for (const Simplex& s : simplices_membership2) K.membership(s);
+  end = std::chrono::system_clock::now();
+  auto c2 = std::chrono::duration_cast<std::chrono::milliseconds>(end - start).count();
+
+  if (c3 > 0)
+    std::cout << c1 << "\t \t" << c2 << "\t \t" << c3 << "\t \t" << K.num_maximal_simplices() << std::endl;
+  else
+    std::cout << c1 << "\t \t" << c2 << "\t \tN/A\t \t" << K.num_maximal_simplices() << std::endl;
+}
+
+int main() {
+  for (int d = 5; d <= 40; d += 5) {
+    std::cout << "d=" << d << " \t  Insertions \t   Membership \t Contractions \t        Size" << std::endl;
+    std::cout << "T Map \t \t";
+    chrono<Toplex_map>(n, d);
+    std::cout << "Lazy \t \t";
+    chrono<Lazy_toplex_map>(n, d);
+    if (d <= 15) {
+      std::cout << "ST \t \t";
+      chrono<ST_wrapper>(n, d);
+    }
+    std::cout << std::endl;
+  }
+}
diff --git a/src/Toplex_map/doc/Intro_Toplex_map.h b/src/Toplex_map/doc/Intro_Toplex_map.h
new file mode 100644
index 00000000..a925dc2b
--- /dev/null
+++ b/src/Toplex_map/doc/Intro_Toplex_map.h
@@ -0,0 +1,61 @@
+/*    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:       François Godi, Vincent Rouvreau
+ *
+ *    Copyright (C) 2017  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 DOC_TOPLEX_MAP_INTRO_TOPLEX_MAP_H_
+#define DOC_TOPLEX_MAP_INTRO_TOPLEX_MAP_H_
+
+// needs namespace for Doxygen to link on classes
+namespace Gudhi {
+
+/**  \defgroup toplex_map Toplex Map
+ * 
+ * \author    Fran&ccedil;ois Godi
+ * @{
+ * 
+ * \section toplexmapdefinition Definition
+ *
+ * A Toplex_map is a data structure to represent and store a simplicial complex. A "toplex" is the contraction of
+ * "top-simplex", also known as a maximal simplex. The plural form of "toplex" will be called "toplices".
+ *
+ * Let's consider a simplicial complex, denote by \f$d\f$ its dimension and by \f$k\f$ its number of maximal simplices.
+ * Furthermore, denote by \f$\Gamma_0\f$ the maximal number of toplices, i.e. maximal simplices, that contain a same
+ * vertex.
+ *
+ * The goal of the Toplex Map is both to represent the complex in optimal O(kd) space and to provide fast standard
+ * operations such as : insertion, removal, contraction of an edge, collapses and membership of a simplex. The time
+ * needed for these operation is linear or quadratic in \f$\Gamma_0\f$ and \f$d\f$.
+ *
+ * Toplex map is composed firstly of a raw storage of toplices and secondly of a map which associate any vertex to a
+ * set of pointers toward all toplices containing this vertex. The data structure is described in
+ * \cite boissonnat_et_al:LIPIcs:2015:5098 (aka. Simplex Array List or SAL).
+ *
+ * \image html map.png
+ *
+ * The performances are a lot better than the `Simplex_tree` as soon you use maximal simplices and not simplices
+ * (cf. \cite DBLP:journals/corr/BoissonnatS16 ).
+ *
+ */
+/** @} */  // end defgroup toplex_map
+
+}  // namespace Gudhi
+
+#endif  // DOC_TOPLEX_MAP_INTRO_TOPLEX_MAP_H_
diff --git a/src/Toplex_map/doc/map.png b/src/Toplex_map/doc/map.png
new file mode 100644
index 00000000..0f9cde2b
--- /dev/null
+++ b/src/Toplex_map/doc/map.png
diff --git a/src/Toplex_map/example/CMakeLists.txt b/src/Toplex_map/example/CMakeLists.txt
new file mode 100644
index 00000000..1d54fe87
--- /dev/null
+++ b/src/Toplex_map/example/CMakeLists.txt
@@ -0,0 +1,4 @@
+project(Toplex_map_examples)
+
+add_executable(Toplex_map_example_simple_toplex_map simple_toplex_map.cpp)
+add_test(NAME Toplex_map_example_simple_toplex_map COMMAND $<TARGET_FILE:Toplex_map_example_simple_toplex_map>)
diff --git a/src/Toplex_map/example/simple_toplex_map.cpp b/src/Toplex_map/example/simple_toplex_map.cpp
new file mode 100644
index 00000000..e1c12ed6
--- /dev/null
+++ b/src/Toplex_map/example/simple_toplex_map.cpp
@@ -0,0 +1,165 @@
+/*    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):       Vincent Rouvreau
+ *
+ *    Copyright (C) 2018
+ *
+ *    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/Toplex_map.h>
+
+#include <iostream>
+#include <utility>  // for pair
+#include <vector>
+#include <cassert>
+
+int main(int argc, char* const argv[]) {
+  using Simplex = Gudhi::Toplex_map::Simplex;
+  Simplex sigma1 = {1, 2, 3};
+  Simplex sigma2 = {2, 3, 4, 5};
+
+  Gudhi::Toplex_map tm;
+  tm.insert_simplex(sigma1);
+  tm.insert_simplex(sigma2);
+
+  /* Simplex is:   */
+  /*    2   4      */
+  /*    o---o      */
+  /*   /X\5/       */
+  /*  o---o        */
+  /*  1   3        */
+
+  std::cout << "num max simplices = " << tm.num_maximal_simplices() << " - num vertices = " << tm.num_vertices()
+            << std::endl;
+
+  // Browse maximal cofaces
+  Simplex sigma3 = {2, 3};
+  std::cout << "Maximal cofaces of {2, 3} are :" << std::endl;
+  for (auto simplex_ptr : tm.maximal_cofaces(sigma3, 2)) {
+    for (auto v : *simplex_ptr) {
+      std::cout << v << ", ";
+    }
+    std::cout << std::endl;
+  }
+
+  // Browse maximal simplices
+  std::cout << "Maximal simplices are :" << std::endl;
+  for (auto simplex_ptr : tm.maximal_simplices()) {
+    for (auto v : *simplex_ptr) {
+      std::cout << v << ", ";
+    }
+    std::cout << std::endl;
+  }
+
+  Simplex sigma4 = {1, 3};
+  assert(tm.membership(sigma4));
+
+  Gudhi::Toplex_map::Vertex v = tm.contraction(1, 3);
+  std::cout << "After contraction(1, 3) - " << v << std::endl;
+  /* Simplex is:   */
+  /*    2   4      */
+  /*    o---o      */
+  /*     \5/       */
+  /*      o        */
+  /*      3        */
+  std::cout << "num max simplices = " << tm.num_maximal_simplices() << " - num vertices = " << tm.num_vertices()
+            << std::endl;
+
+  // Browse maximal simplices
+  std::cout << "Maximal simplices are :" << std::endl;
+  for (auto simplex_ptr : tm.maximal_simplices()) {
+    for (auto v : *simplex_ptr) {
+      std::cout << v << ", ";
+    }
+    std::cout << std::endl;
+  }
+
+  Simplex sigma5 = {3, 4};
+  assert(tm.membership(sigma5));
+
+  v = tm.contraction(3, 4);
+  std::cout << "After contraction(3, 4) - " << v << std::endl;
+  /* Simplex is:   */
+  /*    2   4      */
+  /*    o---o      */
+  /*     \X/       */
+  /*      o        */
+  /*      5        */
+  std::cout << "num max simplices = " << tm.num_maximal_simplices() << " - num vertices = " << tm.num_vertices()
+            << std::endl;
+
+  // Browse maximal simplices
+  std::cout << "Maximal simplices are :" << std::endl;
+  for (auto simplex_ptr : tm.maximal_simplices()) {
+    for (auto v : *simplex_ptr) {
+      std::cout << v << ", ";
+    }
+    std::cout << std::endl;
+  }
+
+  tm.insert_simplex(sigma1);
+  tm.insert_simplex(sigma2);
+  /* Simplex is:   */
+  /*    2   4      */
+  /*    o---o      */
+  /*   /X\5/       */
+  /*  o---o        */
+  /*  1   3        */
+  tm.remove_simplex(sigma1);
+
+  std::cout << "After remove_simplex(1, 2, 3)" << std::endl;
+  /* Simplex is:   */
+  /*    2   4      */
+  /*    o---o      */
+  /*   / \5/       */
+  /*  o---o        */
+  /*  1   3        */
+  std::cout << "num max simplices = " << tm.num_maximal_simplices() << " - num vertices = " << tm.num_vertices()
+            << std::endl;
+
+  // Browse maximal simplices
+  std::cout << "Maximal simplices are :" << std::endl;
+  for (auto simplex_ptr : tm.maximal_simplices()) {
+    for (auto v : *simplex_ptr) {
+      std::cout << v << ", ";
+    }
+    std::cout << std::endl;
+  }
+
+  tm.remove_vertex(1);
+
+  std::cout << "After remove_vertex(1)" << std::endl;
+  /* Simplex is:   */
+  /*    2   4      */
+  /*    o---o      */
+  /*     \5/       */
+  /*      o        */
+  /*      3        */
+  std::cout << "num max simplices = " << tm.num_maximal_simplices() << " - num vertices = " << tm.num_vertices()
+            << std::endl;
+
+  // Browse maximal simplices
+  std::cout << "Maximal simplices are :" << std::endl;
+  for (auto simplex_ptr : tm.maximal_simplices()) {
+    for (auto v : *simplex_ptr) {
+      std::cout << v << ", ";
+    }
+    std::cout << std::endl;
+  }
+
+  return 0;
+}
diff --git a/src/Toplex_map/include/gudhi/Lazy_toplex_map.h b/src/Toplex_map/include/gudhi/Lazy_toplex_map.h
new file mode 100644
index 00000000..b0b3706e
--- /dev/null
+++ b/src/Toplex_map/include/gudhi/Lazy_toplex_map.h
@@ -0,0 +1,271 @@
+/*    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:       François Godi, Vincent Rouvreau
+ *
+ *    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 LAZY_TOPLEX_MAP_H
+#define LAZY_TOPLEX_MAP_H
+
+#include <gudhi/Toplex_map.h>
+#include <boost/heap/fibonacci_heap.hpp>
+
+namespace Gudhi {
+
+/**
+ * \brief Lazy toplex map data structure for representing unfiltered simplicial complexes.
+ *
+ * \details A Toplex_map is an unordered map from vertices to maximal simplices (aka. toplices).
+ * The lazy version is not always up to date as it requires clean operation in order to be.
+ *
+ * \ingroup toplex_map */
+class Lazy_toplex_map {
+ public:
+  /** Vertex is the type of vertices. */
+  using Vertex = Toplex_map::Vertex;
+
+  /** Simplex is the type of simplices. */
+  using Simplex = Toplex_map::Simplex;
+
+  /** The type of the pointers to maximal simplices. */
+  using Simplex_ptr = Toplex_map::Simplex_ptr;
+
+  /** The type of the sets of Simplex_ptr. */
+  using Simplex_ptr_set = Toplex_map::Simplex_ptr_set;
+
+  /** Adds the given simplex to the complex.
+   * The simplex must not be in the complex already, and it must not contain one of the current toplices. */
+  template <typename Input_vertex_range>
+  void insert_independent_simplex(const Input_vertex_range &vertex_range);
+
+  /** \brief Adds the given simplex to the complex.
+   * Nothing happens if the simplex is already in the complex (i.e. it is a face of one of the toplices). */
+  template <typename Input_vertex_range>
+  bool insert_simplex(const Input_vertex_range &vertex_range);
+
+  /** \brief Removes the given simplex and its cofaces from the complex.
+   * Its faces are kept inside. */
+  template <typename Input_vertex_range>
+  void remove_simplex(const Input_vertex_range &vertex_range);
+
+  /** Does a simplex belong to the complex ? */
+  template <typename Input_vertex_range>
+  bool membership(const Input_vertex_range &vertex_range);
+
+  /** Do all the facets of a simplex belong to the complex ? */
+  template <typename Input_vertex_range>
+  bool all_facets_inside(const Input_vertex_range &vertex_range);
+
+  /** Contracts one edge in the complex.
+   * The edge has to verify the link condition if you want to preserve topology.
+   * Returns the remaining vertex. */
+  Vertex contraction(const Vertex x, const Vertex y);
+
+  /** \brief Number of maximal simplices. */
+  std::size_t num_maximal_simplices() const { return size; }
+
+  /** \brief Number of vertices. */
+  std::size_t num_vertices() const { return t0.size(); }
+
+ private:
+  template <typename Input_vertex_range>
+  void erase_max(const Input_vertex_range &vertex_range);
+  template <typename Input_vertex_range>
+  Vertex best_index(const Input_vertex_range &vertex_range);
+  void clean(const Vertex v);
+
+  std::unordered_map<Vertex, std::size_t> gamma0_lbounds;
+
+  std::unordered_map<Vertex, Simplex_ptr_set> t0;
+  bool empty_toplex;  // Is the empty simplex a toplex ?
+
+  typedef boost::heap::fibonacci_heap<std::pair<std::size_t, Vertex>> PriorityQueue;
+  PriorityQueue cleaning_priority;
+  std::unordered_map<Vertex, PriorityQueue::handle_type> cp_handles;
+
+  std::size_t get_gamma0_lbound(const Vertex v) const;
+
+  std::size_t size_lbound = 0;
+  std::size_t size = 0;
+
+  const double ALPHA = 4;  // time
+  const double BETTA = 8;  // memory
+};
+
+template <typename Input_vertex_range>
+void Lazy_toplex_map::insert_independent_simplex(const Input_vertex_range &vertex_range) {
+  for (const Vertex &v : vertex_range)
+    if (!gamma0_lbounds.count(v))
+      gamma0_lbounds.emplace(v, 1);
+    else
+      gamma0_lbounds[v]++;
+  size_lbound++;
+  insert_simplex(vertex_range);
+}
+
+template <typename Input_vertex_range>
+bool Lazy_toplex_map::insert_simplex(const Input_vertex_range &vertex_range) {
+  Simplex sigma(vertex_range.begin(), vertex_range.end());
+  // Check empty face management
+  empty_toplex = (sigma.size() == 0);
+  Simplex_ptr sptr = std::make_shared<Simplex>(sigma);
+  bool inserted = false;
+  for (const Vertex &v : sigma) {
+    if (!t0.count(v)) {
+      t0.emplace(v, Simplex_ptr_set());
+      auto v_handle = cleaning_priority.push(std::make_pair(0, v));
+      cp_handles.emplace(v, v_handle);
+    }
+    inserted = t0.at(v).emplace(sptr).second;
+    cleaning_priority.update(cp_handles.at(v), std::make_pair(t0.at(v).size() - get_gamma0_lbound(v), v));
+  }
+  if (inserted) size++;
+  if (size > (size_lbound + 1) * BETTA) clean(cleaning_priority.top().second);
+  return inserted;
+}
+
+template <typename Input_vertex_range>
+void Lazy_toplex_map::remove_simplex(const Input_vertex_range &vertex_range) {
+  if (vertex_range.begin() == vertex_range.end()) {
+    t0.clear();
+    gamma0_lbounds.clear();
+    cleaning_priority.clear();
+    size_lbound = 0;
+    size = 0;
+    empty_toplex = false;
+  } else {
+    const Vertex &v = best_index(vertex_range);
+    // Copy constructor needed because the set is modified
+    if (t0.count(v))
+      for (const Simplex_ptr &sptr : Simplex_ptr_set(t0.at(v)))
+        if (included(vertex_range, *sptr)) {
+          erase_max(*sptr);
+          for (const Simplex &f : facets(vertex_range)) insert_independent_simplex(f);
+        }
+  }
+}
+
+template <typename Input_vertex_range>
+bool Lazy_toplex_map::membership(const Input_vertex_range &vertex_range) {
+  if (t0.size() == 0 && !empty_toplex) return false;            // empty complex
+  if (vertex_range.begin() == vertex_range.end()) return true;  // empty query simplex
+  Vertex v = best_index(vertex_range);
+  if (!t0.count(v)) return false;
+  for (const Simplex_ptr &sptr : t0.at(v))
+    if (included(vertex_range, *sptr)) return true;
+  return false;
+}
+
+template <typename Input_vertex_range>
+bool Lazy_toplex_map::all_facets_inside(const Input_vertex_range &vertex_range) {
+  Simplex sigma(vertex_range.begin(), vertex_range.end());
+  Vertex v = best_index(sigma);
+  if (!t0.count(v)) return false;
+  Simplex f = sigma;
+  f.erase(v);
+  if (!membership(f)) return false;
+  std::unordered_set<Vertex> facets_inside;
+  for (const Simplex_ptr &sptr : t0.at(v))
+    for (const Vertex &w : sigma) {
+      f = sigma;
+      f.erase(w);
+      if (included(f, *sptr)) facets_inside.insert(w);
+    }
+  return facets_inside.size() == sigma.size() - 1;
+}
+
+/* Returns the remaining vertex */
+Lazy_toplex_map::Vertex Lazy_toplex_map::contraction(const Vertex x, const Vertex y) {
+  if (!t0.count(x)) return y;
+  if (!t0.count(y)) return x;
+  Vertex k, d;
+  if (t0.at(x).size() > t0.at(y).size())
+    k = x, d = y;
+  else
+    k = y, d = x;
+  // Copy constructor needed because the set is modified
+  for (const Simplex_ptr &sptr : Simplex_ptr_set(t0.at(d))) {
+    Simplex sigma(*sptr);
+    erase_max(sigma);
+    sigma.erase(d);
+    sigma.insert(k);
+    insert_simplex(sigma);
+  }
+  t0.erase(d);
+  return k;
+}
+
+/* No facets insert_simplexed */
+template <typename Input_vertex_range>
+inline void Lazy_toplex_map::erase_max(const Input_vertex_range &vertex_range) {
+  Simplex sigma(vertex_range.begin(), vertex_range.end());
+  empty_toplex = false;
+  Simplex_ptr sptr = std::make_shared<Simplex>(sigma);
+  bool erased = false;
+  for (const Vertex &v : sigma) {
+    erased = t0.at(v).erase(sptr) > 0;
+    if (t0.at(v).size() == 0) t0.erase(v);
+  }
+  if (erased) size--;
+}
+
+template <typename Input_vertex_range>
+Lazy_toplex_map::Vertex Lazy_toplex_map::best_index(const Input_vertex_range &vertex_range) {
+  Simplex tau(vertex_range.begin(), vertex_range.end());
+  std::size_t min = std::numeric_limits<size_t>::max();
+  Vertex arg_min = -1;
+  for (const Vertex &v : tau)
+    if (!t0.count(v))
+      return v;
+    else if (t0.at(v).size() < min)
+      min = t0.at(v).size(), arg_min = v;
+  if (min > ALPHA * get_gamma0_lbound(arg_min)) clean(arg_min);
+  return arg_min;
+}
+
+std::size_t Lazy_toplex_map::get_gamma0_lbound(const Vertex v) const {
+  return gamma0_lbounds.count(v) ? gamma0_lbounds.at(v) : 0;
+}
+
+void Lazy_toplex_map::clean(const Vertex v) {
+  Toplex_map toplices;
+  std::unordered_map<int, std::vector<Simplex>> dsorted_simplices;
+  std::size_t max_dim = 0;
+  for (const Simplex_ptr &sptr : Simplex_ptr_set(t0.at(v))) {
+    if (sptr->size() > max_dim) {
+      for (std::size_t d = max_dim + 1; d <= sptr->size(); d++) dsorted_simplices.emplace(d, std::vector<Simplex>());
+      max_dim = sptr->size();
+    }
+    dsorted_simplices[sptr->size()].emplace_back(*sptr);
+    erase_max(*sptr);
+  }
+  for (std::size_t d = max_dim; d >= 1; d--)
+    for (const Simplex &s : dsorted_simplices.at(d))
+      if (!toplices.membership(s)) toplices.insert_independent_simplex(s);
+  Simplex sv;
+  sv.insert(v);
+  auto clean_cofaces = toplices.maximal_cofaces(sv);
+  size_lbound = size_lbound - get_gamma0_lbound(v) + clean_cofaces.size();
+  gamma0_lbounds[v] = clean_cofaces.size();
+  for (const Simplex_ptr &sptr : clean_cofaces) insert_simplex(*sptr);
+}
+
+}  // namespace Gudhi
+
+#endif /* LAZY_TOPLEX_MAP_H */
diff --git a/src/Toplex_map/include/gudhi/Toplex_map.h b/src/Toplex_map/include/gudhi/Toplex_map.h
new file mode 100644
index 00000000..4dc2331f
--- /dev/null
+++ b/src/Toplex_map/include/gudhi/Toplex_map.h
@@ -0,0 +1,338 @@
+/*    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:       François Godi, Vincent Rouvreau
+ *
+ *    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 TOPLEX_MAP_H
+#define TOPLEX_MAP_H
+
+#include <vector>
+#include <set>
+#include <unordered_set>
+#include <unordered_map>
+#include <memory>
+#include <limits>
+
+namespace Gudhi {
+
+/**
+ * \brief Toplex map data structure for representing unfiltered simplicial complexes.
+ *
+ * \details A Toplex_map is an unordered map from vertices to maximal simplices (aka. toplices).
+ *
+ * \ingroup toplex_map */
+class Toplex_map {
+ public:
+  /** Vertex is the type of vertices. */
+  using Vertex = std::size_t;
+
+  /** Simplex is the type of simplices. */
+  using Simplex = std::set<Vertex>;
+
+  /** The type of the pointers to maximal simplices. */
+  using Simplex_ptr = std::shared_ptr<Toplex_map::Simplex>;
+
+  struct Sptr_hash {
+    std::size_t operator()(const Toplex_map::Simplex_ptr& s) const;
+  };
+
+  struct Sptr_equal {
+    std::size_t operator()(const Toplex_map::Simplex_ptr& a, const Toplex_map::Simplex_ptr& b) const;
+  };
+
+  /** The type of the sets of Toplex_map::Simplex_ptr. */
+  using Simplex_ptr_set = std::unordered_set<Toplex_map::Simplex_ptr, Sptr_hash, Sptr_equal>;
+
+  /** \brief Adds the given simplex to the complex.
+   * Nothing happens if the simplex is already in the complex (i.e. it is a face of one of the toplices). */
+  template <typename Input_vertex_range>
+  void insert_simplex(const Input_vertex_range& vertex_range);
+
+  /** \brief Removes the given simplex and its cofaces from the complex.
+   * Its faces are kept inside. */
+  template <typename Input_vertex_range>
+  void remove_simplex(const Input_vertex_range& vertex_range);
+
+  /** Does a simplex belong to the complex ? */
+  template <typename Input_vertex_range>
+  bool membership(const Input_vertex_range& vertex_range) const;
+
+  /** Does a simplex is a toplex ? */
+  template <typename Input_vertex_range>
+  bool maximality(const Input_vertex_range& vertex_range) const;
+
+  /** Gives a set of pointers to the maximal cofaces of a simplex.
+   * Gives all the toplices if given the empty simplex.
+   * Gives not more than max_number maximal cofaces if max_number is strictly positive. */
+  template <typename Input_vertex_range>
+  Toplex_map::Simplex_ptr_set maximal_cofaces(const Input_vertex_range& vertex_range,
+                                              const std::size_t max_number = 0) const;
+
+  /** Gives a set of pointers to the maximal simplices.
+   * Gives not more than max_number maximal cofaces if max_number is strictly positive. */
+  Toplex_map::Simplex_ptr_set maximal_simplices(const std::size_t max_number = 0) const {
+    return maximal_cofaces(Simplex(), max_number);
+  }
+
+  /** Contracts one edge in the complex.
+   * The edge has to verify the link condition if you want to preserve topology.
+   * Returns the remaining vertex. */
+  Vertex contraction(const Vertex x, const Vertex y);
+
+  /** Remove the vertex and all its cofaces from the complex. */
+  void remove_vertex(const Vertex x);
+
+  /** \brief Number of maximal simplices. */
+  std::size_t num_maximal_simplices() const { return maximal_simplices().size(); }
+
+  /** \brief Number of vertices. */
+  std::size_t num_vertices() const { return t0.size(); }
+
+  std::set<Vertex> unitary_collapse(const Vertex k, const Vertex d);
+
+  /** Adds the given simplex to the complex.
+   * The simplex must not be in the complex already, and it must not contain one of the current toplices. */
+  template <typename Input_vertex_range>
+  void insert_independent_simplex(const Input_vertex_range& vertex_range);
+
+ protected:
+  /** \internal Gives an index in order to look for a simplex quickly. */
+  template <typename Input_vertex_range>
+  Vertex best_index(const Input_vertex_range& vertex_range) const;
+
+  /** \internal The map from vertices to toplices */
+  std::unordered_map<Vertex, Toplex_map::Simplex_ptr_set> t0;
+
+  const Vertex VERTEX_UPPER_BOUND = std::numeric_limits<Vertex>::max();
+
+  /** \internal Removes a toplex without adding facets after. */
+  void erase_maximal(const Toplex_map::Simplex_ptr& sptr);
+};
+
+// Pointers are also used as key in the hash sets.
+template <typename Input_vertex_range>
+Toplex_map::Simplex_ptr get_key(const Input_vertex_range& vertex_range);
+
+// Is the first simplex a face of the second ?
+template <typename Input_vertex_range1, typename Input_vertex_range2>
+bool included(const Input_vertex_range1& vertex_range1, const Input_vertex_range2& vertex_range2);
+
+// All the facets of the given simplex.
+template <typename Input_vertex_range>
+std::vector<Toplex_map::Simplex> facets(const Input_vertex_range& vertex_range);
+
+template <typename Input_vertex_range>
+void Toplex_map::insert_simplex(const Input_vertex_range& vertex_range) {
+  if (membership(vertex_range)) return;
+  bool replace_facets = true;
+  for (const Toplex_map::Simplex& facet : facets(vertex_range))
+    if (!maximality(facet)) {
+      replace_facets = false;
+      break;
+    }
+  if (replace_facets)
+    for (const Toplex_map::Simplex& facet : facets(vertex_range)) erase_maximal(get_key(facet));
+  else
+    for (const Vertex& v : vertex_range)
+      if (t0.count(v))
+        for (const Toplex_map::Simplex_ptr& fptr : Simplex_ptr_set(t0.at(v)))
+          // Copy constructor needed because the set is modified
+          if (included(*fptr, vertex_range)) erase_maximal(fptr);
+  // We erase all the maximal faces of the simplex
+  insert_independent_simplex(vertex_range);
+}
+
+template <typename Input_vertex_range>
+void Toplex_map::remove_simplex(const Input_vertex_range& vertex_range) {
+  if (vertex_range.begin() == vertex_range.end()) t0.clear();
+  // Removal of the empty simplex means cleaning everything
+  else {
+    const Vertex& v = best_index(vertex_range);
+    if (t0.count(v))
+      for (const Toplex_map::Simplex_ptr& sptr : Simplex_ptr_set(t0.at(v)))
+        // Copy constructor needed because the set is modified
+        if (included(vertex_range, *sptr)) {
+          erase_maximal(sptr);
+          for (const Toplex_map::Simplex& f : facets(vertex_range))
+            if (!membership(f)) insert_independent_simplex(f);
+          // We add the facets which are new maximal simplices
+        }
+  }
+}
+
+template <typename Input_vertex_range>
+bool Toplex_map::membership(const Input_vertex_range& vertex_range) const {
+  if (t0.size() == 0) return false;
+  const Vertex& v = best_index(vertex_range);
+  if (!t0.count(v)) return false;
+  if (maximality(vertex_range)) return true;
+  for (const Toplex_map::Simplex_ptr& sptr : t0.at(v))
+    if (included(vertex_range, *sptr)) return true;
+  return false;
+}
+
+template <typename Input_vertex_range>
+bool Toplex_map::maximality(const Input_vertex_range& vertex_range) const {
+  const Vertex& v = best_index(vertex_range);
+  if (!t0.count(v)) return false;
+  return t0.at(v).count(get_key(vertex_range));
+}
+
+template <typename Input_vertex_range>
+Toplex_map::Simplex_ptr_set Toplex_map::maximal_cofaces(const Input_vertex_range& vertex_range,
+                                                        const std::size_t max_number) const {
+  Simplex_ptr_set cofaces;
+  if (maximality(vertex_range))
+    cofaces.emplace(get_key(vertex_range));
+  else if (vertex_range.begin() == vertex_range.end())
+    for (const auto& kv : t0)
+      for (const Toplex_map::Simplex_ptr& sptr : kv.second) {
+        // kv.second is a Simplex_ptr_set
+        cofaces.emplace(sptr);
+        if (cofaces.size() == max_number) return cofaces;
+      }
+  else {
+    const Vertex& v = best_index(vertex_range);
+    if (t0.count(v))
+      for (const Toplex_map::Simplex_ptr& sptr : t0.at(v))
+        if (included(vertex_range, *sptr)) {
+          cofaces.emplace(sptr);
+          if (cofaces.size() == max_number) return cofaces;
+        }
+  }
+  return cofaces;
+}
+
+Toplex_map::Vertex Toplex_map::contraction(const Toplex_map::Vertex x, const Toplex_map::Vertex y) {
+  if (!t0.count(x)) return y;
+  if (!t0.count(y)) return x;
+  int k, d;
+  if (t0.at(x).size() > t0.at(y).size())
+    k = x, d = y;
+  else
+    k = y, d = x;
+  for (const Toplex_map::Simplex_ptr& sptr : Simplex_ptr_set(t0.at(d))) {
+    // Copy constructor needed because the set is modified
+    Simplex sigma(*sptr);
+    erase_maximal(sptr);
+    sigma.erase(d);
+    sigma.insert(k);
+    insert_simplex(sigma);
+  }
+  return k;
+}
+
+std::set<Toplex_map::Vertex> Toplex_map::unitary_collapse(const Toplex_map::Vertex k, const Toplex_map::Vertex d) {
+  Toplex_map::Simplex r;
+  for (const Toplex_map::Simplex_ptr& sptr : Simplex_ptr_set(t0.at(d))) {
+    // Copy constructor needed because the set is modified
+    Simplex sigma(*sptr);
+    erase_maximal(sptr);
+    sigma.erase(d);
+    for (const Vertex v : sigma) r.insert(v);
+    sigma.insert(k);
+    insert_simplex(sigma);
+  }
+  return r;
+}
+
+template <typename Input_vertex_range>
+void Toplex_map::insert_independent_simplex(const Input_vertex_range& vertex_range) {
+  auto key = get_key(vertex_range);
+  for (const Vertex& v : vertex_range) {
+    if (!t0.count(v)) t0.emplace(v, Simplex_ptr_set());
+    t0.at(v).emplace(key);
+  }
+}
+
+void Toplex_map::remove_vertex(const Toplex_map::Vertex x) {
+  for (const Toplex_map::Simplex_ptr& sptr : Simplex_ptr_set(t0.at(x))) {
+    Simplex sigma(*sptr);
+    erase_maximal(sptr);
+    sigma.erase(x);
+    insert_simplex(sigma);
+  }
+}
+
+inline void Toplex_map::erase_maximal(const Toplex_map::Simplex_ptr& sptr) {
+  Simplex sigma(*sptr);
+  if (sptr->size() == 0) sigma.insert(VERTEX_UPPER_BOUND);
+  for (const Vertex& v : sigma) {
+    t0.at(v).erase(sptr);
+    if (t0.at(v).size() == 0) t0.erase(v);
+  }
+}
+
+template <typename Input_vertex_range>
+Toplex_map::Vertex Toplex_map::best_index(const Input_vertex_range& vertex_range) const {
+  std::size_t min = std::numeric_limits<size_t>::max();
+  Vertex arg_min = VERTEX_UPPER_BOUND;
+  for (const Vertex& v : vertex_range)
+    if (!t0.count(v))
+      return v;
+    else if (t0.at(v).size() < min)
+      min = t0.at(v).size(), arg_min = v;
+  return arg_min;
+}
+
+std::size_t Toplex_map::Sptr_equal::operator()(const Toplex_map::Simplex_ptr& s1,
+                                               const Toplex_map::Simplex_ptr& s2) const {
+  if (s1->size() != s2->size()) return false;
+  return included(*s1, *s2);
+  // inclusion tests equality for same size simplices
+}
+
+std::size_t Toplex_map::Sptr_hash::operator()(const Toplex_map::Simplex_ptr& s) const {
+  std::hash<double> h_f;
+  // double hash works better than int hash
+  size_t h = 0;
+  for (const Vertex& v : *s) h += h_f(static_cast<double>(v));
+  return h;
+}
+
+template <typename Input_vertex_range>
+Toplex_map::Simplex_ptr get_key(const Input_vertex_range& vertex_range) {
+  Toplex_map::Simplex s(vertex_range.begin(), vertex_range.end());
+  return std::make_shared<Toplex_map::Simplex>(s);
+}
+
+template <typename Input_vertex_range1, typename Input_vertex_range2>
+bool included(const Input_vertex_range1& vertex_range1, const Input_vertex_range2& vertex_range2) {
+  Toplex_map::Simplex s2(vertex_range2.begin(), vertex_range2.end());
+  for (const Toplex_map::Vertex& v : vertex_range1)
+    if (!s2.count(v)) return false;
+  return true;
+}
+
+template <typename Input_vertex_range>
+std::vector<Toplex_map::Simplex> facets(const Input_vertex_range& vertex_range) {
+  std::vector<Toplex_map::Simplex> facets;
+  Toplex_map::Simplex f(vertex_range.begin(), vertex_range.end());
+  for (const Toplex_map::Vertex& v : vertex_range) {
+    f.erase(v);
+    facets.emplace_back(f);
+    f.insert(v);
+  }
+  return facets;
+}
+
+}  // namespace Gudhi
+
+#endif /* TOPLEX_MAP_H */
diff --git a/src/Toplex_map/test/CMakeLists.txt b/src/Toplex_map/test/CMakeLists.txt
new file mode 100644
index 00000000..59517db5
--- /dev/null
+++ b/src/Toplex_map/test/CMakeLists.txt
@@ -0,0 +1,11 @@
+project(Toplex_map_tests)
+
+include(GUDHI_test_coverage)
+
+add_executable( Toplex_map_unit_test toplex_map_unit_test.cpp )
+target_link_libraries(Toplex_map_unit_test ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+gudhi_add_coverage_test(Toplex_map_unit_test)
+
+add_executable( Lazy_toplex_map_unit_test lazy_toplex_map_unit_test.cpp )
+target_link_libraries(Lazy_toplex_map_unit_test ${Boost_UNIT_TEST_FRAMEWORK_LIBRARY})
+gudhi_add_coverage_test(Lazy_toplex_map_unit_test)
diff --git a/src/Toplex_map/test/lazy_toplex_map_unit_test.cpp b/src/Toplex_map/test/lazy_toplex_map_unit_test.cpp
new file mode 100644
index 00000000..a050cc92
--- /dev/null
+++ b/src/Toplex_map/test/lazy_toplex_map_unit_test.cpp
@@ -0,0 +1,170 @@
+/*    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:       François Godi, Vincent Rouvreau
+ *
+ *    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 <iostream>
+#include <vector>
+#include <gudhi/Lazy_toplex_map.h>
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "lazy toplex map"
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_CASE(toplex_map) {
+  using Vertex = Gudhi::Lazy_toplex_map::Vertex;
+
+  Gudhi::Lazy_toplex_map tm;
+  std::cout << "insert_simplex {1, 2, 3, 4}" << std::endl;
+  std::vector<Vertex> sigma1 = {1, 2, 3, 4};
+  tm.insert_simplex(sigma1);
+  std::cout << "insert_simplex {5, 2, 3, 6}" << std::endl;
+  std::vector<Vertex> sigma2 = {5, 2, 3, 6};
+  tm.insert_simplex(sigma2);
+  std::cout << "insert_simplex {5}" << std::endl;
+  std::vector<Vertex> sigma3 = {5};
+  tm.insert_simplex(sigma3);
+  std::cout << "insert_simplex {4, 5, 3}" << std::endl;
+  std::vector<Vertex> sigma6 = {4, 5, 3};
+  tm.insert_simplex(sigma6);
+  std::cout << "insert_simplex {4, 5, 9}" << std::endl;
+  std::vector<Vertex> sigma7 = {4, 5, 9};
+  tm.insert_simplex(sigma7);
+
+  std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl;
+  BOOST_CHECK(tm.num_maximal_simplices() == 5);
+
+  std::vector<Vertex> sigma4 = {5, 2, 3};
+  std::vector<Vertex> sigma5 = {5, 2, 7};
+  BOOST_CHECK(tm.membership(sigma4));
+  BOOST_CHECK(!tm.membership(sigma5));
+  std::cout << "insert_simplex {5, 2, 7}" << std::endl;
+  tm.insert_simplex(sigma5);
+
+  std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl;
+  BOOST_CHECK(tm.num_maximal_simplices() == 6);
+
+  BOOST_CHECK(tm.membership(sigma5));
+
+  std::cout << "contraction(4,5)" << std::endl;
+  auto r = tm.contraction(4, 5);
+  std::cout << "r=" << r << std::endl;
+  BOOST_CHECK(r == 5);
+
+  std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl;
+  BOOST_CHECK(tm.num_maximal_simplices() == 6);
+
+  std::vector<Vertex> sigma8 = {1, 2, 3};
+  std::vector<Vertex> sigma9 = {2, 7};
+
+  sigma8.emplace_back(r);
+  sigma9.emplace_back(r);
+  BOOST_CHECK(!tm.membership(sigma6));
+  BOOST_CHECK(tm.membership(sigma8));
+  BOOST_CHECK(tm.membership(sigma9));
+
+  std::cout << "remove_simplex({2, 7, r = 5})" << std::endl;
+  tm.remove_simplex(sigma9);
+  BOOST_CHECK(!tm.membership(sigma9));
+
+  std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl;
+  BOOST_CHECK(tm.num_maximal_simplices() == 8);
+
+  // {2, 7, 5} is removed, but verify its edges are still there
+  std::vector<Vertex> edge = {2, 7};
+  BOOST_CHECK(tm.membership(edge));
+  edge = {2, 5};
+  BOOST_CHECK(tm.membership(edge));
+  edge = {7, 5};
+  BOOST_CHECK(tm.membership(edge));
+}
+
+BOOST_AUTO_TEST_CASE(toplex_map_empty_toplex) {
+  using Vertex = Gudhi::Lazy_toplex_map::Vertex;
+
+  Gudhi::Lazy_toplex_map tm;
+  std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl;
+  BOOST_CHECK(tm.num_maximal_simplices() == 0);
+  std::cout << "num_vertices = " << tm.num_vertices() << std::endl;
+  BOOST_CHECK(tm.num_vertices() == 0);
+
+  std::cout << "Check an empty simplex is a member." << std::endl;
+  std::vector<Vertex> empty_sigma = {};
+  BOOST_CHECK(tm.membership(empty_sigma));
+
+  std::cout << "Check the edge 2,7 is not a member." << std::endl;
+  std::vector<Vertex> edge = {2, 7};
+  BOOST_CHECK(!tm.membership(edge));
+
+  std::cout << "Insert an empty simplex." << std::endl;
+  tm.insert_simplex(empty_sigma);
+
+  std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl;
+  BOOST_CHECK(tm.num_maximal_simplices() == 0);
+  std::cout << "num_vertices = " << tm.num_vertices() << std::endl;
+  BOOST_CHECK(tm.num_vertices() == 0);
+
+  std::cout << "Check an empty simplex is a member." << std::endl;
+  BOOST_CHECK(tm.membership(empty_sigma));
+  std::cout << "Check the edge 2,7 is not a member." << std::endl;
+  BOOST_CHECK(!tm.membership(edge));
+
+  std::cout << "Insert edge 2,7." << std::endl;
+  tm.insert_simplex(edge);
+
+  std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl;
+  BOOST_CHECK(tm.num_maximal_simplices() == 1);
+  std::cout << "num_vertices = " << tm.num_vertices() << std::endl;
+  BOOST_CHECK(tm.num_vertices() == 2);
+
+  std::cout << "Check an empty simplex is a member." << std::endl;
+  BOOST_CHECK(tm.membership(empty_sigma));
+  std::cout << "Check the edge 2,7 is a member." << std::endl;
+  BOOST_CHECK(tm.membership(edge));
+
+  std::cout << "contraction(2,7)" << std::endl;
+  auto r = tm.contraction(2, 7);
+  std::cout << "r=" << r << std::endl;
+  BOOST_CHECK(r == 7);
+
+  std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl;
+  BOOST_CHECK(tm.num_maximal_simplices() == 1);
+  std::cout << "num_vertices = " << tm.num_vertices() << std::endl;
+  BOOST_CHECK(tm.num_vertices() == 1);
+
+  std::cout << "Check an empty simplex is a member." << std::endl;
+  BOOST_CHECK(tm.membership(empty_sigma));
+  std::cout << "Check the edge 2,7 is not a member." << std::endl;
+  BOOST_CHECK(!tm.membership(edge));
+
+  std::cout << "Remove the vertex 7." << std::endl;
+  std::vector<Vertex> vertex = {7};
+  tm.remove_simplex(vertex);
+
+  std::cout << "num_maximal_simplices = " << tm.num_maximal_simplices() << std::endl;
+  BOOST_CHECK(tm.num_maximal_simplices() == 0);
+  std::cout << "num_vertices = " << tm.num_vertices() << std::endl;
+  BOOST_CHECK(tm.num_vertices() == 0);
+
+  std::cout << "Check an empty simplex is a member." << std::endl;
+  BOOST_CHECK(tm.membership(empty_sigma));
+  std::cout << "Check the edge 2,7 is not a member." << std::endl;
+  BOOST_CHECK(!tm.membership(edge));
+}
diff --git a/src/Toplex_map/test/toplex_map_unit_test.cpp b/src/Toplex_map/test/toplex_map_unit_test.cpp
new file mode 100644
index 00000000..2bd27936
--- /dev/null
+++ b/src/Toplex_map/test/toplex_map_unit_test.cpp
@@ -0,0 +1,138 @@
+/*    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:       François Godi, Vincent Rouvreau
+ *
+ *    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 <iostream>
+#include <vector>
+#include <gudhi/Toplex_map.h>
+
+#define BOOST_TEST_DYN_LINK
+#define BOOST_TEST_MODULE "toplex map"
+#include <boost/test/unit_test.hpp>
+
+BOOST_AUTO_TEST_CASE(toplex_map) {
+  using Vertex = Gudhi::Toplex_map::Vertex;
+
+  Gudhi::Toplex_map tm;
+  std::cout << "insert_simplex {1, 2, 3, 4}" << std::endl;
+  std::vector<Vertex> sigma1 = {1, 2, 3, 4};
+  tm.insert_simplex(sigma1);
+  std::cout << "insert_simplex {5, 2, 3, 6}" << std::endl;
+  std::vector<Vertex> sigma2 = {5, 2, 3, 6};
+  tm.insert_simplex(sigma2);
+  std::cout << "insert_simplex {5}" << std::endl;
+  std::vector<Vertex> sigma3 = {5};
+  tm.insert_simplex(sigma3);
+  std::cout << "insert_simplex {4, 5, 3}" << std::endl;
+  std::vector<Vertex> sigma6 = {4, 5, 3};
+  tm.insert_simplex(sigma6);
+  std::cout << "insert_simplex {4, 5, 9}" << std::endl;
+  std::vector<Vertex> sigma7 = {4, 5, 9};
+  tm.insert_simplex(sigma7);
+
+  std::cout << "num_maximal_simplices" << tm.num_maximal_simplices() << std::endl;
+  BOOST_CHECK(tm.num_maximal_simplices() == 4);
+  // Browse maximal simplices
+  std::cout << "Maximal simplices are :" << std::endl;
+  for (auto simplex_ptr : tm.maximal_simplices()) {
+    for (auto v : *simplex_ptr) {
+      std::cout << v << ", ";
+    }
+    std::cout << std::endl;
+    BOOST_CHECK(tm.maximality(*simplex_ptr));
+  }
+
+  BOOST_CHECK(tm.maximality(sigma1));
+  BOOST_CHECK(tm.maximality(sigma2));
+  BOOST_CHECK(!tm.maximality(sigma3));
+  BOOST_CHECK(tm.maximality(sigma6));
+  BOOST_CHECK(tm.maximality(sigma7));
+
+  std::vector<Vertex> sigma4 = {5, 2, 3};
+  std::vector<Vertex> sigma5 = {5, 2, 7};
+  BOOST_CHECK(tm.membership(sigma4));
+  BOOST_CHECK(!tm.membership(sigma5));
+  std::cout << "insert_simplex {5, 2, 7}" << std::endl;
+  tm.insert_simplex(sigma5);
+
+  std::cout << "num_maximal_simplices" << tm.num_maximal_simplices() << std::endl;
+  BOOST_CHECK(tm.num_maximal_simplices() == 5);
+  // Browse maximal simplices
+  std::cout << "Maximal simplices are :" << std::endl;
+  for (auto simplex_ptr : tm.maximal_simplices()) {
+    for (auto v : *simplex_ptr) {
+      std::cout << v << ", ";
+    }
+    std::cout << std::endl;
+    BOOST_CHECK(tm.maximality(*simplex_ptr));
+  }
+
+  BOOST_CHECK(tm.membership(sigma5));
+
+  std::cout << "contraction(4,5)" << std::endl;
+  auto r = tm.contraction(4, 5);
+  std::cout << "r=" << r << std::endl;
+  BOOST_CHECK(r == 5);
+
+  std::cout << "num_maximal_simplices" << tm.num_maximal_simplices() << std::endl;
+  BOOST_CHECK(tm.num_maximal_simplices() == 4);
+  // Browse maximal simplices
+  std::cout << "Maximal simplices are :" << std::endl;
+  for (auto simplex_ptr : tm.maximal_simplices()) {
+    for (auto v : *simplex_ptr) {
+      std::cout << v << ", ";
+    }
+    std::cout << std::endl;
+    BOOST_CHECK(tm.maximality(*simplex_ptr));
+  }
+
+  std::vector<Vertex> sigma8 = {1, 2, 3};
+  std::vector<Vertex> sigma9 = {2, 7};
+
+  sigma8.emplace_back(r);
+  sigma9.emplace_back(r);
+  BOOST_CHECK(!tm.membership(sigma6));
+  BOOST_CHECK(tm.membership(sigma8));
+  BOOST_CHECK(tm.membership(sigma9));
+
+  std::cout << "remove_simplex({2, 7, r = 5})" << std::endl;
+  tm.remove_simplex(sigma9);
+  BOOST_CHECK(!tm.membership(sigma9));
+
+  std::cout << "num_maximal_simplices" << tm.num_maximal_simplices() << std::endl;
+  BOOST_CHECK(tm.num_maximal_simplices() == 5);
+  // Browse maximal simplices
+  std::cout << "Maximal simplices are :" << std::endl;
+  for (auto simplex_ptr : tm.maximal_simplices()) {
+    for (auto v : *simplex_ptr) {
+      std::cout << v << ", ";
+    }
+    std::cout << std::endl;
+    BOOST_CHECK(tm.maximality(*simplex_ptr));
+  }
+  // {2, 7, 5} is removed, but verify its edges are still there
+  std::vector<Vertex> edge = {2, 7};
+  BOOST_CHECK(tm.membership(edge));
+  edge = {2, 5};
+  BOOST_CHECK(tm.membership(edge));
+  edge = {7, 5};
+  BOOST_CHECK(tm.membership(edge));
+}
diff --git a/src/cmake/modules/GUDHI_third_party_libraries.cmake b/src/cmake/modules/GUDHI_third_party_libraries.cmake
index 7b0d350d..57ea7d14 100644
--- a/src/cmake/modules/GUDHI_third_party_libraries.cmake
+++ b/src/cmake/modules/GUDHI_third_party_libraries.cmake
@@ -148,6 +148,7 @@ if( PYTHONINTERP_FOUND )
   find_python_module("matplotlib")
   find_python_module("numpy")
   find_python_module("scipy")
+  find_python_module("sphinx")
 endif()
 
 if(NOT GUDHI_CYTHON_PATH)
@@ -157,20 +158,16 @@ endif(NOT GUDHI_CYTHON_PATH)
 option(WITH_GUDHI_CYTHON_RUNTIME_LIBRARY_DIRS "Build with setting runtime_library_dirs. Usefull when setting rpath is not allowed" ON)
 
 if(PYTHONINTERP_FOUND AND CYTHON_FOUND)
-  # Default found version 2
-  if(PYTHON_VERSION_MAJOR EQUAL 2)
+  if(SPHINX_FOUND)
     # Documentation generation is available through sphinx
     find_program( SPHINX_PATH sphinx-build )
-  elseif(PYTHON_VERSION_MAJOR EQUAL 3)
-    # No sphinx-build in Pyton3, just hack it
-    set(SPHINX_PATH "${PYTHON_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/${GUDHI_CYTHON_PATH}/doc/python3-sphinx-build.py")
-  else()
-    message(FATAL_ERROR "ERROR: Try to compile the Cython interface. Python version ${PYTHON_VERSION_STRING} is not valid.")
-  endif(PYTHON_VERSION_MAJOR EQUAL 2)
-  # get PYTHON_SITE_PACKAGES relative path from a python command line
-  execute_process(
-    COMMAND "${PYTHON_EXECUTABLE}" -c "from distutils.sysconfig import get_python_lib; print (get_python_lib(prefix='', plat_specific=True))"
-    OUTPUT_VARIABLE PYTHON_SITE_PACKAGES
-    OUTPUT_STRIP_TRAILING_WHITESPACE)
+
+    if(NOT SPHINX_PATH)
+      if(PYTHON_VERSION_MAJOR EQUAL 3)
+        # In Python3, just hack sphinx-build if it does not exist
+        set(SPHINX_PATH "${PYTHON_EXECUTABLE}" "${CMAKE_CURRENT_SOURCE_DIR}/${GUDHI_CYTHON_PATH}/doc/python3-sphinx-build.py")
+      endif(PYTHON_VERSION_MAJOR EQUAL 3)
+    endif(NOT SPHINX_PATH)
+  endif(SPHINX_FOUND)
 endif(PYTHONINTERP_FOUND AND CYTHON_FOUND)
 
diff --git a/src/cmake/modules/GUDHI_user_version_target.cmake b/src/cmake/modules/GUDHI_user_version_target.cmake
index d43a6fa6..2ed48c48 100644
--- a/src/cmake/modules/GUDHI_user_version_target.cmake
+++ b/src/cmake/modules/GUDHI_user_version_target.cmake
@@ -27,7 +27,7 @@ add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E
 add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E
                    copy ${CMAKE_SOURCE_DIR}/Conventions.txt ${GUDHI_USER_VERSION_DIR}/Conventions.txt)
 add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E
-                   copy ${CMAKE_SOURCE_DIR}/README ${GUDHI_USER_VERSION_DIR}/README)
+                   copy ${CMAKE_SOURCE_DIR}/README.md ${GUDHI_USER_VERSION_DIR}/README.md)
 add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E
                    copy ${CMAKE_SOURCE_DIR}/COPYING ${GUDHI_USER_VERSION_DIR}/COPYING)
 add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E
diff --git a/src/common/doc/main_page.h b/src/common/doc/main_page.h
index b33a95a1..afe9b68c 100644
--- a/src/common/doc/main_page.h
+++ b/src/common/doc/main_page.h
@@ -171,6 +171,23 @@
     </td>
  </tr>
 </table>
+ \subsection ToplexMapDataStructure Toplex Map
+ \image html "map.png" "Toplex map representation"
+<table border="0">
+  <tr>
+    <td width="25%">
+      <b>Author:</b> Fran&ccedil;ois Godi<br>
+      <b>Introduced in:</b> GUDHI 2.1.0<br>
+      <b>Copyright:</b> GPL v3<br>
+    </td>
+    <td width="75%">
+    The Toplex map data structure is composed firstly of a raw storage of toplices (the maximal simplices)
+    and secondly of a map which associate any vertex to a set of pointers toward all toplices
+    containing this vertex.
+    <b>User manual:</b> \ref toplex_map - <b>Reference manual:</b> Gudhi::Toplex_map
+    </td>
+ </tr>
+
  \subsection WitnessComplexDataStructure Witness complex
  \image html "Witness_complex_representation.png" "Witness complex representation"
 <table border="0">
diff --git a/src/cython/CMakeLists.txt b/src/cython/CMakeLists.txt
index 97859c98..480332d7 100644
--- a/src/cython/CMakeLists.txt
+++ b/src/cython/CMakeLists.txt
@@ -215,12 +215,7 @@ if(PYTHONINTERP_FOUND)
     add_custom_target(cython ALL DEPENDS gudhi.so
                       COMMENT "Do not forget to add ${CMAKE_CURRENT_BINARY_DIR}/ to your PYTHONPATH before using examples or tests")
 
-    # For installation purpose
-    # TODO(VR) : files matching pattern mechanism is copying all cython directory
-    install(DIRECTORY "${CMAKE_CURRENT_BINARY_DIR}" DESTINATION "${PYTHON_SITE_PACKAGES}/" FILES_MATCHING
-        PATTERN "*.so"
-        PATTERN "*.dylib"
-        PATTERN "*.pyd")
+    install(CODE "execute_process(COMMAND ${PYTHON_EXECUTABLE} ${CMAKE_CURRENT_BINARY_DIR}/setup.py install)")
 
     # Test examples
     if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.8.1)
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/installation.rst b/src/cython/doc/installation.rst
index 040f6b4a..855dea44 100644
--- a/src/cython/doc/installation.rst
+++ b/src/cython/doc/installation.rst
@@ -7,24 +7,23 @@ Installation
 
 Compiling
 *********
-The library uses c++11 and requires `Boost <https://www.boost.org/>`_ ≥ 1.56.0
-and `CMake <https://www.cmake.org/>`_   ≥ 3.1.
+The library uses c++11 and requires `Boost <https://www.boost.org/>`_ ≥ 1.56.0,
+`CMake <https://www.cmake.org/>`_   ≥ 3.1  to generate makefiles, and
+`Cython <https://www.cython.org/>`_ to compile the GUDHI Python module.
 It is a multi-platform library and compiles on Linux, Mac OSX and Visual
 Studio 2015.
-It also requires cmake to generate makefiles, and cython to compile the
-library.
 
 On `Windows <https://wiki.python.org/moin/WindowsCompilers>`_ , only Python
 3.5 and 3.6 are available because of the required Visual Studio version.
 
-On other systems, if you have several Python/cython installed, the version 2.X
+On other systems, if you have several Python/Cython installed, the version 2.X
 will be used by default, but you can force it by adding
 :code:`-DPython_ADDITIONAL_VERSIONS=3` to the cmake command.
 
-GUDHI Cythonization
-===================
+GUDHI Python module compilation
+===============================
 
-To build the GUDHI cython module, run the following commands in a terminal:
+To build the GUDHI Python module, run the following commands in a terminal:
 
 .. code-block:: bash
 
@@ -32,7 +31,28 @@ To build the GUDHI cython module, run the following commands in a terminal:
     mkdir build
     cd build/
     cmake ..
-    make cython
+    cd cython
+    make
+
+GUDHI Python module installation
+================================
+
+Once the compilation succeeds, one can add the GUDHI Python module path to the
+PYTHONPATH:
+
+.. code-block:: bash
+
+    # For windows, you have to set PYTHONPATH environment variable
+    export PYTHONPATH='$PYTHONPATH:/path-to-gudhi/build/cython'
+
+Or install it definitely in your Python packages folder:
+
+.. code-block:: bash
+
+    cd /path-to-gudhi/build/cython
+    # May require sudo or administrator privileges
+    make install
+
 
 Test suites
 ===========
@@ -45,7 +65,7 @@ following command in a terminal:
     cd /path-to-gudhi/build/cython
     # For windows, you have to set PYTHONPATH environment variable
     export PYTHONPATH='$PYTHONPATH:/path-to-gudhi/build/cython'
-    ctest -R py_test
+    make test
 
 Debugging issues
 ================
@@ -54,7 +74,7 @@ If tests fail, please check your PYTHONPATH and try to :code:`import gudhi`
 and check the errors.
 The problem can come from a third-party library bad link or installation.
 
-If :code:`import gudhi` succeeds, please have a look to debug informations:
+If :code:`import gudhi` succeeds, please have a look to debug information:
 
 .. code-block:: python
 
@@ -105,13 +125,17 @@ A complete configuration would be :
 Documentation
 =============
 
-To build the documentation, `sphinx-doc <http://http://www.sphinx-doc.org>`_ is
-required. Please refer to *conf.py* file to see which
-`sphinx-doc <http://http://www.sphinx-doc.org>`_ modules are required to
-generate the documentation. Run the following commands in a terminal:
+To build the documentation, `sphinx-doc <http://www.sphinx-doc.org>`_ and
+`sphinxcontrib-bibtex <https://sphinxcontrib-bibtex.readthedocs.io>`_ are
+required. As the documentation is auto-tested, `CGAL`_, `Eigen3`_,
+`Matplotlib`_, `NumPy`_ and `SciPy`_ are also mandatory to build the
+documentation.
+
+Run the following commands in a terminal:
 
 .. code-block:: bash
 
+    cd /path-to-gudhi/build/cython
     make sphinx
 
 Optional third-party library
@@ -127,7 +151,7 @@ The :doc:`Alpha complex </alpha_complex_user>`,
 C++ library which provides easy access to efficient and reliable geometric
 algorithms.
 
-Having CGAL,  the Computational Geometry Algorithms Library, version 4.7.0 or
+Having CGAL, the Computational Geometry Algorithms Library, version 4.7.0 or
 higher installed is recommended. The procedure to install this library
 according to your operating system is detailed
 `here <http://doc.cgal.org/latest/Manual/installation.html>`_.
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/doc/tangential_complex_user.rst b/src/cython/doc/tangential_complex_user.rst
index 5ce69e86..ebfe1e29 100644
--- a/src/cython/doc/tangential_complex_user.rst
+++ b/src/cython/doc/tangential_complex_user.rst
@@ -23,8 +23,10 @@ What is a Tangential Complex?
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 Let us start with the description of the Tangential complex of a simple
-example, with :math:`k = 1` and :math:`d = 2`. The input data is 4 points
-:math:`P` located on a curve embedded in 2D.
+example, with :math:`k = 1` and :math:`d = 2`. The point set
+:math:`\mathscr P` is located on a closed curve embedded in 2D.
+Only 4 points will be displayed (more are required for PCA) to simplify the
+figures.
 
 .. figure:: ../../doc/Tangential_complex/tc_example_01.png
     :alt: The input
@@ -32,8 +34,7 @@ example, with :math:`k = 1` and :math:`d = 2`. The input data is 4 points
 
     The input
 
-For each point :math:`p`, estimate its tangent subspace :math:`T_p` (e.g.
-using PCA). 
+For each point :math:`P`, estimate its tangent subspace :math:`T_P` using PCA.
 
 .. figure:: ../../doc/Tangential_complex/tc_example_02.png
     :alt: The estimated normals
@@ -43,8 +44,8 @@ using PCA).
 
 
 Let us add the Voronoi diagram of the points in orange. For each point
-:math:`p`, construct its star in the Delaunay triangulation of :math:`P`
-restricted to :math:`T_p`.
+:math:`P`, construct its star in the Delaunay triangulation of
+:math:`\mathscr P` restricted to :math:`T_P`.
 
 .. figure:: ../../doc/Tangential_complex/tc_example_03.png
     :alt: The Voronoi diagram
@@ -72,7 +73,7 @@ Let us take the same example.
 
     Before
 
-Let us slightly move the tangent subspace :math:`T_q`
+Let us slightly move the tangent subspace :math:`T_Q`
 
 .. figure:: ../../doc/Tangential_complex/tc_example_07_after.png
     :alt: After
@@ -128,6 +129,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 +177,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/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<std::vector<double>>& 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<double>::infinity(), threshold);
   }
   void init_matrix_sparse(const std::vector<std::vector<double>>& matrix, double threshold, double epsilon) {
-    sparse_rips_complex_.emplace(matrix, epsilon);
-    threshold_ = threshold;
+    sparse_rips_complex_.emplace(matrix, epsilon, -std::numeric_limits<double>::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<Rips_complex<Simplex_tree_interface<>::Filtration_value>> rips_complex_;
   boost::optional<Sparse_rips_complex<Simplex_tree_interface<>::Filtration_value>> sparse_rips_complex_;
-  double threshold_;
 };
 
 }  // namespace rips_complex
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<Point_d> 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<double> get_point(unsigned vh) {
     std::vector<double> 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