Merge master and resolve conflicts (mainly src/python/CMakeLists.txt)

24 files changed, 634 insertions, 158 deletions

diff --git a/src/python/CMakeLists.txt b/src/python/CMakeLists.txt
index 179117d2..309d85d3 100644
--- a/src/python/CMakeLists.txt
+++ b/src/python/CMakeLists.txt
@@ -41,14 +41,15 @@ function( add_gudhi_debug_info DEBUG_INFO )
 endfunction( add_gudhi_debug_info )
 
 if(PYTHONINTERP_FOUND)
-  if(PYBIND11_FOUND)
+  if(PYBIND11_FOUND AND CYTHON_FOUND)
     add_gudhi_debug_info("Pybind11 version ${PYBIND11_VERSION}")
+    # PyBind11 modules
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'bottleneck', ")
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'hera', ")
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'clustering', ")
     set(GUDHI_PYTHON_MODULES_EXTRA "${GUDHI_PYTHON_MODULES_EXTRA}'datasets', ")
-  endif()
-  if(CYTHON_FOUND)
+
+    # Cython modules
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'off_reader', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'simplex_tree', ")
     set(GUDHI_PYTHON_MODULES "${GUDHI_PYTHON_MODULES}'rips_complex', ")
@@ -162,7 +163,7 @@ if(PYTHONINTERP_FOUND)
     set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'hera/wasserstein', ")
     set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'hera/bottleneck', ")
     if (NOT CGAL_VERSION VERSION_LESS 4.11.0)
-      set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'datasets/generators/points', ")
+      set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'datasets/generators/_points', ")
       set(GUDHI_PYBIND11_MODULES "${GUDHI_PYBIND11_MODULES}'bottleneck', ")
       set(GUDHI_CYTHON_MODULES "${GUDHI_CYTHON_MODULES}'nerve_gic', ")
     endif ()
@@ -278,8 +279,10 @@ if(PYTHONINTERP_FOUND)
     file(COPY "gudhi/hera/__init__.py" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi/hera")
     file(COPY "gudhi/datasets" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/gudhi" FILES_MATCHING PATTERN "*.py")
 
+
     # Some files for pip package
     file(COPY "introduction.rst" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/")
+    file(COPY "pyproject.toml" DESTINATION "${CMAKE_CURRENT_BINARY_DIR}/")
 
     add_custom_command(
         OUTPUT gudhi.so
@@ -300,35 +303,30 @@ if(PYTHONINTERP_FOUND)
             if(SCIPY_FOUND)
               if(SKLEARN_FOUND)
                 if(OT_FOUND)
-                  if(PYBIND11_FOUND)
-                    if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
-                      set (GUDHI_SPHINX_MESSAGE "Generating API documentation with Sphinx in ${CMAKE_CURRENT_BINARY_DIR}/sphinx/")
-                      # User warning - Sphinx is a static pages generator, and configured to work fine with user_version
-                      # Images and biblio warnings because not found on developper version
-                      if (GUDHI_PYTHON_PATH STREQUAL "src/python")
-                        set (GUDHI_SPHINX_MESSAGE "${GUDHI_SPHINX_MESSAGE} \n WARNING : Sphinx is configured for user version, you run it on developper version. Images and biblio will miss")
-                      endif()
-                      # sphinx target requires gudhi.so, because conf.py reads gudhi version from it
-                      add_custom_target(sphinx
-                          WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/doc
-                          COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}"
-                          ${SPHINX_PATH} -b html ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/sphinx
-                          DEPENDS "${CMAKE_CURRENT_BINARY_DIR}/gudhi.so"
-                          COMMENT "${GUDHI_SPHINX_MESSAGE}" VERBATIM)
-                      add_test(NAME sphinx_py_test
-                               WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
-                               COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}"
-                               ${SPHINX_PATH} -b doctest ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/doctest)
-                      # Set missing or not modules
-                      set(GUDHI_MODULES ${GUDHI_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MODULES")
-                    else(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
-                      message("++ Python documentation module will not be compiled because it requires a Eigen3 and CGAL version >= 5.1.0")
-                      set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES")
-                    endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
-                  else(PYBIND11_FOUND)
-                    message("++ Python documentation module will not be compiled because pybind11 was not found")
+                  if(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
+                    set (GUDHI_SPHINX_MESSAGE "Generating API documentation with Sphinx in ${CMAKE_CURRENT_BINARY_DIR}/sphinx/")
+                    # User warning - Sphinx is a static pages generator, and configured to work fine with user_version
+                    # Images and biblio warnings because not found on developper version
+                    if (GUDHI_PYTHON_PATH STREQUAL "src/python")
+                      set (GUDHI_SPHINX_MESSAGE "${GUDHI_SPHINX_MESSAGE} \n WARNING : Sphinx is configured for user version, you run it on developper version. Images and biblio will miss")
+                    endif()
+                    # sphinx target requires gudhi.so, because conf.py reads gudhi version from it
+                    add_custom_target(sphinx
+                        WORKING_DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR}/doc
+                        COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}"
+                        ${SPHINX_PATH} -b html ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/sphinx
+                        DEPENDS "${CMAKE_CURRENT_BINARY_DIR}/gudhi.so"
+                        COMMENT "${GUDHI_SPHINX_MESSAGE}" VERBATIM)
+                    add_test(NAME sphinx_py_test
+                             WORKING_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}
+                             COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}"
+                             ${SPHINX_PATH} -b doctest ${CMAKE_CURRENT_SOURCE_DIR}/doc ${CMAKE_CURRENT_BINARY_DIR}/doctest)
+                    # Set missing or not modules
+                    set(GUDHI_MODULES ${GUDHI_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MODULES")
+                  else(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
+                    message("++ Python documentation module will not be compiled because it requires a Eigen3 and CGAL version >= 4.11.0")
                     set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES")
-                  endif(PYBIND11_FOUND)
+                  endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
                 else(OT_FOUND)
                   message("++ Python documentation module will not be compiled because POT was not found")
                   set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python-documentation" CACHE INTERNAL "GUDHI_MISSING_MODULES")
@@ -371,7 +369,7 @@ if(PYTHONINTERP_FOUND)
                COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}"
                ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/alpha_rips_persistence_bottleneck_distance.py"
                -f ${CMAKE_SOURCE_DIR}/data/points/tore3D_300.off -t 0.15 -d 3)
-    endif(NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
+    endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 5.1.0)
     if (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
       # Tangential
       add_test(NAME tangential_complex_plain_homology_from_off_file_example_py_test
@@ -406,9 +404,7 @@ if(PYTHONINTERP_FOUND)
                COMMAND ${CMAKE_COMMAND} -E env "${GUDHI_PYTHON_PATH_ENV}"
                ${PYTHON_EXECUTABLE} "${CMAKE_CURRENT_SOURCE_DIR}/example/bottleneck_basic_example.py")
 
-      if (PYBIND11_FOUND)
-        add_gudhi_py_test(test_bottleneck_distance)
-      endif()
+      add_gudhi_py_test(test_bottleneck_distance)
 
       # Cover complex
       file(COPY ${CMAKE_SOURCE_DIR}/data/points/human.off DESTINATION ${CMAKE_CURRENT_BINARY_DIR}/)
@@ -464,6 +460,9 @@ if(PYTHONINTERP_FOUND)
       # Euclidean witness
       add_gudhi_py_test(test_euclidean_witness_complex)
 
+      # Datasets generators
+      add_gudhi_py_test(test_datasets_generators) # TODO separate full python datasets generators in another test file independant from CGAL ?
+
     endif (NOT CGAL_WITH_EIGEN3_VERSION VERSION_LESS 4.11.0)
 
     # Cubical
@@ -522,14 +521,14 @@ if(PYTHONINTERP_FOUND)
     add_gudhi_py_test(test_reader_utils)
 
     # Wasserstein
-    if(OT_FOUND AND PYBIND11_FOUND)
+    if(OT_FOUND)
       # EagerPy dependency because of enable_autodiff=True
       if(EAGERPY_FOUND)
         add_gudhi_py_test(test_wasserstein_distance)
       endif()
+
       add_gudhi_py_test(test_wasserstein_barycenter)
-    endif()
-    if(OT_FOUND)
+
       if(TORCH_FOUND AND TENSORFLOW_FOUND AND EAGERPY_FOUND)
         add_gudhi_py_test(test_wasserstein_with_tensors)
       endif()
@@ -550,7 +549,7 @@ if(PYTHONINTERP_FOUND)
     endif()
 
     # Tomato
-    if(SCIPY_FOUND AND SKLEARN_FOUND AND PYBIND11_FOUND)
+    if(SCIPY_FOUND AND SKLEARN_FOUND)
       add_gudhi_py_test(test_tomato)
     endif()
 
@@ -567,11 +566,11 @@ if(PYTHONINTERP_FOUND)
 
     # Set missing or not modules
     set(GUDHI_MODULES ${GUDHI_MODULES} "python" CACHE INTERNAL "GUDHI_MODULES")
-  else(CYTHON_FOUND)
-    message("++ Python module will not be compiled because cython was not found")
+  else(PYBIND11_FOUND AND CYTHON_FOUND)
+    message("++ Python module will not be compiled because cython and/or pybind11 was/were not found")
     set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python" CACHE INTERNAL "GUDHI_MISSING_MODULES")
-  endif(CYTHON_FOUND)
+  endif(PYBIND11_FOUND AND CYTHON_FOUND)
 else(PYTHONINTERP_FOUND)
   message("++ Python module will not be compiled because no Python interpreter was found")
   set(GUDHI_MISSING_MODULES ${GUDHI_MISSING_MODULES} "python" CACHE INTERNAL "GUDHI_MISSING_MODULES")
-endif(PYTHONINTERP_FOUND)
+endif(PYTHONINTERP_FOUND)
\ No newline at end of file
diff --git a/src/python/doc/alpha_complex_user.rst b/src/python/doc/alpha_complex_user.rst
index 84932c18..94c54053 100644
--- a/src/python/doc/alpha_complex_user.rst
+++ b/src/python/doc/alpha_complex_user.rst
@@ -159,7 +159,10 @@ As the squared radii computed by CGAL are an approximation, it might happen that
 :math:`\alpha^2` values do not quite define a proper filtration (i.e. non-decreasing with
 respect to inclusion).
 We fix that up by calling :func:`~gudhi.SimplexTree.make_filtration_non_decreasing` (cf.
-`C++ version <http://gudhi.gforge.inria.fr/doc/latest/index.html>`_).
+`C++ version <https://gudhi.inria.fr/doc/latest/class_gudhi_1_1_simplex__tree.html>`_).
+
+.. note::
+    This is not the case in `exact` version, this is the reason why it is not called in this case.
 
 Prune above given filtration value
 ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
diff --git a/src/python/doc/datasets_generators.inc b/src/python/doc/datasets_generators.inc
new file mode 100644
index 00000000..8d169275
--- /dev/null
+++ b/src/python/doc/datasets_generators.inc
@@ -0,0 +1,14 @@
+.. table::
+   :widths: 30 40 30
+
+   +-----------------------------------+--------------------------------------------+--------------------------------------------------------------------------------------+
+   | .. figure::                       | Datasets generators (points).              | :Authors: Hind Montassif                                                             |
+   |     img/sphere_3d.png             |                                            |                                                                                      |
+   |                                   |                                            | :Since: GUDHI 3.5.0                                                                  |
+   |                                   |                                            |                                                                                      |
+   |                                   |                                            | :License: MIT (`LGPL v3 </licensing/>`_)                                             |
+   |                                   |                                            |                                                                                      |
+   |                                   |                                            | :Requires: `CGAL <installation.html#cgal>`_                                          |
+   +-----------------------------------+--------------------------------------------+--------------------------------------------------------------------------------------+
+   | * :doc:`datasets_generators`                                                                                                                                          |
+   +-----------------------------------+-----------------------------------------------------------------------------------------------------------------------------------+
diff --git a/src/python/doc/datasets_generators.rst b/src/python/doc/datasets_generators.rst
new file mode 100644
index 00000000..6f36bce1
--- /dev/null
+++ b/src/python/doc/datasets_generators.rst
@@ -0,0 +1,105 @@
+
+:orphan:
+
+.. To get rid of WARNING: document isn't included in any toctree
+
+===========================
+Datasets generators manual
+===========================
+
+We provide the generation of different customizable datasets to use as inputs for Gudhi complexes and data structures.
+
+
+Points generators
+------------------
+
+The module **points** enables the generation of random points on a sphere, random points on a torus and as a grid.
+
+Points on sphere
+^^^^^^^^^^^^^^^^
+
+The function **sphere** enables the generation of random i.i.d. points uniformly on a (d-1)-sphere in :math:`R^d`.
+The user should provide the number of points to be generated on the sphere :code:`n_samples` and the ambient dimension :code:`ambient_dim`.
+The :code:`radius` of sphere is optional and is equal to **1** by default.
+Only random points generation is currently available.
+
+The generated points are given as an array of shape :math:`(n\_samples, ambient\_dim)`.
+
+Example
+"""""""
+
+.. code-block:: python
+
+   from gudhi.datasets.generators import points
+   from gudhi import AlphaComplex
+
+   # Generate 50 points on a sphere in R^2
+   gen_points = points.sphere(n_samples = 50, ambient_dim = 2, radius = 1, sample = "random")
+
+   # Create an alpha complex from the generated points
+   alpha_complex = AlphaComplex(points = gen_points)
+   
+.. autofunction:: gudhi.datasets.generators.points.sphere
+
+Points on a flat torus
+^^^^^^^^^^^^^^^^
+
+You can also generate points on a torus.
+
+Two functions are available and give the same output: the first one depends on **CGAL** and the second does not and consists of full python code.
+
+On another hand, two sample types are provided: you can either generate i.i.d. points on a d-torus in :math:`R^{2d}` *randomly* or on a *grid*.
+
+First function: **ctorus**
+"""""""""""""""""""""""""""
+
+The user should provide the number of points to be generated on the torus :code:`n_samples`, and the dimension :code:`dim` of the torus on which points would be generated in :math:`R^{2dim}`.
+The :code:`sample` argument is optional and is set to **'random'** by default.
+In this case, the returned generated points would be an array of shape :math:`(n\_samples, 2*dim)`.
+Otherwise, if set to **'grid'**, the points are generated on a grid and would be given as an array of shape:
+
+.. math::
+
+   ( ⌊n\_samples^{1 \over {dim}}⌋^{dim}, 2*dim )
+
+**Note 1:** The output array first shape is rounded down to the closest perfect :math:`dim^{th}` power.
+
+**Note 2:** This version is recommended when the user wishes to use **'grid'** as sample type, or **'random'** with a relatively small number of samples (~ less than 150).
+
+Example
+"""""""
+.. code-block:: python
+
+   from gudhi.datasets.generators import points
+
+   # Generate 50 points randomly on a torus in R^6
+   gen_points = points.ctorus(n_samples = 50, dim = 3)
+   
+   # Generate 27 points on a torus as a grid in R^6
+   gen_points = points.ctorus(n_samples = 50, dim = 3, sample = 'grid')
+
+.. autofunction:: gudhi.datasets.generators.points.ctorus
+
+Second function: **torus**
+"""""""""""""""""""""""""""
+
+The user should provide the number of points to be generated on the torus :code:`n_samples` and the dimension :code:`dim` of the torus on which points would be generated in :math:`R^{2dim}`.
+The :code:`sample` argument is optional and is set to **'random'** by default.
+The other allowed value of sample type is **'grid'**.
+
+**Note:** This version is recommended when the user wishes to use **'random'** as sample type with a great number of samples and a low dimension.
+
+Example
+"""""""
+.. code-block:: python
+
+   from gudhi.datasets.generators import points
+
+   # Generate 50 points randomly on a torus in R^6
+   gen_points = points.torus(n_samples = 50, dim = 3)
+   
+   # Generate 27 points on a torus as a grid in R^6
+   gen_points = points.torus(n_samples = 50, dim = 3, sample = 'grid')
+
+
+.. autofunction:: gudhi.datasets.generators.points.torus
diff --git a/src/python/doc/examples.rst b/src/python/doc/examples.rst
index 76e5d4c7..1442f185 100644
--- a/src/python/doc/examples.rst
+++ b/src/python/doc/examples.rst
@@ -8,6 +8,7 @@ Examples
 .. only:: builder_html
 
     * :download:`alpha_complex_diagram_persistence_from_off_file_example.py <../example/alpha_complex_diagram_persistence_from_off_file_example.py>`
+    * :download:`alpha_complex_from_generated_points_on_sphere_example.py <../example/alpha_complex_from_generated_points_on_sphere_example.py>`
     * :download:`alpha_complex_from_points_example.py <../example/alpha_complex_from_points_example.py>`
     * :download:`alpha_rips_persistence_bottleneck_distance.py <../example/alpha_rips_persistence_bottleneck_distance.py>`
     * :download:`bottleneck_basic_example.py <../example/bottleneck_basic_example.py>`
diff --git a/src/python/doc/img/sphere_3d.png b/src/python/doc/img/sphere_3d.png
new file mode 100644
index 00000000..70f3184f
--- /dev/null
+++ b/src/python/doc/img/sphere_3d.png
diff --git a/src/python/doc/index.rst b/src/python/doc/index.rst
index 040e57a4..2d7921ae 100644
--- a/src/python/doc/index.rst
+++ b/src/python/doc/index.rst
@@ -91,3 +91,8 @@ Clustering
 **********
 
 .. include:: clustering.inc
+
+Datasets generators
+*******************
+
+.. include:: datasets_generators.inc
diff --git a/src/python/doc/installation.rst b/src/python/doc/installation.rst
index cb7700ce..35c344e3 100644
--- a/src/python/doc/installation.rst
+++ b/src/python/doc/installation.rst
@@ -359,7 +359,7 @@ Python Optimal Transport
 ------------------------
 
 The :doc:`Wasserstein distance </wasserstein_distance_user>`
-module requires `POT <https://pot.readthedocs.io/>`_, a library that provides
+module requires `POT <https://pythonot.github.io/>`_, a library that provides
 several solvers for optimization problems related to Optimal Transport.
 
 PyTorch
diff --git a/src/python/example/alpha_complex_from_generated_points_on_sphere_example.py b/src/python/example/alpha_complex_from_generated_points_on_sphere_example.py
index 267e6436..3558077e 100644
--- a/src/python/example/alpha_complex_from_generated_points_on_sphere_example.py
+++ b/src/python/example/alpha_complex_from_generated_points_on_sphere_example.py
@@ -1,6 +1,6 @@
 #!/usr/bin/env python
 
-from gudhi.datasets.generators import points
+from gudhi.datasets.generators import _points
 from gudhi import AlphaComplex
 
 
@@ -22,7 +22,7 @@ print("#####################################################################")
 print("AlphaComplex creation from generated points on sphere")
 
 
-gen_points = points.sphere(n_samples = 50, ambient_dim = 2, radius = 1, sample = "random")
+gen_points = _points.sphere(n_samples = 50, ambient_dim = 2, radius = 1, sample = "random")
 
 # Create an alpha complex
 alpha_complex = AlphaComplex(points = gen_points)
diff --git a/src/python/gudhi/cubical_complex.pyx b/src/python/gudhi/cubical_complex.pyx
index 97c69a2d..8e244bb8 100644
--- a/src/python/gudhi/cubical_complex.pyx
+++ b/src/python/gudhi/cubical_complex.pyx
@@ -281,4 +281,8 @@ cdef class CubicalComplex:
             launched first.
         """
         assert self.pcohptr != NULL, "compute_persistence() must be called before persistence_intervals_in_dimension()"
-        return np.array(self.pcohptr.intervals_in_dimension(dimension))
+        piid = np.array(self.pcohptr.intervals_in_dimension(dimension))
+        # Workaround https://github.com/GUDHI/gudhi-devel/issues/507
+        if len(piid) == 0:
+            return np.empty(shape = [0, 2])
+        return piid
diff --git a/src/python/gudhi/datasets/generators/_points.cc b/src/python/gudhi/datasets/generators/_points.cc
new file mode 100644
index 00000000..70ce4925
--- /dev/null
+++ b/src/python/gudhi/datasets/generators/_points.cc
@@ -0,0 +1,120 @@
+/*    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
+ */
+
+#include <pybind11/pybind11.h>
+#include <pybind11/numpy.h>
+
+#include <gudhi/random_point_generators.h>
+#include <gudhi/Debug_utils.h>
+
+#include <CGAL/Epick_d.h>
+
+namespace py = pybind11;
+
+
+typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kern;
+
+py::array_t<double> generate_points_on_sphere(size_t n_samples, int ambient_dim, double radius, std::string sample) {
+
+    if (sample != "random") {
+        throw pybind11::value_error("This sample type is not supported");
+    }
+
+    py::array_t<double> points({n_samples, (size_t)ambient_dim});
+
+    py::buffer_info buf = points.request();
+    double *ptr = static_cast<double *>(buf.ptr);
+
+    GUDHI_CHECK(n_samples == buf.shape[0], "Py array first dimension not matching n_samples on sphere");
+    GUDHI_CHECK(ambient_dim == buf.shape[1], "Py array second dimension not matching the ambient space dimension");
+
+
+    std::vector<typename Kern::Point_d> points_generated;
+
+    {
+        py::gil_scoped_release release;
+        points_generated = Gudhi::generate_points_on_sphere_d<Kern>(n_samples, ambient_dim, radius);
+    }
+
+    for (size_t i = 0; i < n_samples; i++)
+        for (int j = 0; j < ambient_dim; j++)
+            ptr[i*ambient_dim+j] = points_generated[i][j];
+
+    return points;
+}
+
+py::array_t<double> generate_points_on_torus(size_t n_samples, int dim, std::string sample) {
+
+    if ( (sample != "random") && (sample != "grid")) {
+        throw pybind11::value_error("This sample type is not supported");
+    }
+
+    std::vector<typename Kern::Point_d> points_generated;
+
+    {
+        py::gil_scoped_release release;
+        points_generated = Gudhi::generate_points_on_torus_d<Kern>(n_samples, dim, sample);
+    }
+
+    size_t npoints = points_generated.size();
+
+    GUDHI_CHECK(2*dim == points_generated[0].size(), "Py array second dimension not matching the double torus dimension");
+
+    py::array_t<double> points({npoints, (size_t)2*dim});
+
+    py::buffer_info buf = points.request();
+    double *ptr = static_cast<double *>(buf.ptr);
+
+    for (size_t i = 0; i < npoints; i++)
+        for (int j = 0; j < 2*dim; j++)
+            ptr[i*(2*dim)+j] = points_generated[i][j];
+
+    return points;
+}
+
+PYBIND11_MODULE(_points, m) {
+    m.attr("__license__") = "LGPL v3";
+
+    m.def("sphere", &generate_points_on_sphere,
+          py::arg("n_samples"), py::arg("ambient_dim"),
+          py::arg("radius") = 1., py::arg("sample") = "random",
+          R"pbdoc(
+          Generate random i.i.d. points uniformly on a (d-1)-sphere in R^d
+
+          :param n_samples: The number of points to be generated.
+          :type n_samples: integer
+          :param ambient_dim: The ambient dimension d.
+          :type ambient_dim: integer
+          :param radius: The radius. Default value is `1.`.
+          :type radius: float
+          :param sample: The sample type. Default and only available value is `"random"`.
+          :type sample: string
+          :rtype: numpy array of float
+          :returns: the generated points on a sphere.
+          )pbdoc");
+
+    m.def("ctorus", &generate_points_on_torus,
+          py::arg("n_samples"), py::arg("dim"), py::arg("sample") = "random",
+          R"pbdoc(
+          Generate random i.i.d. points on a d-torus in R^2d or as a grid
+
+          :param n_samples: The number of points to be generated.
+          :type n_samples: integer
+          :param dim: The dimension of the torus on which points would be generated in R^2*dim.
+          :type dim: integer
+          :param sample: The sample type. Available values are: `"random"` and `"grid"`. Default value is `"random"`.
+          :type sample: string
+          :rtype: numpy array of float.
+          The shape of returned numpy array is :
+              if sample is 'random' : (n_samples, 2*dim).
+              if sample is 'grid' : (⌊n_samples**(1./dim)⌋**dim, 2*dim), where shape[0] is rounded down to the closest perfect 'dim'th power.
+          :returns: the generated points on a torus.
+          )pbdoc");
+}
diff --git a/src/python/gudhi/datasets/generators/points.cc b/src/python/gudhi/datasets/generators/points.cc
deleted file mode 100644
index d658946b..00000000
--- a/src/python/gudhi/datasets/generators/points.cc
+++ /dev/null
@@ -1,68 +0,0 @@
-/*    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
- */
-
-#include <pybind11/pybind11.h>
-#include <pybind11/numpy.h>
-
-#include <gudhi/random_point_generators.h>
-#include <gudhi/Debug_utils.h>
-
-#include <CGAL/Epick_d.h>
-
-namespace py = pybind11;
-
-
-typedef CGAL::Epick_d< CGAL::Dynamic_dimension_tag > Kern;
-
-py::array_t<double> generate_points_on_sphere(size_t n_samples, int ambient_dim, double radius, std::string sample) {
-    
-    if (sample != "random") {
-        throw pybind11::value_error("This sample type is not supported");
-    }
-
-    py::array_t<double> points({n_samples, (size_t)ambient_dim});
-
-    py::buffer_info buf = points.request();
-    double *ptr = static_cast<double *>(buf.ptr);
-
-    GUDHI_CHECK(n_samples == buf.shape[0], "Py array first dimension not matching n_samples on sphere");
-    GUDHI_CHECK(ambient_dim == buf.shape[1], "Py array second dimension not matching the ambient space dimension");
-
-
-    py::gil_scoped_release release;
-    auto points_generated = Gudhi::generate_points_on_sphere_d<Kern>(n_samples, ambient_dim, radius);
-
-    for (size_t i = 0; i < n_samples; i++)
-        for (int j = 0; j < ambient_dim; j++)
-            ptr[i*ambient_dim+j] = points_generated[i][j];
-
-    return points;
-}
-
-PYBIND11_MODULE(points, m) {
-      m.attr("__license__") = "LGPL v3";
-      m.def("sphere", &generate_points_on_sphere,
-          py::arg("n_samples"), py::arg("ambient_dim"),
-          py::arg("radius") = 1., py::arg("sample") = "random",
-          R"pbdoc(
-    Generate random i.i.d. points uniformly on a (d-1)-sphere in R^d
-
-    :param n_samples: The number of points to be generated.
-    :type n_samples: integer
-    :param ambient_dim: The ambient dimension d.
-    :type ambient_dim: integer
-    :param radius: The radius. Default value is `1.`.
-    :type radius: float
-    :param sample: The sample type. Default and only available value is `"random"`.
-    :type sample: string
-    :rtype: numpy array of float
-    :returns: the generated points on a sphere.
-    )pbdoc");
-}
diff --git a/src/python/gudhi/datasets/generators/points.py b/src/python/gudhi/datasets/generators/points.py
new file mode 100644
index 00000000..cf97777d
--- /dev/null
+++ b/src/python/gudhi/datasets/generators/points.py
@@ -0,0 +1,56 @@
+# 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
+
+import numpy as np
+
+from ._points import ctorus
+from ._points import sphere
+
+def _generate_random_points_on_torus(n_samples, dim):
+
+    # Generate random angles of size n_samples*dim
+    alpha = 2*np.pi*np.random.rand(n_samples*dim)
+
+    # Based on angles, construct points of size n_samples*dim on a circle and reshape the result in a n_samples*2*dim array
+    array_points = np.column_stack([np.cos(alpha), np.sin(alpha)]).reshape(-1, 2*dim)
+    
+    return array_points
+
+def _generate_grid_points_on_torus(n_samples, dim):
+    
+    # Generate points on a dim-torus as a grid
+    n_samples_grid = int((n_samples+.5)**(1./dim)) # add .5 to avoid rounding down with numerical approximations
+    alpha = np.linspace(0, 2*np.pi, n_samples_grid, endpoint=False)
+    
+    array_points = np.column_stack([np.cos(alpha), np.sin(alpha)])
+    array_points_idx = np.empty([n_samples_grid]*dim + [dim], dtype=int)
+    for i, x in enumerate(np.ix_(*([np.arange(n_samples_grid)]*dim))):
+        array_points_idx[...,i] = x
+    return array_points[array_points_idx].reshape(-1, 2*dim)
+
+def torus(n_samples, dim, sample='random'):
+    """ 
+    Generate points on a flat dim-torus in R^2dim either randomly or on a grid
+    
+    :param n_samples: The number of points to be generated.
+    :param dim: The dimension of the torus on which points would be generated in R^2*dim.
+    :param sample: The sample type of the generated points. Can be 'random' or 'grid'.
+    :returns: numpy array containing the generated points on a torus.
+        The shape of returned numpy array is:
+        if sample is 'random' : (n_samples, 2*dim).
+        if sample is 'grid' : (⌊n_samples**(1./dim)⌋**dim, 2*dim), where shape[0] is rounded down to the closest perfect 'dim'th power.
+    """
+    if sample == 'random':
+        # Generate points randomly
+        return _generate_random_points_on_torus(n_samples, dim)
+    elif sample == 'grid':
+        # Generate points on a grid
+        return _generate_grid_points_on_torus(n_samples, dim)
+    else:
+        raise ValueError("Sample type '{}' is not supported".format(sample))
diff --git a/src/python/gudhi/periodic_cubical_complex.pyx b/src/python/gudhi/periodic_cubical_complex.pyx
index ef1d3080..6c21e902 100644
--- a/src/python/gudhi/periodic_cubical_complex.pyx
+++ b/src/python/gudhi/periodic_cubical_complex.pyx
@@ -280,4 +280,8 @@ cdef class PeriodicCubicalComplex:
             launched first.
         """
         assert self.pcohptr != NULL, "compute_persistence() must be called before persistence_intervals_in_dimension()"
-        return np.array(self.pcohptr.intervals_in_dimension(dimension))
+        piid = np.array(self.pcohptr.intervals_in_dimension(dimension))
+        # Workaround https://github.com/GUDHI/gudhi-devel/issues/507
+        if len(piid) == 0:
+            return np.empty(shape = [0, 2])
+        return piid
diff --git a/src/python/gudhi/point_cloud/knn.py b/src/python/gudhi/point_cloud/knn.py
index 829bf1bf..de5844f9 100644
--- a/src/python/gudhi/point_cloud/knn.py
+++ b/src/python/gudhi/point_cloud/knn.py
@@ -8,6 +8,7 @@
 #   - YYYY/MM Author: Description of the modification
 
 import numpy
+import warnings
 
 # TODO: https://github.com/facebookresearch/faiss
 
@@ -257,6 +258,9 @@ class KNearestNeighbors:
             if ef is not None:
                 self.graph.set_ef(ef)
             neighbors, distances = self.graph.knn_query(X, k, num_threads=self.params["num_threads"])
+            with warnings.catch_warnings():
+                if not(numpy.all(numpy.isfinite(distances))):
+                    warnings.warn("Overflow/infinite value encountered while computing 'distances'", RuntimeWarning)
             # The k nearest neighbors are always sorted. I couldn't find it in the doc, but the code calls searchKnn,
             # which returns a priority_queue, and then fills the return array backwards with top/pop on the queue.
             if self.return_index:
@@ -290,6 +294,9 @@ class KNearestNeighbors:
             if self.return_index:
                 if self.return_distance:
                     distances, neighbors = mat.Kmin_argKmin(k, dim=1)
+                    with warnings.catch_warnings():
+                        if not(torch.isfinite(distances).all()):
+                            warnings.warn("Overflow/infinite value encountered while computing 'distances'", RuntimeWarning)
                     if p != numpy.inf:
                         distances = distances ** (1.0 / p)
                     return neighbors, distances
@@ -298,6 +305,9 @@ class KNearestNeighbors:
                     return neighbors
             if self.return_distance:
                 distances = mat.Kmin(k, dim=1)
+                with warnings.catch_warnings():
+                    if not(torch.isfinite(distances).all()):
+                        warnings.warn("Overflow/infinite value encountered while computing 'distances'", RuntimeWarning)
                 if p != numpy.inf:
                     distances = distances ** (1.0 / p)
                 return distances
diff --git a/src/python/gudhi/representations/vector_methods.py b/src/python/gudhi/representations/vector_methods.py
index 84bc99a2..e883b5dd 100644
--- a/src/python/gudhi/representations/vector_methods.py
+++ b/src/python/gudhi/representations/vector_methods.py
@@ -6,6 +6,7 @@
 #
 # Modification(s):
 #   - 2020/06 Martin: ATOL integration
+#   - 2021/11 Vincent Rouvreau: factorize _automatic_sample_range
 
 import numpy as np
 from sklearn.base          import BaseEstimator, TransformerMixin
@@ -45,10 +46,14 @@ class PersistenceImage(BaseEstimator, TransformerMixin):
             y (n x 1 array): persistence diagram labels (unused).
         """
         if np.isnan(np.array(self.im_range)).any():
-            new_X = BirthPersistenceTransform().fit_transform(X)
-            pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(new_X,y)
-            [mx,my],[Mx,My] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]], [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]]
-            self.im_range = np.where(np.isnan(np.array(self.im_range)), np.array([mx, Mx, my, My]), np.array(self.im_range))
+            try:
+                new_X = BirthPersistenceTransform().fit_transform(X)
+                pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(new_X,y)
+                [mx,my],[Mx,My] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]], [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]]
+                self.im_range = np.where(np.isnan(np.array(self.im_range)), np.array([mx, Mx, my, My]), np.array(self.im_range))
+            except ValueError:
+                # Empty persistence diagram case - https://github.com/GUDHI/gudhi-devel/issues/507
+                pass
         return self
 
     def transform(self, X):
@@ -94,6 +99,28 @@ class PersistenceImage(BaseEstimator, TransformerMixin):
         """
         return self.fit_transform([diag])[0,:]
 
+def _automatic_sample_range(sample_range, X, y):
+        """
+        Compute and returns sample range from the persistence diagrams if one of the sample_range values is numpy.nan.
+
+        Parameters:
+            sample_range (a numpy array of 2 float): minimum and maximum of all piecewise-linear function domains, of
+                the form [x_min, x_max].
+            X (list of n x 2 numpy arrays): input persistence diagrams.
+            y (n x 1 array): persistence diagram labels (unused).
+        """
+        nan_in_range = np.isnan(sample_range)
+        if nan_in_range.any():
+            try:
+                pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(X,y)
+                [mx,my] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]]
+                [Mx,My] = [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]]
+                return np.where(nan_in_range, np.array([mx, My]), sample_range)
+            except ValueError:
+                # Empty persistence diagram case - https://github.com/GUDHI/gudhi-devel/issues/507
+                pass
+        return sample_range
+
 class Landscape(BaseEstimator, TransformerMixin):
     """
     This is a class for computing persistence landscapes from a list of persistence diagrams. A persistence landscape is a collection of 1D piecewise-linear functions computed from the rank function associated to the persistence diagram. These piecewise-linear functions are then sampled evenly on a given range and the corresponding vectors of samples are concatenated and returned. See http://jmlr.org/papers/v16/bubenik15a.html for more details.
@@ -119,10 +146,7 @@ class Landscape(BaseEstimator, TransformerMixin):
             X (list of n x 2 numpy arrays): input persistence diagrams.
             y (n x 1 array): persistence diagram labels (unused).
         """
-        if self.nan_in_range.any():
-            pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(X,y)
-            [mx,my],[Mx,My] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]], [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]]
-            self.sample_range = np.where(self.nan_in_range, np.array([mx, My]), np.array(self.sample_range))
+        self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y)
         return self
 
     def transform(self, X):
@@ -218,10 +242,7 @@ class Silhouette(BaseEstimator, TransformerMixin):
             X (list of n x 2 numpy arrays): input persistence diagrams.
             y (n x 1 array): persistence diagram labels (unused).
         """
-        if np.isnan(np.array(self.sample_range)).any():
-            pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(X,y)
-            [mx,my],[Mx,My] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]], [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]]
-            self.sample_range = np.where(np.isnan(np.array(self.sample_range)), np.array([mx, My]), np.array(self.sample_range))
+        self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y)
         return self
 
     def transform(self, X):
@@ -307,10 +328,7 @@ class BettiCurve(BaseEstimator, TransformerMixin):
             X (list of n x 2 numpy arrays): input persistence diagrams.
             y (n x 1 array): persistence diagram labels (unused).
         """
-        if np.isnan(np.array(self.sample_range)).any():
-            pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(X,y)
-            [mx,my],[Mx,My] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]], [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]]
-            self.sample_range = np.where(np.isnan(np.array(self.sample_range)), np.array([mx, My]), np.array(self.sample_range))
+        self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y)
         return self
 
     def transform(self, X):
@@ -374,10 +392,7 @@ class Entropy(BaseEstimator, TransformerMixin):
             X (list of n x 2 numpy arrays): input persistence diagrams.
             y (n x 1 array): persistence diagram labels (unused).
         """
-        if np.isnan(np.array(self.sample_range)).any():
-            pre = DiagramScaler(use=True, scalers=[([0], MinMaxScaler()), ([1], MinMaxScaler())]).fit(X,y)
-            [mx,my],[Mx,My] = [pre.scalers[0][1].data_min_[0], pre.scalers[1][1].data_min_[0]], [pre.scalers[0][1].data_max_[0], pre.scalers[1][1].data_max_[0]]
-            self.sample_range = np.where(np.isnan(np.array(self.sample_range)), np.array([mx, My]), np.array(self.sample_range))
+        self.sample_range = _automatic_sample_range(np.array(self.sample_range), X, y)
         return self
 
     def transform(self, X):
@@ -396,9 +411,13 @@ 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]
-            new_diagram = DiagramScaler(use=True, scalers=[([1], MaxAbsScaler())]).fit_transform([diagram])[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])
 
             if self.mode == "scalar":
                 ent = - np.sum( np.multiply(new_diagram[:,1], np.log(new_diagram[:,1])) )
@@ -412,12 +431,11 @@ class Entropy(BaseEstimator, TransformerMixin):
                     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])
-                    if self.normalized:
-                        ent = ent / np.linalg.norm(ent, ord=1)
-                    Xfit.append(np.reshape(ent,[1,-1]))
-
-        Xfit = np.concatenate(Xfit, 0)
+                if self.normalized:
+                    ent = ent / np.linalg.norm(ent, ord=1)
+                Xfit.append(np.reshape(ent,[1,-1]))
 
+        Xfit = np.concatenate(Xfit, axis=0)
         return Xfit
 
     def __call__(self, diag):
@@ -478,7 +496,13 @@ class TopologicalVector(BaseEstimator, TransformerMixin):
             diagram, num_pts_in_diag = X[i], X[i].shape[0]
             pers = 0.5 * (diagram[:,1]-diagram[:,0])
             min_pers = np.minimum(pers,np.transpose(pers))
-            distances = DistanceMetric.get_metric("chebyshev").pairwise(diagram)
+            # Works fine with sklearn 1.0, but an ValueError exception is thrown on past versions
+            try:
+                distances = DistanceMetric.get_metric("chebyshev").pairwise(diagram)
+            except ValueError:
+                # Empty persistence diagram case - https://github.com/GUDHI/gudhi-devel/issues/507
+                assert len(diagram) == 0
+                distances = np.empty(shape = [0, 0])
             vect = np.flip(np.sort(np.triu(np.minimum(distances, min_pers)), axis=None), 0)
             dim = min(len(vect), thresh)
             Xfit[i, :dim] = vect[:dim]
diff --git a/src/python/gudhi/simplex_tree.pyx b/src/python/gudhi/simplex_tree.pyx
index 9c51cb46..c3720936 100644
--- a/src/python/gudhi/simplex_tree.pyx
+++ b/src/python/gudhi/simplex_tree.pyx
@@ -9,8 +9,7 @@
 
 from cython.operator import dereference, preincrement
 from libc.stdint cimport intptr_t
-import numpy
-from numpy import array as np_array
+import numpy as np
 cimport gudhi.simplex_tree
 
 __author__ = "Vincent Rouvreau"
@@ -542,7 +541,11 @@ cdef class SimplexTree:
             function to be launched first.
         """
         assert self.pcohptr != NULL, "compute_persistence() must be called before persistence_intervals_in_dimension()"
-        return np_array(self.pcohptr.intervals_in_dimension(dimension))
+        piid = np.array(self.pcohptr.intervals_in_dimension(dimension))
+        # Workaround https://github.com/GUDHI/gudhi-devel/issues/507
+        if len(piid) == 0:
+            return np.empty(shape = [0, 2])
+        return piid
 
     def persistence_pairs(self):
         """This function returns a list of persistence birth and death simplices pairs.
@@ -583,8 +586,8 @@ cdef class SimplexTree:
         """
         assert self.pcohptr != NULL, "lower_star_persistence_generators() requires that persistence() be called first."
         gen = self.pcohptr.lower_star_generators()
-        normal = [np_array(d).reshape(-1,2) for d in gen.first]
-        infinite = [np_array(d) for d in gen.second]
+        normal = [np.array(d).reshape(-1,2) for d in gen.first]
+        infinite = [np.array(d) for d in gen.second]
         return (normal, infinite)
 
     def flag_persistence_generators(self):
@@ -602,19 +605,19 @@ cdef class SimplexTree:
         assert self.pcohptr != NULL, "flag_persistence_generators() requires that persistence() be called first."
         gen = self.pcohptr.flag_generators()
         if len(gen.first) == 0:
-            normal0 = numpy.empty((0,3))
+            normal0 = np.empty((0,3))
             normals = []
         else:
             l = iter(gen.first)
-            normal0 = np_array(next(l)).reshape(-1,3)
-            normals = [np_array(d).reshape(-1,4) for d in l]
+            normal0 = np.array(next(l)).reshape(-1,3)
+            normals = [np.array(d).reshape(-1,4) for d in l]
         if len(gen.second) == 0:
-            infinite0 = numpy.empty(0)
+            infinite0 = np.empty(0)
             infinites = []
         else:
             l = iter(gen.second)
-            infinite0 = np_array(next(l))
-            infinites = [np_array(d).reshape(-1,2) for d in l]
+            infinite0 = np.array(next(l))
+            infinites = [np.array(d).reshape(-1,2) for d in l]
         return (normal0, normals, infinite0, infinites)
 
     def collapse_edges(self, nb_iterations = 1):
diff --git a/src/python/pyproject.toml b/src/python/pyproject.toml
new file mode 100644
index 00000000..a9fb4985
--- /dev/null
+++ b/src/python/pyproject.toml
@@ -0,0 +1,3 @@
+[build-system]
+requires = ["setuptools", "wheel", "numpy>=1.15.0", "cython", "pybind11"]
+build-backend = "setuptools.build_meta"
diff --git a/src/python/setup.py.in b/src/python/setup.py.in
index c400b601..23746998 100644
--- a/src/python/setup.py.in
+++ b/src/python/setup.py.in
@@ -82,10 +82,10 @@ setup(
         },
     description='The Gudhi library is an open source library for ' \
         'Computational Topology and Topological Data Analysis (TDA).',
+    data_files=[('.', ['./introduction.rst'])],
     long_description_content_type='text/x-rst',
     long_description=long_description,
     ext_modules = ext_modules,
     install_requires = ['numpy >= 1.15.0',],
-    setup_requires = ['cython','numpy >= 1.15.0','pybind11',],
     package_data={"": ["*.dll"], },
 )
diff --git a/src/python/test/test_cubical_complex.py b/src/python/test/test_cubical_complex.py
index d0e4e9e8..29d559b3 100755
--- a/src/python/test/test_cubical_complex.py
+++ b/src/python/test/test_cubical_complex.py
@@ -174,3 +174,28 @@ def test_periodic_cofaces_of_persistence_pairs_when_pd_has_no_paired_birth_and_d
     assert np.array_equal(pairs[1][0], np.array([0]))
     assert np.array_equal(pairs[1][1], np.array([0, 1]))
     assert np.array_equal(pairs[1][2], np.array([1]))
+
+def test_cubical_persistence_intervals_in_dimension():
+    cub = CubicalComplex(
+        dimensions=[3, 3],
+        top_dimensional_cells=[1, 2, 3, 4, 5, 6, 7, 8, 9],
+    )
+    cub.compute_persistence()
+    H0 = cub.persistence_intervals_in_dimension(0)
+    assert np.array_equal(H0, np.array([[ 1., float("inf")]]))
+    assert cub.persistence_intervals_in_dimension(1).shape == (0, 2)
+
+def test_periodic_cubical_persistence_intervals_in_dimension():
+    cub = PeriodicCubicalComplex(
+        dimensions=[3, 3],
+        top_dimensional_cells=[1, 2, 3, 4, 5, 6, 7, 8, 9],
+        periodic_dimensions = [True, True]
+    )
+    cub.compute_persistence()
+    H0 = cub.persistence_intervals_in_dimension(0)
+    assert np.array_equal(H0, np.array([[ 1., float("inf")]]))
+    H1 = cub.persistence_intervals_in_dimension(1)
+    assert np.array_equal(H1, np.array([[ 3., float("inf")], [ 7., float("inf")]]))
+    H2 = cub.persistence_intervals_in_dimension(2)
+    assert np.array_equal(H2, np.array([[ 9., float("inf")]]))
+    assert cub.persistence_intervals_in_dimension(3).shape == (0, 2)
diff --git a/src/python/test/test_datasets_generators.py b/src/python/test/test_datasets_generators.py
new file mode 100755
index 00000000..91ec4a65
--- /dev/null
+++ b/src/python/test/test_datasets_generators.py
@@ -0,0 +1,39 @@
+""" 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 gudhi.datasets.generators import points
+
+import pytest
+
+def test_sphere():
+    assert points.sphere(n_samples = 10, ambient_dim = 2, radius = 1., sample = 'random').shape == (10, 2)
+
+    with pytest.raises(ValueError):
+        points.sphere(n_samples = 10, ambient_dim = 2, radius = 1., sample = 'other')
+
+def _basic_torus(impl):
+    assert impl(n_samples = 64, dim = 3, sample = 'random').shape == (64, 6)
+    assert impl(n_samples = 64, dim = 3, sample = 'grid').shape == (64, 6)
+
+    assert impl(n_samples = 10, dim = 4, sample = 'random').shape == (10, 8)
+
+    # Here 1**dim < n_samples < 2**dim, the output shape is therefore (1, 2*dim) = (1, 8), where shape[0] is rounded down to the closest perfect 'dim'th power
+    assert impl(n_samples = 10, dim = 4, sample = 'grid').shape == (1, 8)
+
+    with pytest.raises(ValueError):
+        impl(n_samples = 10, dim = 4, sample = 'other')
+
+def test_torus():
+    for torus_impl in [points.torus, points.ctorus]:
+        _basic_torus(torus_impl)
+    # Check that the two versions (torus and ctorus) generate the same output
+    assert points.ctorus(n_samples = 64, dim = 3, sample = 'random').all() == points.torus(n_samples = 64, dim = 3, sample = 'random').all()
+    assert points.ctorus(n_samples = 64, dim = 3, sample = 'grid').all() == points.torus(n_samples = 64, dim = 3, sample = 'grid').all()
+    assert points.ctorus(n_samples = 10, dim = 3, sample = 'grid').all() == points.torus(n_samples = 10, dim = 3, sample = 'grid').all()
diff --git a/src/python/test/test_dtm.py b/src/python/test/test_dtm.py
index 0a52279e..c29471cf 100755
--- a/src/python/test/test_dtm.py
+++ b/src/python/test/test_dtm.py
@@ -13,6 +13,7 @@ import numpy
 import pytest
 import torch
 import math
+import warnings
 
 
 def test_dtm_compare_euclidean():
@@ -87,3 +88,16 @@ def test_density():
     assert density == pytest.approx(expected)
     density = DTMDensity(weights=[0.5, 0.5], metric="neighbors", dim=1).fit_transform(distances)
     assert density == pytest.approx(expected)
+
+def test_dtm_overflow_warnings():
+    pts = numpy.array([[10., 100000000000000000000000000000.], [1000., 100000000000000000000000000.]])
+    impl_warn = ["keops", "hnsw"]
+
+    with warnings.catch_warnings(record=True) as w:
+        for impl in impl_warn:
+            dtm = DistanceToMeasure(2, q=10000, implementation=impl)
+            r = dtm.fit_transform(pts)
+        assert len(w) == 2
+        for i in range(len(w)):
+            assert issubclass(w[i].category, RuntimeWarning)
+            assert "Overflow" in str(w[i].message)
diff --git a/src/python/test/test_representations.py b/src/python/test/test_representations.py
index cda1a15b..93461f1e 100755
--- a/src/python/test/test_representations.py
+++ b/src/python/test/test_representations.py
@@ -3,9 +3,23 @@ import sys
 import matplotlib.pyplot as plt
 import numpy as np
 import pytest
+import random
 
 from sklearn.cluster import KMeans
 
+# Vectorization
+from gudhi.representations import (Landscape, Silhouette, BettiCurve, ComplexPolynomial,\
+  TopologicalVector, PersistenceImage, Entropy)
+
+# Preprocessing
+from gudhi.representations import (BirthPersistenceTransform, Clamping, DiagramScaler, Padding, ProminentPoints, \
+  DiagramSelector)
+
+# Kernel
+from gudhi.representations import (PersistenceWeightedGaussianKernel, \
+  PersistenceScaleSpaceKernel, SlicedWassersteinDistance,\
+  SlicedWassersteinKernel, PersistenceFisherKernel, WassersteinDistance)
+
 
 def test_representations_examples():
     # Disable graphics for testing purposes
@@ -98,3 +112,60 @@ def test_infinity():
     assert c[1] == 0
     assert c[7] == 3
     assert c[9] == 2
+
+
+def test_preprocessing_empty_diagrams():
+    empty_diag = np.empty(shape = [0, 2])
+    assert not np.any(BirthPersistenceTransform()(empty_diag))
+    assert not np.any(Clamping().fit_transform(empty_diag))
+    assert not np.any(DiagramScaler()(empty_diag))
+    assert not np.any(Padding()(empty_diag))
+    assert not np.any(ProminentPoints()(empty_diag))
+    assert not np.any(DiagramSelector()(empty_diag))
+
+def pow(n):
+  return lambda x: np.power(x[1]-x[0],n)
+
+def test_vectorization_empty_diagrams():
+    empty_diag = np.empty(shape = [0, 2])
+    random_resolution = random.randint(50,100)*10 # between 500 and 1000
+    print("resolution = ", random_resolution)
+    lsc = Landscape(resolution=random_resolution)(empty_diag)
+    assert not np.any(lsc)
+    assert lsc.shape[0]%random_resolution == 0
+    slt = Silhouette(resolution=random_resolution, weight=pow(2))(empty_diag)
+    assert not np.any(slt)
+    assert slt.shape[0] == random_resolution
+    btc = BettiCurve(resolution=random_resolution)(empty_diag)
+    assert not np.any(btc)
+    assert btc.shape[0] == random_resolution
+    cpp = ComplexPolynomial(threshold=random_resolution, polynomial_type="T")(empty_diag)
+    assert not np.any(cpp)
+    assert cpp.shape[0] == random_resolution
+    tpv = TopologicalVector(threshold=random_resolution)(empty_diag)
+    assert tpv.shape[0] == random_resolution
+    assert not np.any(tpv)
+    prmg = PersistenceImage(resolution=[random_resolution,random_resolution])(empty_diag)
+    assert not np.any(prmg)
+    assert prmg.shape[0] == random_resolution * random_resolution
+    sce = Entropy(mode="scalar", resolution=random_resolution)(empty_diag)
+    assert not np.any(sce)
+    assert sce.shape[0] == 1
+    scv = Entropy(mode="vector", normalized=False, resolution=random_resolution)(empty_diag)
+    assert not np.any(scv)
+    assert scv.shape[0] == random_resolution
+
+def test_kernel_empty_diagrams():
+    empty_diag = np.empty(shape = [0, 2])
+    assert SlicedWassersteinDistance(num_directions=100)(empty_diag, empty_diag) == 0.
+    assert SlicedWassersteinKernel(num_directions=100, bandwidth=1.)(empty_diag, empty_diag) == 1.
+    assert WassersteinDistance(mode="hera", delta=0.0001)(empty_diag, empty_diag) == 0.
+    assert WassersteinDistance(mode="pot")(empty_diag, empty_diag) == 0.
+    assert BottleneckDistance(epsilon=.001)(empty_diag, empty_diag) == 0.
+    assert BottleneckDistance()(empty_diag, empty_diag) == 0.
+#    PersistenceWeightedGaussianKernel(bandwidth=1., kernel_approx=None, weight=arctan(1.,1.))(empty_diag, empty_diag)
+#    PersistenceWeightedGaussianKernel(kernel_approx=RBFSampler(gamma=1./2, n_components=100000).fit(np.ones([1,2])), weight=arctan(1.,1.))(empty_diag, empty_diag)
+#    PersistenceScaleSpaceKernel(bandwidth=1.)(empty_diag, empty_diag)
+#    PersistenceScaleSpaceKernel(kernel_approx=RBFSampler(gamma=1./2, n_components=100000).fit(np.ones([1,2])))(empty_diag, empty_diag)
+#    PersistenceFisherKernel(bandwidth_fisher=1., bandwidth=1.)(empty_diag, empty_diag)
+#    PersistenceFisherKernel(bandwidth_fisher=1., bandwidth=1., kernel_approx=RBFSampler(gamma=1./2, n_components=100000).fit(np.ones([1,2])))(empty_diag, empty_diag)
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index a3eacaa9..31c46213 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -9,6 +9,7 @@
 """
 
 from gudhi import SimplexTree, __GUDHI_USE_EIGEN3
+import numpy as np
 import pytest
 
 __author__ = "Vincent Rouvreau"
@@ -404,3 +405,46 @@ def test_boundaries_iterator():
 
     with pytest.raises(RuntimeError):
         list(st.get_boundaries([6])) # (6) does not exist
+
+def test_persistence_intervals_in_dimension():
+    # Here is our triangulation of a 2-torus - taken from https://dioscuri-tda.org/Paris_TDA_Tutorial_2021.html
+    #   0-----3-----4-----0
+    #   | \   | \   | \   | \   |
+    #   |   \ |   \ |    \|   \ | 
+    #   1-----8-----7-----1
+    #   | \   | \   | \   | \   |
+    #   |   \ |   \ |   \ |   \ |
+    #   2-----5-----6-----2
+    #   | \   | \   | \   | \   |
+    #   |   \ |   \ |   \ |   \ |
+    #   0-----3-----4-----0
+    st = SimplexTree()
+    st.insert([0,1,8])
+    st.insert([0,3,8])
+    st.insert([3,7,8])
+    st.insert([3,4,7])
+    st.insert([1,4,7])
+    st.insert([0,1,4])
+    st.insert([1,2,5])
+    st.insert([1,5,8])
+    st.insert([5,6,8])
+    st.insert([6,7,8])
+    st.insert([2,6,7])
+    st.insert([1,2,7])
+    st.insert([0,2,3])
+    st.insert([2,3,5])
+    st.insert([3,4,5])
+    st.insert([4,5,6])
+    st.insert([0,4,6])
+    st.insert([0,2,6])
+    st.compute_persistence(persistence_dim_max=True)
+    
+    H0 = st.persistence_intervals_in_dimension(0)
+    assert np.array_equal(H0, np.array([[ 0., float("inf")]]))
+    H1 = st.persistence_intervals_in_dimension(1)
+    assert np.array_equal(H1, np.array([[ 0., float("inf")], [ 0., float("inf")]]))
+    H2 = st.persistence_intervals_in_dimension(2)
+    assert np.array_equal(H2, np.array([[ 0., float("inf")]]))
+    # Test empty case
+    assert st.persistence_intervals_in_dimension(3).shape == (0, 2)
+    
\ No newline at end of file