Merge branch 'master' into toplex_map

20 files changed, 315 insertions, 108 deletions

diff --git a/src/cython/CMakeLists.txt b/src/cython/CMakeLists.txt
index dc2e9278..480332d7 100644
--- a/src/cython/CMakeLists.txt
+++ b/src/cython/CMakeLists.txt
@@ -198,9 +198,9 @@ if(PYTHONINTERP_FOUND)
       set(GUDHI_CYTHON_INCLUDE_DIRS "${GUDHI_CYTHON_INCLUDE_DIRS}'${TBB_INCLUDE_DIRS}', ")
     endif()
 
-    if(UNIX)
+    if(UNIX AND WITH_GUDHI_CYTHON_RUNTIME_LIBRARY_DIRS)
       set( GUDHI_CYTHON_RUNTIME_LIBRARY_DIRS "${GUDHI_CYTHON_LIBRARY_DIRS}")
-    endif(UNIX)
+    endif(UNIX AND WITH_GUDHI_CYTHON_RUNTIME_LIBRARY_DIRS)
 
     # Generate setup.py file to cythonize Gudhi - This file must be named setup.py by convention
     configure_file(setup.py.in "${CMAKE_CURRENT_BINARY_DIR}/setup.py" @ONLY)
@@ -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/rips_complex.pyx b/src/cython/cython/rips_complex.pyx
index 30ca4443..7c83241c 100644
--- a/src/cython/cython/rips_complex.pyx
+++ b/src/cython/cython/rips_complex.pyx
@@ -33,7 +33,11 @@ __license__ = "GPL v3"
 
 cdef extern from "Rips_complex_interface.h" namespace "Gudhi":
     cdef cppclass Rips_complex_interface "Gudhi::rips_complex::Rips_complex_interface":
-        Rips_complex_interface(vector[vector[double]] values, double threshold, bool euclidean)
+        Rips_complex_interface()
+        void init_points(vector[vector[double]] values, double threshold)
+        void init_matrix(vector[vector[double]] values, double threshold)
+        void init_points_sparse(vector[vector[double]] values, double threshold, double sparse)
+        void init_matrix_sparse(vector[vector[double]] values, double threshold, double sparse)
         void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, int dim_max)
 
 # RipsComplex python interface
@@ -44,10 +48,11 @@ cdef class RipsComplex:
     function, or a distance matrix.
     """
 
-    cdef Rips_complex_interface * thisptr
+    cdef Rips_complex_interface thisref
 
     # Fake constructor that does nothing but documenting the constructor
-    def __init__(self, points=None, distance_matrix=None, max_edge_length=float('inf')):
+    def __init__(self, points=None, distance_matrix=None,
+                 max_edge_length=float('inf'), sparse=None):
         """RipsComplex constructor.
 
         :param max_edge_length: Rips value.
@@ -59,29 +64,38 @@ cdef class RipsComplex:
         Or
 
         :param distance_matrix: A distance matrix (full square or lower
-        triangular).
+            triangular).
         :type points: list of list of double
+
+        And in both cases
+
+        :param sparse: If this is not None, it switches to building a sparse
+            Rips and represents the approximation parameter epsilon.
+        :type sparse: float
         """
 
     # The real cython constructor
-    def __cinit__(self, points=None, distance_matrix=None, max_edge_length=float('inf')):
-        if distance_matrix is not None:
-            self.thisptr = new Rips_complex_interface(distance_matrix, max_edge_length, False)
+    def __cinit__(self, points=None, distance_matrix=None,
+                  max_edge_length=float('inf'), sparse=None):
+        if sparse is not None:
+          if distance_matrix is not None:
+              self.thisref.init_matrix_sparse(distance_matrix,
+                                              max_edge_length,
+                                              sparse)
+          else:
+              if points is None:
+                  # Empty Rips construction
+                  points=[]
+              self.thisref.init_points_sparse(points, max_edge_length, sparse)
         else:
-            if points is None:
-                # Empty Rips construction
-                points=[]
-            self.thisptr = new Rips_complex_interface(points, max_edge_length, True)
-
-
-    def __dealloc__(self):
-        if self.thisptr != NULL:
-            del self.thisptr
+          if distance_matrix is not None:
+              self.thisref.init_matrix(distance_matrix, max_edge_length)
+          else:
+              if points is None:
+                  # Empty Rips construction
+                  points=[]
+              self.thisref.init_points(points, max_edge_length)
 
-    def __is_defined(self):
-        """Returns true if RipsComplex pointer is not NULL.
-         """
-        return self.thisptr != NULL
 
     def create_simplex_tree(self, max_dimension=1):
         """
@@ -92,5 +106,5 @@ cdef class RipsComplex:
         :rtype: SimplexTree
         """
         simplex_tree = SimplexTree()
-        self.thisptr.create_simplex_tree(simplex_tree.thisptr, max_dimension)
+        self.thisref.create_simplex_tree(simplex_tree.thisptr, max_dimension)
         return simplex_tree
diff --git a/src/cython/cython/simplex_tree.pyx b/src/cython/cython/simplex_tree.pyx
index 8397d9d9..0ab97f80 100644
--- a/src/cython/cython/simplex_tree.pyx
+++ b/src/cython/cython/simplex_tree.pyx
@@ -44,8 +44,8 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         void set_dimension(int dimension)
         int dimension()
         int upper_bound_dimension()
-        bint find_simplex(vector[int] simplex)
-        bint insert_simplex_and_subfaces(vector[int] simplex,
+        bool find_simplex(vector[int] simplex)
+        bool insert_simplex_and_subfaces(vector[int] simplex,
                                          double filtration)
         vector[pair[vector[int], double]] get_filtration()
         vector[pair[vector[int], double]] get_skeleton(int dimension)
@@ -355,7 +355,7 @@ cdef class SimplexTree:
         :param filtration: Maximum threshold value.
         :type filtration: float.
         :returns: The filtration modification information.
-        :rtype: bint
+        :rtype: bool
 
 
         .. note::
@@ -406,7 +406,7 @@ cdef class SimplexTree:
         value than its faces by increasing the filtration values.
 
         :returns: The filtration modification information.
-        :rtype: bint
+        :rtype: bool
 
 
         .. note::
@@ -432,6 +432,10 @@ cdef class SimplexTree:
             0.0.
             Sets min_persistence to -1.0 to see all values.
         :type min_persistence: float.
+        :param persistence_dim_max: If true, the persistent homology for the
+            maximal dimension in the complex is computed. If false, it is
+            ignored. Default is false.
+        :type persistence_dim_max: bool
         :returns: The persistence of the simplicial complex.
         :rtype:  list of pairs(dimension, pair(birth, death))
         """
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/cubical_complex_user.rst b/src/cython/doc/cubical_complex_user.rst
index 320bd79b..19120360 100644
--- a/src/cython/doc/cubical_complex_user.rst
+++ b/src/cython/doc/cubical_complex_user.rst
@@ -83,9 +83,15 @@ Input Format.
 In the current implantation, filtration is given at the maximal cubes, and it is then extended by the lower star
 filtration to all cubes. There are a number of constructors that can be used to construct cubical complex by users
 who want to use the code directly. They can be found in the :doc:`cubical_complex_ref`.
-Currently one input from a text file is used. It uses a format used already in
+Currently one input from a text file is used. It uses a format inspired from the Perseus software
 `Perseus software <http://www.sas.upenn.edu/~vnanda/perseus/>`_ by Vidit Nanda.
-The file format is described here: :doc:`Perseus <fileformats>`.
+
+.. note::
+    While Perseus assume the filtration of all maximal cubes to be non-negative, over here we do not enforce this and
+    we allow any filtration values. As a consequence one cannot use ``-1``'s to indicate missing cubes. If you have
+    missing cubes in your complex, please set their filtration to :math:`+\infty` (aka. ``inf`` in the file).
+
+The file format is described in details in :ref:`Perseus file format` file format section.
 
 .. testcode::
 
diff --git a/src/cython/doc/examples.rst b/src/cython/doc/examples.rst
index 1f02f8a2..edbc2f72 100644
--- a/src/cython/doc/examples.rst
+++ b/src/cython/doc/examples.rst
@@ -22,6 +22,7 @@ Examples
     * :download:`rips_complex_diagram_persistence_from_off_file_example.py <../example/rips_complex_diagram_persistence_from_off_file_example.py>`
     * :download:`rips_complex_diagram_persistence_from_distance_matrix_file_example.py <../example/rips_complex_diagram_persistence_from_distance_matrix_file_example.py>`
     * :download:`rips_persistence_diagram.py <../example/rips_persistence_diagram.py>`
+    * :download:`sparse_rips_persistence_diagram.py <../example/sparse_rips_persistence_diagram.py>`
     * :download:`random_cubical_complex_persistence_example.py <../example/random_cubical_complex_persistence_example.py>`
     * :download:`coordinate_graph_induced_complex.py <../example/coordinate_graph_induced_complex.py>`
     * :download:`functional_graph_induced_complex.py <../example/functional_graph_induced_complex.py>`
diff --git a/src/cython/doc/fileformats.rst b/src/cython/doc/fileformats.rst
index ff20f26e..e205cc8b 100644
--- a/src/cython/doc/fileformats.rst
+++ b/src/cython/doc/fileformats.rst
@@ -51,10 +51,12 @@ Here is a simple sample file in the 3D case::
     1. 1. 1.
 
 
+.. _Perseus file format:
+
 Perseus
 *******
 
-This file format is the format used by the
+This file format is a format inspired from the
 `Perseus software <http://www.sas.upenn.edu/~vnanda/perseus/>`_ by Vidit Nanda.
 The first line contains a number d begin the dimension of the bitmap (2 in the
 example below). Next d lines are the numbers of top dimensional cubes in each
@@ -88,3 +90,9 @@ Indicate that we have imposed periodic boundary conditions in the direction x,
 but not in the direction y.
 
 Other sample files can be found in the `data/bitmap` folder.
+
+.. note::
+    Unlike in Perseus format the filtration on the maximal cubes can be any
+    double precision number. Consequently one cannot mark the cubes that are
+    not present with ``-1``'s. To do that please set their filtration value to
+    :math:`+\infty` (aka. ``inf`` in the file).
\ No newline at end of file
diff --git a/src/cython/doc/installation.rst b/src/cython/doc/installation.rst
index ef2f7af2..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.48.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/nerve_gic_complex_ref.rst b/src/cython/doc/nerve_gic_complex_ref.rst
index e24e01fc..abde2e8c 100644
--- a/src/cython/doc/nerve_gic_complex_ref.rst
+++ b/src/cython/doc/nerve_gic_complex_ref.rst
@@ -1,3 +1,7 @@
+:orphan:
+
+.. To get rid of WARNING: document isn't included in any toctree
+
 ================================
 Cover complexes reference manual
 ================================
diff --git a/src/cython/doc/nerve_gic_complex_user.rst b/src/cython/doc/nerve_gic_complex_user.rst
index d774827e..44f30e1a 100644
--- a/src/cython/doc/nerve_gic_complex_user.rst
+++ b/src/cython/doc/nerve_gic_complex_user.rst
@@ -1,3 +1,7 @@
+:orphan:
+
+.. To get rid of WARNING: document isn't included in any toctree
+
 Cover complexes user manual
 ===========================
 Definition
diff --git a/src/cython/doc/persistent_cohomology_user.rst b/src/cython/doc/persistent_cohomology_user.rst
index ce7fc685..de83cda1 100644
--- a/src/cython/doc/persistent_cohomology_user.rst
+++ b/src/cython/doc/persistent_cohomology_user.rst
@@ -10,12 +10,14 @@ Definition
 :Author: Clément Maria                 :Introduced in: GUDHI PYTHON 2.0.0     :Copyright: GPL v3
 =====================================  =====================================  =====================================
 
-+---------------------------------------------+----------------------------------------------------------------------+
-|  :doc:`persistent_cohomology_user`          | Please refer to each data structure that contains persistence        |
-|                                             | feature for reference:                                               |
-|                                             |                                                                      |
-|                                             | * :doc:`simplex_tree_ref`                                            |
-+---------------------------------------------+----------------------------------------------------------------------+
++-----------------------------------------------------------------+-----------------------------------------------------------------------+
+|  :doc:`persistent_cohomology_user`                              | Please refer to each data structure that contains persistence         |
+|                                                                 | feature for reference:                                                |
+|                                                                 |                                                                       |
+|                                                                 | * :doc:`simplex_tree_ref`                                             |
+|                                                                 | * :doc:`cubical_complex_ref`                                          |
+|                                                                 | * :doc:`periodic_cubical_complex_ref`                                 |
++-----------------------------------------------------------------+-----------------------------------------------------------------------+
 
 
 Computation of persistent cohomology using the algorithm of :cite:`DBLP:journals/dcg/SilvaMV11` and
diff --git a/src/cython/doc/rips_complex_sum.inc b/src/cython/doc/rips_complex_sum.inc
index 5616bfa9..ea26769a 100644
--- a/src/cython/doc/rips_complex_sum.inc
+++ b/src/cython/doc/rips_complex_sum.inc
@@ -1,6 +1,6 @@
-=================================================================  ===================================  ===================================
-:Author: Clément Maria, Pawel Dlotko, Vincent Rouvreau             :Introduced in: GUDHI 2.0.0          :Copyright: GPL v3
-=================================================================  ===================================  ===================================
+=====================================================================  ===========================  ===================================
+:Author: Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse    :Introduced in: GUDHI 2.0.0          :Copyright: GPL v3
+=====================================================================  ===========================  ===================================
 
 +----------------------------------------------------------------+------------------------------------------------------------------------+
 | .. figure::                                                    | Rips complex is a simplicial complex constructed from a one skeleton   |
diff --git a/src/cython/doc/rips_complex_user.rst b/src/cython/doc/rips_complex_user.rst
index a8c06cf9..e814b4c3 100644
--- a/src/cython/doc/rips_complex_user.rst
+++ b/src/cython/doc/rips_complex_user.rst
@@ -7,27 +7,27 @@ Rips complex user manual
 Definition
 ----------
 
-=======================================================  =====================================  =====================================
-:Authors: Clément Maria, Pawel Dlotko, Vincent Rouvreau  :Introduced in: GUDHI 2.0.0            :Copyright: GPL v3
-=======================================================  =====================================  =====================================
+====================================================================  ================================  ======================
+:Authors: Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse  :Introduced in: GUDHI 2.0.0       :Copyright: GPL v3
+====================================================================  ================================  ======================
 
 +-------------------------------------------+----------------------------------------------------------------------+
 | :doc:`rips_complex_user`                  | :doc:`rips_complex_ref`                                              |
 +-------------------------------------------+----------------------------------------------------------------------+
 
-`Rips complex <https://en.wikipedia.org/wiki/Vietoris%E2%80%93Rips_complex>`_ is a one skeleton graph that allows to
-construct a simplicial complex from it. The input can be a point cloud with a given distance function, or a distance
-matrix.
+The `Rips complex <https://en.wikipedia.org/wiki/Vietoris%E2%80%93Rips_complex>`_ is a simplicial complex that
+generalizes proximity (:math:`\varepsilon`-ball) graphs to higher dimensions. The vertices correspond to the input
+points, and a simplex is present if and only if its diameter is smaller than some parameter α.  Considering all
+parameters α defines a filtered simplicial complex, where the filtration value of a simplex is its diameter.
+The filtration can be restricted to values α smaller than some threshold, to reduce its size.
 
-The filtration value of each edge is computed from a user-given distance function, or directly from the distance
-matrix.
+The input discrete metric space can be provided as a point cloud plus a distance function, or as a distance matrix.
 
-All edges that have a filtration value strictly greater than a given threshold value are not inserted into the complex.
+When creating a simplicial complex from the graph, :doc:`RipsComplex <rips_complex_ref>` first builds the graph and
+inserts it into the data structure. It then expands the simplicial complex (adds the simplices corresponding to cliques)
+when required. The expansion can be stopped at dimension `max_dimension`, by default 1.
 
-When creating a simplicial complex from this one skeleton graph, Rips inserts the one skeleton graph into the data
-structure, and then expands the simplicial complex when required.
-
-Vertex name correspond to the index of the point in the given range (aka. the point cloud).
+A vertex name corresponds to the index of the point in the given range (aka. the point cloud).
 
 .. figure::
     ../../doc/Rips_complex/rips_complex_representation.png
@@ -38,8 +38,27 @@ Vertex name correspond to the index of the point in the given range (aka. the po
 On this example, as edges (4,5), (4,6) and (5,6) are in the complex, simplex (4,5,6) is added with the filtration value
 set with :math:`max(filtration(4,5), filtration(4,6), filtration(5,6))`. And so on for simplex (0,1,2,3).
 
-If the Rips_complex interfaces are not detailed enough for your need, please refer to rips_persistence_step_by_step.cpp
-example, where the graph construction over the Simplex_tree is more detailed.
+If the `RipsComplex` interfaces are not detailed enough for your need, please refer to rips_persistence_step_by_step.cpp
+C++ example, where the graph construction over the Simplex_tree is more detailed.
+
+A Rips complex can easily become huge, even if we limit the length of the edges
+and the dimension of the simplices. One easy trick, before building a Rips
+complex on a point cloud, is to call `sparsify_point_set` which removes points
+that are too close to each other. This does not change its persistence diagram
+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
+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
+construction of a `RipsComplex` object asks it to build a sparse Rips with
+parameter :math:`\varepsilon=0.3`, while the default `sparse=None` builds the
+regular Rips complex.
+
 
 Point cloud
 -----------
@@ -47,7 +66,7 @@ Point cloud
 Example from a point cloud
 ^^^^^^^^^^^^^^^^^^^^^^^^^^
 
-This example builds the one skeleton graph from the given points, and max_edge_length value.
+This example builds the neighborhood graph from the given points, up to max_edge_length.
 Then it creates a :doc:`Simplex_tree <simplex_tree_ref>` with it.
 
 Finally, it is asked to display information about the simplicial complex.
@@ -56,7 +75,7 @@ Finally, it is asked to display information about the simplicial complex.
 
     import gudhi
     rips_complex = gudhi.RipsComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]],
-        max_edge_length=12.0)
+                                     max_edge_length=12.0)
 
     simplex_tree = rips_complex.create_simplex_tree(max_dimension=1)
     result_str = 'Rips complex is of dimension ' + repr(simplex_tree.dimension()) + ' - ' + \
@@ -92,10 +111,20 @@ until dimension 1 - one skeleton graph in other words), the output is:
     [4, 6] -> 9.49
     [3, 6] -> 11.00
 
+Notice that if we use
+
+.. code-block:: python
+
+    rips_complex = gudhi.RipsComplex(points=[[1, 1], [7, 0], [4, 6], [9, 6], [0, 14], [2, 19], [9, 17]],
+                                     max_edge_length=12.0, sparse=2)
+
+asking for a very sparse version (theory only gives some guarantee on the meaning of the output if `sparse<1`),
+2 to 5 edges disappear, depending on the random vertex used to start the sparsification.
+
 Example from OFF file
 ^^^^^^^^^^^^^^^^^^^^^
 
-This example builds the :doc:`Rips_complex <rips_complex_ref>` from the given
+This example builds the :doc:`RipsComplex <rips_complex_ref>` from the given
 points in an OFF file, and max_edge_length value.
 Then it creates a :doc:`Simplex_tree <simplex_tree_ref>` with it.
 
@@ -200,7 +229,7 @@ until dimension 1 - one skeleton graph in other words), the output is:
 Example from csv file
 ^^^^^^^^^^^^^^^^^^^^^
 
-This example builds the :doc:`Rips_complex <rips_complex_ref>` from the given
+This example builds the :doc:`RipsComplex <rips_complex_ref>` from the given
 distance matrix in a csv file, and max_edge_length value.
 Then it creates a :doc:`Simplex_tree <simplex_tree_ref>` with it.
 
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/sparse_rips_persistence_diagram.py b/src/cython/example/sparse_rips_persistence_diagram.py
new file mode 100755
index 00000000..d58c244c
--- /dev/null
+++ b/src/cython/example/sparse_rips_persistence_diagram.py
@@ -0,0 +1,43 @@
+#!/usr/bin/env python
+
+import gudhi
+
+"""This file is part of the Gudhi Library. The Gudhi library
+   (Geometric Understanding in Higher Dimensions) is a generic C++
+   library for computational topology.
+
+   Author(s):       Marc Glisse
+
+   Copyright (C) 2018 Inria
+
+   This program is free software: you can redistribute it and/or modify
+   it under the terms of the GNU General Public License as published by
+   the Free Software Foundation, either version 3 of the License, or
+   (at your option) any later version.
+
+   This program is distributed in the hope that it will be useful,
+   but WITHOUT ANY WARRANTY; without even the implied warranty of
+   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
+   GNU General Public License for more details.
+
+   You should have received a copy of the GNU General Public License
+   along with this program.  If not, see <http://www.gnu.org/licenses/>.
+"""
+
+__author__ = "Marc Glisse"
+__copyright__ = "Copyright (C) 2018 Inria"
+__license__ = "GPL v3"
+
+print("#####################################################################")
+print("Sparse RipsComplex creation from points")
+rips = gudhi.RipsComplex(points=[[0, 0], [0, 0.1], [1, 0], [0, 1], [1, 1]],
+                         max_edge_length=42, sparse=.5)
+
+simplex_tree = rips.create_simplex_tree(max_dimension=2)
+
+
+diag = simplex_tree.persistence(homology_coeff_field=2, min_persistence=0)
+print("diag=", diag)
+
+pplot = gudhi.plot_persistence_diagram(diag)
+pplot.show()
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 8b6c9c35..1a6e2477 100644
--- a/src/cython/include/Rips_complex_interface.h
+++ b/src/cython/include/Rips_complex_interface.h
@@ -25,8 +25,11 @@
 
 #include <gudhi/Simplex_tree.h>
 #include <gudhi/Rips_complex.h>
+#include <gudhi/Sparse_rips_complex.h>
 #include <gudhi/distance_functions.h>
 
+#include <boost/optional.hpp>
+
 #include "Simplex_tree_interface.h"
 
 #include <iostream>
@@ -43,28 +46,40 @@ class Rips_complex_interface {
   using Distance_matrix = std::vector<std::vector<Simplex_tree_interface<>::Filtration_value>>;
 
  public:
-  Rips_complex_interface(const std::vector<std::vector<double>>& values, double threshold, bool euclidean) {
-    if (euclidean) {
-      // Rips construction where values is a vector of points
-      rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(values, threshold,
-                                                                                   Gudhi::Euclidean_distance());
-    } else {
-      // Rips construction where values is a distance matrix
-      rips_complex_ = new Rips_complex<Simplex_tree_interface<>::Filtration_value>(values, threshold);
-    }
+  void init_points(const std::vector<std::vector<double>>& points, double threshold) {
+    rips_complex_.emplace(points, threshold, Gudhi::Euclidean_distance());
+  }
+  void init_matrix(const std::vector<std::vector<double>>& matrix, double threshold) {
+    rips_complex_.emplace(matrix, threshold);
   }
 
-  ~Rips_complex_interface() {
-    delete rips_complex_;
+  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;
+  }
+  void init_matrix_sparse(const std::vector<std::vector<double>>& matrix, double threshold, double epsilon) {
+    sparse_rips_complex_.emplace(matrix, epsilon);
+    threshold_ = threshold;
   }
 
   void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, int dim_max) {
-    rips_complex_->create_complex(*simplex_tree, dim_max);
+    if (rips_complex_)
+      rips_complex_->create_complex(*simplex_tree, dim_max);
+    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();
   }
 
  private:
-  Rips_complex<Simplex_tree_interface<>::Filtration_value>* rips_complex_;
+  // std::variant would work, but we don't require C++17 yet, and boost::variant is not super convenient.
+  // 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_rips_complex.py b/src/cython/test/test_rips_complex.py
index c37b5400..05dfcaf7 100755
--- a/src/cython/test/test_rips_complex.py
+++ b/src/cython/test/test_rips_complex.py
@@ -30,7 +30,6 @@ __license__ = "GPL v3"
 
 def test_empty_rips():
     rips_complex = RipsComplex()
-    assert rips_complex.__is_defined() == True
 
 def test_rips_from_points():
     point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
@@ -68,12 +67,26 @@ def test_filtered_rips_from_points():
     assert simplex_tree.num_simplices() == 8
     assert simplex_tree.num_vertices() == 4
 
+def test_sparse_filtered_rips_from_points():
+    point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
+    filtered_rips = RipsComplex(points=point_list, max_edge_length=1.0,
+                                sparse=.001)
+
+    simplex_tree = filtered_rips.create_simplex_tree(max_dimension=1)
+
+    assert simplex_tree.__is_defined() == True
+    assert simplex_tree.__is_persistence_defined() == False
+
+    assert simplex_tree.num_simplices() == 8
+    assert simplex_tree.num_vertices() == 4
+
 def test_rips_from_distance_matrix():
     distance_matrix = [[0],
                        [1, 0],
                        [1, sqrt(2), 0],
                        [sqrt(2), 1, 1, 0]]
-    rips_complex = RipsComplex(distance_matrix=distance_matrix, max_edge_length=42)
+    rips_complex = RipsComplex(distance_matrix=distance_matrix,
+                               max_edge_length=42)
 
     simplex_tree = rips_complex.create_simplex_tree(max_dimension=1)
 
@@ -100,7 +113,8 @@ def test_filtered_rips_from_distance_matrix():
                        [1, 0],
                        [1, sqrt(2), 0],
                        [sqrt(2), 1, 1, 0]]
-    filtered_rips = RipsComplex(distance_matrix=distance_matrix, max_edge_length=1.0)
+    filtered_rips = RipsComplex(distance_matrix=distance_matrix,
+                                max_edge_length=1.0)
 
     simplex_tree = filtered_rips.create_simplex_tree(max_dimension=1)
 
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