Merge tag 'tags/gudhi-release-3.6.0' into dfsg/latest

17 files changed, 1112 insertions, 278 deletions

diff --git a/src/python/gudhi/__init__.py.in b/src/python/gudhi/__init__.py.in
index 3043201a..79e12fbc 100644
--- a/src/python/gudhi/__init__.py.in
+++ b/src/python/gudhi/__init__.py.in
@@ -23,10 +23,6 @@ __all__ = [@GUDHI_PYTHON_MODULES@ @GUDHI_PYTHON_MODULES_EXTRA@]
 __available_modules = ''
 __missing_modules = ''
 
-# For unitary tests purpose
-# could use "if 'collapse_edges' in gudhi.__all__" when collapse edges will have a python module
-__GUDHI_USE_EIGEN3 = @GUDHI_USE_EIGEN3@
-
 # Try to import * from gudhi.__module_name for default modules.
 # Extra modules require an explicit import by the user (mostly because of
 # unusual dependencies, but also to avoid cluttering namespace gudhi and
diff --git a/src/python/gudhi/alpha_complex.pyx b/src/python/gudhi/alpha_complex.pyx
index ea128743..375e1561 100644
--- a/src/python/gudhi/alpha_complex.pyx
+++ b/src/python/gudhi/alpha_complex.pyx
@@ -16,7 +16,7 @@ from libcpp.utility cimport pair
 from libcpp.string cimport string
 from libcpp cimport bool
 from libc.stdint cimport intptr_t
-import os
+import warnings
 
 from gudhi.simplex_tree cimport *
 from gudhi.simplex_tree import SimplexTree
@@ -28,66 +28,76 @@ __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, bool fast_version, bool exact_version) nogil except +
+        Alpha_complex_interface(vector[vector[double]] points, vector[double] weights, bool fast_version, bool exact_version) 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, bool default_filtration_value) nogil except +
+        @staticmethod
+        void set_float_relative_precision(double precision) nogil
+        @staticmethod
+        double get_float_relative_precision() nogil
 
 # AlphaComplex python interface
 cdef class AlphaComplex:
-    """AlphaComplex is a simplicial complex constructed from the finite cells
-    of a Delaunay Triangulation.
+    """AlphaComplex is a simplicial complex constructed from the finite cells of a Delaunay Triangulation.
 
-    The filtration value of each simplex is computed as the square of the
-    circumradius of the simplex if the circumsphere is empty (the simplex is
-    then said to be Gabriel), and as the minimum of the filtration values of
-    the codimension 1 cofaces that make it not Gabriel otherwise.
+    The filtration value of each simplex is computed as the square of the circumradius of the simplex if the
+    circumsphere is empty (the simplex is then said to be Gabriel), and as the minimum of the filtration values of the
+    codimension 1 cofaces that make it not Gabriel otherwise.
 
-    All simplices that have a filtration value strictly greater than a given
-    alpha squared value are not inserted into the complex.
+    All simplices that have a filtration value strictly greater than a given alpha squared value are not inserted into
+    the complex.
 
     .. note::
 
-        When Alpha_complex is constructed with an infinite value of alpha, the
-        complex is a Delaunay complex.
-
+        When Alpha_complex is constructed with an infinite value of alpha, the complex is a Delaunay complex.
     """
 
     cdef Alpha_complex_interface * this_ptr
 
     # Fake constructor that does nothing but documenting the constructor
-    def __init__(self, points=None, off_file='', precision='safe'):
+    def __init__(self, points=[], off_file='', weights=None, precision='safe'):
         """AlphaComplex constructor.
 
         :param points: A list of points in d-Dimension.
-        :type points: list of list of double
-
-        Or
+        :type points: Iterable[Iterable[float]]
 
-        :param off_file: An OFF file style name.
+        :param off_file: **[deprecated]** An `OFF file style <fileformats.html#off-file-format>`_ name.
+            If an `off_file` is given with `points` as arguments, only points from the file are taken into account.
         :type off_file: string
 
+        :param weights: A list of weights. If set, the number of weights must correspond to the number of points.
+        :type weights: Iterable[float]
+
         :param precision: Alpha complex precision can be 'fast', 'safe' or 'exact'. Default is 'safe'.
         :type precision: string
+
+        :raises FileNotFoundError: **[deprecated]** If `off_file` is set but not found.
+        :raises ValueError: In case of inconsistency between the number of points and weights.
         """
 
     # The real cython constructor
-    def __cinit__(self, points = None, off_file = '', precision = 'safe'):
+    def __cinit__(self, points = [], off_file = '', weights=None, precision = 'safe'):
         assert precision in ['fast', 'safe', 'exact'], "Alpha complex precision can only be 'fast', 'safe' or 'exact'"
         cdef bool fast = precision == 'fast'
         cdef bool exact = precision == 'exact'
 
-        cdef vector[vector[double]] pts
         if off_file:
-            if os.path.isfile(off_file):
-                points = read_points_from_off_file(off_file = off_file)
-            else:
-                print("file " + off_file + " not found.")
-        if points is None:
-            # Empty Alpha construction
-            points=[]
+            warnings.warn("off_file is a deprecated parameter, please consider using gudhi.read_points_from_off_file",
+                          DeprecationWarning)
+            points = read_points_from_off_file(off_file = off_file)
+
+        # weights are set but is inconsistent with the number of points
+        if weights != None and len(weights) != len(points):
+            raise ValueError("Inconsistency between the number of points and weights")
+
+        # need to copy the points to use them without the gil
+        cdef vector[vector[double]] pts
+        cdef vector[double] wgts
         pts = points
+        if weights != None:
+            wgts = weights
         with nogil:
-            self.this_ptr = new Alpha_complex_interface(pts, fast, exact)
+            self.this_ptr = new Alpha_complex_interface(pts, wgts, fast, exact)
 
     def __dealloc__(self):
         if self.this_ptr != NULL:
@@ -127,3 +137,28 @@ cdef class AlphaComplex:
             self.this_ptr.create_simplex_tree(<Simplex_tree_interface_full_featured*>stree_int_ptr,
                                               mas, compute_filtration)
         return stree
+
+    @staticmethod
+    def set_float_relative_precision(precision):
+        """
+        :param precision: When the AlphaComplex is constructed with :code:`precision = 'safe'` (the default),
+            one can set the float relative precision of filtration values computed in
+            :func:`~gudhi.AlphaComplex.create_simplex_tree`.
+            Default is :code:`1e-5` (cf. :func:`~gudhi.AlphaComplex.get_float_relative_precision`).
+            For more details, please refer to
+            `CGAL::Lazy_exact_nt<NT>::set_relative_precision_of_to_double <https://doc.cgal.org/latest/Number_types/classCGAL_1_1Lazy__exact__nt.html>`_
+        :type precision: float
+        """
+        if precision <=0. or precision >= 1.:
+            raise ValueError("Relative precision value must be strictly greater than 0 and strictly lower than 1")
+        Alpha_complex_interface.set_float_relative_precision(precision)
+    
+    @staticmethod
+    def get_float_relative_precision():
+        """
+        :returns: The float relative precision of filtration values computation in
+            :func:`~gudhi.AlphaComplex.create_simplex_tree` when the AlphaComplex is constructed with
+            :code:`precision = 'safe'` (the default).
+        :rtype: float
+        """
+        return Alpha_complex_interface.get_float_relative_precision()
diff --git a/src/python/gudhi/datasets/remote.py b/src/python/gudhi/datasets/remote.py
new file mode 100644
index 00000000..f6d3fe56
--- /dev/null
+++ b/src/python/gudhi/datasets/remote.py
@@ -0,0 +1,223 @@
+# 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):       Hind Montassif
+#
+# Copyright (C) 2021 Inria
+#
+# Modification(s):
+#   - YYYY/MM Author: Description of the modification
+
+from os.path import join, split, exists, expanduser
+from os import makedirs, remove, environ
+
+from urllib.request import urlretrieve
+import hashlib
+import shutil
+
+import numpy as np
+
+def _get_data_home(data_home = None):
+    """
+    Return the path of the remote datasets directory.
+    This folder is used to store remotely fetched datasets.
+    By default the datasets directory is set to a folder named 'gudhi_data' in the user home folder.
+    Alternatively, it can be set by the 'GUDHI_DATA' environment variable.
+    The '~' symbol is expanded to the user home folder.
+    If the folder does not already exist, it is automatically created.
+
+    Parameters
+    ----------
+    data_home : string
+        The path to remote datasets directory.
+        Default is `None`, meaning that the data home directory will be set to "~/gudhi_data",
+        if the 'GUDHI_DATA' environment variable does not exist.
+
+    Returns
+    -------
+    data_home: string
+        The path to remote datasets directory.
+    """
+    if data_home is None:
+        data_home = environ.get("GUDHI_DATA", join("~", "gudhi_data"))
+    data_home = expanduser(data_home)
+    makedirs(data_home, exist_ok=True)
+    return data_home
+
+
+def clear_data_home(data_home = None):
+    """
+    Delete the data home cache directory and all its content.
+
+    Parameters
+    ----------
+    data_home : string, default is None.
+        The path to remote datasets directory.
+        If `None` and the 'GUDHI_DATA' environment variable does not exist,
+        the default directory to be removed is set to "~/gudhi_data".
+    """
+    data_home = _get_data_home(data_home)
+    shutil.rmtree(data_home)
+
+def _checksum_sha256(file_path):
+    """
+    Compute the file checksum using sha256.
+
+    Parameters
+    ----------
+    file_path: string
+        Full path of the created file including filename.
+
+    Returns
+    -------
+        The hex digest of file_path.
+    """
+    sha256_hash = hashlib.sha256()
+    chunk_size = 4096
+    with open(file_path,"rb") as f:
+        # Read and update hash string value in blocks of 4K
+        while True:
+            buffer = f.read(chunk_size)
+            if not buffer:
+                break
+            sha256_hash.update(buffer)
+    return sha256_hash.hexdigest()
+
+def _fetch_remote(url, file_path, file_checksum = None):
+    """
+    Fetch the wanted dataset from the given url and save it in file_path.
+
+    Parameters
+    ----------
+    url : string
+        The url to fetch the dataset from.
+    file_path : string
+        Full path of the downloaded file including filename.
+    file_checksum : string
+        The file checksum using sha256 to check against the one computed on the downloaded file.
+        Default is 'None', which means the checksum is not checked.
+
+    Raises
+    ------
+    IOError
+        If the computed SHA256 checksum of file does not match the one given by the user.
+    """
+
+    # Get the file
+    urlretrieve(url, file_path)
+
+    if file_checksum is not None:
+        checksum = _checksum_sha256(file_path)
+        if file_checksum != checksum:
+            # Remove file and raise error
+            remove(file_path)
+            raise IOError("{} has a SHA256 checksum : {}, "
+                        "different from expected : {}."
+                        "The file may be corrupted or the given url may be wrong !".format(file_path, checksum, file_checksum))
+
+def _get_archive_path(file_path, label):
+    """
+    Get archive path based on file_path given by user and label.
+
+    Parameters
+    ----------
+    file_path: string
+        Full path of the file to get including filename, or None.
+    label: string
+        Label used along with 'data_home' to get archive path, in case 'file_path' is None.
+
+    Returns
+    -------
+        Full path of archive including filename.
+    """
+    if file_path is None:
+        archive_path = join(_get_data_home(), label)
+        dirname = split(archive_path)[0]
+        makedirs(dirname, exist_ok=True)
+    else:
+        archive_path = file_path
+        dirname = split(archive_path)[0]
+        makedirs(dirname, exist_ok=True)
+
+    return archive_path
+
+def fetch_spiral_2d(file_path = None):
+    """
+    Load the spiral_2d dataset.
+
+    Note that if the dataset already exists in the target location, it is not downloaded again,
+    and the corresponding array is returned from cache.
+
+    Parameters
+    ----------
+    file_path : string
+        Full path of the downloaded file including filename.
+
+        Default is None, meaning that it's set to "data_home/points/spiral_2d/spiral_2d.npy".
+
+        The "data_home" directory is set by default to "~/gudhi_data",
+        unless the 'GUDHI_DATA' environment variable is set.
+
+    Returns
+    -------
+    points: numpy array
+        Array of shape (114562, 2).
+    """
+    file_url = "https://raw.githubusercontent.com/GUDHI/gudhi-data/main/points/spiral_2d/spiral_2d.npy"
+    file_checksum = '2226024da76c073dd2f24b884baefbfd14928b52296df41ad2d9b9dc170f2401'
+
+    archive_path = _get_archive_path(file_path, "points/spiral_2d/spiral_2d.npy")
+
+    if not exists(archive_path):
+        _fetch_remote(file_url, archive_path, file_checksum)
+
+    return np.load(archive_path, mmap_mode='r')
+
+def fetch_bunny(file_path = None, accept_license = False):
+    """
+    Load the Stanford bunny dataset.
+
+    This dataset contains 35947 vertices.
+
+    Note that if the dataset already exists in the target location, it is not downloaded again,
+    and the corresponding array is returned from cache.
+
+    Parameters
+    ----------
+    file_path : string
+        Full path of the downloaded file including filename.
+
+        Default is None, meaning that it's set to "data_home/points/bunny/bunny.npy".
+        In this case, the LICENSE file would be downloaded as "data_home/points/bunny/bunny.LICENSE".
+
+        The "data_home" directory is set by default to "~/gudhi_data",
+        unless the 'GUDHI_DATA' environment variable is set.
+
+    accept_license : boolean
+        Flag to specify if user accepts the file LICENSE and prevents from printing the corresponding license terms.
+
+        Default is False.
+
+    Returns
+    -------
+    points: numpy array
+        Array of shape (35947, 3).
+    """
+
+    file_url = "https://raw.githubusercontent.com/GUDHI/gudhi-data/main/points/bunny/bunny.npy"
+    file_checksum = 'f382482fd89df8d6444152dc8fd454444fe597581b193fd139725a85af4a6c6e'
+    license_url = "https://raw.githubusercontent.com/GUDHI/gudhi-data/main/points/bunny/bunny.LICENSE"
+    license_checksum = 'b763dbe1b2fc6015d05cbf7bcc686412a2eb100a1f2220296e3b4a644c69633a'
+
+    archive_path = _get_archive_path(file_path, "points/bunny/bunny.npy")
+
+    if not exists(archive_path):
+        _fetch_remote(file_url, archive_path, file_checksum)
+        license_path = join(split(archive_path)[0], "bunny.LICENSE")
+        _fetch_remote(license_url, license_path, license_checksum)
+        # Print license terms unless accept_license is set to True
+        if not accept_license:
+            if exists(license_path):
+                with open(license_path, 'r') as f:
+                    print(f.read())
+
+    return np.load(archive_path, mmap_mode='r')
diff --git a/src/python/gudhi/hera/wasserstein.cc b/src/python/gudhi/hera/wasserstein.cc
index 1a21f02f..fa0cf8aa 100644
--- a/src/python/gudhi/hera/wasserstein.cc
+++ b/src/python/gudhi/hera/wasserstein.cc
@@ -29,7 +29,7 @@ double wasserstein_distance(
   if(std::isinf(internal_p)) internal_p = hera::get_infinity<double>();
   params.internal_p = internal_p;
   params.delta = delta;
-  // The extra parameters are purposedly not exposed for now.
+  // The extra parameters are purposely not exposed for now.
   return hera::wasserstein_dist(diag1, diag2, params);
 }
 
diff --git a/src/python/gudhi/persistence_graphical_tools.py b/src/python/gudhi/persistence_graphical_tools.py
index 848dc03e..21275cdd 100644
--- a/src/python/gudhi/persistence_graphical_tools.py
+++ b/src/python/gudhi/persistence_graphical_tools.py
@@ -12,6 +12,9 @@ from os import path
 from math import isfinite
 import numpy as np
 from functools import lru_cache
+import warnings
+import errno
+import os
 
 from gudhi.reader_utils import read_persistence_intervals_in_dimension
 from gudhi.reader_utils import read_persistence_intervals_grouped_by_dimension
@@ -22,6 +25,7 @@ __license__ = "MIT"
 
 _gudhi_matplotlib_use_tex = True
 
+
 def __min_birth_max_death(persistence, band=0.0):
     """This function returns (min_birth, max_death) from the persistence.
 
@@ -44,20 +48,46 @@ def __min_birth_max_death(persistence, band=0.0):
             min_birth = float(interval[1][0])
     if band > 0.0:
         max_death += band
+    # can happen if only points at inf death
+    if min_birth == max_death:
+        max_death = max_death + 1.0
     return (min_birth, max_death)
 
 
 def _array_handler(a):
-    '''
+    """
     :param a: if array, assumes it is a (n x 2) np.array and return a
                 persistence-compatible list (padding with 0), so that the
                 plot can be performed seamlessly.
-    '''
-    if isinstance(a[0][1], np.float64) or isinstance(a[0][1], float):
+    """
+    if isinstance(a[0][1], (np.floating, float)):
         return [[0, x] for x in a]
     else:
         return a
 
+
+def _limit_to_max_intervals(persistence, max_intervals, key):
+    """This function returns truncated persistence if length is bigger than max_intervals.
+    :param persistence: Persistence intervals values list. Can be grouped by dimension or not.
+    :type persistence: an array of (dimension, array of (birth, death)) or an array of (birth, death).
+    :param max_intervals: maximal number of intervals to display.
+        Selected intervals are those with the longest life time. Set it
+        to 0 to see all. Default value is 1000.
+    :type max_intervals: int.
+    :param key: key function for sort algorithm.
+    :type key: function or lambda.
+    """
+    if max_intervals > 0 and max_intervals < len(persistence):
+        warnings.warn(
+            "There are %s intervals given as input, whereas max_intervals is set to %s."
+            % (len(persistence), max_intervals)
+        )
+        # Sort by life time, then takes only the max_intervals elements
+        return sorted(persistence, key=key, reverse=True)[:max_intervals]
+    else:
+        return persistence
+
+
 @lru_cache(maxsize=1)
 def _matplotlib_can_use_tex():
     """This function returns True if matplotlib can deal with LaTeX, False otherwise.
@@ -65,17 +95,17 @@ def _matplotlib_can_use_tex():
     """
     try:
         from matplotlib import checkdep_usetex
+
         return checkdep_usetex(True)
-    except ImportError:
-        print("This function is not available, you may be missing matplotlib.")
+    except ImportError as import_error:
+        warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.")
 
 
 def plot_persistence_barcode(
     persistence=[],
     persistence_file="",
     alpha=0.6,
-    max_intervals=1000,
-    max_barcodes=1000,
+    max_intervals=20000,
     inf_delta=0.1,
     legend=False,
     colormap=None,
@@ -97,7 +127,7 @@ def plot_persistence_barcode(
     :type alpha: float.
     :param max_intervals: maximal number of intervals to display.
         Selected intervals are those with the longest life time. Set it
-        to 0 to see all. Default value is 1000.
+        to 0 to see all. Default value is 20000.
     :type max_intervals: int.
     :param inf_delta: Infinity is placed at :code:`((max_death - min_birth) x
         inf_delta)` above :code:`max_death` value. A reasonable value is
@@ -119,99 +149,68 @@ def plot_persistence_barcode(
         import matplotlib.pyplot as plt
         import matplotlib.patches as mpatches
         from matplotlib import rc
+
         if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex():
-            plt.rc('text', usetex=True)
-            plt.rc('font', family='serif')
+            plt.rc("text", usetex=True)
+            plt.rc("font", family="serif")
         else:
-            plt.rc('text', usetex=False)
-            plt.rc('font', family='DejaVu Sans')
+            plt.rc("text", usetex=False)
+            plt.rc("font", family="DejaVu Sans")
 
         if persistence_file != "":
             if path.isfile(persistence_file):
                 # Reset persistence
                 persistence = []
-                diag = read_persistence_intervals_grouped_by_dimension(
-                    persistence_file=persistence_file
-                )
+                diag = read_persistence_intervals_grouped_by_dimension(persistence_file=persistence_file)
                 for key in diag.keys():
                     for persistence_interval in diag[key]:
                         persistence.append((key, persistence_interval))
             else:
-                print("file " + persistence_file + " not found.")
-                return None
-
-        persistence = _array_handler(persistence)
-
-        if max_barcodes != 1000:
-            print("Deprecated parameter. It has been replaced by max_intervals")
-            max_intervals = max_barcodes
-
-        if max_intervals > 0 and max_intervals < len(persistence):
-            # Sort by life time, then takes only the max_intervals elements
-            persistence = sorted(
-                persistence,
-                key=lambda life_time: life_time[1][1] - life_time[1][0],
-                reverse=True,
-            )[:max_intervals]
-        
-        if colormap == None:
-            colormap = plt.cm.Set1.colors
-        if axes == None:
-            fig, axes = plt.subplots(1, 1)
+                raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), persistence_file)
 
-        persistence = sorted(persistence, key=lambda birth: birth[1][0])
+        try:
+            persistence = _array_handler(persistence)
+            persistence = _limit_to_max_intervals(
+                persistence, max_intervals, key=lambda life_time: life_time[1][1] - life_time[1][0]
+            )
+            (min_birth, max_death) = __min_birth_max_death(persistence)
+            persistence = sorted(persistence, key=lambda birth: birth[1][0])
+        except IndexError:
+            min_birth, max_death = 0.0, 1.0
+            pass
 
-        (min_birth, max_death) = __min_birth_max_death(persistence)
-        ind = 0
         delta = (max_death - min_birth) * inf_delta
         # Replace infinity values with max_death + delta for bar code to be more
         # readable
         infinity = max_death + delta
         axis_start = min_birth - delta
-        # Draw horizontal bars in loop
-        for interval in reversed(persistence):
-            if float(interval[1][1]) != float("inf"):
-                # Finite death case
-                axes.barh(
-                    ind,
-                    (interval[1][1] - interval[1][0]),
-                    height=0.8,
-                    left=interval[1][0],
-                    alpha=alpha,
-                    color=colormap[interval[0]],
-                    linewidth=0,
-                )
-            else:
-                # Infinite death case for diagram to be nicer
-                axes.barh(
-                    ind,
-                    (infinity - interval[1][0]),
-                    height=0.8,
-                    left=interval[1][0],
-                    alpha=alpha,
-                    color=colormap[interval[0]],
-                    linewidth=0,
-                )
-            ind = ind + 1
+
+        if axes == None:
+            _, axes = plt.subplots(1, 1)
+        if colormap == None:
+            colormap = plt.cm.Set1.colors
+
+        x=[birth for (dim,(birth,death)) in persistence]
+        y=[(death - birth) if death != float("inf") else (infinity - birth) for (dim,(birth,death)) in persistence]
+        c=[colormap[dim] for (dim,(birth,death)) in persistence]
+
+        axes.barh(list(reversed(range(len(x)))), y, height=0.8, left=x, alpha=alpha, color=c, linewidth=0)
 
         if legend:
             dimensions = list(set(item[0] for item in persistence))
             axes.legend(
-                handles=[
-                    mpatches.Patch(color=colormap[dim], label=str(dim))
-                    for dim in dimensions
-                ],
-                loc="lower right",
+                handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions], loc="lower right",
             )
 
         axes.set_title("Persistence barcode", fontsize=fontsize)
 
         # Ends plot on infinity value and starts a little bit before min_birth
-        axes.axis([axis_start, infinity, 0, ind])
+        if len(x) != 0:
+            axes.axis([axis_start, infinity, 0, len(x)])
         return axes
 
-    except ImportError:
-        print("This function is not available, you may be missing matplotlib.")
+    except ImportError as import_error:
+        warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.")
 
 
 def plot_persistence_diagram(
@@ -219,14 +218,13 @@ def plot_persistence_diagram(
     persistence_file="",
     alpha=0.6,
     band=0.0,
-    max_intervals=1000,
-    max_plots=1000,
+    max_intervals=1000000,
     inf_delta=0.1,
     legend=False,
     colormap=None,
     axes=None,
     fontsize=16,
-    greyblock=True
+    greyblock=True,
 ):
     """This function plots the persistence diagram from persistence values
     list, a np.array of shape (N x 2) representing a diagram in a single
@@ -244,7 +242,7 @@ def plot_persistence_diagram(
     :type band: float.
     :param max_intervals: maximal number of intervals to display.
         Selected intervals are those with the longest life time. Set it
-        to 0 to see all. Default value is 1000.
+        to 0 to see all. Default value is 1000000.
     :type max_intervals: int.
     :param inf_delta: Infinity is placed at :code:`((max_death - min_birth) x
         inf_delta)` above :code:`max_death` value. A reasonable value is
@@ -268,47 +266,35 @@ def plot_persistence_diagram(
         import matplotlib.pyplot as plt
         import matplotlib.patches as mpatches
         from matplotlib import rc
+
         if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex():
-            plt.rc('text', usetex=True)
-            plt.rc('font', family='serif')
+            plt.rc("text", usetex=True)
+            plt.rc("font", family="serif")
         else:
-            plt.rc('text', usetex=False)
-            plt.rc('font', family='DejaVu Sans')
+            plt.rc("text", usetex=False)
+            plt.rc("font", family="DejaVu Sans")
 
         if persistence_file != "":
             if path.isfile(persistence_file):
                 # Reset persistence
                 persistence = []
-                diag = read_persistence_intervals_grouped_by_dimension(
-                    persistence_file=persistence_file
-                )
+                diag = read_persistence_intervals_grouped_by_dimension(persistence_file=persistence_file)
                 for key in diag.keys():
                     for persistence_interval in diag[key]:
                         persistence.append((key, persistence_interval))
             else:
-                print("file " + persistence_file + " not found.")
-                return None
-
-        persistence = _array_handler(persistence)
+                raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), persistence_file)
 
-        if max_plots != 1000:
-            print("Deprecated parameter. It has been replaced by max_intervals")
-            max_intervals = max_plots
-
-        if max_intervals > 0 and max_intervals < len(persistence):
-            # Sort by life time, then takes only the max_intervals elements
-            persistence = sorted(
-                persistence,
-                key=lambda life_time: life_time[1][1] - life_time[1][0],
-                reverse=True,
-            )[:max_intervals]
-
-        if colormap == None:
-            colormap = plt.cm.Set1.colors
-        if axes == None:
-            fig, axes = plt.subplots(1, 1)
+        try:
+            persistence = _array_handler(persistence)
+            persistence = _limit_to_max_intervals(
+                persistence, max_intervals, key=lambda life_time: life_time[1][1] - life_time[1][0]
+            )
+            min_birth, max_death = __min_birth_max_death(persistence, band)
+        except IndexError:
+            min_birth, max_death = 0.0, 1.0
+            pass
 
-        (min_birth, max_death) = __min_birth_max_death(persistence, band)
         delta = (max_death - min_birth) * inf_delta
         # Replace infinity values with max_death + delta for diagram to be more
         # readable
@@ -316,61 +302,56 @@ def plot_persistence_diagram(
         axis_end = max_death + delta / 2
         axis_start = min_birth - delta
 
+        if axes == None:
+            _, axes = plt.subplots(1, 1)
+        if colormap == None:
+            colormap = plt.cm.Set1.colors
         # bootstrap band
         if band > 0.0:
             x = np.linspace(axis_start, infinity, 1000)
             axes.fill_between(x, x, x + band, alpha=alpha, facecolor="red")
         # lower diag patch
         if greyblock:
-            axes.add_patch(mpatches.Polygon([[axis_start, axis_start], [axis_end, axis_start], [axis_end, axis_end]], fill=True, color='lightgrey'))
-        # Draw points in loop
-        pts_at_infty = False  # Records presence of pts at infty
-        for interval in reversed(persistence):
-            if float(interval[1][1]) != float("inf"):
-                # Finite death case
-                axes.scatter(
-                    interval[1][0],
-                    interval[1][1],
-                    alpha=alpha,
-                    color=colormap[interval[0]],
+            axes.add_patch(
+                mpatches.Polygon(
+                    [[axis_start, axis_start], [axis_end, axis_start], [axis_end, axis_end]],
+                    fill=True,
+                    color="lightgrey",
                 )
-            else:
-                pts_at_infty = True
-                # Infinite death case for diagram to be nicer
-                axes.scatter(
-                    interval[1][0], infinity, alpha=alpha, color=colormap[interval[0]]
-                )
-        if pts_at_infty:
+            )
+        # line display of equation : birth = death
+        axes.plot([axis_start, axis_end], [axis_start, axis_end], linewidth=1.0, color="k")
+
+        x=[birth for (dim,(birth,death)) in persistence]
+        y=[death if death != float("inf") else infinity for (dim,(birth,death)) in persistence]
+        c=[colormap[dim] for (dim,(birth,death)) in persistence]
+
+        axes.scatter(x,y,alpha=alpha,color=c)
+        if float("inf") in (death for (dim,(birth,death)) in persistence):
             # infinity line and text
-            axes.plot([axis_start, axis_end], [axis_start, axis_end], linewidth=1.0, color="k")
             axes.plot([axis_start, axis_end], [infinity, infinity], linewidth=1.0, color="k", alpha=alpha)
             # Infinity label
             yt = axes.get_yticks()
-            yt = yt[np.where(yt < axis_end)] # to avoid ploting ticklabel higher than infinity
+            yt = yt[np.where(yt < axis_end)]  # to avoid plotting ticklabel higher than infinity
             yt = np.append(yt, infinity)
             ytl = ["%.3f" % e for e in yt]  # to avoid float precision error
-            ytl[-1] = r'$+\infty$'
+            ytl[-1] = r"$+\infty$"
             axes.set_yticks(yt)
             axes.set_yticklabels(ytl)
 
         if legend:
             dimensions = list(set(item[0] for item in persistence))
-            axes.legend(
-                handles=[
-                    mpatches.Patch(color=colormap[dim], label=str(dim))
-                    for dim in dimensions
-                ]
-            )
+            axes.legend(handles=[mpatches.Patch(color=colormap[dim], label=str(dim)) for dim in dimensions])
 
         axes.set_xlabel("Birth", fontsize=fontsize)
         axes.set_ylabel("Death", fontsize=fontsize)
         axes.set_title("Persistence diagram", fontsize=fontsize)
         # Ends plot on infinity value and starts a little bit before min_birth
-        axes.axis([axis_start, axis_end, axis_start, infinity + delta/2])
+        axes.axis([axis_start, axis_end, axis_start, infinity + delta / 2])
         return axes
 
-    except ImportError:
-        print("This function is not available, you may be missing matplotlib.")
+    except ImportError as import_error:
+        warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.")
 
 
 def plot_persistence_density(
@@ -384,7 +365,7 @@ def plot_persistence_density(
     legend=False,
     axes=None,
     fontsize=16,
-    greyblock=False
+    greyblock=False,
 ):
     """This function plots the persistence density from persistence
     values list, np.array of shape (N x 2) representing a diagram
@@ -444,12 +425,13 @@ def plot_persistence_density(
         import matplotlib.patches as mpatches
         from scipy.stats import kde
         from matplotlib import rc
+
         if _gudhi_matplotlib_use_tex and _matplotlib_can_use_tex():
-            plt.rc('text', usetex=True)
-            plt.rc('font', family='serif')
+            plt.rc("text", usetex=True)
+            plt.rc("font", family="serif")
         else:
-            plt.rc('text', usetex=False)
-            plt.rc('font', family='DejaVu Sans')
+            plt.rc("text", usetex=False)
+            plt.rc("font", family="DejaVu Sans")
 
         if persistence_file != "":
             if dimension is None:
@@ -460,10 +442,16 @@ def plot_persistence_density(
                     persistence_file=persistence_file, only_this_dim=dimension
                 )
             else:
-                print("file " + persistence_file + " not found.")
-                return None
+                raise FileNotFoundError(errno.ENOENT, os.strerror(errno.ENOENT), persistence_file)
+
+        # default cmap value cannot be done at argument definition level as matplotlib is not yet defined.
+        if cmap is None:
+            cmap = plt.cm.hot_r
+        if axes == None:
+            _, axes = plt.subplots(1, 1)
 
-        if len(persistence) > 0:
+        try:
+            # if not read from file but given by an argument
             persistence = _array_handler(persistence)
             persistence_dim = np.array(
                 [
@@ -472,47 +460,54 @@ def plot_persistence_density(
                     if (dim_interval[0] == dimension) or (dimension is None)
                 ]
             )
-
-        persistence_dim = persistence_dim[np.isfinite(persistence_dim[:, 1])]
-        if max_intervals > 0 and max_intervals < len(persistence_dim):
-            # Sort by life time, then takes only the max_intervals elements
+            persistence_dim = persistence_dim[np.isfinite(persistence_dim[:, 1])]
             persistence_dim = np.array(
-                sorted(
-                    persistence_dim,
-                    key=lambda life_time: life_time[1] - life_time[0],
-                    reverse=True,
-                )[:max_intervals]
+                _limit_to_max_intervals(
+                    persistence_dim, max_intervals, key=lambda life_time: life_time[1] - life_time[0]
+                )
             )
 
-        # Set as numpy array birth and death (remove undefined values - inf and NaN)
-        birth = persistence_dim[:, 0]
-        death = persistence_dim[:, 1]
-
-        # default cmap value cannot be done at argument definition level as matplotlib is not yet defined.
-        if cmap is None:
-            cmap = plt.cm.hot_r
-        if axes == None:
-            fig, axes = plt.subplots(1, 1)
+            # Set as numpy array birth and death (remove undefined values - inf and NaN)
+            birth = persistence_dim[:, 0]
+            death = persistence_dim[:, 1]
+            birth_min = birth.min()
+            birth_max = birth.max()
+            death_min = death.min()
+            death_max = death.max()
+
+            # Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents
+            k = kde.gaussian_kde([birth, death], bw_method=bw_method)
+            xi, yi = np.mgrid[
+                birth_min : birth_max : nbins * 1j, death_min : death_max : nbins * 1j,
+            ]
+            zi = k(np.vstack([xi.flatten(), yi.flatten()]))
+            # Make the plot
+            img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap, shading="auto")
+            plot_success = True
+
+        # IndexError on empty diagrams, ValueError on only inf death values
+        except (IndexError, ValueError):
+            birth_min = 0.0
+            birth_max = 1.0
+            death_min = 0.0
+            death_max = 1.0
+            plot_success = False
+            pass
 
         # line display of equation : birth = death
-        x = np.linspace(death.min(), birth.max(), 1000)
+        x = np.linspace(death_min, birth_max, 1000)
         axes.plot(x, x, color="k", linewidth=1.0)
 
-        # Evaluate a gaussian kde on a regular grid of nbins x nbins over data extents
-        k = kde.gaussian_kde([birth, death], bw_method=bw_method)
-        xi, yi = np.mgrid[
-            birth.min() : birth.max() : nbins * 1j,
-            death.min() : death.max() : nbins * 1j,
-        ]
-        zi = k(np.vstack([xi.flatten(), yi.flatten()]))
-
-        # Make the plot
-        img = axes.pcolormesh(xi, yi, zi.reshape(xi.shape), cmap=cmap)
-
         if greyblock:
-            axes.add_patch(mpatches.Polygon([[birth.min(), birth.min()], [death.max(), birth.min()], [death.max(), death.max()]], fill=True, color='lightgrey'))
+            axes.add_patch(
+                mpatches.Polygon(
+                    [[birth_min, birth_min], [death_max, birth_min], [death_max, death_max]],
+                    fill=True,
+                    color="lightgrey",
+                )
+            )
 
-        if legend:
+        if plot_success and legend:
             plt.colorbar(img, ax=axes)
 
         axes.set_xlabel("Birth", fontsize=fontsize)
@@ -521,7 +516,5 @@ def plot_persistence_density(
 
         return axes
 
-    except ImportError:
-        print(
-            "This function is not available, you may be missing matplotlib and/or scipy."
-        )
+    except ImportError as import_error:
+        warnings.warn(f"This function is not available.\nModuleNotFoundError: No module named '{import_error.name}'.")
diff --git a/src/python/gudhi/representations/preprocessing.py b/src/python/gudhi/representations/preprocessing.py
index a8545349..8722e162 100644
--- a/src/python/gudhi/representations/preprocessing.py
+++ b/src/python/gudhi/representations/preprocessing.py
@@ -1,10 +1,11 @@
 # 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):       Mathieu Carrière
+# Author(s):       Mathieu Carrière, Vincent Rouvreau
 #
 # Copyright (C) 2018-2019 Inria
 #
 # Modification(s):
+#   - 2021/10 Vincent Rouvreau: Add DimensionSelector
 #   - YYYY/MM Author: Description of the modification
 
 import numpy as np
@@ -75,7 +76,7 @@ class Clamping(BaseEstimator, TransformerMixin):
         Constructor for the Clamping class.
 
         Parameters:
-            limit (double): clamping value (default np.inf).
+            limit (float): clamping value (default np.inf).
         """
         self.minimum = minimum
         self.maximum = maximum
@@ -234,7 +235,7 @@ class ProminentPoints(BaseEstimator, TransformerMixin):
             use (bool): whether to use the class or not (default False).
             location (string): either "upper" or "lower" (default "upper"). Whether to keep the points that are far away ("upper") or close ("lower") to the diagonal.
             num_pts (int): cardinality threshold (default 10). If location == "upper", keep the top **num_pts** points that are the farthest away from the diagonal. If location == "lower", keep the top **num_pts** points that are the closest to the diagonal. 
-            threshold (double): distance-to-diagonal threshold (default -1). If location == "upper", keep the points that are at least at a distance **threshold** from the diagonal. If location == "lower", keep the points that are at most at a distance **threshold** from the diagonal. 
+            threshold (float): distance-to-diagonal threshold (default -1). If location == "upper", keep the points that are at least at a distance **threshold** from the diagonal. If location == "lower", keep the points that are at most at a distance **threshold** from the diagonal. 
         """
         self.num_pts    = num_pts
         self.threshold  = threshold
@@ -317,7 +318,7 @@ class DiagramSelector(BaseEstimator, TransformerMixin):
 
         Parameters:
             use (bool): whether to use the class or not (default False).
-            limit (double): second coordinate value that is the criterion for being an essential point (default numpy.inf).
+            limit (float): second coordinate value that is the criterion for being an essential point (default numpy.inf).
             point_type (string): either "finite" or "essential". The type of the points that are going to be extracted.
         """
         self.use, self.limit, self.point_type = use, limit, point_type
@@ -363,3 +364,51 @@ class DiagramSelector(BaseEstimator, TransformerMixin):
             n x 2 numpy array: extracted persistence diagram.
         """
         return self.fit_transform([diag])[0]
+
+
+# Mermaid sequence diagram - https://mermaid-js.github.io/mermaid-live-editor/
+# sequenceDiagram
+#     USER->>DimensionSelector: fit_transform(<br/>[[array( Hi(X0) ), array( Hj(X0) ), ...],<br/> [array( Hi(X1) ), array( Hj(X1) ), ...],<br/> ...])
+#     DimensionSelector->>thread1: _transform([array( Hi(X0) ), array( Hj(X0) )], ...)
+#     DimensionSelector->>thread2: _transform([array( Hi(X1) ), array( Hj(X1) )], ...)
+#     Note right of DimensionSelector: ...
+#     thread1->>DimensionSelector: array( Hn(X0) )
+#     thread2->>DimensionSelector: array( Hn(X1) )
+#     Note right of DimensionSelector: ...
+#     DimensionSelector->>USER: [array( Hn(X0) ), <br/> array( Hn(X1) ), <br/> ...]
+
+class DimensionSelector(BaseEstimator, TransformerMixin):
+    """
+    This is a class to select persistence diagrams in a specific dimension from its index.
+    """
+
+    def __init__(self, index=0):
+        """
+        Constructor for the DimensionSelector class.
+
+        Parameters:
+            index (int): The returned persistence diagrams dimension index. Default value is `0`.
+        """
+        self.index = index
+
+    def fit(self, X, Y=None):
+        """
+        Nothing to be done, but useful when included in a scikit-learn Pipeline.
+        """
+        return self
+
+    def transform(self, X, Y=None):
+        """
+        Select persistence diagrams from its dimension.
+
+        Parameters:
+            X (list of list of tuple): List of list of persistence pairs, i.e.
+                `[[array( Hi(X0) ), array( Hj(X0) ), ...], [array( Hi(X1) ), array( Hj(X1) ), ...], ...]` 
+
+        Returns:
+            list of tuple:
+            Persistence diagrams in a specific dimension. i.e. if `index` was set to `m` and `Hn` is at index `m` of
+            the input, it returns `[array( Hn(X0) ), array( Hn(X1), ...]`
+        """
+
+        return [persistence[self.index] for persistence in X]
diff --git a/src/python/gudhi/representations/vector_methods.py b/src/python/gudhi/representations/vector_methods.py
index e883b5dd..69ff5e1e 100644
--- a/src/python/gudhi/representations/vector_methods.py
+++ b/src/python/gudhi/representations/vector_methods.py
@@ -1,15 +1,17 @@
 # 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):       Mathieu Carrière, Martin Royer
+# Author(s):       Mathieu Carrière, Martin Royer, Gard Spreemann
 #
 # Copyright (C) 2018-2020 Inria
 #
 # Modification(s):
 #   - 2020/06 Martin: ATOL integration
+#   - 2020/12 Gard: A more flexible Betti curve class capable of computing exact curves.
 #   - 2021/11 Vincent Rouvreau: factorize _automatic_sample_range
 
 import numpy as np
 from sklearn.base          import BaseEstimator, TransformerMixin
+from sklearn.exceptions    import NotFittedError
 from sklearn.preprocessing import MinMaxScaler, MaxAbsScaler
 from sklearn.neighbors     import DistanceMetric
 from sklearn.metrics       import pairwise
@@ -306,67 +308,162 @@ class Silhouette(BaseEstimator, TransformerMixin):
         """
         return self.fit_transform([diag])[0,:]
 
+
 class BettiCurve(BaseEstimator, TransformerMixin):
     """
-    This is a class for computing Betti curves from a list of persistence diagrams. A Betti curve is a 1D piecewise-constant function obtained from the rank function. It is sampled evenly on a given range and the vector of samples is returned. See https://www.researchgate.net/publication/316604237_Time_Series_Classification_via_Topological_Data_Analysis for more details.
+    Compute Betti curves from persistence diagrams. There are several modes of operation: with a given resolution (with or without a sample_range), with a predefined grid, and with none of the previous. With a predefined grid, the class computes the Betti numbers at those grid points. Without a predefined grid, if the resolution is set to None, it can be fit to a list of persistence diagrams and produce a grid that consists of (at least) the filtration values at which at least one of those persistence diagrams changes Betti numbers, and then compute the Betti numbers at those grid points. In the latter mode, the exact Betti curve is computed for the entire real line. Otherwise, if the resolution is given, the Betti curve is obtained by sampling evenly using either the given sample_range or based on the persistence diagrams.
     """
-    def __init__(self, resolution=100, sample_range=[np.nan, np.nan]):
+
+    def __init__(self, resolution=100, sample_range=[np.nan, np.nan], predefined_grid=None):
         """
         Constructor for the BettiCurve class.
 
         Parameters:
             resolution (int): number of sample for the piecewise-constant function (default 100).
             sample_range ([double, double]): minimum and maximum of the piecewise-constant function domain, of the form [x_min, x_max] (default [numpy.nan, numpy.nan]). It is the interval on which samples will be drawn evenly. If one of the values is numpy.nan, it can be computed from the persistence diagrams with the fit() method.
+            predefined_grid (1d array or None, default=None): Predefined filtration grid points at which to compute the Betti curves. Must be strictly ordered. Infinities are ok. If None (default), and resolution is given, the grid will be uniform from x_min to x_max in 'resolution' steps, otherwise a grid will be computed that captures all changes in Betti numbers in the provided data.
+
+        Attributes:
+            grid_ (1d array): The grid on which the Betti numbers are computed. If predefined_grid was specified, `grid_` will always be that grid, independently of data. If not, the grid is fitted to capture all filtration values at which the Betti numbers change.
+
+        Examples
+        --------
+        If pd is a persistence diagram and xs is a nonempty grid of finite values such that xs[0] >= pd.min(), then the results of:
+
+        >>> bc = BettiCurve(predefined_grid=xs) # doctest: +SKIP
+        >>> result = bc(pd) # doctest: +SKIP
+
+        and
+
+        >>> from scipy.interpolate import interp1d # doctest: +SKIP
+        >>> bc = BettiCurve(resolution=None, predefined_grid=None) # doctest: +SKIP
+        >>> bettis = bc.fit_transform([pd]) # doctest: +SKIP
+        >>> interp = interp1d(bc.grid_, bettis[0, :], kind="previous", fill_value="extrapolate") # doctest: +SKIP
+        >>> result = np.array(interp(xs), dtype=int) # doctest: +SKIP
+
+        are the same.
         """
-        self.resolution, self.sample_range = resolution, sample_range
 
-    def fit(self, X, y=None):
+        if (predefined_grid is not None) and (not isinstance(predefined_grid, np.ndarray)):
+            raise ValueError("Expected predefined_grid as array or None.")
+
+        self.predefined_grid = predefined_grid
+        self.resolution = resolution
+        self.sample_range = sample_range
+
+    def is_fitted(self):
+        return hasattr(self, "grid_")
+
+    def fit(self, X, y = None):
         """
-        Fit the BettiCurve class on a list of persistence diagrams: if any of the values in **sample_range** is numpy.nan, replace it with the corresponding value computed on the given list of persistence diagrams.
+        Fit the BettiCurve class on a list of persistence diagrams: if any of the values in **sample_range** is numpy.nan, replace it with the corresponding value computed on the given list of persistence diagrams. When no predefined grid is provided and resolution set to None, compute a filtration grid that captures all changes in Betti numbers for all the given persistence diagrams.
 
         Parameters:
-            X (list of n x 2 numpy arrays): input persistence diagrams.
-            y (n x 1 array): persistence diagram labels (unused).
+            X (list of 2d arrays): Persistence diagrams.
+            y (None): Ignored.
         """
-        self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y)
+
+        if self.predefined_grid is None:
+            if self.resolution is None: # Flexible/exact version
+                events = np.unique(np.concatenate([pd.flatten() for pd in X] + [[-np.inf]], axis=0))
+                self.grid_ = np.array(events)
+            else:
+                self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y)
+                self.grid_ = np.linspace(self.sample_range[0], self.sample_range[1], self.resolution)
+        else:
+            self.grid_ = self.predefined_grid # Get the predefined grid from user
+
         return self
 
     def transform(self, X):
         """
-        Compute the Betti curve for each persistence diagram individually and concatenate the results.
+        Compute Betti curves.
 
         Parameters:
-            X (list of n x 2 numpy arrays): input persistence diagrams.
-    
+            X (list of 2d arrays): Persistence diagrams.
+
         Returns:
-            numpy array with shape (number of diagrams) x (**resolution**): output Betti curves.
+            `len(X).len(self.grid_)` array of ints: Betti numbers of the given persistence diagrams at the grid points given in `self.grid_`
         """
-        Xfit = []
-        x_values = np.linspace(self.sample_range[0], self.sample_range[1], self.resolution)
-        step_x = x_values[1] - x_values[0]
 
-        for diagram in X:
-            diagram_int = np.clip(np.ceil((diagram[:,:2] - self.sample_range[0]) / step_x), 0, self.resolution).astype(int)
-            bc =  np.zeros(self.resolution)
-            for interval in diagram_int:
-                bc[interval[0]:interval[1]] += 1
-            Xfit.append(np.reshape(bc,[1,-1]))
+        if not self.is_fitted():
+            raise NotFittedError("Not fitted.")
 
-        Xfit = np.concatenate(Xfit, 0)
+        if not X:
+            X = [np.zeros((0, 2))]
+        
+        N = len(X)
 
-        return Xfit
+        events = np.concatenate([pd.flatten(order="F") for pd in X], axis=0)
+        sorting = np.argsort(events)
+        offsets = np.zeros(1 + N, dtype=int)
+        for i in range(0, N):
+            offsets[i+1] = offsets[i] + 2*X[i].shape[0]
+        starts = offsets[0:N]
+        ends = offsets[1:N + 1] - 1
 
-    def __call__(self, diag):
+        bettis = [[0] for i in range(0, N)]
+
+        i = 0
+        for x in self.grid_:
+            while i < len(sorting) and events[sorting[i]] <= x:
+                j = np.searchsorted(ends, sorting[i])
+                delta = 1 if sorting[i] - starts[j] < len(X[j]) else -1
+                bettis[j][-1] += delta
+                i += 1
+            for k in range(0, N):
+                bettis[k].append(bettis[k][-1])
+
+        return np.array(bettis, dtype=int)[:, 0:-1]
+
+    def fit_transform(self, X):
+        """
+        The result is the same as fit(X) followed by transform(X), but potentially faster.
         """
-        Apply BettiCurve on a single persistence diagram and outputs the result.
 
-        Parameters:
-            diag (n x 2 numpy array): input persistence diagram.
+        if self.predefined_grid is None and self.resolution is None:
+            if not X:
+                X = [np.zeros((0, 2))]
 
-        Returns:
-            numpy array with shape (**resolution**): output Betti curve.
+            N = len(X)
+
+            events = np.concatenate([pd.flatten(order="F") for pd in X], axis=0)
+            sorting = np.argsort(events)
+            offsets = np.zeros(1 + N, dtype=int)
+            for i in range(0, N):
+                offsets[i+1] = offsets[i] + 2*X[i].shape[0]
+            starts = offsets[0:N]
+            ends = offsets[1:N + 1] - 1
+
+            xs = [-np.inf]
+            bettis = [[0] for i in range(0, N)]
+
+            for i in sorting:
+                j = np.searchsorted(ends, i)
+                delta = 1 if i - starts[j] < len(X[j]) else -1
+                if events[i] == xs[-1]:
+                    bettis[j][-1] += delta
+                else:
+                    xs.append(events[i])
+                    for k in range(0, j):
+                        bettis[k].append(bettis[k][-1])
+                    bettis[j].append(bettis[j][-1] + delta)
+                    for k in range(j+1, N):
+                        bettis[k].append(bettis[k][-1])
+
+            self.grid_ = np.array(xs)
+            return np.array(bettis, dtype=int)
+
+        else:
+            return self.fit(X).transform(X)
+
+    def __call__(self, diag):
         """
-        return self.fit_transform([diag])[0,:]
+        Shorthand for transform on a single persistence diagram.
+        """
+        return self.fit_transform([diag])[0, :]
+
+
 
 class Entropy(BaseEstimator, TransformerMixin):
     """
@@ -411,26 +508,20 @@ class Entropy(BaseEstimator, TransformerMixin):
         new_X = BirthPersistenceTransform().fit_transform(X)        
 
         for i in range(num_diag):
-            orig_diagram, diagram, num_pts_in_diag = X[i], new_X[i], X[i].shape[0]
-            try:
-                new_diagram = DiagramScaler(use=True, scalers=[([1], MaxAbsScaler())]).fit_transform([diagram])[0]
-            except ValueError:
-                # Empty persistence diagram case - https://github.com/GUDHI/gudhi-devel/issues/507
-                assert len(diagram) == 0
-                new_diagram = np.empty(shape = [0, 2])
-
+            orig_diagram, new_diagram, num_pts_in_diag = X[i], new_X[i], X[i].shape[0]
+                
+            p = new_diagram[:,1]
+            p = p/np.sum(p)
             if self.mode == "scalar":
-                ent = - np.sum( np.multiply(new_diagram[:,1], np.log(new_diagram[:,1])) )
+                ent = -np.dot(p, np.log(p))
                 Xfit.append(np.array([[ent]]))
-
             else:
                 ent = np.zeros(self.resolution)
                 for j in range(num_pts_in_diag):
                     [px,py] = orig_diagram[j,:2]
                     min_idx = np.clip(np.ceil((px - self.sample_range[0]) / step_x).astype(int), 0, self.resolution)
                     max_idx = np.clip(np.ceil((py - self.sample_range[0]) / step_x).astype(int), 0, self.resolution)
-                    for k in range(min_idx, max_idx):
-                        ent[k] += (-1) * new_diagram[j,1] * np.log(new_diagram[j,1])
+                    ent[min_idx:max_idx]-=p[j]*np.log(p[j])
                 if self.normalized:
                     ent = ent / np.linalg.norm(ent, ord=1)
                 Xfit.append(np.reshape(ent,[1,-1]))
diff --git a/src/python/gudhi/simplex_tree.pxd b/src/python/gudhi/simplex_tree.pxd
index 006a24ed..5642f82d 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -45,6 +45,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
 
     cdef cppclass Simplex_tree_interface_full_featured "Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>":
         Simplex_tree_interface_full_featured() nogil
+        Simplex_tree_interface_full_featured(Simplex_tree_interface_full_featured&) nogil
         double simplex_filtration(vector[int] simplex) nogil
         void assign_simplex_filtration(vector[int] simplex, double filtration) nogil
         void initialize_filtration() nogil
@@ -62,9 +63,9 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         bool prune_above_filtration(double filtration) nogil
         bool make_filtration_non_decreasing() nogil
         void compute_extended_filtration() nogil
-        vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(vector[pair[int, pair[double, double]]] dgm, double min_persistence) nogil
         Simplex_tree_interface_full_featured* collapse_edges(int nb_collapse_iteration) nogil except +
         void reset_filtration(double filtration, int dimension) nogil
+        bint operator==(Simplex_tree_interface_full_featured) nogil
         # Iterators over Simplex tree
         pair[vector[int], double] get_simplex_and_filtration(Simplex_tree_simplex_handle f_simplex) nogil
         Simplex_tree_simplices_iterator get_simplices_iterator_begin() nogil
@@ -74,9 +75,12 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         Simplex_tree_skeleton_iterator get_skeleton_iterator_begin(int dimension) nogil
         Simplex_tree_skeleton_iterator get_skeleton_iterator_end(int dimension) nogil
         pair[Simplex_tree_boundary_iterator, Simplex_tree_boundary_iterator] get_boundary_iterators(vector[int] simplex) nogil except +
+        # Expansion with blockers
+        ctypedef bool (*blocker_func_t)(vector[int], void *user_data)
+        void expansion_with_blockers_callback(int dimension, blocker_func_t user_func, void *user_data)
 
 cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
-    cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Simplex_tree<Gudhi::Simplex_tree_options_full_featured>>":
+    cdef cppclass Simplex_tree_persistence_interface "Gudhi::Persistent_cohomology_interface<Gudhi::Simplex_tree_interface<Gudhi::Simplex_tree_options_full_featured>>":
         Simplex_tree_persistence_interface(Simplex_tree_interface_full_featured * st, bool persistence_dim_max) nogil
         void compute_persistence(int homology_coeff_field, double min_persistence) nogil except +
         vector[pair[int, pair[double, double]]] get_persistence() nogil
@@ -87,3 +91,4 @@ cdef extern from "Persistent_cohomology_interface.h" namespace "Gudhi":
         vector[pair[vector[int], vector[int]]] persistence_pairs() nogil
         pair[vector[vector[int]], vector[vector[int]]] lower_star_generators() nogil
         pair[vector[vector[int]], vector[vector[int]]] flag_generators() nogil
+        vector[vector[pair[int, pair[double, double]]]] compute_extended_persistence_subdiagrams(double min_persistence) nogil
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index c3720936..05bfe22e 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -16,6 +16,9 @@ __author__ = "Vincent Rouvreau"
 __copyright__ = "Copyright (C) 2016 Inria"
 __license__ = "MIT"
 
+cdef bool callback(vector[int] simplex, void *blocker_func):
+    return (<object>blocker_func)(simplex)
+
 # SimplexTree python interface
 cdef class SimplexTree:
     """The simplex tree is an efficient and flexible data structure for
@@ -38,13 +41,29 @@ cdef class SimplexTree:
     cdef Simplex_tree_persistence_interface * pcohptr
 
     # Fake constructor that does nothing but documenting the constructor
-    def __init__(self):
+    def __init__(self, other = None):
         """SimplexTree constructor.
+
+        :param other: If `other` is `None` (default value), an empty `SimplexTree` is created.
+            If `other` is a `SimplexTree`, the `SimplexTree` is constructed from a deep copy of `other`.
+        :type other: SimplexTree (Optional)
+        :returns: An empty or a copy simplex tree.
+        :rtype: SimplexTree
+
+        :raises TypeError: In case `other` is neither `None`, nor a `SimplexTree`.
+        :note: If the `SimplexTree` is a copy, the persistence information is not copied. If you need it in the clone,
+            you have to call :func:`compute_persistence` on it even if you had already computed it in the original.
         """
 
     # The real cython constructor
-    def __cinit__(self):
-        self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured())
+    def __cinit__(self, other = None):
+        if other:
+            if isinstance(other, SimplexTree):
+                self.thisptr = _get_copy_intptr(other)
+            else:
+                raise TypeError("`other` argument requires to be of type `SimplexTree`, or `None`.")
+        else:
+            self.thisptr = <intptr_t>(new Simplex_tree_interface_full_featured())
 
     def __dealloc__(self):
         cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr()
@@ -63,6 +82,21 @@ cdef class SimplexTree:
          """
         return self.pcohptr != NULL
 
+    def copy(self):
+        """ 
+        :returns: A simplex tree that is a deep copy of itself.
+        :rtype: SimplexTree
+
+        :note: The persistence information is not copied. If you need it in the clone, you have to call
+            :func:`compute_persistence` on it even if you had already computed it in the original.
+        """
+        stree = SimplexTree()
+        stree.thisptr = _get_copy_intptr(self)
+        return stree
+
+    def __deepcopy__(self):
+        return self.copy()
+
     def filtration(self, simplex):
         """This function returns the filtration value for a given N-simplex in
         this simplicial complex, or +infinity if it is not in the complex.
@@ -440,9 +474,29 @@ cdef class SimplexTree:
             del self.pcohptr
         self.pcohptr = new Simplex_tree_persistence_interface(self.get_ptr(), False)
         self.pcohptr.compute_persistence(homology_coeff_field, -1.)
-        persistence_result = self.pcohptr.get_persistence()
-        return self.get_ptr().compute_extended_persistence_subdiagrams(persistence_result, min_persistence)
+        return self.pcohptr.compute_extended_persistence_subdiagrams(min_persistence)
+
+    def expansion_with_blocker(self, max_dim, blocker_func):
+        """Expands the Simplex_tree containing only a graph. Simplices corresponding to cliques in the graph are added
+        incrementally, faces before cofaces, unless the simplex has dimension larger than `max_dim` or `blocker_func`
+        returns `True` for this simplex.
+
+        The function identifies a candidate simplex whose faces are all already in the complex, inserts it with a
+        filtration value corresponding to the maximum of the filtration values of the faces, then calls `blocker_func`
+        with this new simplex (represented as a list of int). If `blocker_func` returns `True`, the simplex is removed,
+        otherwise it is kept. The algorithm then proceeds with the next candidate.
 
+        .. warning::
+            Several candidates of the same dimension may be inserted simultaneously before calling `blocker_func`, so
+            if you examine the complex in `blocker_func`, you may hit a few simplices of the same dimension that have
+            not been vetted by `blocker_func` yet, or have already been rejected but not yet removed.
+
+        :param max_dim: Expansion maximal dimension value.
+        :type max_dim: int
+        :param blocker_func: Blocker oracle.
+        :type blocker_func: Callable[[List[int]], bool]
+        """
+        self.get_ptr().expansion_with_blockers_callback(max_dim, callback, <void*>blocker_func)
 
     def persistence(self, homology_coeff_field=11, min_persistence=0, persistence_dim_max = False):
         """This function computes and returns the persistence of the simplicial complex.
@@ -621,18 +675,17 @@ cdef class SimplexTree:
         return (normal0, normals, infinite0, infinites)
 
     def collapse_edges(self, nb_iterations = 1):
-        """Assuming the simplex tree is a 1-skeleton graph, this method collapse edges (simplices of higher dimension
-        are ignored) and resets the simplex tree from the remaining edges.
-        A good candidate is to build a simplex tree on top of a :class:`~gudhi.RipsComplex` of dimension 1 before
-        collapsing edges
+        """Assuming the complex is a graph (simplices of higher dimension are ignored), this method implicitly
+        interprets it as the 1-skeleton of a flag complex, and replaces it with another (smaller) graph whose
+        expansion has the same persistent homology, using a technique known as edge collapses
+        (see :cite:`edgecollapsearxiv`).
+
+        A natural application is to get a simplex tree of dimension 1 from :class:`~gudhi.RipsComplex`,
+        then collapse edges, perform :meth:`expansion()` and finally compute persistence
         (cf. :download:`rips_complex_edge_collapse_example.py <../example/rips_complex_edge_collapse_example.py>`).
-        For implementation details, please refer to :cite:`edgecollapsesocg2020`.
 
         :param nb_iterations: The number of edge collapse iterations to perform. Default is 1.
         :type nb_iterations: int
-
-        :note: collapse_edges method requires `Eigen <installation.html#eigen>`_ >= 3.1.0 and an exception is thrown
-            if this method is not available.
         """
         # Backup old pointer
         cdef Simplex_tree_interface_full_featured* ptr = self.get_ptr()
@@ -642,3 +695,13 @@ cdef class SimplexTree:
             self.thisptr = <intptr_t>(ptr.collapse_edges(nb_iter))
             # Delete old pointer
             del ptr
+
+    def __eq__(self, other:SimplexTree):
+        """Test for structural equality
+        :returns: True if the 2 simplex trees are equal, False otherwise.
+        :rtype: bool
+        """
+        return dereference(self.get_ptr()) == dereference(other.get_ptr())
+
+cdef intptr_t _get_copy_intptr(SimplexTree stree) nogil:
+    return <intptr_t>(new Simplex_tree_interface_full_featured(dereference(stree.get_ptr())))
diff --git a/src/python/gudhi/sklearn/__init__.py b/src/python/gudhi/sklearn/__init__.py
new file mode 100644
index 00000000..e69de29b
--- /dev/null
+++ b/src/python/gudhi/sklearn/__init__.py
diff --git a/src/python/gudhi/sklearn/cubical_persistence.py b/src/python/gudhi/sklearn/cubical_persistence.py
new file mode 100644
index 00000000..672af278
--- /dev/null
+++ b/src/python/gudhi/sklearn/cubical_persistence.py
@@ -0,0 +1,110 @@
+# 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):       Vincent Rouvreau
+#
+# Copyright (C) 2021 Inria
+#
+# Modification(s):
+#   - YYYY/MM Author: Description of the modification
+
+from .. import CubicalComplex
+from sklearn.base import BaseEstimator, TransformerMixin
+
+import numpy as np
+# joblib is required by scikit-learn
+from joblib import Parallel, delayed
+
+# Mermaid sequence diagram - https://mermaid-js.github.io/mermaid-live-editor/
+# sequenceDiagram
+#     USER->>CubicalPersistence: fit_transform(X)
+#     CubicalPersistence->>thread1: _tranform(X[0])
+#     CubicalPersistence->>thread2: _tranform(X[1])
+#     Note right of CubicalPersistence: ...
+#     thread1->>CubicalPersistence: [array( H0(X[0]) ), array( H1(X[0]) )]
+#     thread2->>CubicalPersistence: [array( H0(X[1]) ), array( H1(X[1]) )]
+#     Note right of CubicalPersistence: ...
+#     CubicalPersistence->>USER: [[array( H0(X[0]) ), array( H1(X[0]) )],<br/> [array( H0(X[1]) ), array( H1(X[1]) )],<br/> ...]
+
+
+class CubicalPersistence(BaseEstimator, TransformerMixin):
+    """
+    This is a class for computing the persistence diagrams from a cubical complex.
+    """
+
+    def __init__(
+        self,
+        homology_dimensions,
+        newshape=None,
+        homology_coeff_field=11,
+        min_persistence=0.0,
+        n_jobs=None,
+    ):
+        """
+        Constructor for the CubicalPersistence class.
+
+        Parameters:
+            homology_dimensions (int or list of int): The returned persistence diagrams dimension(s).
+                Short circuit the use of :class:`~gudhi.representations.preprocessing.DimensionSelector` when only one
+                dimension matters (in other words, when `homology_dimensions` is an int).
+            newshape (tuple of ints): If cells filtration values require to be reshaped
+                (cf. :func:`~gudhi.sklearn.cubical_persistence.CubicalPersistence.transform`), set `newshape`
+                to perform `numpy.reshape(X, newshape, order='C')` in
+                :func:`~gudhi.sklearn.cubical_persistence.CubicalPersistence.transform` method.
+            homology_coeff_field (int): The homology coefficient field. Must be a prime number. Default value is 11.
+            min_persistence (float): The minimum persistence value to take into account (strictly greater than
+                `min_persistence`). Default value is `0.0`. Set `min_persistence` to `-1.0` to see all values.
+            n_jobs (int): cf. https://joblib.readthedocs.io/en/latest/generated/joblib.Parallel.html
+        """
+        self.homology_dimensions = homology_dimensions
+        self.newshape = newshape
+        self.homology_coeff_field = homology_coeff_field
+        self.min_persistence = min_persistence
+        self.n_jobs = n_jobs
+
+    def fit(self, X, Y=None):
+        """
+        Nothing to be done, but useful when included in a scikit-learn Pipeline.
+        """
+        return self
+
+    def __transform(self, cells):
+        cubical_complex = CubicalComplex(top_dimensional_cells=cells)
+        cubical_complex.compute_persistence(
+            homology_coeff_field=self.homology_coeff_field, min_persistence=self.min_persistence
+        )
+        return [
+            cubical_complex.persistence_intervals_in_dimension(dim) for dim in self.homology_dimensions
+        ]
+
+    def __transform_only_this_dim(self, cells):
+        cubical_complex = CubicalComplex(top_dimensional_cells=cells)
+        cubical_complex.compute_persistence(
+            homology_coeff_field=self.homology_coeff_field, min_persistence=self.min_persistence
+        )
+        return cubical_complex.persistence_intervals_in_dimension(self.homology_dimensions)
+
+    def transform(self, X, Y=None):
+        """Compute all the cubical complexes and their associated persistence diagrams.
+
+        :param X: List of cells filtration values (`numpy.reshape(X, newshape, order='C'` if `newshape` is set with a tuple of ints).
+        :type X: list of list of float OR list of numpy.ndarray
+
+        :return: Persistence diagrams in the format:
+
+              - If `homology_dimensions` was set to `n`: `[array( Hn(X[0]) ), array( Hn(X[1]) ), ...]` 
+              - If `homology_dimensions` was set to `[i, j]`: `[[array( Hi(X[0]) ), array( Hj(X[0]) )], [array( Hi(X[1]) ), array( Hj(X[1]) )], ...]`
+        :rtype: list of (,2) array_like or list of list of (,2) array_like
+        """
+        if self.newshape is not None:
+            X = np.reshape(X, self.newshape, order='C')
+        
+        # Depends on homology_dimensions is an integer or a list of integer (else case)
+        if isinstance(self.homology_dimensions, int):
+            # threads is preferred as cubical construction and persistence computation releases the GIL
+            return Parallel(n_jobs=self.n_jobs, prefer="threads")(
+                delayed(self.__transform_only_this_dim)(cells) for cells in X
+            )
+        else:
+            # threads is preferred as cubical construction and persistence computation releases the GIL
+            return Parallel(n_jobs=self.n_jobs, prefer="threads")(delayed(self.__transform)(cells) for cells in X)
+
diff --git a/src/python/gudhi/tensorflow/__init__.py b/src/python/gudhi/tensorflow/__init__.py
new file mode 100644
index 00000000..1599cf52
--- /dev/null
+++ b/src/python/gudhi/tensorflow/__init__.py
@@ -0,0 +1,5 @@
+from .cubical_layer import CubicalLayer
+from .lower_star_simplex_tree_layer import LowerStarSimplexTreeLayer
+from .rips_layer import RipsLayer
+
+__all__ = ["LowerStarSimplexTreeLayer", "RipsLayer", "CubicalLayer"]
diff --git a/src/python/gudhi/tensorflow/cubical_layer.py b/src/python/gudhi/tensorflow/cubical_layer.py
new file mode 100644
index 00000000..3304e719
--- /dev/null
+++ b/src/python/gudhi/tensorflow/cubical_layer.py
@@ -0,0 +1,82 @@
+import numpy               as np
+import tensorflow          as tf
+from ..cubical_complex  import CubicalComplex
+
+######################
+# Cubical filtration #
+######################
+
+# The parameters of the model are the pixel values.
+
+def _Cubical(Xflat, Xdim, dimensions, homology_coeff_field):
+    # Parameters: Xflat (flattened image),
+    #             Xdim (shape of non-flattened image)
+    #             dimensions (homology dimensions)
+
+    # Compute the persistence pairs with Gudhi
+    # We reverse the dimensions because CubicalComplex uses Fortran ordering
+    cc = CubicalComplex(dimensions=Xdim[::-1], top_dimensional_cells=Xflat)
+    cc.compute_persistence(homology_coeff_field=homology_coeff_field)
+
+    # Retrieve and ouput image indices/pixels corresponding to positive and negative simplices    
+    cof_pp = cc.cofaces_of_persistence_pairs()
+    
+    L_cofs = []
+    for dim in dimensions:
+
+        try:
+            cof = cof_pp[0][dim]
+        except IndexError:
+            cof = np.array([])
+
+        L_cofs.append(np.array(cof, dtype=np.int32))
+
+    return L_cofs
+
+class CubicalLayer(tf.keras.layers.Layer):
+    """
+    TensorFlow layer for computing the persistent homology of a cubical complex
+    """
+    def __init__(self, homology_dimensions, min_persistence=None, homology_coeff_field=11, **kwargs):
+        """
+        Constructor for the CubicalLayer class
+
+        Parameters:
+            homology_dimensions (List[int]): list of homology dimensions
+            min_persistence (List[float]): minimum distance-to-diagonal of the points in the output persistence diagrams (default None, in which case 0. is used for all dimensions)
+            homology_coeff_field (int): homology field coefficient. Must be a prime number. Default value is 11. Max is 46337.
+        """
+        super().__init__(dynamic=True, **kwargs)
+        self.dimensions = homology_dimensions
+        self.min_persistence = min_persistence if min_persistence != None else [0.] * len(self.dimensions)
+        self.hcf = homology_coeff_field
+        assert len(self.min_persistence) == len(self.dimensions)
+
+    def call(self, X):
+        """
+        Compute persistence diagram associated to a cubical complex filtered by some pixel values 
+
+        Parameters:
+            X (TensorFlow variable): pixel values of the cubical complex
+
+        Returns:
+            List[Tuple[tf.Tensor,tf.Tensor]]: List of cubical persistence diagrams. The length of this list is the same than that of dimensions, i.e., there is one persistence diagram per homology dimension provided in the input list dimensions. Moreover, the finite and essential parts of the persistence diagrams are provided separately: each element of this list is a tuple of size two that contains the finite and essential parts of the corresponding persistence diagram, of shapes [num_finite_points, 2] and [num_essential_points, 1] respectively. Note that the essential part is always empty in cubical persistence diagrams, except in homology dimension zero, where the essential part always contains a single point, with abscissa equal to the smallest value in the complex, and infinite ordinate
+        """
+        # Compute pixels associated to positive and negative simplices 
+        # Don't compute gradient for this operation
+        Xflat = tf.reshape(X, [-1])
+        Xdim, Xflat_numpy = X.shape, Xflat.numpy()
+        indices_list = _Cubical(Xflat_numpy, Xdim, self.dimensions, self.hcf)
+        index_essential = np.argmin(Xflat_numpy) # index of minimum pixel value for essential persistence diagram
+        # Get persistence diagram by simply picking the corresponding entries in the image
+        self.dgms = []
+        for idx_dim, dimension in enumerate(self.dimensions):
+            finite_dgm = tf.reshape(tf.gather(Xflat, indices_list[idx_dim]), [-1,2])
+            essential_dgm = tf.reshape(tf.gather(Xflat, index_essential), [-1,1]) if dimension == 0 else tf.zeros([0, 1])
+            min_pers = self.min_persistence[idx_dim]
+            if min_pers >= 0:
+                persistent_indices = tf.where(tf.math.abs(finite_dgm[:,1]-finite_dgm[:,0]) > min_pers)
+                self.dgms.append((tf.reshape(tf.gather(finite_dgm, indices=persistent_indices), [-1,2]), essential_dgm))
+            else:
+                self.dgms.append((finite_dgm, essential_dgm))
+        return self.dgms
diff --git a/src/python/gudhi/tensorflow/lower_star_simplex_tree_layer.py b/src/python/gudhi/tensorflow/lower_star_simplex_tree_layer.py
new file mode 100644
index 00000000..5a8e5b75
--- /dev/null
+++ b/src/python/gudhi/tensorflow/lower_star_simplex_tree_layer.py
@@ -0,0 +1,87 @@
+import numpy               as np
+import tensorflow          as tf
+
+#########################################
+# Lower star filtration on simplex tree #
+#########################################
+
+# The parameters of the model are the vertex function values of the simplex tree.
+
+def _LowerStarSimplexTree(simplextree, filtration, dimensions, homology_coeff_field):
+    # Parameters: simplextree (simplex tree on which to compute persistence)
+    #             filtration (function values on the vertices of st),
+    #             dimensions (homology dimensions),
+    #             homology_coeff_field (homology field coefficient)
+    
+    simplextree.reset_filtration(-np.inf, 0)
+
+    # Assign new filtration values
+    for i in range(simplextree.num_vertices()):
+        simplextree.assign_filtration([i], filtration[i])
+    simplextree.make_filtration_non_decreasing()
+    
+    # Compute persistence diagram
+    simplextree.compute_persistence(homology_coeff_field=homology_coeff_field)
+    
+    # Get vertex pairs for optimization. First, get all simplex pairs
+    pairs = simplextree.lower_star_persistence_generators()
+    
+    L_indices = []
+    for dimension in dimensions:
+    
+        finite_pairs = pairs[0][dimension] if len(pairs[0]) >= dimension+1 else np.empty(shape=[0,2])
+        essential_pairs = pairs[1][dimension] if len(pairs[1]) >= dimension+1 else np.empty(shape=[0,1])
+        
+        finite_indices = np.array(finite_pairs.flatten(), dtype=np.int32)
+        essential_indices = np.array(essential_pairs.flatten(), dtype=np.int32)
+
+        L_indices.append((finite_indices, essential_indices))
+
+    return L_indices
+
+class LowerStarSimplexTreeLayer(tf.keras.layers.Layer):
+    """
+    TensorFlow layer for computing lower-star persistence out of a simplex tree
+    """
+    def __init__(self, simplextree, homology_dimensions, min_persistence=None, homology_coeff_field=11, **kwargs):
+        """
+        Constructor for the LowerStarSimplexTreeLayer class
+  
+        Parameters:
+            simplextree (gudhi.SimplexTree): underlying simplex tree. Its vertices MUST be named with integers from 0 to n-1, where n is its number of vertices. Note that its filtration values are modified in each call of the class.
+            homology_dimensions (List[int]): list of homology dimensions
+            min_persistence (List[float]): minimum distance-to-diagonal of the points in the output persistence diagrams (default None, in which case 0. is used for all dimensions)
+            homology_coeff_field (int): homology field coefficient. Must be a prime number. Default value is 11. Max is 46337.
+        """
+        super().__init__(dynamic=True, **kwargs)
+        self.dimensions  = homology_dimensions
+        self.simplextree = simplextree
+        self.min_persistence = min_persistence if min_persistence != None else [0. for _ in range(len(self.dimensions))]
+        self.hcf = homology_coeff_field
+        assert len(self.min_persistence) == len(self.dimensions)
+
+    def call(self, filtration):
+        """
+        Compute lower-star persistence diagram associated to a function defined on the vertices of the simplex tree
+
+        Parameters:
+            F (TensorFlow variable): filter function values over the vertices of the simplex tree. The ith entry of F corresponds to vertex i in self.simplextree
+
+        Returns:
+            List[Tuple[tf.Tensor,tf.Tensor]]: List of lower-star persistence diagrams. The length of this list is the same than that of dimensions, i.e., there is one persistence diagram per homology dimension provided in the input list dimensions. Moreover, the finite and essential parts of the persistence diagrams are provided separately: each element of this list is a tuple of size two that contains the finite and essential parts of the corresponding persistence diagram, of shapes [num_finite_points, 2] and [num_essential_points, 1] respectively
+        """
+        # Don't try to compute gradients for the vertex pairs
+        indices = _LowerStarSimplexTree(self.simplextree, filtration.numpy(), self.dimensions, self.hcf)
+        # Get persistence diagrams
+        self.dgms = []
+        for idx_dim, dimension in enumerate(self.dimensions):
+            finite_dgm = tf.reshape(tf.gather(filtration, indices[idx_dim][0]), [-1,2])
+            essential_dgm = tf.reshape(tf.gather(filtration, indices[idx_dim][1]), [-1,1])
+            min_pers = self.min_persistence[idx_dim]
+            if min_pers >= 0:
+                persistent_indices = tf.where(tf.math.abs(finite_dgm[:,1]-finite_dgm[:,0]) > min_pers)
+                self.dgms.append((tf.reshape(tf.gather(finite_dgm, indices=persistent_indices),[-1,2]), essential_dgm))
+            else:
+                self.dgms.append((finite_dgm, essential_dgm))
+        return self.dgms
+
diff --git a/src/python/gudhi/tensorflow/rips_layer.py b/src/python/gudhi/tensorflow/rips_layer.py
new file mode 100644
index 00000000..2a73472c
--- /dev/null
+++ b/src/python/gudhi/tensorflow/rips_layer.py
@@ -0,0 +1,93 @@
+import numpy               as np
+import tensorflow          as tf
+from ..rips_complex     import RipsComplex
+
+############################
+# Vietoris-Rips filtration #
+############################
+
+# The parameters of the model are the point coordinates.
+
+def _Rips(DX, max_edge, dimensions, homology_coeff_field):
+    # Parameters: DX (distance matrix), 
+    #             max_edge (maximum edge length for Rips filtration), 
+    #             dimensions (homology dimensions)
+
+    # Compute the persistence pairs with Gudhi
+    rc = RipsComplex(distance_matrix=DX, max_edge_length=max_edge)
+    st = rc.create_simplex_tree(max_dimension=max(dimensions)+1)
+    st.compute_persistence(homology_coeff_field=homology_coeff_field)
+    pairs = st.flag_persistence_generators()
+
+    L_indices = []
+    for dimension in dimensions:
+
+        if dimension == 0:
+            finite_pairs = pairs[0]
+            essential_pairs = pairs[2]
+        else:
+            finite_pairs = pairs[1][dimension-1] if len(pairs[1]) >= dimension else np.empty(shape=[0,4])
+            essential_pairs = pairs[3][dimension-1] if len(pairs[3]) >= dimension else np.empty(shape=[0,2])
+        
+        finite_indices = np.array(finite_pairs.flatten(), dtype=np.int32)
+        essential_indices = np.array(essential_pairs.flatten(), dtype=np.int32)
+
+        L_indices.append((finite_indices, essential_indices))
+
+    return L_indices
+
+class RipsLayer(tf.keras.layers.Layer):
+    """
+    TensorFlow layer for computing Rips persistence out of a point cloud
+    """
+    def __init__(self, homology_dimensions, maximum_edge_length=np.inf, min_persistence=None, homology_coeff_field=11, **kwargs):
+        """
+        Constructor for the RipsLayer class
+
+        Parameters:
+            maximum_edge_length (float): maximum edge length for the Rips complex 
+            homology_dimensions (List[int]): list of homology dimensions
+            min_persistence (List[float]): minimum distance-to-diagonal of the points in the output persistence diagrams (default None, in which case 0. is used for all dimensions)
+            homology_coeff_field (int): homology field coefficient. Must be a prime number. Default value is 11. Max is 46337.
+        """
+        super().__init__(dynamic=True, **kwargs)
+        self.max_edge = maximum_edge_length
+        self.dimensions = homology_dimensions
+        self.min_persistence = min_persistence if min_persistence != None else [0. for _ in range(len(self.dimensions))]
+        self.hcf = homology_coeff_field
+        assert len(self.min_persistence) == len(self.dimensions)
+        
+    def call(self, X):
+        """
+        Compute Rips persistence diagram associated to a point cloud
+
+        Parameters:   
+            X (TensorFlow variable): point cloud of shape [number of points, number of dimensions]
+
+        Returns:
+            List[Tuple[tf.Tensor,tf.Tensor]]: List of Rips persistence diagrams. The length of this list is the same than that of dimensions, i.e., there is one persistence diagram per homology dimension provided in the input list dimensions. Moreover, the finite and essential parts of the persistence diagrams are provided separately: each element of this list is a tuple of size two that contains the finite and essential parts of the corresponding persistence diagram, of shapes [num_finite_points, 2] and [num_essential_points, 1] respectively
+        """    
+        # Compute distance matrix
+        DX = tf.norm(tf.expand_dims(X, 1)-tf.expand_dims(X, 0), axis=2)
+        # Compute vertices associated to positive and negative simplices 
+        # Don't compute gradient for this operation
+        indices = _Rips(DX.numpy(), self.max_edge, self.dimensions, self.hcf)
+        # Get persistence diagrams by simply picking the corresponding entries in the distance matrix
+        self.dgms = []
+        for idx_dim, dimension in enumerate(self.dimensions):
+            cur_idx = indices[idx_dim]
+            if dimension > 0:
+                finite_dgm = tf.reshape(tf.gather_nd(DX, tf.reshape(cur_idx[0], [-1,2])), [-1,2])
+                essential_dgm = tf.reshape(tf.gather_nd(DX, tf.reshape(cur_idx[1], [-1,2])), [-1,1])
+            else:
+                reshaped_cur_idx = tf.reshape(cur_idx[0], [-1,3])
+                finite_dgm = tf.concat([tf.zeros([reshaped_cur_idx.shape[0],1]), tf.reshape(tf.gather_nd(DX, reshaped_cur_idx[:,1:]), [-1,1])], axis=1)
+                essential_dgm = tf.zeros([cur_idx[1].shape[0],1])
+            min_pers = self.min_persistence[idx_dim]
+            if min_pers >= 0:
+                persistent_indices = tf.where(tf.math.abs(finite_dgm[:,1]-finite_dgm[:,0]) > min_pers)
+                self.dgms.append((tf.reshape(tf.gather(finite_dgm, indices=persistent_indices),[-1,2]), essential_dgm))
+            else:
+                self.dgms.append((finite_dgm, essential_dgm))
+        return self.dgms
+
diff --git a/src/python/gudhi/wasserstein/barycenter.py b/src/python/gudhi/wasserstein/barycenter.py
index d67bcde7..bb6e641e 100644
--- a/src/python/gudhi/wasserstein/barycenter.py
+++ b/src/python/gudhi/wasserstein/barycenter.py
@@ -37,7 +37,7 @@ def lagrangian_barycenter(pdiagset, init=None, verbose=False):
     :param init: The initial value for barycenter estimate.
         If ``None``, init is made on a random diagram from the dataset.
         Otherwise, it can be an ``int`` (then initialization is made on ``pdiagset[init]``)
-        or a `(n x 2)` ``numpy.array`` enconding a persistence diagram with `n` points.
+        or a `(n x 2)` ``numpy.array`` encoding a persistence diagram with `n` points.
     :type init: ``int``, or (n x 2) ``np.array``
     :param verbose: if ``True``, returns additional information about the barycenter.
     :type verbose: boolean
@@ -45,7 +45,7 @@ def lagrangian_barycenter(pdiagset, init=None, verbose=False):
         (local minimum of the energy function).
         If ``pdiagset`` is empty, returns ``None``.
         If verbose, returns a couple ``(Y, log)`` where ``Y`` is the barycenter estimate,
-        and ``log`` is a ``dict`` that contains additional informations:
+        and ``log`` is a ``dict`` that contains additional information:
 
         - `"groupings"`, a list of list of pairs ``(i,j)``. Namely, ``G[k] = [...(i, j)...]``, where ``(i,j)`` indicates that `pdiagset[k][i]`` is matched to ``Y[j]`` if ``i = -1`` or ``j = -1``, it means they represent the diagonal.
 
@@ -73,7 +73,7 @@ def lagrangian_barycenter(pdiagset, init=None, verbose=False):
 
     nb_iter = 0
 
-    converged = False  # stoping criterion
+    converged = False  # stopping criterion
     while not converged:
         nb_iter += 1
         K = len(Y)  # current nb of points in Y (some might be on diagonal)
diff --git a/src/python/gudhi/weighted_rips_complex.py b/src/python/gudhi/weighted_rips_complex.py
index 0541572b..16f63c3d 100644
--- a/src/python/gudhi/weighted_rips_complex.py
+++ b/src/python/gudhi/weighted_rips_complex.py
@@ -12,9 +12,11 @@ from gudhi import SimplexTree
 class WeightedRipsComplex:
     """
     Class to generate a weighted Rips complex from a distance matrix and weights on vertices, 
-    in the way described in :cite:`dtmfiltrations`.
+    in the way described in :cite:`dtmfiltrations` with `p=1`. The filtration value of vertex `i` is `2*weights[i]`,
+    and the filtration value of edge `ij` is `distance_matrix[i][j]+weights[i]+weights[j]`,
+    or the maximum of the filtrations of its extremities, whichever is largest.
     Remark that all the filtration values are doubled compared to the definition in the paper 
-    for the consistency with RipsComplex.
+    for consistency with RipsComplex.
     """
     def __init__(self, 
                 distance_matrix,