Merge branch 'master' into insert

16 files changed, 876 insertions, 217 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 a4888914..375e1561 100644
--- a/src/python/gudhi/alpha_complex.pyx
+++ b/src/python/gudhi/alpha_complex.pyx
@@ -31,6 +31,10 @@ cdef extern from "Alpha_complex_interface.h" namespace "Gudhi":
         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:
@@ -133,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 f8078d03..69ff5e1e 100644
--- a/src/python/gudhi/representations/vector_methods.py
+++ b/src/python/gudhi/representations/vector_methods.py
@@ -508,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 284daa96..f86f1232 100644
--- a/src/python/gudhi/simplex_tree.pxd
+++ b/src/python/gudhi/simplex_tree.pxd
@@ -64,7 +64,6 @@ 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
@@ -82,7 +81,7 @@ cdef extern from "Simplex_tree_interface.h" namespace "Gudhi":
         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
@@ -93,3 +92,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 1ac03afa..6b1b5c00 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -558,8 +558,7 @@ 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
@@ -572,9 +571,9 @@ cdef class SimplexTree:
         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 `block_simplex`, so
-            if you examine the complex in `block_simplex`, you may hit a few simplices of the same dimension that have
-            not been vetted by `block_simplex` yet, or have already been rejected but not yet removed.
+            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
@@ -760,18 +759,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()
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..5df2c370
--- /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 output 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)