Merge branch 'master' into fix342

41 files changed, 1371 insertions, 165 deletions

diff --git a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
index 7c898dfd..e17831d9 100644
--- a/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
+++ b/src/Alpha_complex/utilities/alpha_complex_persistence.cpp
@@ -11,6 +11,7 @@
 #include <boost/program_options.hpp>
 
 #include <CGAL/Epick_d.h>
+#include <CGAL/Epeck_d.h>
 
 #include <gudhi/Alpha_complex.h>
 #include <gudhi/Persistent_cohomology.h>
diff --git a/src/cmake/modules/GUDHI_third_party_libraries.cmake b/src/cmake/modules/GUDHI_third_party_libraries.cmake
index 0abe66b7..a56a2756 100644
--- a/src/cmake/modules/GUDHI_third_party_libraries.cmake
+++ b/src/cmake/modules/GUDHI_third_party_libraries.cmake
@@ -68,7 +68,10 @@ if(CGAL_FOUND)
 endif()
 
 # For those who dislike bundled dependencies, this indicates where to find a preinstalled Hera.
-set(HERA_WASSERSTEIN_INCLUDE_DIR ${CMAKE_SOURCE_DIR}/ext/hera/wasserstein/include CACHE PATH "Directory where one can find Hera's wasserstein.h")
+set(HERA_WASSERSTEIN_INTERNAL_INCLUDE_DIR ${CMAKE_SOURCE_DIR}/ext/hera/wasserstein/include)
+set(HERA_WASSERSTEIN_INCLUDE_DIR ${HERA_WASSERSTEIN_INTERNAL_INCLUDE_DIR} CACHE PATH "Directory where one can find Hera's wasserstein.h")
+set(HERA_BOTTLENECK_INTERNAL_INCLUDE_DIR ${CMAKE_SOURCE_DIR}/ext/hera/bottleneck/include)
+set(HERA_BOTTLENECK_INCLUDE_DIR ${HERA_BOTTLENECK_INTERNAL_INCLUDE_DIR} CACHE PATH "Directory where one can find Hera's bottleneck.h")
 
 option(WITH_GUDHI_USE_TBB "Build with Intel TBB parallelization" ON)
 
@@ -199,7 +202,7 @@ if(PYTHONINTERP_FOUND AND CYTHON_FOUND)
     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_PYTHON_PATH}/doc/python3-sphinx-build.py")
+        set(SPHINX_PATH "${PYTHON_EXECUTABLE}" "-m" "sphinx.cmd.build")
       endif(PYTHON_VERSION_MAJOR EQUAL 3)
     endif(NOT SPHINX_PATH)
   endif(SPHINX_FOUND)
diff --git a/src/cmake/modules/GUDHI_user_version_target.cmake b/src/cmake/modules/GUDHI_user_version_target.cmake
index e99bb42d..e4f39aae 100644
--- a/src/cmake/modules/GUDHI_user_version_target.cmake
+++ b/src/cmake/modules/GUDHI_user_version_target.cmake
@@ -67,8 +67,11 @@ 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_directory ${CMAKE_SOURCE_DIR}/src/GudhUI ${GUDHI_USER_VERSION_DIR}/GudhUI)
 
-add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E
-                   copy_directory ${CMAKE_SOURCE_DIR}/ext/hera/wasserstein/include ${GUDHI_USER_VERSION_DIR}/ext/hera/wasserstein/include)
+if(HERA_WASSERSTEIN_INCLUDE_DIR STREQUAL HERA_WASSERSTEIN_INTERNAL_INCLUDE_DIR OR
+    HERA_BOTTLENECK_INCLUDE_DIR STREQUAL HERA_BOTTLENECK_INTERNAL_INCLUDE_DIR)
+  add_custom_command(TARGET user_version PRE_BUILD COMMAND ${CMAKE_COMMAND} -E
+                     copy_directory ${CMAKE_SOURCE_DIR}/ext/hera ${GUDHI_USER_VERSION_DIR}/ext/hera)
+endif()
 
 set(GUDHI_DIRECTORIES "doc;example;concept;utilities")
 
diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index ab08cd6d..04f64ce3 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -36,6 +36,7 @@ if(PYTHONINTERP_FOUND)
     add_gudhi_debug_info("Pybind11 version ${PYBIND11_VERSION}")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'bottleneck', ")
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'hera', ")
+    set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'clustering', ")
   endif()
   if(CYTHON_FOUND)
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'off_reader', ")
@@ -58,6 +59,7 @@ if(PYTHONINTERP_FOUND)
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'wasserstein', ")
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'point_cloud', ")
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'weighted_rips_complex', ")
+    set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'dtm_rips_complex', ")
 
     add_gudhi_debug_info("Python version ${PYTHON_VERSION_STRING}")
     add_gudhi_debug_info("Cython version ${CYTHON_VERSION}")
@@ -129,7 +131,9 @@ if(PYTHONINTERP_FOUND)
     set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'reader_utils', ")
     set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'witness_complex', ")
     set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'strong_witness_complex', ")
-    set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'hera', ")
+    set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'clustering/_tomato', ")
+    set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'hera/wasserstein', ")
+    set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'hera/bottleneck', ")
     if (NOT CGAL_VERSION VERSION_LESS 4.11.0)
       set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'bottleneck', ")
       set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'nerve_gic', ")
@@ -233,7 +237,13 @@ if(PYTHONINTERP_FOUND)
     file(COPY "gudhi/representations" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/")
     file(COPY "gudhi/wasserstein" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
     file(COPY "gudhi/point_cloud" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
+    file(COPY "gudhi/clustering" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi" FILES_MATCHING PATTERN "*.py")
     file(COPY "gudhi/weighted_rips_complex.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
+    file(COPY "gudhi/dtm_rips_complex.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi")
+    file(COPY "gudhi/hera/__init__.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/hera")
+
+    # Some files for pip package
+    file(COPY "introduction.rst" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/")
 
     add_custom_command(
         OUTPUT gudhi.so
@@ -353,7 +363,9 @@ if(PYTHONINTERP_FOUND)
                COMMAND ${CMAKE_COMMAND} -E env "PYTHONPATH=${CMAKE_CURRENT_BINARY_DIR}"
                ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/bottleneck_basic_example.py")
 
-      add_gudhi_py_test(test_bottleneck_distance)
+      if (PYBIND11_FOUND)
+        add_gudhi_py_test(test_bottleneck_distance)
+      endif()
 
       # Cover complex
       file(COPY ${CMAKE_SOURCE_DIR}/data/points/human.off DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
@@ -490,11 +502,22 @@ if(PYTHONINTERP_FOUND)
       add_gudhi_py_test(test_dtm)
     endif()
 
+    # Tomato
+    if(SCIPY_FOUND AND SKLEARN_FOUND AND PYBIND11_FOUND)
+      add_gudhi_py_test(test_tomato)
+    endif()
+
     # Weighted Rips
     if(SCIPY_FOUND)
-      add_gudhi_py_test(test_weighted_rips)
+      add_gudhi_py_test(test_weighted_rips_complex)
     endif()
 
+    # DTM Rips
+    if(SCIPY_FOUND)
+      add_gudhi_py_test(test_dtm_rips_complex)
+    endif()
+
+
     # Set missing or not modules
     set(GUDHI_MODULES ${GUDHI_MODULES} "python" CACHE INTERNAL "GUDHI_MODULES")
   else(CYTHON_FOUND)
diff --git a/src/python/doc/alpha_complex_sum.inc b/src/python/doc/alpha_complex_sum.inc
index 3aba0d71..aeab493f 100644
--- a/src/python/doc/alpha_complex_sum.inc
+++ b/src/python/doc/alpha_complex_sum.inc
@@ -3,15 +3,13 @@
 
    +----------------------------------------------------------------+-------------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------------+
    | .. figure::                                                    | Alpha complex is a simplicial complex constructed from the finite       | :Author: Vincent Rouvreau                                                                                                 |
-   |      ../../doc/Alpha_complex/alpha_complex_representation.png  | cells of a Delaunay Triangulation.                                      |                                                                                                                           |
-   |      :alt: Alpha complex representation                        |                                                                         | :Since: GUDHI 2.0.0                                                                                                       |
-   |      :figclass: align-center                                   | The filtration value of each simplex is computed as the **square** of   |                                                                                                                           |
-   |                                                                | the circumradius of the simplex if the circumsphere is empty (the       | :License: MIT (`GPL v3 </licensing/>`_)                                                                                   |
-   |                                                                | simplex is then said to be Gabriel), and as the minimum of the          |                                                                                                                           |
-   |                                                                | filtration values of the codimension 1 cofaces that make it not         | :Requires: `Eigen <installation.html#eigen>`_ :math:`\geq` 3.1.0 and `CGAL <installation.html#cgal>`_ :math:`\geq` 4.11.0 |
-   |                                                                | Gabriel otherwise.                                                      |                                                                                                                           |
+   |      ../../doc/Alpha_complex/alpha_complex_representation.png  | cells of a Delaunay Triangulation. It has the same persistent homology  |                                                                                                                           |
+   |      :alt: Alpha complex representation                        | as the Čech complex and is significantly smaller.                       | :Since: GUDHI 2.0.0                                                                                                       |
+   |      :figclass: align-center                                   |                                                                         |                                                                                                                           |
+   |                                                                |                                                                         | :License: MIT (`GPL v3 </licensing/>`_)                                                                                   |
+   |                                                                |                                                                         |                                                                                                                           |
+   |                                                                |                                                                         | :Requires: `Eigen <installation.html#eigen>`_ :math:`\geq` 3.1.0 and `CGAL <installation.html#cgal>`_ :math:`\geq` 4.11.0 |
    |                                                                |                                                                         |                                                                                                                           |
-   |                                                                | For performances reasons, it is advised to use CGAL :math:`\geq` 5.0.0. |                                                                                                                           |
    +----------------------------------------------------------------+-------------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------------+
    | * :doc:`alpha_complex_user`                                    | * :doc:`alpha_complex_ref`                                                                                                                                                                          |
    +----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------------+
diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst
index d49f45b4..097470c1 100644
--- a/src/python/doc/alpha_complex_user.rst
+++ b/src/python/doc/alpha_complex_user.rst
@@ -9,15 +9,33 @@ Definition
 
 .. include:: alpha_complex_sum.inc
 
-`AlphaComplex` is constructing a :doc:`SimplexTree <simplex_tree_ref>` using
+:doc:`AlphaComplex <alpha_complex_ref>` is constructing a :doc:`SimplexTree <simplex_tree_ref>` using
 `Delaunay Triangulation  <http://doc.cgal.org/latest/Triangulation/index.html#Chapter_Triangulations>`_
 :cite:`cgal:hdj-t-19b` from the `Computational Geometry Algorithms Library <http://www.cgal.org/>`_
 :cite:`cgal:eb-19b`.
 
 Remarks
 ^^^^^^^
-When an :math:`\alpha`-complex is constructed with an infinite value of :math:`\alpha^2`,
-the complex is a Delaunay complex (with special filtration values).
+* When an :math:`\alpha`-complex is constructed with an infinite value of :math:`\alpha^2`, the complex is a Delaunay
+  complex (with special filtration values). The Delaunay complex without filtration values is also available by
+  passing :code:`default_filtration_value = True` to :func:`~gudhi.AlphaComplex.create_simplex_tree`.
+* For people only interested in the topology of the Alpha complex (for instance persistence), Alpha complex is
+  equivalent to the `Čech complex <https://gudhi.inria.fr/doc/latest/group__cech__complex.html>`_ and much smaller if
+  you do not bound the radii. `Čech complex <https://gudhi.inria.fr/doc/latest/group__cech__complex.html>`_ can still
+  make sense in higher dimension precisely because you can bound the radii.
+* Using the default :code:`precision = 'safe'` makes the construction safe.
+  If you pass :code:`precision = 'exact'` to :func:`~gudhi.AlphaComplex.__init__`, the filtration values are the exact
+  ones converted to float. This can be very slow.
+  If you pass :code:`precision = 'safe'` (the default), the filtration values are only
+  guaranteed to have a small multiplicative error compared to the exact value.
+  A drawback, when computing persistence, is that an empty exact interval [10^12,10^12] may become a
+  non-empty approximate interval [10^12,10^12+10^6].
+  Using :code:`precision = 'fast'` makes the computations slightly faster, and the combinatorics are still exact, but
+  the computation of filtration values can exceptionally be arbitrarily bad. In all cases, we still guarantee that the
+  output is a valid filtration (faces have a filtration value no larger than their cofaces).
+* For performances reasons, it is advised to use Alpha_complex with `CGAL <installation.html#cgal>`_ :math:`\geq` 5.0.0.
+
+For performances reasons, it is advised to use CGAL :math:`\geq` 5.0.0.
 
 Example from points
 -------------------
diff --git a/src/python/doc/bottleneck_distance_user.rst b/src/python/doc/bottleneck_distance_user.rst
index 89da89d3..6c6e08d9 100644
--- a/src/python/doc/bottleneck_distance_user.rst
+++ b/src/python/doc/bottleneck_distance_user.rst
@@ -9,14 +9,23 @@ Definition
 
 .. include:: bottleneck_distance_sum.inc
 
-This implementation is based on ideas from "Geometry Helps in Bottleneck Matching and Related Problems"
-:cite:`DBLP:journals/algorithmica/EfratIK01`. Another relevant publication, although it was not used is
-"Geometry Helps to Compare Persistence Diagrams" :cite:`Kerber:2017:GHC:3047249.3064175`.
+This implementation by François Godi is based on ideas from "Geometry Helps in Bottleneck Matching and Related Problems"
+:cite:`DBLP:journals/algorithmica/EfratIK01` and requires `CGAL <installation.html#cgal>`_ (`GPL v3 </licensing/>`_).
 
-Function
---------
 .. autofunction:: gudhi.bottleneck_distance
 
+This other implementation comes from `Hera
+<https://bitbucket.org/grey_narn/hera/src/master/>`_ (BSD-3-Clause) which is
+based on "Geometry Helps to Compare Persistence Diagrams"
+:cite:`Kerber:2017:GHC:3047249.3064175` by Michael Kerber, Dmitriy
+Morozov, and Arnur Nigmetov.
+
+.. warning::
+   Beware that its approximation allows for a multiplicative error, while the function above uses an additive error.
+
+.. autofunction:: gudhi.hera.bottleneck_distance
+
+
 Distance computation
 --------------------
 
diff --git a/src/python/doc/clustering.inc b/src/python/doc/clustering.inc
new file mode 100644
index 00000000..2d07ae88
--- /dev/null
+++ b/src/python/doc/clustering.inc
@@ -0,0 +1,12 @@
+.. table::
+   :widths: 30 40 30
+
+   +--------------------------+-------------------------------------------------------+---------------------------------+
+   | .. figure::              | Clustering tools.                                     | :Author: Marc Glisse            |
+   |     img/spiral-color.png |                                                       |                                 |
+   |                          |                                                       | :Since: GUDHI 3.3.0             |
+   |                          |                                                       |                                 |
+   |                          |                                                       | :License: MIT                   |
+   +--------------------------+-------------------------------------------------------+---------------------------------+
+   | * :doc:`clustering`                                                                                                |
+   +--------------------------+-----------------------------------------------------------------------------------------+
diff --git a/src/python/doc/clustering.rst b/src/python/doc/clustering.rst
new file mode 100644
index 00000000..c5a57d3c
--- /dev/null
+++ b/src/python/doc/clustering.rst
@@ -0,0 +1,73 @@
+:orphan:
+
+.. To get rid of WARNING: document isn't included in any toctree
+
+=================
+Clustering manual
+=================
+
+We provide an implementation of ToMATo :cite:`tomato`, a persistence-based clustering algorithm. In short, this algorithm uses a density estimator and a neighborhood graph, starts with a mode-seeking phase (naive hill-climbing) to build initial clusters, and finishes by merging clusters based on their prominence.
+
+The merging phase depends on a parameter, which is the minimum prominence a cluster needs to avoid getting merged into another, bigger cluster. This parameter determines the number of clusters, and for convenience we allow you to choose instead the number of clusters. Decreasing the prominence threshold defines a hierarchy of clusters: if 2 points are in separate clusters when we have k clusters, they are still in different clusters for k+1 clusters.
+
+As a by-product, we produce the persistence diagram of the merge tree of the initial clusters. This is a convenient graphical tool to help decide how many clusters we want.
+
+.. plot::
+   :context:
+   :include-source:
+
+   import gudhi
+   f = open(gudhi.__root_source_dir__ + '/data/points/spiral_2d.csv', 'r')
+   import numpy as np
+   data = np.loadtxt(f)
+   import matplotlib.pyplot as plt
+   plt.scatter(data[:,0],data[:,1],marker='.',s=1)
+   plt.show()
+
+.. plot::
+   :context: close-figs
+   :include-source:
+
+   from gudhi.clustering.tomato import Tomato
+   t = Tomato()
+   t.fit(data)
+   t.plot_diagram()
+
+As one can see in `t.n_clusters_`, the algorithm found 6316 initial clusters. The diagram shows their prominence as their distance to the diagonal. There is always one point infinitely far: there is at least one cluster. Among the others, one point seems significantly farther from the diagonal than the others, which indicates that splitting the points into 2 clusters may be a sensible idea.
+
+.. plot::
+   :context: close-figs
+   :include-source:
+
+   t.n_clusters_=2
+   plt.scatter(data[:,0],data[:,1],marker='.',s=1,c=t.labels_)
+   plt.show()
+
+Of course this is just the result for one set of parameters. We can ask for a different density estimator and a different neighborhood graph computed with different parameters.
+
+.. plot::
+   :context: close-figs
+   :include-source:
+
+   t = Tomato(density_type='DTM', k=100)
+   t.fit(data)
+   t.plot_diagram()
+
+Makes the number of clusters clearer, and changes a bit the shape of the clusters.
+
+A quick look at the corresponding density estimate
+
+.. plot::
+   :context: close-figs
+   :include-source:
+
+   plt.scatter(data[:,0],data[:,1],marker='.',s=1,c=t.weights_)
+   plt.show()
+
+The code provides a few density estimators and graph constructions for convenience when first experimenting, but it is actually expected that advanced users provide their own graph and density estimates instead of point coordinates.
+
+Since the algorithm essentially computes basins of attraction, it is also encouraged to use it on functions that do not represent densities at all.
+
+.. autoclass:: gudhi.clustering.tomato.Tomato
+   :members:
+   :special-members: __init__
diff --git a/src/python/doc/cubical_complex_sum.inc b/src/python/doc/cubical_complex_sum.inc
index 28bf8e94..87db184d 100644
--- a/src/python/doc/cubical_complex_sum.inc
+++ b/src/python/doc/cubical_complex_sum.inc
@@ -2,9 +2,9 @@
    :widths: 30 40 30
 
    +--------------------------------------------------------------------------+----------------------------------------------------------------------+-----------------------------+
-   | .. figure::                                                              | The cubical complex is an example of a structured complex useful in  | :Author: Pawel Dlotko       |
-   |      ../../doc/Bitmap_cubical_complex/Cubical_complex_representation.png | computational mathematics (specially rigorous numerics) and image    |                             |
-   |      :alt: Cubical complex representation                                | analysis.                                                            | :Since: GUDHI 2.0.0         |
+   | .. figure::                                                              | The cubical complex represents a grid as a cell complex with         | :Author: Pawel Dlotko       |
+   |      ../../doc/Bitmap_cubical_complex/Cubical_complex_representation.png | cells of all dimensions.                                             |                             |
+   |      :alt: Cubical complex representation                                |                                                                      | :Since: GUDHI 2.0.0         |
    |      :figclass: align-center                                             |                                                                      |                             |
    |                                                                          |                                                                      | :License: MIT               |
    |                                                                          |                                                                      |                             |
diff --git a/src/python/doc/img/spiral-color.png b/src/python/doc/img/spiral-color.png
new file mode 100644
index 00000000..21b62dfc
--- /dev/null
+++ b/src/python/doc/img/spiral-color.png
diff --git a/src/python/doc/index.rst b/src/python/doc/index.rst
index 13e51047..040e57a4 100644
--- a/src/python/doc/index.rst
+++ b/src/python/doc/index.rst
@@ -53,8 +53,8 @@ Tangential complex
 Topological descriptors computation
 ***********************************
 
-Persistence cohomology
-======================
+Persistent cohomology
+=====================
 
 .. include:: persistent_cohomology_sum.inc
 
@@ -86,3 +86,8 @@ Point cloud utilities
 *********************
 
 .. include:: point_cloud_sum.inc
+
+Clustering
+**********
+
+.. include:: clustering.inc
diff --git a/src/python/doc/installation.rst b/src/python/doc/installation.rst
index de09c5b3..525ca84e 100644
--- a/src/python/doc/installation.rst
+++ b/src/python/doc/installation.rst
@@ -5,20 +5,47 @@
 Installation
 ############
 
-Conda
-*****
-The easiest way to install the Python version of GUDHI is using
-`conda <https://gudhi.inria.fr/conda/>`_.
+Packages
+********
+The easiest way to install the Python version of GUDHI is using pre-built packages.
+We recommend `conda <https://gudhi.inria.fr/conda/>`_
+
+.. code-block:: bash
+
+    conda install -c conda-forge gudhi
+
+Gudhi is also available on `PyPI <https://pypi.org/project/gudhi/>`_
+
+.. code-block:: bash
+
+    pip install gudhi
+
+Third party packages are also available, for instance on Debian or Ubuntu
+
+.. code-block:: bash
+
+    apt install python3-gudhi
+
+In all cases, you may still want to install some of the optional `run time dependencies`_.
 
 Compiling
 *********
+These instructions are for people who want to compile gudhi from source, they are
+unnecessary if you installed a binary package of Gudhi as above. They assume that
+you have downloaded a `release <https://github.com/GUDHI/gudhi-devel/releases>`_,
+with a name like `gudhi.3.2.0.tar.gz`, then run `tar xf gudhi.3.2.0.tar.gz`, which
+created a directory `gudhi.3.2.0`, hereinafter referred to as `/path-to-gudhi/`.
+If you are instead using a git checkout, beware that the paths are a bit
+different, and in particular the `python/` subdirectory is actually `src/python/`
+there.
+
 The library uses c++14 and requires `Boost <https://www.boost.org/>`_ :math:`\geq` 1.56.0,
 `CMake <https://www.cmake.org/>`_ :math:`\geq` 3.1  to generate makefiles,
 `NumPy <http://numpy.org>`_, `Cython <https://www.cython.org/>`_ and
 `pybind11 <https://github.com/pybind/pybind11>`_ to compile
 the GUDHI Python module.
 It is a multi-platform library and compiles on Linux, Mac OSX and Visual
-Studio 2017.
+Studio 2017 or later.
 
 On `Windows <https://wiki.python.org/moin/WindowsCompilers>`_ , only Python
 :math:`\geq` 3.5 are available because of the required Visual Studio version.
@@ -260,6 +287,13 @@ a flag `enable_autodiff=True` is used). In order to reduce code duplication, we
 use `EagerPy <https://eagerpy.jonasrauber.de/>`_ which wraps arrays from
 PyTorch, TensorFlow and JAX in a common interface.
 
+Joblib
+------
+
+`Joblib <https://joblib.readthedocs.io/>`_ is used both as a dependency of `Scikit-learn`_,
+and directly for parallelism in some modules (:class:`~gudhi.point_cloud.knn.KNearestNeighbors`,
+:func:`~gudhi.representations.metrics.pairwise_persistence_diagram_distances`).
+
 Hnswlib
 -------
 
diff --git a/src/python/doc/nerve_gic_complex_user.rst b/src/python/doc/nerve_gic_complex_user.rst
index 0e67fc78..0b820abf 100644
--- a/src/python/doc/nerve_gic_complex_user.rst
+++ b/src/python/doc/nerve_gic_complex_user.rst
@@ -50,7 +50,7 @@ The cover C comes from the preimages of intervals (10 intervals with gain 0.3)
 covering the height function (coordinate 2),
 which are then refined into their connected components using the triangulation of the .OFF file.
 
-.. testcode::
+.. code-block:: python
 
     import gudhi
     nerve_complex = gudhi.CoverComplex()
@@ -99,9 +99,6 @@ the program output is:
       [-0.171433, 0.367393]
       [-0.909111, 0.745853]
     0 interval(s) in dimension 1:
-
-.. testoutput::
-
     Nerve is of dimension 1 - 41 simplices - 21 vertices.
     [0]
     [1]
diff --git a/src/python/doc/persistent_cohomology_sum.inc b/src/python/doc/persistent_cohomology_sum.inc
index a1ff2eee..58e44b8a 100644
--- a/src/python/doc/persistent_cohomology_sum.inc
+++ b/src/python/doc/persistent_cohomology_sum.inc
@@ -6,9 +6,7 @@
    |      ../../doc/Persistent_cohomology/3DTorus_poch.png           | a sequence of (homology) groups, capturing global topological         |                                               |
    |      :figclass: align-center                                    | features like connected components, holes, cavities, etc. Persistent  | :Since: GUDHI 2.0.0                           |
    |                                                                 | homology studies the evolution -- birth, life and death -- of these   |                                               |
-   |      Rips Persistent Cohomology on a 3D                         | features when the topological space is changing. Consequently, the    | :License: MIT                                 |
-   |      Torus                                                      | theory is essentially composed of three elements: topological spaces, |                                               |
-   |                                                                 | their homology groups and an evolution scheme.                        |                                               |
+   |      Rips Persistent Cohomology on a 3D Torus                   | features when the topological space is changing.                      | :License: MIT                                 |
    |                                                                 |                                                                       |                                               |
    |                                                                 | Computation of persistent cohomology using the algorithm of           |                                               |
    |                                                                 | :cite:`DBLP:journals/dcg/SilvaMV11` and                               |                                               |
diff --git a/src/python/doc/python3-sphinx-build.py b/src/python/doc/python3-sphinx-build.py
deleted file mode 100755
index 84d158cf..00000000
--- a/src/python/doc/python3-sphinx-build.py
+++ /dev/null
@@ -1,11 +0,0 @@
-#!/usr/bin/env python3
-
-"""
-Emulate sphinx-build for python3
-"""
-
-from sys import exit, argv
-from sphinx import main
-
-if __name__ == '__main__':
-    exit(main(argv))
diff --git a/src/python/doc/rips_complex_ref.rst b/src/python/doc/rips_complex_ref.rst
index 5f3e46c1..f0582d5c 100644
--- a/src/python/doc/rips_complex_ref.rst
+++ b/src/python/doc/rips_complex_ref.rst
@@ -13,8 +13,6 @@ Rips complex reference manual
 
    .. automethod:: gudhi.RipsComplex.__init__
 
-.. _weighted-rips-complex-reference-manual:
-
 ======================================
 Weighted Rips complex reference manual
 ======================================
@@ -25,3 +23,14 @@ Weighted Rips complex reference manual
    :show-inheritance:
 
    .. automethod:: gudhi.weighted_rips_complex.WeightedRipsComplex.__init__
+
+=================================
+DTM Rips complex reference manual
+=================================
+
+.. autoclass:: gudhi.dtm_rips_complex.DTMRipsComplex
+   :members:
+   :undoc-members:
+   :show-inheritance:
+
+   .. automethod:: gudhi.dtm_rips_complex.DTMRipsComplex.__init__
\ No newline at end of file
diff --git a/src/python/doc/rips_complex_sum.inc b/src/python/doc/rips_complex_sum.inc
index f7580714..c123ea2a 100644
--- a/src/python/doc/rips_complex_sum.inc
+++ b/src/python/doc/rips_complex_sum.inc
@@ -2,18 +2,13 @@
    :widths: 30 40 30
 
    +----------------------------------------------------------------+------------------------------------------------------------------------+----------------------------------------------------------------------+
-   | .. figure::                                                    | Rips complex is a simplicial complex constructed from a one skeleton   | :Authors: Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse |
-   |      ../../doc/Rips_complex/rips_complex_representation.png    | graph.                                                                 |                                                                      |
+   | .. figure::                                                    | The Vietoris-Rips complex is a simplicial complex built as the         | :Authors: Clément Maria, Pawel Dlotko, Vincent Rouvreau, Marc Glisse |
+   |      ../../doc/Rips_complex/rips_complex_representation.png    | clique-complex of a proximity graph.                                   |                                                                      |
    |      :figclass: align-center                                   |                                                                        | :Since: GUDHI 2.0.0                                                  |
-   |                                                                | The filtration value of each edge is computed from a user-given        |                                                                      |
-   |                                                                | distance function and is inserted until a user-given threshold         | :License: MIT                                                        |
-   |                                                                | value.                                                                 |                                                                      |
+   |                                                                | We also provide sparse approximations, to speed-up the computation     |                                                                      |
+   |                                                                | of persistent homology, and weighted versions, which are more robust   | :License: MIT                                                        |
+   |                                                                | to outliers.                                                           |                                                                      |
    |                                                                |                                                                        |                                                                      |
-   |                                                                | This complex can be built from a point cloud and a distance function,  |                                                                      |
-   |                                                                | or from a distance matrix.                                             |                                                                      |
-   |                                                                |                                                                        |                                                                      |
-   |                                                                | Weighted Rips complex constructs a simplicial complex from a distance  |                                                                      |
-   |                                                                | matrix and weights on vertices.                                        |                                                                      |
    +----------------------------------------------------------------+------------------------------------------------------------------------+----------------------------------------------------------------------+
    | * :doc:`rips_complex_user`                                     | * :doc:`rips_complex_ref`                                                                                                                     |
    +----------------------------------------------------------------+-----------------------------------------------------------------------------------------------------------------------------------------------+
diff --git a/src/python/doc/rips_complex_user.rst b/src/python/doc/rips_complex_user.rst
index 819568be..6048cc4e 100644
--- a/src/python/doc/rips_complex_user.rst
+++ b/src/python/doc/rips_complex_user.rst
@@ -378,6 +378,7 @@ Example from a point cloud combined with DistanceToMeasure
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
 
 Combining with DistanceToMeasure, one can compute the DTM-filtration of a point set, as in `this notebook <https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-DTM-filtrations.ipynb>`_. 
+Remark that `DTMRipsComplex <rips_complex_user.html#dtm-rips-complex>`_ class provides exactly this function.
 
 .. testcode::
 
@@ -398,3 +399,24 @@ The output is:
 .. testoutput::
 
     [(0, (3.1622776601683795, inf)), (0, (3.1622776601683795, 5.39834563766817)), (0, (3.1622776601683795, 5.39834563766817))]
+
+DTM Rips Complex
+----------------
+
+:class:`~gudhi.dtm_rips_complex.DTMRipsComplex` builds a simplicial complex from a point set or a full distance matrix (in the form of ndarray), as described in the above example.
+This class constructs a weighted Rips complex giving larger weights to outliers, which reduces their impact on the persistence diagram. See `this notebook <https://github.com/GUDHI/TDA-tutorial/blob/master/Tuto-GUDHI-DTM-filtrations.ipynb>`_ for some experiments.
+
+.. testcode::
+
+    import numpy as np
+    from gudhi.dtm_rips_complex import DTMRipsComplex
+    pts = np.array([[2.0, 2.0], [0.0, 1.0], [3.0, 4.0]])
+    dtm_rips = DTMRipsComplex(points=pts, k=2)
+    st = dtm_rips.create_simplex_tree(max_dimension=2)
+    print(st.persistence())
+
+The output is:
+
+.. testoutput::
+
+    [(0, (3.1622776601683795, inf)), (0, (3.1622776601683795, 5.39834563766817)), (0, (3.1622776601683795, 5.39834563766817))]
diff --git a/src/python/doc/tangential_complex_sum.inc b/src/python/doc/tangential_complex_sum.inc
index 22314a2d..2f330a07 100644
--- a/src/python/doc/tangential_complex_sum.inc
+++ b/src/python/doc/tangential_complex_sum.inc
@@ -3,10 +3,10 @@
 
    +----------------------------------------------------------------+------------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------------+
    | .. figure::                                                    | A Tangential Delaunay complex is a simplicial complex designed to      | :Author:    Clément Jamin                                                                                                 |
-   |      ../../doc/Tangential_complex/tc_examples.png              | reconstruct a :math:`k`-dimensional manifold embedded in :math:`d`-    |                                                                                                                           |
-   |      :figclass: align-center                                   | dimensional Euclidean space. The input is a point sample coming from   | :Since:    GUDHI 2.0.0                                                                                                    |
-   |                                                                | an unknown manifold. The running time depends only linearly on the     |                                                                                                                           |
-   |                                                                | extrinsic dimension :math:`d` and exponentially on the intrinsic       | :License: MIT (`GPL v3 </licensing/>`_)                                                                                   |
+   |      ../../doc/Tangential_complex/tc_examples.png              | reconstruct a :math:`k`-dimensional manifold embedded in               |                                                                                                                           |
+   |      :figclass: align-center                                   | :math:`d`-dimensional Euclidean space. The input is a point sample     | :Since:    GUDHI 2.0.0                                                                                                    |
+   |                                                                | coming from an unknown manifold. The running time depends only linearly|                                                                                                                           |
+   |                                                                | on the extrinsic dimension :math:`d` and exponentially on the intrinsic| :License: MIT (`GPL v3 </licensing/>`_)                                                                                   |
    |                                                                | dimension :math:`k`.                                                   |                                                                                                                           |
    |                                                                |                                                                        | :Requires: `Eigen <installation.html#eigen>`_ :math:`\geq` 3.1.0 and `CGAL <installation.html#cgal>`_ :math:`\geq` 4.11.0 |
    +----------------------------------------------------------------+------------------------------------------------------------------------+---------------------------------------------------------------------------------------------------------------------------+
diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx
index d75e374a..a356384d 100644
--- a/src/python/gudhi/alpha_complex.pyx
+++ b/src/python/gudhi/alpha_complex.pyx
@@ -27,11 +27,11 @@ __license__ = "GPL v3"
 
 cdef extern from "Alpha_complex_interface.h" namespace "Gudhi":
     cdef cppclass Alpha_complex_interface "Gudhi::alpha_complex::Alpha_complex_interface":
-        Alpha_complex_interface(vector[vector[double]] points) nogil except +
+        Alpha_complex_interface(vector[vector[double]] points, bool fast_version) nogil except +
         # bool from_file is a workaround for cython to find the correct signature
-        Alpha_complex_interface(string off_file, bool from_file) nogil except +
+        Alpha_complex_interface(string off_file, bool fast_version, bool from_file) nogil except +
         vector[double] get_point(int vertex) nogil except +
-        void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square) nogil except +
+        void create_simplex_tree(Simplex_tree_interface_full_featured* simplex_tree, double max_alpha_square, bool exact_version, bool default_filtration_value) nogil except +
 
 # AlphaComplex python interface
 cdef class AlphaComplex:
@@ -53,10 +53,11 @@ cdef class AlphaComplex:
 
     """
 
-    cdef Alpha_complex_interface * thisptr
+    cdef Alpha_complex_interface * this_ptr
+    cdef bool exact
 
     # Fake constructor that does nothing but documenting the constructor
-    def __init__(self, points=None, off_file=''):
+    def __init__(self, points=None, off_file='', precision='safe'):
         """AlphaComplex constructor.
 
         :param points: A list of points in d-Dimension.
@@ -66,15 +67,21 @@ cdef class AlphaComplex:
 
         :param off_file: An OFF file style name.
         :type off_file: string
+
+        :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'.
+        :type precision: string
         """
 
     # The real cython constructor
-    def __cinit__(self, points = None, off_file = ''):
+    def __cinit__(self, points = None, off_file = '', precision = 'safe'):
+        assert precision in ['fast', 'safe', 'exact'], "Alpha complex precision can only be 'fast', 'safe' or 'exact'"
+        cdef bool fast = precision == 'fast'
+        self.exact = precision == 'exact'
+
         cdef vector[vector[double]] pts
         if off_file:
             if os.path.isfile(off_file):
-                self.thisptr = new Alpha_complex_interface(
-                    off_file.encode('utf-8'), True)
+                self.this_ptr = new Alpha_complex_interface(off_file.encode('utf-8'), fast, True)
             else:
                 print("file " + off_file + " not found.")
         else:
@@ -83,17 +90,16 @@ cdef class AlphaComplex:
                 points=[]
             pts = points
             with nogil:
-                self.thisptr = new Alpha_complex_interface(pts)
-                
+                self.this_ptr = new Alpha_complex_interface(pts, fast)
 
     def __dealloc__(self):
-        if self.thisptr != NULL:
-            del self.thisptr
+        if self.this_ptr != NULL:
+            del self.this_ptr
 
     def __is_defined(self):
         """Returns true if AlphaComplex pointer is not NULL.
          """
-        return self.thisptr != NULL
+        return self.this_ptr != NULL
 
     def get_point(self, vertex):
         """This function returns the point corresponding to a given vertex.
@@ -103,21 +109,24 @@ cdef class AlphaComplex:
         :rtype: list of float
         :returns: the point.
         """
-        return self.thisptr.get_point(vertex)
+        return self.this_ptr.get_point(vertex)
 
-    def create_simplex_tree(self, max_alpha_square = float('inf')):
+    def create_simplex_tree(self, max_alpha_square = float('inf'), default_filtration_value = False):
         """
-        :param max_alpha_square: The maximum alpha square threshold the
-            simplices shall not exceed. Default is set to infinity, and
-            there is very little point using anything else since it does
-            not save time.
+        :param max_alpha_square: The maximum alpha square threshold the simplices shall not exceed. Default is set to
+            infinity, and there is very little point using anything else since it does not save time.
         :type max_alpha_square: float
+        :param default_filtration_value: Set this value to `True` if filtration values are not needed to be computed
+            (will be set to `NaN`). Default value is `False` (which means compute the filtration values).
+        :type default_filtration_value: bool
         :returns: A simplex tree created from the Delaunay Triangulation.
         :rtype: SimplexTree
         """
         stree = SimplexTree()
         cdef double mas = max_alpha_square
         cdef intptr_t stree_int_ptr=stree.thisptr
+        cdef bool compute_filtration = default_filtration_value == True
         with nogil:
-            self.thisptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr, mas)
+            self.this_ptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr,
+                                              mas, self.exact, compute_filtration)
         return stree
diff --git a/src/python/gudhi/bottleneck.cc b/src/python/gudhi/bottleneck.cc
index 732cb9a8..8a3d669a 100644
--- a/src/python/gudhi/bottleneck.cc
+++ b/src/python/gudhi/bottleneck.cc
@@ -12,24 +12,29 @@
 
 #include <pybind11_diagram_utils.h>
 
-double bottleneck(Dgm d1, Dgm d2, double epsilon)
+// For compatibility with older versions, we want to support e=None.
+// In C++17, the recommended way is std::optional<double>.
+double bottleneck(Dgm d1, Dgm d2, py::object epsilon)
 {
+  double e = (std::numeric_limits<double>::min)();
+  if (!epsilon.is_none()) e = epsilon.cast<double>();
   // I *think* the call to request() has to be before releasing the GIL.
   auto diag1 = numpy_to_range_of_pairs(d1);
   auto diag2 = numpy_to_range_of_pairs(d2);
 
   py::gil_scoped_release release;
 
-  return Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, epsilon);
+  return Gudhi::persistence_diagram::bottleneck_distance(diag1, diag2, e);
 }
 
 PYBIND11_MODULE(bottleneck, m) {
       m.attr("__license__") = "GPL v3";
       m.def("bottleneck_distance", &bottleneck,
           py::arg("diagram_1"), py::arg("diagram_2"),
-          py::arg("e") = (std::numeric_limits<double>::min)(),
+          py::arg("e") = py::none(),
           R"pbdoc(
-    This function returns the point corresponding to a given vertex.
+    Compute the Bottleneck distance between two diagrams.
+    Points at infinity and on the diagonal are supported.
 
     :param diagram_1: The first diagram.
     :type diagram_1: numpy array of shape (m,2)
@@ -42,7 +47,7 @@ PYBIND11_MODULE(bottleneck, m) {
         bits of the mantissa may be wrong). This version of the algorithm takes
         advantage of the limited precision of `double` and is usually a lot
         faster to compute, whatever the value of `e`.
-        Thus, by default, `e` is the smallest positive double.
+        Thus, by default (`e=None`), `e` is the smallest positive double.
     :type e: float
     :rtype: float
     :returns: the bottleneck distance.
diff --git a/src/python/gudhi/clustering/__init__.py b/src/python/gudhi/clustering/__init__.py
new file mode 100644
index 00000000..e69de29b
--- /dev/null
+++ b/src/python/gudhi/clustering/__init__.py
diff --git a/src/python/gudhi/clustering/_tomato.cc b/src/python/gudhi/clustering/_tomato.cc
new file mode 100644
index 00000000..a76a2c3a
--- /dev/null
+++ b/src/python/gudhi/clustering/_tomato.cc
@@ -0,0 +1,277 @@
+/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ *    Author(s):       Marc Glisse
+ *
+ *    Copyright (C) 2020 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#include <boost/container/flat_map.hpp>
+#include <boost/pending/disjoint_sets.hpp>
+#include <boost/property_map/property_map.hpp>
+#include <boost/property_map/transform_value_property_map.hpp>
+#include <boost/property_map/vector_property_map.hpp>
+#include <boost/property_map/function_property_map.hpp>
+#include <boost/iterator/counting_iterator.hpp>
+#include <boost/range/irange.hpp>
+#include <boost/range/adaptor/transformed.hpp>
+#include <vector>
+#include <unordered_map>
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+#include <iostream>
+
+namespace py = pybind11;
+
+template <class T, class = std::enable_if_t<std::is_integral<T>::value>>
+int getint(int i) {
+  return i;
+}
+// Gcc-8 has a bug that breaks this version, fixed in gcc-9
+// template<class T, class=decltype(std::declval<T>().template cast<int>())>
+// int getint(T i){return i.template cast<int>();}
+template <class T>
+auto getint(T i) -> decltype(i.template cast<int>()) {
+  return i.template cast<int>();
+}
+
+// Raw clusters are clusters obtained through single-linkage, no merging.
+
+typedef int Point_index;
+typedef int Cluster_index;
+struct Merge {
+  Cluster_index first, second;
+  double persist;
+};
+
+template <class Neighbors, class Density, class Order, class ROrder>
+auto tomato(Point_index num_points, Neighbors const& neighbors, Density const& density, Order const& order,
+            ROrder const& rorder) {
+  // point index --> index of raw cluster it belongs to
+  std::vector<Cluster_index> raw_cluster;
+  raw_cluster.reserve(num_points);
+  // cluster index --> index of top point in the cluster
+  Cluster_index n_raw_clusters = 0;  // current number of raw clusters seen
+  //
+  std::vector<Merge> merges;
+  struct Data {
+    Cluster_index parent;
+    int rank;
+    Point_index max;
+  };  // information on a cluster
+  std::vector<Data> ds_base;
+  // boost::vector_property_map does resize(size+1) for every new element, don't use it
+  auto ds_data =
+      boost::make_function_property_map<std::size_t>([&ds_base](std::size_t n) -> Data& { return ds_base[n]; });
+  auto ds_parent =
+      boost::make_transform_value_property_map([](auto& p) -> Cluster_index& { return p.parent; }, ds_data);
+  auto ds_rank = boost::make_transform_value_property_map([](auto& p) -> int& { return p.rank; }, ds_data);
+  boost::disjoint_sets<decltype(ds_rank), decltype(ds_parent)> ds(
+      ds_rank, ds_parent);                         // on the clusters, not directly the points
+  std::vector<std::array<double, 2>> persistence;  // diagram (finite points)
+  boost::container::flat_map<Cluster_index, Cluster_index>
+      adj_clusters;  // first: the merged cluster, second: the raw cluster
+  // we only care about the raw cluster, we could use a vector to store the second, store first into a set, and only
+  // insert in the vector if merged is absent from the set
+
+  for (Point_index i = 0; i < num_points; ++i) {
+    // auto&& ngb = neighbors[order[i]];
+    // TODO: have a specialization where we don't need python types and py::cast
+    // TODO: move py::cast and getint into Neighbors
+    py::object ngb = neighbors[py::cast(order[i])];  // auto&& also seems to work
+    adj_clusters.clear();
+    Point_index j = i;  // highest neighbor
+    // for(Point_index k : ngb)
+    for (auto k_py : ngb) {
+      Point_index k = rorder[getint(k_py)];
+      if (k >= i || k < 0)  // ???
+        continue;
+      if (k < j) j = k;
+      Cluster_index rk = raw_cluster[k];
+      adj_clusters.emplace(ds.find_set(rk), rk);
+      // does not insert if ck=ds.find_set(rk) already seen
+      // which rk we keep from those with the same ck is arbitrary
+    }
+    assert((Point_index)raw_cluster.size() == i);
+    if (i == j) {  // local maximum -> new cluster
+      Cluster_index c = n_raw_clusters++;
+      ds_base.emplace_back();  // could be integrated in ds_data, but then we would check the size for every access
+      ds.make_set(c);
+      ds_base[c].max = i;  // max
+      raw_cluster.push_back(c);
+      continue;
+    } else {  // add i to the cluster of j
+      assert(j < i);
+      raw_cluster.push_back(raw_cluster[j]);
+      // FIXME: we are adding point i to the raw cluster of j, but that one might not be in adj_clusters, so we may
+      // merge clusters A and B through a point of C. It is strange, but I don't know if it can really cause problems.
+      // We could just not set j at all and use arbitrarily the first element of adj_clusters.
+    }
+    // possibly merge clusters
+    // we could sort, in case there are several merges, but that doesn't seem so useful
+    Cluster_index rj = raw_cluster[j];
+    Cluster_index cj = ds.find_set(rj);
+    Cluster_index orig_cj = cj;
+    for (auto ckk : adj_clusters) {
+      Cluster_index rk = ckk.second;
+      Cluster_index ck = ckk.first;
+      if (ck == orig_cj) continue;
+      assert(ck == ds.find_set(rk));
+      Point_index j = ds_base[cj].max;
+      Point_index k = ds_base[ck].max;
+      Point_index young = std::max(j, k);
+      Point_index old = std::min(j, k);
+      auto d_young = density[order[young]];
+      auto d_i = density[order[i]];
+      assert(d_young >= d_i);
+      // Always merge (the non-hierarchical algorithm would only conditionally merge here
+      persistence.push_back({d_young, d_i});
+      assert(ds.find_set(rj) != ds.find_set(rk));
+      ds.link(cj, ck);
+      cj = ds.find_set(cj);
+      ds_base[cj].max = old;  // just one parent, no need for find_set
+      // record the raw clusters, we don't know what will have already been merged.
+      merges.push_back({rj, rk, d_young - d_i});
+    }
+  }
+  {
+    boost::counting_iterator<int> b(0), e(ds_base.size());
+    ds.compress_sets(b, e);
+    // Now we stop using find_sets and look at the parent directly
+    // rank is reused to rename clusters contiguously 0, 1, etc
+  }
+  // Maximum for each connected component
+  std::vector<double> max_cc;
+  for (Cluster_index i = 0; i < n_raw_clusters; ++i) {
+    if (ds_base[i].parent == i) max_cc.push_back(density[order[ds_base[i].max]]);
+  }
+  assert((Cluster_index)(merges.size() + max_cc.size()) == n_raw_clusters);
+
+  // TODO: create a "noise" cluster, merging all those not prominent enough?
+
+  // Replay the merges, in increasing order of prominence, to build the hierarchy
+  std::sort(merges.begin(), merges.end(), [](Merge const& a, Merge const& b) { return a.persist < b.persist; });
+  std::vector<std::array<Cluster_index, 2>> children;
+  children.reserve(merges.size());
+  {
+    struct Dat {
+      Cluster_index parent;
+      int rank;
+      Cluster_index name;
+    };
+    std::vector<Dat> ds_bas(2 * n_raw_clusters - 1);
+    Cluster_index i;
+    auto ds_dat =
+        boost::make_function_property_map<std::size_t>([&ds_bas](std::size_t n) -> Dat& { return ds_bas[n]; });
+    auto ds_par = boost::make_transform_value_property_map([](auto& p) -> Cluster_index& { return p.parent; }, ds_dat);
+    auto ds_ran = boost::make_transform_value_property_map([](auto& p) -> int& { return p.rank; }, ds_dat);
+    boost::disjoint_sets<decltype(ds_ran), decltype(ds_par)> ds(ds_ran, ds_par);
+    for (i = 0; i < n_raw_clusters; ++i) {
+      ds.make_set(i);
+      ds_bas[i].name = i;
+    }
+    for (Merge const& m : merges) {
+      Cluster_index j = ds.find_set(m.first);
+      Cluster_index k = ds.find_set(m.second);
+      assert(j != k);
+      children.push_back({ds_bas[j].name, ds_bas[k].name});
+      ds.make_set(i);
+      ds.link(i, j);
+      ds.link(i, k);
+      ds_bas[ds.find_set(i)].name = i;
+      ++i;
+    }
+  }
+
+  std::vector<Cluster_index> raw_cluster_ordered(num_points);
+  for (int i = 0; i < num_points; ++i) raw_cluster_ordered[i] = raw_cluster[rorder[i]];
+  // return raw_cluster, children, persistence
+  // TODO avoid copies: https://github.com/pybind/pybind11/issues/1042
+  return py::make_tuple(py::array(raw_cluster_ordered.size(), raw_cluster_ordered.data()),
+                        py::array(children.size(), children.data()), py::array(persistence.size(), persistence.data()),
+                        py::array(max_cc.size(), max_cc.data()));
+}
+
+auto merge(py::array_t<Cluster_index, py::array::c_style> children, Cluster_index n_leaves, Cluster_index n_final) {
+  if (n_final > n_leaves) {
+    std::cerr << "The number of clusters required " << n_final << " is larger than the number of mini-clusters " << n_leaves << '\n';
+    n_final = n_leaves; // or return something special and let Tomato use leaf_labels_?
+  }
+  py::buffer_info cbuf = children.request();
+  if ((cbuf.ndim != 2 || cbuf.shape[1] != 2) && (cbuf.ndim != 1 || cbuf.shape[0] != 0))
+    throw std::runtime_error("internal error: children have to be (n,2) or empty");
+  const int n_merges = cbuf.shape[0];
+  Cluster_index* d = (Cluster_index*)cbuf.ptr;
+  if (n_merges + n_final < n_leaves) {
+    std::cerr << "The number of clusters required " << n_final << " is smaller than the number of connected components " << n_leaves - n_merges << '\n';
+    n_final = n_leaves - n_merges;
+  }
+  struct Dat {
+    Cluster_index parent;
+    int rank;
+    int name;
+  };
+  std::vector<Dat> ds_bas(2 * n_leaves - 1);
+  auto ds_dat = boost::make_function_property_map<std::size_t>([&ds_bas](std::size_t n) -> Dat& { return ds_bas[n]; });
+  auto ds_par = boost::make_transform_value_property_map([](auto& p) -> Cluster_index& { return p.parent; }, ds_dat);
+  auto ds_ran = boost::make_transform_value_property_map([](auto& p) -> int& { return p.rank; }, ds_dat);
+  boost::disjoint_sets<decltype(ds_ran), decltype(ds_par)> ds(ds_ran, ds_par);
+  Cluster_index i;
+  for (i = 0; i < n_leaves; ++i) {
+    ds.make_set(i);
+    ds_bas[i].name = -1;
+  }
+  for (Cluster_index m = 0; m < n_leaves - n_final; ++m) {
+    Cluster_index j = ds.find_set(d[2 * m]);
+    Cluster_index k = ds.find_set(d[2 * m + 1]);
+    assert(j != k);
+    ds.make_set(i);
+    ds.link(i, j);
+    ds.link(i, k);
+    ++i;
+  }
+  Cluster_index next_cluster_name = 0;
+  std::vector<Cluster_index> ret;
+  ret.reserve(n_leaves);
+  for (Cluster_index j = 0; j < n_leaves; ++j) {
+    Cluster_index k = ds.find_set(j);
+    if (ds_bas[k].name == -1) ds_bas[k].name = next_cluster_name++;
+    ret.push_back(ds_bas[k].name);
+  }
+  return py::array(ret.size(), ret.data());
+}
+
+// TODO: Do a special version when ngb is a numpy array, where we can cast to int[k][n] ?
+// py::isinstance<py::array_t<std::int32_t>> (ou py::isinstance<py::array> et tester dtype) et flags&c_style
+// ou overload (en virant forcecast?)
+// aussi le faire au cas où on n'aurait pas un tableau, mais où chaque liste de voisins serait un tableau ?
+auto hierarchy(py::handle ngb, py::array_t<double, py::array::c_style | py::array::forcecast> density) {
+  // used to be py::iterable ngb, but that's inconvenient if it doesn't come pre-sorted
+  // use py::handle and check if [] (aka __getitem__) works? But then we need to build an object to pass it to []
+  // (I _think_ handle is ok and we don't need object here)
+  py::buffer_info wbuf = density.request();
+  if (wbuf.ndim != 1) throw std::runtime_error("density must be 1D");
+  const int n = wbuf.shape[0];
+  double* d = (double*)wbuf.ptr;
+  // Vector { 0, 1, ..., n-1 }
+  std::vector<Point_index> order(boost::counting_iterator<Point_index>(0), boost::counting_iterator<Point_index>(n));
+  // Permutation of the indices to get points in decreasing order of density
+  std::sort(std::begin(order), std::end(order), [=](Point_index i, Point_index j) { return d[i] > d[j]; });
+  // Inverse permutation
+  std::vector<Point_index> rorder(n);
+  for (Point_index i : boost::irange(0, n)) rorder[order[i]] = i;
+  // Used as:
+  // order[i] is the index of the point with i-th largest density
+  // rorder[i] is the rank of the i-th point in order of decreasing density
+  // TODO: put a wrapper on ngb and d so we don't need to pass (r)order (there is still the issue of reordering the
+  // output)
+  return tomato(n, ngb, d, order, rorder);
+}
+
+PYBIND11_MODULE(_tomato, m) {
+  m.doc() = "Internals of tomato clustering";
+  m.def("hierarchy", &hierarchy, "does the clustering");
+  m.def("merge", &merge, "merge clusters");
+}
diff --git a/src/python/gudhi/clustering/tomato.py b/src/python/gudhi/clustering/tomato.py
new file mode 100644
index 00000000..fbba3cc8
--- /dev/null
+++ b/src/python/gudhi/clustering/tomato.py
@@ -0,0 +1,321 @@
+# This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+# See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+# Author(s):       Marc Glisse
+#
+# Copyright (C) 2020 Inria
+#
+# Modification(s):
+#   - YYYY/MM Author: Description of the modification
+
+import numpy
+from ..point_cloud.knn import KNearestNeighbors
+from ..point_cloud.dtm import DTMDensity
+from ._tomato import *
+
+# The fit/predict interface is not so well suited...
+
+
+class Tomato:
+    """
+    This clustering algorithm needs a neighborhood graph on the points, and an estimation of the density at each point.
+    A few possible graph constructions and density estimators are provided for convenience, but it is perfectly natural
+    to provide your own.
+
+    :Requires: `SciPy <installation.html#scipy>`_, `Scikit-learn <installation.html#scikit-learn>`_ or others
+        (see :class:`~gudhi.point_cloud.knn.KNearestNeighbors`) in function of the options.
+
+    Attributes
+    ----------
+    n_clusters_: int
+        The number of clusters. Writing to it automatically adjusts `labels_`.
+    merge_threshold_: float
+        minimum prominence of a cluster so it doesn't get merged. Writing to it automatically adjusts `labels_`.
+    n_leaves_: int
+        number of leaves (unstable clusters) in the hierarchical tree
+    leaf_labels_: ndarray of shape (n_samples,)
+        cluster labels for each point, at the very bottom of the hierarchy
+    labels_: ndarray of shape (n_samples,)
+        cluster labels for each point, after merging
+    diagram_: ndarray of shape (`n_leaves_`, 2)
+        persistence diagram (only the finite points)
+    max_weight_per_cc_: ndarray of shape (n_connected_components,)
+        maximum of the density function on each connected component. This corresponds to the abscissa of infinite
+        points in the diagram
+    children_: ndarray of shape (`n_leaves_`-n_connected_components, 2)
+        The children of each non-leaf node. Values less than `n_leaves_` correspond to leaves of the tree.
+        A node i greater than or equal to `n_leaves_` is a non-leaf node and has children children_[i - `n_leaves_`].
+        Alternatively at the i-th iteration, children[i][0] and children[i][1] are merged to form node `n_leaves_` + i
+    weights_: ndarray of shape (n_samples,)
+        weights of the points, as computed by the density estimator or provided by the user
+    params_: dict
+        Parameters like metric, etc
+    """
+
+    def __init__(
+        self,
+        graph_type="knn",
+        density_type="logDTM",
+        n_clusters=None,
+        merge_threshold=None,
+        #       eliminate_threshold=None,
+        #           eliminate_threshold (float): minimum max weight of a cluster so it doesn't get eliminated
+        **params
+    ):
+        """
+        Args:
+            graph_type (str): 'manual', 'knn' or 'radius'. Default is 'knn'.
+            density_type (str): 'manual', 'DTM', 'logDTM', 'KDE' or 'logKDE'. When you have many points,
+                'KDE' and 'logKDE' tend to be slower. Default is 'logDTM'.
+            metric (str|Callable): metric used when calculating the distance between instances in a feature array.
+                Defaults to Minkowski of parameter p.
+            kde_params (dict): if density_type is 'KDE' or 'logKDE', additional parameters passed directly to
+                sklearn.neighbors.KernelDensity.
+            k (int): number of neighbors for a knn graph (including the vertex itself). Defaults to 10.
+            k_DTM (int): number of neighbors for the DTM density estimation (including the vertex itself).
+                Defaults to k.
+            r (float): size of a neighborhood if graph_type is 'radius'. Also used as default bandwidth in kde_params.
+            eps (float): (1+eps) approximation factor when computing distances (ignored in many cases).
+            n_clusters (int): number of clusters requested. Defaults to None, i.e. no merging occurs and we get
+                the maximal number of clusters.
+            merge_threshold (float): minimum prominence of a cluster so it doesn't get merged.
+            symmetrize_graph (bool): whether we should add edges to make the neighborhood graph symmetric.
+                This can be useful with k-NN for small k. Defaults to false.
+            p (float): norm L^p on input points. Defaults to 2.
+            q (float): order used to compute the distance to measure. Defaults to dim.
+                Beware that when the dimension is large, this can easily cause overflows.
+            dim (float): final exponent in DTM density estimation, representing the dimension. Defaults to the
+                dimension, or 2 when the dimension cannot be read from the input (metric is "precomputed").
+            n_jobs (int): Number of jobs to schedule for parallel processing on the CPU.
+                If -1 is given all processors are used. Default: 1.
+            params: extra parameters are passed to :class:`~gudhi.point_cloud.knn.KNearestNeighbors` and
+                :class:`~gudhi.point_cloud.dtm.DTMDensity`.
+        """
+        # Should metric='precomputed' mean input_type='distance_matrix'?
+        # Should we be able to pass metric='minkowski' (what None does currently)?
+        self.graph_type_ = graph_type
+        self.density_type_ = density_type
+        self.params_ = params
+        self.__n_clusters = n_clusters
+        self.__merge_threshold = merge_threshold
+        # self.eliminate_threshold_ = eliminate_threshold
+        if n_clusters and merge_threshold:
+            raise ValueError("Cannot specify both a merge threshold and a number of clusters")
+
+    def fit(self, X, y=None, weights=None):
+        """
+        Args:
+            X ((n,d)-array of float|(n,n)-array of float|Sequence[Iterable[int]]): coordinates of the points,
+                or distance matrix (full, not just a triangle) if metric is "precomputed", or list of neighbors
+                for each point (points are represented by their index, starting from 0) if graph_type is "manual".
+                The number of points is currently limited to about 2 billion.
+            weights (ndarray of shape (n_samples)): if density_type is 'manual', a density estimate at each point
+            y: Not used, present here for API consistency with scikit-learn by convention.
+        """
+        # TODO: First detect if this is a new call with the same data (only threshold changed?)
+        # TODO: less code duplication (subroutines?), less spaghetti, but don't compute neighbors twice if not needed. Clear error message for missing or contradictory parameters.
+        if weights is not None:
+            density_type = "manual"
+        else:
+            density_type = self.density_type_
+            if density_type == "manual":
+                raise ValueError("If density_type is 'manual', you must provide weights to fit()")
+
+        if self.graph_type_ == "manual":
+            self.neighbors_ = X
+            # FIXME: uniformize "message 'option'" vs 'message "option"'
+            assert density_type == "manual", 'If graph_type is "manual", density_type must be as well'
+        else:
+            metric = self.params_.get("metric", "minkowski")
+            if metric != "precomputed":
+                self.points_ = X
+
+        # Slight complication to avoid computing knn twice.
+        need_knn = 0
+        need_knn_ngb = False
+        need_knn_dist = False
+        if self.graph_type_ == "knn":
+            k_graph = self.params_.get("k", 10)
+            # If X has fewer than k points...
+            if k_graph > len(X):
+                k_graph = len(X)
+            need_knn = k_graph
+            need_knn_ngb = True
+        if self.density_type_ in ["DTM", "logDTM"]:
+            k = self.params_.get("k", 10)
+            k_DTM = self.params_.get("k_DTM", k)
+            # If X has fewer than k points...
+            if k_DTM > len(X):
+                k_DTM = len(X)
+            need_knn = max(need_knn, k_DTM)
+            need_knn_dist = True
+            # if we ask for more neighbors for the graph than the DTM, getting the distances is a slight waste,
+            # but it looks negligible
+        if need_knn > 0:
+            knn_args = dict(self.params_)
+            knn_args["k"] = need_knn
+            knn = KNearestNeighbors(return_index=need_knn_ngb, return_distance=need_knn_dist, **knn_args).fit_transform(
+                X
+            )
+            if need_knn_ngb:
+                if need_knn_dist:
+                    self.neighbors_ = knn[0][:, 0:k_graph]
+                    knn_dist = knn[1]
+                else:
+                    self.neighbors_ = knn
+            elif need_knn_dist:
+                knn_dist = knn
+        if self.density_type_ in ["DTM", "logDTM"]:
+            dim = self.params_.get("dim")
+            if dim is None:
+                dim = len(X[0]) if metric != "precomputed" else 2
+            q = self.params_.get("q", dim)
+            weights = DTMDensity(k=k_DTM, metric="neighbors", dim=dim, q=q).fit_transform(knn_dist)
+            if self.density_type_ == "logDTM":
+                weights = numpy.log(weights)
+
+        if self.graph_type_ == "radius":
+            if metric in ["minkowski", "euclidean", "manhattan", "chebyshev"]:
+                from scipy.spatial import cKDTree
+
+                tree = cKDTree(X)
+                # TODO: handle "l1" and "l2" aliases?
+                p = self.params_.get("p")
+                if metric == "euclidean":
+                    assert p is None or p == 2, "p=" + str(p) + " is not consistent with metric='euclidean'"
+                    p = 2
+                elif metric == "manhattan":
+                    assert p is None or p == 1, "p=" + str(p) + " is not consistent with metric='manhattan'"
+                    p = 1
+                elif metric == "chebyshev":
+                    assert p is None or p == numpy.inf, "p=" + str(p) + " is not consistent with metric='chebyshev'"
+                    p = numpy.inf
+                elif p is None:
+                    p = 2  # the default
+                eps = self.params_.get("eps", 0)
+                self.neighbors_ = tree.query_ball_tree(tree, r=self.params_["r"], p=p, eps=eps)
+
+            # TODO: sklearn's NearestNeighbors.radius_neighbors can handle more metrics efficiently via its BallTree
+            # (don't bother with the _graph variant, it just calls radius_neighbors).
+            elif metric != "precomputed":
+                from sklearn.metrics import pairwise_distances
+
+                X = pairwise_distances(X, metric=metric, n_jobs=self.params_.get("n_jobs"))
+                metric = "precomputed"
+
+            if metric == "precomputed":
+                # TODO: parallelize? May not be worth it.
+                X = numpy.asarray(X)
+                r = self.params_["r"]
+                self.neighbors_ = [numpy.flatnonzero(l <= r) for l in X]
+
+        if self.density_type_ in {"KDE", "logKDE"}:
+            # Slow...
+            assert (
+                self.graph_type_ != "manual" and metric != "precomputed"
+            ), "Scikit-learn's KernelDensity requires point coordinates"
+            kde_params = dict(self.params_.get("kde_params", dict()))
+            kde_params.setdefault("metric", metric)
+            r = self.params_.get("r")
+            if r is not None:
+                kde_params.setdefault("bandwidth", r)
+            # Should we default rtol to eps?
+            from sklearn.neighbors import KernelDensity
+
+            weights = KernelDensity(**kde_params).fit(self.points_).score_samples(self.points_)
+            if self.density_type_ == "KDE":
+                weights = numpy.exp(weights)
+
+        # TODO: do it at the C++ level and/or in parallel if this is too slow?
+        if self.params_.get("symmetrize_graph"):
+            self.neighbors_ = [set(line) for line in self.neighbors_]
+            for i, line in enumerate(self.neighbors_):
+                line.discard(i)
+                for j in line:
+                    self.neighbors_[j].add(i)
+
+        self.weights_ = weights
+        # This is where the main computation happens
+        self.leaf_labels_, self.children_, self.diagram_, self.max_weight_per_cc_ = hierarchy(self.neighbors_, weights)
+        self.n_leaves_ = len(self.max_weight_per_cc_) + len(self.children_)
+        assert self.leaf_labels_.max() + 1 == len(self.max_weight_per_cc_) + len(self.children_)
+        # TODO: deduplicate this code with the setters below
+        if self.__merge_threshold:
+            assert not self.__n_clusters
+            self.__n_clusters = numpy.count_nonzero(
+                self.diagram_[:, 0] - self.diagram_[:, 1] > self.__merge_threshold
+            ) + len(self.max_weight_per_cc_)
+        if self.__n_clusters:
+            # TODO: set corresponding merge_threshold?
+            renaming = merge(self.children_, self.n_leaves_, self.__n_clusters)
+            self.labels_ = renaming[self.leaf_labels_]
+            # In case the user asked for something impossible.
+            # TODO: check for impossible situations before calling merge.
+            self.__n_clusters = self.labels_.max() + 1
+        else:
+            self.labels_ = self.leaf_labels_
+            self.__n_clusters = self.n_leaves_
+        return self
+
+    def fit_predict(self, X, y=None, weights=None):
+        """
+        Equivalent to fit(), and returns the `labels_`.
+        """
+        return self.fit(X, y, weights).labels_
+
+    # TODO: add argument k or threshold? Have a version where you can click and it shows the line and the corresponding k?
+    def plot_diagram(self):
+        """
+        """
+        import matplotlib.pyplot as plt
+
+        l = self.max_weight_per_cc_.min()
+        r = self.max_weight_per_cc_.max()
+        if self.diagram_.size > 0:
+            plt.plot(self.diagram_[:, 0], self.diagram_[:, 1], "ro")
+            l = min(l, self.diagram_[:, 1].min())
+            r = max(r, self.diagram_[:, 0].max())
+        if l == r:
+            if l > 0:
+                l, r = 0.9 * l, 1.1 * r
+            elif l < 0:
+                l, r = 1.1 * l, 0.9 * r
+            else:
+                l, r = -1.0, 1.0
+        plt.plot([l, r], [l, r])
+        plt.plot(
+            self.max_weight_per_cc_, numpy.full(self.max_weight_per_cc_.shape, 1.1 * l - 0.1 * r), "ro", color="green"
+        )
+        plt.show()
+
+    # Use set_params instead?
+    @property
+    def n_clusters_(self):
+        return self.__n_clusters
+
+    @n_clusters_.setter
+    def n_clusters_(self, n_clusters):
+        if n_clusters == self.__n_clusters:
+            return
+        self.__n_clusters = n_clusters
+        self.__merge_threshold = None
+        if hasattr(self, "leaf_labels_"):
+            renaming = merge(self.children_, self.n_leaves_, self.__n_clusters)
+            self.labels_ = renaming[self.leaf_labels_]
+            # In case the user asked for something impossible
+            self.__n_clusters = self.labels_.max() + 1
+
+    @property
+    def merge_threshold_(self):
+        return self.__merge_threshold
+
+    @merge_threshold_.setter
+    def merge_threshold_(self, merge_threshold):
+        if merge_threshold == self.__merge_threshold:
+            return
+        if hasattr(self, "leaf_labels_"):
+            self.n_clusters_ = numpy.count_nonzero(self.diagram_[:, 0] - self.diagram_[:, 1] > merge_threshold) + len(
+                self.max_weight_per_cc_
+            )
+        else:
+            self.__n_clusters = None
+        self.__merge_threshold = merge_threshold
diff --git a/src/python/gudhi/dtm_rips_complex.py b/src/python/gudhi/dtm_rips_complex.py
new file mode 100644
index 00000000..63c9b138
--- /dev/null
+++ b/src/python/gudhi/dtm_rips_complex.py
@@ -0,0 +1,51 @@
+# This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+# See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+# Author(s):       Yuichi Ike, Raphaël Tinarrage
+#
+# Copyright (C) 2020 Inria, Copyright (C) 2020 FUjitsu Laboratories Ltd.
+#
+# Modification(s):
+#   - YYYY/MM Author: Description of the modification
+
+
+from gudhi.weighted_rips_complex import WeightedRipsComplex
+from gudhi.point_cloud.dtm import DistanceToMeasure
+from scipy.spatial.distance import cdist
+
+class DTMRipsComplex(WeightedRipsComplex):
+    """
+    Class to generate a DTM Rips complex from a distance matrix or a point set, 
+    in the way described in :cite:`dtmfiltrations`.
+    Remark that all the filtration values are doubled compared to the definition in the paper 
+    for the consistency with RipsComplex.
+    :Requires: `SciPy <installation.html#scipy>`_
+    """
+    def __init__(self, 
+                 points=None, 
+                 distance_matrix=None, 
+                 k=1, 
+                 q=2,
+                 max_filtration=float('inf')):
+        """
+        Args:
+            points (numpy.ndarray): array of points.
+            distance_matrix (numpy.ndarray): full distance matrix.
+            k (int): number of neighbors for the computation of DTM. Defaults to 1, which is equivalent to the usual Rips complex.
+            q (float): order used to compute the distance to measure. Defaults to 2.
+            max_filtration (float): specifies the maximal filtration value to be considered.      
+        """
+        if distance_matrix is None:
+            if points is None:
+                # Empty Rips construction
+                points=[]
+            distance_matrix = cdist(points,points)
+        self.distance_matrix = distance_matrix
+
+        # TODO: address the error when k is too large 
+        if k <= 1:
+            self.weights = [0] * len(distance_matrix)
+        else:
+            dtm = DistanceToMeasure(k, q=q, metric="precomputed")        
+            self.weights = dtm.fit_transform(distance_matrix)
+        self.max_filtration = max_filtration
+        
diff --git a/src/python/gudhi/hera/__init__.py b/src/python/gudhi/hera/__init__.py
new file mode 100644
index 00000000..f70b92b9
--- /dev/null
+++ b/src/python/gudhi/hera/__init__.py
@@ -0,0 +1,7 @@
+from .wasserstein import wasserstein_distance
+from .bottleneck import bottleneck_distance
+
+
+__author__ = "Marc Glisse"
+__copyright__ = "Copyright (C) 2020 Inria"
+__license__ = "MIT"
diff --git a/src/python/gudhi/hera/bottleneck.cc b/src/python/gudhi/hera/bottleneck.cc
new file mode 100644
index 00000000..0cb562ce
--- /dev/null
+++ b/src/python/gudhi/hera/bottleneck.cc
@@ -0,0 +1,54 @@
+/*    This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+ *    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+ *    Author(s):       Marc Glisse
+ *
+ *    Copyright (C) 2020 Inria
+ *
+ *    Modification(s):
+ *      - YYYY/MM Author: Description of the modification
+ */
+
+#include <pybind11_diagram_utils.h>
+
+#ifdef _MSC_VER
+// https://github.com/grey-narn/hera/issues/3
+// ssize_t is a non-standard type (well, posix)
+using py::ssize_t;
+#endif
+
+#include <bottleneck.h> // Hera
+
+double bottleneck_distance(Dgm d1, Dgm d2, double delta)
+{
+  // I *think* the call to request() has to be before releasing the GIL.
+  auto diag1 = numpy_to_range_of_pairs(d1);
+  auto diag2 = numpy_to_range_of_pairs(d2);
+
+  py::gil_scoped_release release;
+
+  if (delta == 0)
+    return hera::bottleneckDistExact(diag1, diag2);
+  else
+    return hera::bottleneckDistApprox(diag1, diag2, delta);
+}
+
+PYBIND11_MODULE(bottleneck, m) {
+      m.def("bottleneck_distance", &bottleneck_distance,
+          py::arg("X"), py::arg("Y"),
+          py::arg("delta") = .01,
+          R"pbdoc(
+        Compute the Bottleneck distance between two diagrams.
+        Points at infinity are supported.
+
+        .. note::
+           Points on the diagonal are not supported and must be filtered out before calling this function.
+
+        Parameters:
+            X (n x 2 numpy array): First diagram
+            Y (n x 2 numpy array): Second diagram
+            delta (float): Relative error 1+delta
+
+        Returns:
+            float: (approximate) bottleneck distance d_B(X,Y)
+    )pbdoc");
+}
diff --git a/src/python/gudhi/hera.cc b/src/python/gudhi/hera/wasserstein.cc
index ea80a9a8..1a21f02f 100644
--- a/src/python/gudhi/hera.cc
+++ b/src/python/gudhi/hera/wasserstein.cc
@@ -33,7 +33,7 @@ double wasserstein_distance(
   return hera::wasserstein_dist(diag1, diag2, params);
 }
 
-PYBIND11_MODULE(hera, m) {
+PYBIND11_MODULE(wasserstein, m) {
       m.def("wasserstein_distance", &wasserstein_distance,
           py::arg("X"), py::arg("Y"),
           py::arg("order") = 1,
diff --git a/src/python/gudhi/representations/kernel_methods.py b/src/python/gudhi/representations/kernel_methods.py
index 596f4f07..23fd23c7 100644
--- a/src/python/gudhi/representations/kernel_methods.py
+++ b/src/python/gudhi/representations/kernel_methods.py
@@ -10,7 +10,7 @@
 import numpy as np
 from sklearn.base import BaseEstimator, TransformerMixin
 from sklearn.metrics import pairwise_distances, pairwise_kernels
-from .metrics import SlicedWassersteinDistance, PersistenceFisherDistance, _sklearn_wrapper, pairwise_persistence_diagram_distances, _sliced_wasserstein_distance, _persistence_fisher_distance
+from .metrics import SlicedWassersteinDistance, PersistenceFisherDistance, _sklearn_wrapper, _pairwise, pairwise_persistence_diagram_distances, _sliced_wasserstein_distance, _persistence_fisher_distance
 from .preprocessing import Padding
 
 #############################################
@@ -60,7 +60,7 @@ def _persistence_scale_space_kernel(D1, D2, kernel_approx=None, bandwidth=1.):
     weight_pss = lambda x: 1 if x[1] >= x[0] else -1
     return 0.5 * _persistence_weighted_gaussian_kernel(DD1, DD2, weight=weight_pss, kernel_approx=kernel_approx, bandwidth=bandwidth)
 
-def pairwise_persistence_diagram_kernels(X, Y=None, kernel="sliced_wasserstein", **kwargs):
+def pairwise_persistence_diagram_kernels(X, Y=None, kernel="sliced_wasserstein", n_jobs=None, **kwargs):
     """
     This function computes the kernel matrix between two lists of persistence diagrams given as numpy arrays of shape (nx2).
 
@@ -68,38 +68,41 @@ def pairwise_persistence_diagram_kernels(X, Y=None, kernel="sliced_wasserstein",
         X (list of n numpy arrays of shape (numx2)): first list of persistence diagrams. 
         Y (list of m numpy arrays of shape (numx2)): second list of persistence diagrams (optional). If None, pairwise kernel values are computed from the first list only.
         kernel: kernel to use. It can be either a string ("sliced_wasserstein", "persistence_scale_space", "persistence_weighted_gaussian", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs. If it is a function, make sure that it is symmetric.
+        n_jobs (int): number of jobs to use for the computation. This uses joblib.Parallel(prefer="threads"), so kernels that do not release the GIL may not scale unless run inside a `joblib.parallel_backend <https://joblib.readthedocs.io/en/latest/parallel.html#joblib.parallel_backend>`_ block.
         **kwargs: optional keyword parameters. Any further parameters are passed directly to the kernel function. See the docs of the various kernel classes in this module.
 
     Returns: 
         numpy array of shape (nxm): kernel matrix.
     """    
     XX = np.reshape(np.arange(len(X)), [-1,1])
-    YY = None if Y is None else np.reshape(np.arange(len(Y)), [-1,1])
+    YY = None if Y is None or Y is X else np.reshape(np.arange(len(Y)), [-1,1])
     if kernel == "sliced_wasserstein":
-        return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="sliced_wasserstein", num_directions=kwargs["num_directions"]) / kwargs["bandwidth"])
+        return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="sliced_wasserstein", num_directions=kwargs["num_directions"], n_jobs=n_jobs) / kwargs["bandwidth"])
     elif kernel == "persistence_fisher":
-        return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="persistence_fisher", kernel_approx=kwargs["kernel_approx"], bandwidth=kwargs["bandwidth"]) / kwargs["bandwidth_fisher"])
+        return np.exp(-pairwise_persistence_diagram_distances(X, Y, metric="persistence_fisher", kernel_approx=kwargs["kernel_approx"], bandwidth=kwargs["bandwidth"], n_jobs=n_jobs) / kwargs["bandwidth_fisher"])
     elif kernel == "persistence_scale_space":
-        return pairwise_kernels(XX, YY, metric=_sklearn_wrapper(_persistence_scale_space_kernel, X, Y, **kwargs))
+        return _pairwise(pairwise_kernels, False, XX, YY, metric=_sklearn_wrapper(_persistence_scale_space_kernel, X, Y, **kwargs), n_jobs=n_jobs)
     elif kernel == "persistence_weighted_gaussian":
-        return pairwise_kernels(XX, YY, metric=_sklearn_wrapper(_persistence_weighted_gaussian_kernel, X, Y, **kwargs))
+        return _pairwise(pairwise_kernels, False, XX, YY, metric=_sklearn_wrapper(_persistence_weighted_gaussian_kernel, X, Y, **kwargs), n_jobs=n_jobs)
     else:
-        return pairwise_kernels(XX, YY, metric=_sklearn_wrapper(metric, **kwargs))
+        return _pairwise(pairwise_kernels, False, XX, YY, metric=_sklearn_wrapper(metric, **kwargs), n_jobs=n_jobs)
 
 class SlicedWassersteinKernel(BaseEstimator, TransformerMixin):
     """
     This is a class for computing the sliced Wasserstein kernel matrix from a list of persistence diagrams. The sliced Wasserstein kernel is computed by exponentiating the corresponding sliced Wasserstein distance with a Gaussian kernel. See http://proceedings.mlr.press/v70/carriere17a.html for more details. 
     """
-    def __init__(self, num_directions=10, bandwidth=1.0):
+    def __init__(self, num_directions=10, bandwidth=1.0, n_jobs=None):
         """
         Constructor for the SlicedWassersteinKernel class.
 
         Parameters:
             bandwidth (double): bandwidth of the Gaussian kernel applied to the sliced Wasserstein distance (default 1.).
             num_directions (int): number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the kernel computation (default 10).
+            n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_kernels` for details.
         """
         self.bandwidth = bandwidth
         self.num_directions = num_directions
+        self.n_jobs = n_jobs
 
     def fit(self, X, y=None):
         """
@@ -122,7 +125,7 @@ class SlicedWassersteinKernel(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise sliced Wasserstein kernel values.
         """
-        return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="sliced_wasserstein", bandwidth=self.bandwidth, num_directions=self.num_directions) 
+        return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="sliced_wasserstein", bandwidth=self.bandwidth, num_directions=self.num_directions, n_jobs=self.n_jobs)
 
     def __call__(self, diag1, diag2):
         """
@@ -141,7 +144,7 @@ class PersistenceWeightedGaussianKernel(BaseEstimator, TransformerMixin):
     """
     This is a class for computing the persistence weighted Gaussian kernel matrix from a list of persistence diagrams. The persistence weighted Gaussian kernel is computed by convolving the persistence diagram points with weighted Gaussian kernels. See http://proceedings.mlr.press/v48/kusano16.html for more details. 
     """
-    def __init__(self, bandwidth=1., weight=lambda x: 1, kernel_approx=None):
+    def __init__(self, bandwidth=1., weight=lambda x: 1, kernel_approx=None, n_jobs=None):
         """
         Constructor for the PersistenceWeightedGaussianKernel class.
   
@@ -149,9 +152,11 @@ class PersistenceWeightedGaussianKernel(BaseEstimator, TransformerMixin):
             bandwidth (double): bandwidth of the Gaussian kernel with which persistence diagrams will be convolved (default 1.)
             weight (function): weight function for the persistence diagram points (default constant function, ie lambda x: 1). This function must be defined on 2D points, ie lists or numpy arrays of the form [p_x,p_y].
             kernel_approx (class): kernel approximation class used to speed up computation (default None). Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
+            n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_kernels` for details.
         """
         self.bandwidth, self.weight = bandwidth, weight
         self.kernel_approx = kernel_approx
+        self.n_jobs = n_jobs
 
     def fit(self, X, y=None):
         """
@@ -174,7 +179,7 @@ class PersistenceWeightedGaussianKernel(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence weighted Gaussian kernel values.
         """
-        return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_weighted_gaussian", bandwidth=self.bandwidth, weight=self.weight, kernel_approx=self.kernel_approx) 
+        return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_weighted_gaussian", bandwidth=self.bandwidth, weight=self.weight, kernel_approx=self.kernel_approx, n_jobs=self.n_jobs)
 
     def __call__(self, diag1, diag2):
         """
@@ -193,15 +198,17 @@ class PersistenceScaleSpaceKernel(BaseEstimator, TransformerMixin):
     """
     This is a class for computing the persistence scale space kernel matrix from a list of persistence diagrams. The persistence scale space kernel is computed by adding the symmetric to the diagonal of each point in each persistence diagram, with negative weight, and then convolving the points with a Gaussian kernel. See https://www.cv-foundation.org/openaccess/content_cvpr_2015/papers/Reininghaus_A_Stable_Multi-Scale_2015_CVPR_paper.pdf for more details. 
     """
-    def __init__(self, bandwidth=1., kernel_approx=None):
+    def __init__(self, bandwidth=1., kernel_approx=None, n_jobs=None):
         """
         Constructor for the PersistenceScaleSpaceKernel class.
   
         Parameters:
             bandwidth (double): bandwidth of the Gaussian kernel with which persistence diagrams will be convolved (default 1.)
             kernel_approx (class): kernel approximation class used to speed up computation (default None). Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
+            n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_kernels` for details.
         """
         self.bandwidth, self.kernel_approx = bandwidth, kernel_approx
+        self.n_jobs = n_jobs
 
     def fit(self, X, y=None):
         """
@@ -224,7 +231,7 @@ class PersistenceScaleSpaceKernel(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence scale space kernel values.
         """
-        return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_scale_space", bandwidth=self.bandwidth, kernel_approx=self.kernel_approx)
+        return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_scale_space", bandwidth=self.bandwidth, kernel_approx=self.kernel_approx, n_jobs=self.n_jobs)
 
     def __call__(self, diag1, diag2):
         """
@@ -243,7 +250,7 @@ class PersistenceFisherKernel(BaseEstimator, TransformerMixin):
     """
     This is a class for computing the persistence Fisher kernel matrix from a list of persistence diagrams. The persistence Fisher kernel is computed by exponentiating the corresponding persistence Fisher distance with a Gaussian kernel. See papers.nips.cc/paper/8205-persistence-fisher-kernel-a-riemannian-manifold-kernel-for-persistence-diagrams for more details. 
     """
-    def __init__(self, bandwidth_fisher=1., bandwidth=1., kernel_approx=None):
+    def __init__(self, bandwidth_fisher=1., bandwidth=1., kernel_approx=None, n_jobs=None):
         """
         Constructor for the PersistenceFisherKernel class.
 
@@ -251,9 +258,11 @@ class PersistenceFisherKernel(BaseEstimator, TransformerMixin):
             bandwidth (double): bandwidth of the Gaussian kernel applied to the persistence Fisher distance (default 1.).
             bandwidth_fisher (double): bandwidth of the Gaussian kernel used to turn persistence diagrams into probability distributions by PersistenceFisherDistance class (default 1.).
             kernel_approx (class): kernel approximation class used to speed up computation (default None). Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).
+            n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_kernels` for details.
         """
         self.bandwidth = bandwidth
         self.bandwidth_fisher, self.kernel_approx = bandwidth_fisher, kernel_approx
+        self.n_jobs = n_jobs
 
     def fit(self, X, y=None):
         """
@@ -276,7 +285,7 @@ class PersistenceFisherKernel(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence Fisher kernel values.
         """
-        return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_fisher", bandwidth=self.bandwidth, bandwidth_fisher=self.bandwidth_fisher, kernel_approx=self.kernel_approx)
+        return pairwise_persistence_diagram_kernels(X, self.diagrams_, kernel="persistence_fisher", bandwidth=self.bandwidth, bandwidth_fisher=self.bandwidth_fisher, kernel_approx=self.kernel_approx, n_jobs=self.n_jobs)
 
     def __call__(self, diag1, diag2):
         """
diff --git a/src/python/gudhi/representations/metrics.py b/src/python/gudhi/representations/metrics.py
index 3e8487c4..142ddef1 100644
--- a/src/python/gudhi/representations/metrics.py
+++ b/src/python/gudhi/representations/metrics.py
@@ -12,6 +12,7 @@ from sklearn.base import BaseEstimator, TransformerMixin
 from sklearn.metrics import pairwise_distances
 from gudhi.hera import wasserstein_distance as hera_wasserstein_distance
 from .preprocessing import Padding
+from joblib import Parallel, delayed
 
 #############################################
 # Metrics ###################################
@@ -116,6 +117,20 @@ def _persistence_fisher_distance(D1, D2, kernel_approx=None, bandwidth=1.):
             vectorj = vectorj/vectorj_sum
         return np.arccos(  min(np.dot(np.sqrt(vectori), np.sqrt(vectorj)), 1.)  )
 
+def _pairwise(fallback, skipdiag, X, Y, metric, n_jobs):
+    if Y is not None:
+        return fallback(X, Y, metric=metric, n_jobs=n_jobs)
+    triu = np.triu_indices(len(X), k=skipdiag)
+    tril = (triu[1], triu[0])
+    par = Parallel(n_jobs=n_jobs, prefer="threads")
+    d = par(delayed(metric)([triu[0][i]], [triu[1][i]]) for i in range(len(triu[0])))
+    m = np.empty((len(X), len(X)))
+    m[triu] = d
+    m[tril] = d
+    if skipdiag:
+        np.fill_diagonal(m, 0)
+    return m
+
 def _sklearn_wrapper(metric, X, Y, **kwargs):
     """
     This function is a wrapper for any metric between two persistence diagrams that takes two numpy arrays of shapes (nx2) and (mx2) as arguments.
@@ -134,7 +149,7 @@ PAIRWISE_DISTANCE_FUNCTIONS = {
     "persistence_fisher": _persistence_fisher_distance,
 }
 
-def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwargs):
+def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", n_jobs=None, **kwargs):
     """
     This function computes the distance matrix between two lists of persistence diagrams given as numpy arrays of shape (nx2).
 
@@ -142,48 +157,51 @@ def pairwise_persistence_diagram_distances(X, Y=None, metric="bottleneck", **kwa
         X (list of n numpy arrays of shape (numx2)): first list of persistence diagrams. 
         Y (list of m numpy arrays of shape (numx2)): second list of persistence diagrams (optional). If None, pairwise distances are computed from the first list only.
         metric: distance to use. It can be either a string ("sliced_wasserstein", "wasserstein", "hera_wasserstein" (Wasserstein distance computed with Hera---note that Hera is also used for the default option "wasserstein"), "pot_wasserstein" (Wasserstein distance computed with POT), "bottleneck", "persistence_fisher") or a function taking two numpy arrays of shape (nx2) and (mx2) as inputs. If it is a function, make sure that it is symmetric and that it outputs 0 if called on the same two arrays. 
+        n_jobs (int): number of jobs to use for the computation. This uses joblib.Parallel(prefer="threads"), so metrics that do not release the GIL may not scale unless run inside a `joblib.parallel_backend <https://joblib.readthedocs.io/en/latest/parallel.html#joblib.parallel_backend>`_ block.
         **kwargs: optional keyword parameters. Any further parameters are passed directly to the distance function. See the docs of the various distance classes in this module.
 
     Returns: 
         numpy array of shape (nxm): distance matrix
     """
     XX = np.reshape(np.arange(len(X)), [-1,1])
-    YY = None if Y is None else np.reshape(np.arange(len(Y)), [-1,1]) 
+    YY = None if Y is None or Y is X else np.reshape(np.arange(len(Y)), [-1,1])
     if metric == "bottleneck":
         try: 
             from .. import bottleneck_distance
-            return pairwise_distances(XX, YY, metric=_sklearn_wrapper(bottleneck_distance, X, Y, **kwargs))
+            return _pairwise(pairwise_distances, True, XX, YY, metric=_sklearn_wrapper(bottleneck_distance, X, Y, **kwargs), n_jobs=n_jobs)
         except ImportError:
             print("Gudhi built without CGAL")
             raise
     elif metric == "pot_wasserstein":
         try:
             from gudhi.wasserstein import wasserstein_distance as pot_wasserstein_distance
-            return pairwise_distances(XX, YY, metric=_sklearn_wrapper(pot_wasserstein_distance,  X, Y, **kwargs))
+            return _pairwise(pairwise_distances, True, XX, YY, metric=_sklearn_wrapper(pot_wasserstein_distance,  X, Y, **kwargs), n_jobs=n_jobs)
         except ImportError:
             print("POT (Python Optimal Transport) is not installed. Please install POT or use metric='wasserstein' or metric='hera_wasserstein'")
             raise
     elif metric == "sliced_wasserstein":
         Xproj = _compute_persistence_diagram_projections(X, **kwargs)
         Yproj = None if Y is None else _compute_persistence_diagram_projections(Y, **kwargs)
-        return pairwise_distances(XX, YY, metric=_sklearn_wrapper(_sliced_wasserstein_distance_on_projections, Xproj, Yproj))
+        return _pairwise(pairwise_distances, True, XX, YY, metric=_sklearn_wrapper(_sliced_wasserstein_distance_on_projections, Xproj, Yproj), n_jobs=n_jobs)
     elif type(metric) == str:
-        return pairwise_distances(XX, YY, metric=_sklearn_wrapper(PAIRWISE_DISTANCE_FUNCTIONS[metric], X, Y, **kwargs))
+        return _pairwise(pairwise_distances, True, XX, YY, metric=_sklearn_wrapper(PAIRWISE_DISTANCE_FUNCTIONS[metric], X, Y, **kwargs), n_jobs=n_jobs)
     else:
-        return pairwise_distances(XX, YY, metric=_sklearn_wrapper(metric, X, Y, **kwargs))
+        return _pairwise(pairwise_distances, True, XX, YY, metric=_sklearn_wrapper(metric, X, Y, **kwargs), n_jobs=n_jobs)
 
 class SlicedWassersteinDistance(BaseEstimator, TransformerMixin):
     """
     This is a class for computing the sliced Wasserstein distance matrix from a list of persistence diagrams. The Sliced Wasserstein distance is computed by projecting the persistence diagrams onto lines, comparing the projections with the 1-norm, and finally integrating over all possible lines. See http://proceedings.mlr.press/v70/carriere17a.html for more details. 
     """
-    def __init__(self, num_directions=10):
+    def __init__(self, num_directions=10, n_jobs=None):
         """
         Constructor for the SlicedWassersteinDistance class.
 
         Parameters:
             num_directions (int): number of lines evenly sampled from [-pi/2,pi/2] in order to approximate and speed up the distance computation (default 10). 
+            n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_distances` for details.
         """
         self.num_directions = num_directions
+        self.n_jobs = n_jobs
 
     def fit(self, X, y=None):
         """
@@ -206,7 +224,7 @@ class SlicedWassersteinDistance(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise sliced Wasserstein distances.
         """
-        return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="sliced_wasserstein", num_directions=self.num_directions)
+        return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="sliced_wasserstein", num_directions=self.num_directions, n_jobs=self.n_jobs)
 
     def __call__(self, diag1, diag2):
         """
@@ -227,14 +245,16 @@ class BottleneckDistance(BaseEstimator, TransformerMixin):
 
     :Requires: `CGAL <installation.html#cgal>`_ :math:`\geq` 4.11.0
     """
-    def __init__(self, epsilon=None):
+    def __init__(self, epsilon=None, n_jobs=None):
         """
         Constructor for the BottleneckDistance class.
 
         Parameters:
             epsilon (double): absolute (additive) error tolerated on the distance (default is the smallest positive float), see :func:`gudhi.bottleneck_distance`.
+            n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_distances` for details.
         """
         self.epsilon = epsilon
+        self.n_jobs = n_jobs
 
     def fit(self, X, y=None):
         """
@@ -257,7 +277,7 @@ class BottleneckDistance(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise bottleneck distances.
         """
-        Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric="bottleneck", e=self.epsilon)
+        Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric="bottleneck", e=self.epsilon, n_jobs=self.n_jobs)
         return Xfit
 
     def __call__(self, diag1, diag2):
@@ -282,15 +302,17 @@ class PersistenceFisherDistance(BaseEstimator, TransformerMixin):
     """
     This is a class for computing the persistence Fisher distance matrix from a list of persistence diagrams. The persistence Fisher distance is obtained by computing the original Fisher distance between the probability distributions associated to the persistence diagrams given by convolving them with a Gaussian kernel. See http://papers.nips.cc/paper/8205-persistence-fisher-kernel-a-riemannian-manifold-kernel-for-persistence-diagrams for more details. 
     """
-    def __init__(self, bandwidth=1., kernel_approx=None):
+    def __init__(self, bandwidth=1., kernel_approx=None, n_jobs=None):
         """
         Constructor for the PersistenceFisherDistance class.
 
         Parameters:
             bandwidth (double): bandwidth of the Gaussian kernel used to turn persistence diagrams into probability distributions (default 1.).
             kernel_approx (class): kernel approximation class used to speed up computation (default None). Common kernel approximations classes can be found in the scikit-learn library (such as RBFSampler for instance).   
+            n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_distances` for details.
         """
         self.bandwidth, self.kernel_approx = bandwidth, kernel_approx
+        self.n_jobs = n_jobs
 
     def fit(self, X, y=None):
         """
@@ -313,7 +335,7 @@ class PersistenceFisherDistance(BaseEstimator, TransformerMixin):
         Returns:
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise persistence Fisher distances.
         """
-        return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="persistence_fisher", bandwidth=self.bandwidth, kernel_approx=self.kernel_approx)
+        return pairwise_persistence_diagram_distances(X, self.diagrams_, metric="persistence_fisher", bandwidth=self.bandwidth, kernel_approx=self.kernel_approx, n_jobs=self.n_jobs)
 
     def __call__(self, diag1, diag2):
         """
@@ -333,7 +355,8 @@ class WassersteinDistance(BaseEstimator, TransformerMixin):
     """
     This is a class for computing the Wasserstein distance matrix from a list of persistence diagrams. 
     """
-    def __init__(self, order=1, internal_p=np.inf, mode="hera", delta=0.01):
+
+    def __init__(self, order=1, internal_p=np.inf, mode="hera", delta=0.01, n_jobs=None):
         """
         Constructor for the WassersteinDistance class.
 
@@ -342,6 +365,7 @@ class WassersteinDistance(BaseEstimator, TransformerMixin):
             internal_p (int): ground metric on the (upper-half) plane (i.e. norm l_p in R^2), default value is `np.inf`, see :func:`gudhi.wasserstein.wasserstein_distance`.
             mode (str): method for computing Wasserstein distance. Either "pot" or "hera". Default set to "hera".
             delta (float): relative error 1+delta. Used only if mode == "hera".
+            n_jobs (int): number of jobs to use for the computation. See :func:`pairwise_persistence_diagram_distances` for details.
         """
         self.order, self.internal_p, self.mode = order, internal_p, mode
         if mode == "pot":
@@ -351,6 +375,7 @@ class WassersteinDistance(BaseEstimator, TransformerMixin):
         else:
             raise NameError("Unknown mode. Current available values for mode are 'hera' and 'pot'")
         self.delta = delta
+        self.n_jobs = n_jobs
 
     def fit(self, X, y=None):
         """
@@ -374,9 +399,9 @@ class WassersteinDistance(BaseEstimator, TransformerMixin):
             numpy array of shape (number of diagrams in **diagrams**) x (number of diagrams in X): matrix of pairwise Wasserstein distances.
         """
         if self.metric == "hera_wasserstein":
-            Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric=self.metric, order=self.order, internal_p=self.internal_p, delta=self.delta)
+            Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric=self.metric, order=self.order, internal_p=self.internal_p, delta=self.delta, n_jobs=self.n_jobs)
         else:
-            Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric=self.metric, order=self.order, internal_p=self.internal_p, matching=False)
+            Xfit = pairwise_persistence_diagram_distances(X, self.diagrams_, metric=self.metric, order=self.order, internal_p=self.internal_p, matching=False, n_jobs=self.n_jobs)
         return Xfit
 
     def __call__(self, diag1, diag2):
diff --git a/src/python/include/Alpha_complex_interface.h b/src/python/include/Alpha_complex_interface.h
index 40de88f3..3ac5db1f 100644
--- a/src/python/include/Alpha_complex_interface.h
+++ b/src/python/include/Alpha_complex_interface.h
@@ -23,45 +23,74 @@
 #include <iostream>
 #include <vector>
 #include <string>
+#include <memory>  // for std::unique_ptr
 
 namespace Gudhi {
 
 namespace alpha_complex {
 
 class Alpha_complex_interface {
-  using Dynamic_kernel = CGAL::Epeck_d< CGAL::Dynamic_dimension_tag >;
-  using Point_d = Dynamic_kernel::Point_d;
+ private:
+  using Exact_kernel = CGAL::Epeck_d<CGAL::Dynamic_dimension_tag>;
+  using Inexact_kernel = CGAL::Epick_d<CGAL::Dynamic_dimension_tag>;
+  using Point_exact_kernel = typename Exact_kernel::Point_d;
+  using Point_inexact_kernel = typename Inexact_kernel::Point_d;
 
- public:
-  Alpha_complex_interface(const std::vector<std::vector<double>>& points) {
-    auto mkpt = [](std::vector<double> const& vec){
-      return Point_d(vec.size(), vec.begin(), vec.end());
-    };
-    alpha_complex_ = new Alpha_complex<Dynamic_kernel>(boost::adaptors::transform(points, mkpt));
+  template <typename CgalPointType>
+  std::vector<double> pt_cgal_to_cython(CgalPointType& point) {
+    std::vector<double> vd;
+    for (auto coord = point.cartesian_begin(); coord != point.cartesian_end(); coord++)
+      vd.push_back(CGAL::to_double(*coord));
+    return vd;
   }
 
-  Alpha_complex_interface(const std::string& off_file_name, bool from_file = true) {
-    alpha_complex_ = new Alpha_complex<Dynamic_kernel>(off_file_name);
+  template <typename CgalPointType>
+  static CgalPointType pt_cython_to_cgal(std::vector<double> const& vec) {
+    return CgalPointType(vec.size(), vec.begin(), vec.end());
   }
 
-  ~Alpha_complex_interface() {
-    delete alpha_complex_;
+ public:
+  Alpha_complex_interface(const std::vector<std::vector<double>>& points, bool fast_version)
+      : fast_version_(fast_version) {
+    if (fast_version_) {
+      ac_inexact_ptr_ = std::make_unique<Alpha_complex<Inexact_kernel>>(
+          boost::adaptors::transform(points, pt_cython_to_cgal<Point_inexact_kernel>));
+    } else {
+      ac_exact_ptr_ = std::make_unique<Alpha_complex<Exact_kernel>>(
+          boost::adaptors::transform(points, pt_cython_to_cgal<Point_exact_kernel>));
+    }
+  }
+
+  Alpha_complex_interface(const std::string& off_file_name, bool fast_version, bool from_file = true)
+      : fast_version_(fast_version) {
+    if (fast_version_)
+      ac_inexact_ptr_ = std::make_unique<Alpha_complex<Inexact_kernel>>(off_file_name);
+    else
+      ac_exact_ptr_ = std::make_unique<Alpha_complex<Exact_kernel>>(off_file_name);
   }
 
   std::vector<double> get_point(int vh) {
-    std::vector<double> vd;
-    Point_d const& ph = alpha_complex_->get_point(vh);
-    for (auto coord = ph.cartesian_begin(); coord != ph.cartesian_end(); coord++)
-      vd.push_back(CGAL::to_double(*coord));
-    return vd;
+    if (fast_version_) {
+      Point_inexact_kernel const& point = ac_inexact_ptr_->get_point(vh);
+      return pt_cgal_to_cython(point);
+    } else {
+      Point_exact_kernel const& point = ac_exact_ptr_->get_point(vh);
+      return pt_cgal_to_cython(point);
+    }
   }
 
-  void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square) {
-    alpha_complex_->create_complex(*simplex_tree, max_alpha_square);
+  void create_simplex_tree(Simplex_tree_interface<>* simplex_tree, double max_alpha_square, bool exact_version,
+                           bool default_filtration_value) {
+    if (fast_version_)
+      ac_inexact_ptr_->create_complex(*simplex_tree, max_alpha_square, exact_version, default_filtration_value);
+    else
+      ac_exact_ptr_->create_complex(*simplex_tree, max_alpha_square, exact_version, default_filtration_value);
   }
 
  private:
-  Alpha_complex<Dynamic_kernel>* alpha_complex_;
+  bool fast_version_;
+  std::unique_ptr<Alpha_complex<Exact_kernel>> ac_exact_ptr_;
+  std::unique_ptr<Alpha_complex<Inexact_kernel>> ac_inexact_ptr_;
 };
 
 }  // namespace alpha_complex
diff --git a/src/python/include/pybind11_diagram_utils.h b/src/python/include/pybind11_diagram_utils.h
index d9627258..2d5194f4 100644
--- a/src/python/include/pybind11_diagram_utils.h
+++ b/src/python/include/pybind11_diagram_utils.h
@@ -18,8 +18,8 @@ namespace py = pybind11;
 typedef py::array_t<double> Dgm;
 
 // Get m[i,0] and m[i,1] as a pair
-static auto pairify(void* p, ssize_t h, ssize_t w) {
-  return [=](ssize_t i){
+static auto pairify(void* p, py::ssize_t h, py::ssize_t w) {
+  return [=](py::ssize_t i){
     char* birth = (char*)p + i * h;
     char* death = birth + w;
     return std::make_pair(*(double*)birth, *(double*)death);
@@ -32,8 +32,8 @@ inline auto numpy_to_range_of_pairs(py::array_t<double> dgm) {
   if((buf.ndim!=2 || buf.shape[1]!=2) && (buf.ndim!=1 || buf.shape[0]!=0))
     throw std::runtime_error("Diagram must be an array of size n x 2");
   // In the case of shape (0), avoid reading non-existing strides[1] even if we won't use it.
-  ssize_t stride1 = buf.ndim == 2 ? buf.strides[1] : 0;
-  auto cnt = boost::counting_range<ssize_t>(0, buf.shape[0]);
+  py::ssize_t stride1 = buf.ndim == 2 ? buf.strides[1] : 0;
+  auto cnt = boost::counting_range<py::ssize_t>(0, buf.shape[0]);
   return boost::adaptors::transform(cnt, pairify(buf.ptr, buf.strides[0], stride1));
   // Be careful that the returned range cannot contain references to dead temporaries.
 }
diff --git a/src/python/introduction.rst b/src/python/introduction.rst
new file mode 100644
index 00000000..11c06ac5
--- /dev/null
+++ b/src/python/introduction.rst
@@ -0,0 +1,24 @@
+The Gudhi library is an open source library for Computational Topology and
+Topological Data Analysis (TDA). It offers state-of-the-art algorithms
+to construct various types of simplicial complexes, data structures to
+represent them, and algorithms to compute geometric approximations of shapes
+and persistent homology.
+
+The GUDHI library offers the following interoperable modules:
+
+* Complexes:
+   * Cubical
+   * Simplicial: Rips, Witness, Alpha and Čech complexes
+   * Cover: Nerve and Graph induced complexes
+* Data structures and basic operations:
+   * Simplex tree, Skeleton blockers and Toplex map
+   * Construction, update, filtration and simplification
+* Topological descriptors computation
+* Manifold reconstruction
+* Topological descriptors tools:
+   * Bottleneck distance
+   * Statistical tools
+   * Persistence diagram and barcode
+
+For more information about Topological Data Analysis and its workflow, please
+refer to the `Wikipedia TDA dedicated page <https://en.wikipedia.org/wiki/Topological_data_analysis>`_.
diff --git a/src/python/setup.py.in b/src/python/setup.py.in
index b9f4e3f0..98d058fc 100644
--- a/src/python/setup.py.in
+++ b/src/python/setup.py.in
@@ -48,10 +48,12 @@ ext_modules = cythonize(ext_modules)
 
 for module in pybind11_modules:
     my_include_dirs = include_dirs + [pybind11.get_include(False), pybind11.get_include(True)]
-    if module == 'hera':
+    if module == 'hera/wasserstein':
         my_include_dirs = ['@HERA_WASSERSTEIN_INCLUDE_DIR@'] + my_include_dirs
+    elif module == 'hera/bottleneck':
+        my_include_dirs = ['@HERA_BOTTLENECK_INCLUDE_DIR@'] + my_include_dirs
     ext_modules.append(Extension(
-        'gudhi.' + module,
+        'gudhi.' + module.replace('/', '.'),
         sources = [source_dir + module + '.cc'],
         language = 'c++',
         include_dirs = my_include_dirs,
@@ -62,14 +64,29 @@ for module in pybind11_modules:
         runtime_library_dirs=runtime_library_dirs,
         ))
 
+# read the contents of introduction.rst
+with open("introduction.rst", "r") as fh:
+    long_description = fh.read()
+
 setup(
     name = 'gudhi',
     packages=find_packages(), # find_namespace_packages(include=["gudhi*"])
     author='GUDHI Editorial Board',
     author_email='gudhi-contact@lists.gforge.inria.fr',
     version='@GUDHI_VERSION@',
-    url='http://gudhi.gforge.inria.fr/',
+    url='https://gudhi.inria.fr/',
+    project_urls={
+        'Bug Tracker': 'https://github.com/GUDHI/gudhi-devel/issues',
+        'Documentation': 'https://gudhi.inria.fr/python/latest/',
+        'Source Code': 'https://github.com/GUDHI/gudhi-devel',
+        'License': 'https://gudhi.inria.fr/licensing/'
+        },
+    description='The Gudhi library is an open source library for ' \
+        'Computational Topology and Topological Data Analysis (TDA).',
+    long_description_content_type='text/x-rst',
+    long_description=long_description,
     ext_modules = ext_modules,
-    install_requires = ['cython','numpy >= 1.9',],
-    setup_requires = ['numpy >= 1.9','pybind11',],
+    install_requires = ['numpy >= 1.9',],
+    setup_requires = ['cython','numpy >= 1.9','pybind11',],
+    package_data={"": ["*.dll"], },
 )
diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py
index 77121302..943ad2c4 100755
--- a/src/python/test/test_alpha_complex.py
+++ b/src/python/test/test_alpha_complex.py
@@ -24,14 +24,18 @@ __copyright__ = "Copyright (C) 2016 Inria"
 __license__ = "MIT"
 
 
-def test_empty_alpha():
-    alpha_complex = AlphaComplex(points=[[0, 0]])
+def _empty_alpha(precision):
+    alpha_complex = AlphaComplex(points=[[0, 0]], precision = precision)
     assert alpha_complex.__is_defined() == True
 
+def test_empty_alpha():
+    _empty_alpha('fast')
+    _empty_alpha('safe')
+    _empty_alpha('exact')
 
-def test_infinite_alpha():
+def _infinite_alpha(precision):
     point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
-    alpha_complex = AlphaComplex(points=point_list)
+    alpha_complex = AlphaComplex(points=point_list, precision = precision)
     assert alpha_complex.__is_defined() == True
 
     simplex_tree = alpha_complex.create_simplex_tree()
@@ -79,10 +83,14 @@ def test_infinite_alpha():
     else:
         assert False
 
+def test_infinite_alpha():
+    _infinite_alpha('fast')
+    _infinite_alpha('safe')
+    _infinite_alpha('exact')
 
-def test_filtered_alpha():
+def _filtered_alpha(precision):
     point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
-    filtered_alpha = AlphaComplex(points=point_list)
+    filtered_alpha = AlphaComplex(points=point_list, precision = precision)
 
     simplex_tree = filtered_alpha.create_simplex_tree(max_alpha_square=0.25)
 
@@ -119,7 +127,12 @@ def test_filtered_alpha():
     assert simplex_tree.get_star([0]) == [([0], 0.0), ([0, 1], 0.25), ([0, 2], 0.25)]
     assert simplex_tree.get_cofaces([0], 1) == [([0, 1], 0.25), ([0, 2], 0.25)]
 
-def test_safe_alpha_persistence_comparison():
+def test_filtered_alpha():
+    _filtered_alpha('fast')
+    _filtered_alpha('safe')
+    _filtered_alpha('exact')
+
+def _safe_alpha_persistence_comparison(precision):
     #generate periodic signal
     time = np.arange(0, 10, 1)
     signal = [math.sin(x) for x in time]
@@ -131,10 +144,10 @@ def test_safe_alpha_persistence_comparison():
     embedding2 = [[signal[i], delayed[i]] for i in range(len(time))]
     
     #build alpha complex and simplex tree
-    alpha_complex1 = AlphaComplex(points=embedding1)
+    alpha_complex1 = AlphaComplex(points=embedding1, precision = precision)
     simplex_tree1 = alpha_complex1.create_simplex_tree()
     
-    alpha_complex2 = AlphaComplex(points=embedding2)
+    alpha_complex2 = AlphaComplex(points=embedding2, precision = precision)
     simplex_tree2 = alpha_complex2.create_simplex_tree()
     
     diag1 = simplex_tree1.persistence()
@@ -143,3 +156,47 @@ def test_safe_alpha_persistence_comparison():
     for (first_p, second_p) in zip_longest(diag1, diag2):
         assert first_p[0] == pytest.approx(second_p[0])
         assert first_p[1] == pytest.approx(second_p[1])
+
+
+def test_safe_alpha_persistence_comparison():
+    # Won't work for 'fast' version
+    _safe_alpha_persistence_comparison('safe')
+    _safe_alpha_persistence_comparison('exact')
+
+def _delaunay_complex(precision):
+    point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
+    filtered_alpha = AlphaComplex(points=point_list, precision = precision)
+
+    simplex_tree = filtered_alpha.create_simplex_tree(default_filtration_value = True)
+
+    assert simplex_tree.num_simplices() == 11
+    assert simplex_tree.num_vertices() == 4
+
+    assert point_list[0] == filtered_alpha.get_point(0)
+    assert point_list[1] == filtered_alpha.get_point(1)
+    assert point_list[2] == filtered_alpha.get_point(2)
+    assert point_list[3] == filtered_alpha.get_point(3)
+    try:
+        filtered_alpha.get_point(4) == []
+    except IndexError:
+        pass
+    else:
+        assert False
+    try:
+        filtered_alpha.get_point(125) == []
+    except IndexError:
+        pass
+    else:
+        assert False
+
+    for filtered_value in simplex_tree.get_filtration():
+        assert math.isnan(filtered_value[1])
+    for filtered_value in simplex_tree.get_star([0]):
+        assert math.isnan(filtered_value[1])
+    for filtered_value in simplex_tree.get_cofaces([0], 1):
+        assert math.isnan(filtered_value[1])
+
+def test_delaunay_complex():
+    _delaunay_complex('fast')
+    _delaunay_complex('safe')
+    _delaunay_complex('exact')
diff --git a/src/python/test/test_bottleneck_distance.py b/src/python/test/test_bottleneck_distance.py
index 70b2abad..6915bea8 100755
--- a/src/python/test/test_bottleneck_distance.py
+++ b/src/python/test/test_bottleneck_distance.py
@@ -9,6 +9,8 @@
 """
 
 import gudhi
+import gudhi.hera
+import pytest
 
 __author__ = "Vincent Rouvreau"
 __copyright__ = "Copyright (C) 2016 Inria"
@@ -19,5 +21,7 @@ def test_basic_bottleneck():
     diag1 = [[2.7, 3.7], [9.6, 14.0], [34.2, 34.974], [3.0, float("Inf")]]
     diag2 = [[2.8, 4.45], [9.5, 14.1], [3.2, float("Inf")]]
 
-    assert gudhi.bottleneck_distance(diag1, diag2, 0.1) == 0.8081763781405569
     assert gudhi.bottleneck_distance(diag1, diag2) == 0.75
+    assert gudhi.bottleneck_distance(diag1, diag2, 0.1) == pytest.approx(0.75, abs=0.1)
+    assert gudhi.hera.bottleneck_distance(diag1, diag2, 0) == 0.75
+    assert gudhi.hera.bottleneck_distance(diag1, diag2, 0.1) == pytest.approx(0.75, rel=0.1)
diff --git a/src/python/test/test_dtm_rips_complex.py b/src/python/test/test_dtm_rips_complex.py
new file mode 100644
index 00000000..e1c0ee44
--- /dev/null
+++ b/src/python/test/test_dtm_rips_complex.py
@@ -0,0 +1,32 @@
+""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+    Author(s):       Yuichi Ike
+
+    Copyright (C) 2020 Inria, Copyright (C) 2020 FUjitsu Laboratories Ltd.
+
+    Modification(s):
+      - YYYY/MM Author: Description of the modification
+"""
+
+from gudhi.dtm_rips_complex import DTMRipsComplex
+from gudhi import RipsComplex
+import numpy as np
+from math import sqrt
+import pytest
+
+def test_dtm_rips_complex():
+    pts = np.array([[2.0, 2.0], [0.0, 1.0], [3.0, 4.0]])
+    dtm_rips = DTMRipsComplex(points=pts, k=2)
+    st = dtm_rips.create_simplex_tree(max_dimension=2)
+    st.persistence()
+    persistence_intervals0 = st.persistence_intervals_in_dimension(0)
+    assert persistence_intervals0 == pytest.approx(np.array([[3.16227766, 5.39834564],[3.16227766, 5.39834564], [3.16227766, float("inf")]]))
+    
+def test_compatibility_with_rips():
+    distance_matrix = np.array([[0, 1, 1, sqrt(2)], [1, 0, sqrt(2), 1], [1, sqrt(2), 0, 1], [sqrt(2), 1, 1, 0]])
+    dtm_rips = DTMRipsComplex(distance_matrix=distance_matrix, max_filtration=42)
+    st = dtm_rips.create_simplex_tree(max_dimension=1)
+    rips_complex = RipsComplex(distance_matrix=distance_matrix, max_edge_length=42)
+    st_from_rips = rips_complex.create_simplex_tree(max_dimension=1)
+    assert list(st.get_filtration()) == list(st_from_rips.get_filtration())
+
diff --git a/src/python/test/test_representations.py b/src/python/test/test_representations.py
index dba7f952..589cee00 100755
--- a/src/python/test/test_representations.py
+++ b/src/python/test/test_representations.py
@@ -1,12 +1,43 @@
 import os
 import sys
 import matplotlib.pyplot as plt
+import numpy as np
+import pytest
+
 
 def test_representations_examples():
     # Disable graphics for testing purposes
-    plt.show = lambda:None
+    plt.show = lambda: None
     here = os.path.dirname(os.path.realpath(__file__))
     sys.path.append(here + "/../example")
     import diagram_vectorizations_distances_kernels
 
     return None
+
+
+from gudhi.representations.metrics import *
+from gudhi.representations.kernel_methods import *
+
+
+def _n_diags(n):
+    l = []
+    for _ in range(n):
+        a = np.random.rand(50, 2)
+        a[:, 1] += a[:, 0]  # So that y >= x
+        l.append(a)
+    return l
+
+
+def test_multiple():
+    l1 = _n_diags(9)
+    l2 = _n_diags(11)
+    l1b = l1.copy()
+    d1 = pairwise_persistence_diagram_distances(l1, e=0.00001, n_jobs=4)
+    d2 = BottleneckDistance(epsilon=0.00001).fit_transform(l1)
+    d3 = pairwise_persistence_diagram_distances(l1, l1b, e=0.00001, n_jobs=4)
+    assert d1 == pytest.approx(d2)
+    assert d3 == pytest.approx(d2, abs=1e-5) # Because of 0 entries (on the diagonal)
+    d1 = pairwise_persistence_diagram_distances(l1, l2, metric="wasserstein", order=2, internal_p=2)
+    d2 = WassersteinDistance(order=2, internal_p=2, n_jobs=4).fit(l2).transform(l1)
+    print(d1.shape, d2.shape)
+    assert d1 == pytest.approx(d2, rel=.02)
diff --git a/src/python/test/test_tomato.py b/src/python/test/test_tomato.py
new file mode 100755
index 00000000..ecab03c4
--- /dev/null
+++ b/src/python/test/test_tomato.py
@@ -0,0 +1,65 @@
+""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+    Author(s):       Marc Glisse
+
+    Copyright (C) 2020 Inria
+
+    Modification(s):
+      - YYYY/MM Author: Description of the modification
+"""
+
+from gudhi.clustering.tomato import Tomato
+import numpy as np
+import pytest
+import matplotlib.pyplot as plt
+
+# Disable graphics for testing purposes
+plt.show = lambda: None
+
+
+def test_tomato_1():
+    a = [(1, 2), (1.1, 1.9), (0.9, 1.8), (10, 0), (10.1, 0.05), (10.2, -0.1), (5.4, 0)]
+    t = Tomato(metric="euclidean", n_clusters=2, k=4, n_jobs=-1, eps=0.05)
+    assert np.array_equal(t.fit_predict(a), [1, 1, 1, 0, 0, 0, 0])  # or with swapped 0 and 1
+    assert np.array_equal(t.children_, [[0, 1]])
+
+    t = Tomato(density_type="KDE", r=1, k=4)
+    t.fit(a)
+    assert np.array_equal(t.leaf_labels_, [1, 1, 1, 0, 0, 0, 0])  # or with swapped 0 and 1
+    assert t.n_clusters_ == 2
+    t.merge_threshold_ = 10
+    assert t.n_clusters_ == 1
+    assert (t.labels_ == 0).all()
+
+    t = Tomato(graph_type="radius", r=0.1, metric="cosine", k=3)
+    assert np.array_equal(t.fit_predict(a), [1, 1, 1, 0, 0, 0, 0])  # or with swapped 0 and 1
+
+    t = Tomato(metric="euclidean", graph_type="radius", r=4.7, k=4)
+    t.fit(a)
+    assert t.max_weight_per_cc_.size == 2
+    assert np.array_equal(t.neighbors_, [[0, 1, 2], [0, 1, 2], [0, 1, 2], [3, 4, 5, 6], [3, 4, 5], [3, 4, 5], [3, 6]])
+    t.plot_diagram()
+
+    t = Tomato(graph_type="radius", r=4.7, k=4, symmetrize_graph=True)
+    t.fit(a)
+    assert t.max_weight_per_cc_.size == 2
+    assert [set(i) for i in t.neighbors_] == [{1, 2}, {0, 2}, {0, 1}, {4, 5, 6}, {3, 5}, {3, 4}, {3}]
+
+    t = Tomato(n_clusters=2, k=4, symmetrize_graph=True)
+    t.fit(a)
+    assert [set(i) for i in t.neighbors_] == [
+        {1, 2, 6},
+        {0, 2, 6},
+        {0, 1, 6},
+        {4, 5, 6},
+        {3, 5, 6},
+        {3, 4, 6},
+        {0, 1, 2, 3, 4, 5},
+    ]
+    t.plot_diagram()
+
+    t = Tomato(k=6, metric="manhattan")
+    t.fit(a)
+    assert t.diagram_.size == 0
+    assert t.max_weight_per_cc_.size == 1
+    t.plot_diagram()
diff --git a/src/python/test/test_weighted_rips.py b/src/python/test/test_weighted_rips_complex.py
index 7ef48333..7ef48333 100644
--- a/src/python/test/test_weighted_rips.py
+++ b/src/python/test/test_weighted_rips_complex.py