Merge remote-tracking branch 'upstream/master' into cech_optimization

9 files changed, 464 insertions, 83 deletions

diff --git a/src/python/test/test_alpha_complex.py b/src/python/test/test_alpha_complex.py
index 814f8289..f15284f3 100755
--- a/src/python/test/test_alpha_complex.py
+++ b/src/python/test/test_alpha_complex.py
@@ -8,10 +8,12 @@
       - YYYY/MM Author: Description of the modification
 """
 
-import gudhi as gd
+from gudhi import AlphaComplex
 import math
 import numpy as np
 import pytest
+import warnings
+
 try:
     # python3
     from itertools import zip_longest
@@ -19,22 +21,24 @@ except ImportError:
     # python2
     from itertools import izip_longest as zip_longest
 
-__author__ = "Vincent Rouvreau"
-__copyright__ = "Copyright (C) 2016 Inria"
-__license__ = "MIT"
 
 
 def _empty_alpha(precision):
-    alpha_complex = gd.AlphaComplex(points=[[0, 0]], precision = precision)
+    alpha_complex = AlphaComplex(precision = precision)
+    assert alpha_complex.__is_defined() == True
+
+def _one_2d_point_alpha(precision):
+    alpha_complex = AlphaComplex(points=[[0, 0]], precision = precision)
     assert alpha_complex.__is_defined() == True
 
 def test_empty_alpha():
     for precision in ['fast', 'safe', 'exact']:
         _empty_alpha(precision)
+        _one_2d_point_alpha(precision)
 
 def _infinite_alpha(precision):
     point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
-    alpha_complex = gd.AlphaComplex(points=point_list, precision = precision)
+    alpha_complex = AlphaComplex(points=point_list, precision = precision)
     assert alpha_complex.__is_defined() == True
 
     simplex_tree = alpha_complex.create_simplex_tree()
@@ -69,18 +73,9 @@ def _infinite_alpha(precision):
     assert point_list[1] == alpha_complex.get_point(1)
     assert point_list[2] == alpha_complex.get_point(2)
     assert point_list[3] == alpha_complex.get_point(3)
-    try:
-        alpha_complex.get_point(4) == []
-    except IndexError:
-        pass
-    else:
-        assert False
-    try:
-        alpha_complex.get_point(125) == []
-    except IndexError:
-        pass
-    else:
-        assert False
+
+    with pytest.raises(IndexError):
+        alpha_complex.get_point(len(point_list))
 
 def test_infinite_alpha():
     for precision in ['fast', 'safe', 'exact']:
@@ -88,7 +83,7 @@ def test_infinite_alpha():
 
 def _filtered_alpha(precision):
     point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
-    filtered_alpha = gd.AlphaComplex(points=point_list, precision = precision)
+    filtered_alpha = AlphaComplex(points=point_list, precision = precision)
 
     simplex_tree = filtered_alpha.create_simplex_tree(max_alpha_square=0.25)
 
@@ -99,18 +94,9 @@ def _filtered_alpha(precision):
     assert point_list[1] == filtered_alpha.get_point(1)
     assert point_list[2] == filtered_alpha.get_point(2)
     assert point_list[3] == filtered_alpha.get_point(3)
-    try:
-        filtered_alpha.get_point(4) == []
-    except IndexError:
-        pass
-    else:
-        assert False
-    try:
-        filtered_alpha.get_point(125) == []
-    except IndexError:
-        pass
-    else:
-        assert False
+
+    with pytest.raises(IndexError):
+        filtered_alpha.get_point(len(point_list))
 
     assert list(simplex_tree.get_filtration()) == [
         ([0], 0.0),
@@ -141,10 +127,10 @@ def _safe_alpha_persistence_comparison(precision):
     embedding2 = [[signal[i], delayed[i]] for i in range(len(time))]
     
     #build alpha complex and simplex tree
-    alpha_complex1 = gd.AlphaComplex(points=embedding1, precision = precision)
+    alpha_complex1 = AlphaComplex(points=embedding1, precision = precision)
     simplex_tree1 = alpha_complex1.create_simplex_tree()
     
-    alpha_complex2 = gd.AlphaComplex(points=embedding2, precision = precision)
+    alpha_complex2 = AlphaComplex(points=embedding2, precision = precision)
     simplex_tree2 = alpha_complex2.create_simplex_tree()
     
     diag1 = simplex_tree1.persistence()
@@ -162,7 +148,7 @@ def test_safe_alpha_persistence_comparison():
 
 def _delaunay_complex(precision):
     point_list = [[0, 0], [1, 0], [0, 1], [1, 1]]
-    filtered_alpha = gd.AlphaComplex(points=point_list, precision = precision)
+    filtered_alpha = AlphaComplex(points=point_list, precision = precision)
 
     simplex_tree = filtered_alpha.create_simplex_tree(default_filtration_value = True)
 
@@ -173,18 +159,11 @@ def _delaunay_complex(precision):
     assert point_list[1] == filtered_alpha.get_point(1)
     assert point_list[2] == filtered_alpha.get_point(2)
     assert point_list[3] == filtered_alpha.get_point(3)
-    try:
-        filtered_alpha.get_point(4) == []
-    except IndexError:
-        pass
-    else:
-        assert False
-    try:
-        filtered_alpha.get_point(125) == []
-    except IndexError:
-        pass
-    else:
-        assert False
+
+    with pytest.raises(IndexError):
+        filtered_alpha.get_point(4)
+    with pytest.raises(IndexError):
+        filtered_alpha.get_point(125)
 
     for filtered_value in simplex_tree.get_filtration():
         assert math.isnan(filtered_value[1])
@@ -198,7 +177,13 @@ def test_delaunay_complex():
         _delaunay_complex(precision)
 
 def _3d_points_on_a_plane(precision, default_filtration_value):
-    alpha = gd.AlphaComplex(off_file='alphacomplexdoc.off', precision = precision)
+    alpha = AlphaComplex(points = [[1.0, 1.0 , 0.0],
+                                      [7.0, 0.0 , 0.0],
+                                      [4.0, 6.0 , 0.0],
+                                      [9.0, 6.0 , 0.0],
+                                      [0.0, 14.0, 0.0],
+                                      [2.0, 19.0, 0.0],
+                                      [9.0, 17.0, 0.0]], precision = precision)
 
     simplex_tree = alpha.create_simplex_tree(default_filtration_value = default_filtration_value)
     assert simplex_tree.dimension() == 2
@@ -206,28 +191,16 @@ def _3d_points_on_a_plane(precision, default_filtration_value):
     assert simplex_tree.num_simplices() == 25
 
 def test_3d_points_on_a_plane():
-    off_file = open("alphacomplexdoc.off", "w")
-    off_file.write("OFF         \n" \
-                   "7 0 0       \n" \
-                   "1.0 1.0  0.0\n" \
-                   "7.0 0.0  0.0\n" \
-                   "4.0 6.0  0.0\n" \
-                   "9.0 6.0  0.0\n" \
-                   "0.0 14.0 0.0\n" \
-                   "2.0 19.0 0.0\n" \
-                   "9.0 17.0 0.0\n"  )
-    off_file.close()
-
     for default_filtration_value in [True, False]:
         for precision in ['fast', 'safe', 'exact']:
             _3d_points_on_a_plane(precision, default_filtration_value)
 
 def _3d_tetrahedrons(precision):
     points = 10*np.random.rand(10, 3)
-    alpha = gd.AlphaComplex(points=points, precision = precision)
+    alpha = AlphaComplex(points = points, precision = precision)
     st_alpha = alpha.create_simplex_tree(default_filtration_value = False)
     # New AlphaComplex for get_point to work
-    delaunay = gd.AlphaComplex(points=points, precision = precision)
+    delaunay = AlphaComplex(points = points, precision = precision)
     st_delaunay = delaunay.create_simplex_tree(default_filtration_value = True)
 
     delaunay_tetra = []
@@ -256,3 +229,60 @@ def _3d_tetrahedrons(precision):
 def test_3d_tetrahedrons():
     for precision in ['fast', 'safe', 'exact']:
         _3d_tetrahedrons(precision)
+
+def test_off_file_deprecation_warning():
+    off_file = open("alphacomplexdoc.off", "w")
+    off_file.write("OFF         \n" \
+                   "7 0 0       \n" \
+                   "1.0 1.0  0.0\n" \
+                   "7.0 0.0  0.0\n" \
+                   "4.0 6.0  0.0\n" \
+                   "9.0 6.0  0.0\n" \
+                   "0.0 14.0 0.0\n" \
+                   "2.0 19.0 0.0\n" \
+                   "9.0 17.0 0.0\n"  )
+    off_file.close()
+
+    with pytest.warns(DeprecationWarning):
+        alpha = AlphaComplex(off_file="alphacomplexdoc.off")
+
+def test_non_existing_off_file():
+    with pytest.warns(DeprecationWarning):
+        with pytest.raises(FileNotFoundError):
+            alpha = AlphaComplex(off_file="pouetpouettralala.toubiloubabdou")
+
+def test_inconsistency_points_and_weights():
+    points = [[1.0, 1.0 , 0.0],
+              [7.0, 0.0 , 0.0],
+              [4.0, 6.0 , 0.0],
+              [9.0, 6.0 , 0.0],
+              [0.0, 14.0, 0.0],
+              [2.0, 19.0, 0.0],
+              [9.0, 17.0, 0.0]]
+    with pytest.raises(ValueError):
+        # 7 points, 8 weights, on purpose
+        alpha = AlphaComplex(points = points,
+                                weights = [1., 2., 3., 4., 5., 6., 7., 8.])
+
+    with pytest.raises(ValueError):
+        # 7 points, 6 weights, on purpose
+        alpha = AlphaComplex(points = points,
+                                weights = [1., 2., 3., 4., 5., 6.])
+
+def _weighted_doc_example(precision):
+    stree = AlphaComplex(points=[[ 1., -1., -1.],
+                                    [-1.,  1., -1.],
+                                    [-1., -1.,  1.],
+                                    [ 1.,  1.,  1.],
+                                    [ 2.,  2.,  2.]],
+                            weights = [4., 4., 4., 4., 1.],
+                            precision = precision).create_simplex_tree()
+
+    assert stree.filtration([0, 1, 2, 3]) == pytest.approx(-1.)
+    assert stree.filtration([0, 1, 3, 4]) == pytest.approx(95.)
+    assert stree.filtration([0, 2, 3, 4]) == pytest.approx(95.)
+    assert stree.filtration([1, 2, 3, 4]) == pytest.approx(95.)
+
+def test_weighted_doc_example():
+    for precision in ['fast', 'safe', 'exact']:
+        _weighted_doc_example(precision)
diff --git a/src/python/test/test_betti_curve_representations.py b/src/python/test/test_betti_curve_representations.py
new file mode 100755
index 00000000..6a45da4d
--- /dev/null
+++ b/src/python/test/test_betti_curve_representations.py
@@ -0,0 +1,59 @@
+import numpy as np
+import scipy.interpolate
+import pytest
+
+from gudhi.representations.vector_methods import BettiCurve
+
+def test_betti_curve_is_irregular_betti_curve_followed_by_interpolation():
+    m = 10
+    n = 1000
+    pinf = 0.05
+    pzero = 0.05
+    res = 100
+    
+    pds = []
+    for i in range(0, m):
+        pd = np.zeros((n, 2))
+        pd[:, 0] = np.random.uniform(0, 10, n)
+        pd[:, 1] = np.random.uniform(pd[:, 0], 10, n)
+        pd[np.random.uniform(0, 1, n) < pzero, 0] = 0
+        pd[np.random.uniform(0, 1, n) < pinf, 1] = np.inf
+        pds.append(pd)
+    
+    bc = BettiCurve(resolution=None, predefined_grid=None)
+    bc.fit(pds)
+    bettis = bc.transform(pds)
+
+    bc2 = BettiCurve(resolution=None, predefined_grid=None)
+    bettis2 = bc2.fit_transform(pds)
+    assert((bc2.grid_ == bc.grid_).all())
+    assert((bettis2 == bettis).all())
+
+    for i in range(0, m):
+        grid = np.linspace(pds[i][np.isfinite(pds[i])].min(), pds[i][np.isfinite(pds[i])].max() + 1, res)
+        bc_gridded = BettiCurve(predefined_grid=grid)
+        bc_gridded.fit([])
+        bettis_gridded = bc_gridded(pds[i])
+
+        interp = scipy.interpolate.interp1d(bc.grid_, bettis[i, :], kind="previous", fill_value="extrapolate")
+        bettis_interp = np.array(interp(grid), dtype=int)
+        assert((bettis_interp == bettis_gridded).all())
+
+
+def test_empty_with_predefined_grid():
+    random_grid = np.sort(np.random.uniform(0, 1, 100))
+    bc = BettiCurve(predefined_grid=random_grid)
+    bettis = bc.fit_transform([])
+    assert((bc.grid_ == random_grid).all())
+    assert((bettis == 0).all())
+
+    
+def test_empty():
+    bc = BettiCurve(resolution=None, predefined_grid=None)
+    bettis = bc.fit_transform([])
+    assert(bc.grid_ == [-np.inf])
+    assert((bettis == 0).all())
+
+def test_wrong_value_of_predefined_grid():
+    with pytest.raises(ValueError):
+        BettiCurve(predefined_grid=[1, 2, 3])
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..e46d616c 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,14 @@ 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.]])
+
+    with warnings.catch_warnings(record=True) as w:
+        # TODO Test "keops" implementation as well when next version of pykeops (current is 1.5) is released (should fix the problem (cf. issue #543))
+        dtm = DistanceToMeasure(2, implementation="hnsw")
+        r = dtm.fit_transform(pts)
+        assert len(w) == 1
+        assert issubclass(w[0].category, RuntimeWarning)
+        assert "Overflow" in str(w[0].message)
diff --git a/src/python/test/test_reader_utils.py b/src/python/test/test_reader_utils.py
index 90da6651..fdfddc4b 100755
--- a/src/python/test/test_reader_utils.py
+++ b/src/python/test/test_reader_utils.py
@@ -8,8 +8,9 @@
       - YYYY/MM Author: Description of the modification
 """
 
-import gudhi
+import gudhi as gd
 import numpy as np
+from pytest import raises
 
 __author__ = "Vincent Rouvreau"
 __copyright__ = "Copyright (C) 2017 Inria"
@@ -18,7 +19,7 @@ __license__ = "MIT"
 
 def test_non_existing_csv_file():
     # Try to open a non existing file
-    matrix = gudhi.read_lower_triangular_matrix_from_csv_file(
+    matrix = gd.read_lower_triangular_matrix_from_csv_file(
         csv_file="pouetpouettralala.toubiloubabdou"
     )
     assert matrix == []
@@ -29,8 +30,8 @@ def test_full_square_distance_matrix_csv_file():
     test_file = open("full_square_distance_matrix.csv", "w")
     test_file.write("0;1;2;3;\n1;0;4;5;\n2;4;0;6;\n3;5;6;0;")
     test_file.close()
-    matrix = gudhi.read_lower_triangular_matrix_from_csv_file(
-        csv_file="full_square_distance_matrix.csv"
+    matrix = gd.read_lower_triangular_matrix_from_csv_file(
+        csv_file="full_square_distance_matrix.csv", separator=";"
     )
     assert matrix == [[], [1.0], [2.0, 4.0], [3.0, 5.0, 6.0]]
 
@@ -40,7 +41,7 @@ def test_lower_triangular_distance_matrix_csv_file():
     test_file = open("lower_triangular_distance_matrix.csv", "w")
     test_file.write("\n1,\n2,3,\n4,5,6,\n7,8,9,10,")
     test_file.close()
-    matrix = gudhi.read_lower_triangular_matrix_from_csv_file(
+    matrix = gd.read_lower_triangular_matrix_from_csv_file(
         csv_file="lower_triangular_distance_matrix.csv", separator=","
     )
     assert matrix == [[], [1.0], [2.0, 3.0], [4.0, 5.0, 6.0], [7.0, 8.0, 9.0, 10.0]]
@@ -48,11 +49,11 @@ def test_lower_triangular_distance_matrix_csv_file():
 
 def test_non_existing_persistence_file():
     # Try to open a non existing file
-    persistence = gudhi.read_persistence_intervals_grouped_by_dimension(
+    persistence = gd.read_persistence_intervals_grouped_by_dimension(
         persistence_file="pouetpouettralala.toubiloubabdou"
     )
     assert persistence == []
-    persistence = gudhi.read_persistence_intervals_in_dimension(
+    persistence = gd.read_persistence_intervals_in_dimension(
         persistence_file="pouetpouettralala.toubiloubabdou", only_this_dim=1
     )
     np.testing.assert_array_equal(persistence, [])
@@ -65,21 +66,21 @@ def test_read_persistence_intervals_without_dimension():
         "# Simple persistence diagram without dimension\n2.7 3.7\n9.6 14.\n34.2 34.974\n3. inf"
     )
     test_file.close()
-    persistence = gudhi.read_persistence_intervals_in_dimension(
+    persistence = gd.read_persistence_intervals_in_dimension(
         persistence_file="persistence_intervals_without_dimension.pers"
     )
     np.testing.assert_array_equal(
         persistence, [(2.7, 3.7), (9.6, 14.0), (34.2, 34.974), (3.0, float("Inf"))]
     )
-    persistence = gudhi.read_persistence_intervals_in_dimension(
+    persistence = gd.read_persistence_intervals_in_dimension(
         persistence_file="persistence_intervals_without_dimension.pers", only_this_dim=0
     )
     np.testing.assert_array_equal(persistence, [])
-    persistence = gudhi.read_persistence_intervals_in_dimension(
+    persistence = gd.read_persistence_intervals_in_dimension(
         persistence_file="persistence_intervals_without_dimension.pers", only_this_dim=1
     )
     np.testing.assert_array_equal(persistence, [])
-    persistence = gudhi.read_persistence_intervals_grouped_by_dimension(
+    persistence = gd.read_persistence_intervals_grouped_by_dimension(
         persistence_file="persistence_intervals_without_dimension.pers"
     )
     assert persistence == {
@@ -94,29 +95,29 @@ def test_read_persistence_intervals_with_dimension():
         "# Simple persistence diagram with dimension\n0 2.7 3.7\n1 9.6 14.\n3 34.2 34.974\n1 3. inf"
     )
     test_file.close()
-    persistence = gudhi.read_persistence_intervals_in_dimension(
+    persistence = gd.read_persistence_intervals_in_dimension(
         persistence_file="persistence_intervals_with_dimension.pers"
     )
     np.testing.assert_array_equal(
         persistence, [(2.7, 3.7), (9.6, 14.0), (34.2, 34.974), (3.0, float("Inf"))]
     )
-    persistence = gudhi.read_persistence_intervals_in_dimension(
+    persistence = gd.read_persistence_intervals_in_dimension(
         persistence_file="persistence_intervals_with_dimension.pers", only_this_dim=0
     )
     np.testing.assert_array_equal(persistence, [(2.7, 3.7)])
-    persistence = gudhi.read_persistence_intervals_in_dimension(
+    persistence = gd.read_persistence_intervals_in_dimension(
         persistence_file="persistence_intervals_with_dimension.pers", only_this_dim=1
     )
     np.testing.assert_array_equal(persistence, [(9.6, 14.0), (3.0, float("Inf"))])
-    persistence = gudhi.read_persistence_intervals_in_dimension(
+    persistence = gd.read_persistence_intervals_in_dimension(
         persistence_file="persistence_intervals_with_dimension.pers", only_this_dim=2
     )
     np.testing.assert_array_equal(persistence, [])
-    persistence = gudhi.read_persistence_intervals_in_dimension(
+    persistence = gd.read_persistence_intervals_in_dimension(
         persistence_file="persistence_intervals_with_dimension.pers", only_this_dim=3
     )
     np.testing.assert_array_equal(persistence, [(34.2, 34.974)])
-    persistence = gudhi.read_persistence_intervals_grouped_by_dimension(
+    persistence = gd.read_persistence_intervals_grouped_by_dimension(
         persistence_file="persistence_intervals_with_dimension.pers"
     )
     assert persistence == {
diff --git a/src/python/test/test_representations.py b/src/python/test/test_representations.py
index cda1a15b..d219ce7a 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
@@ -91,10 +105,66 @@ def test_dummy_atol():
 
 from gudhi.representations.vector_methods import BettiCurve
 
-
 def test_infinity():
     a = np.array([[1.0, 8.0], [2.0, np.inf], [3.0, 4.0]])
     c = BettiCurve(20, [0.0, 10.0])(a)
     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
diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py
index e3b521d6..3a004d77 100755
--- a/src/python/test/test_wasserstein_distance.py
+++ b/src/python/test/test_wasserstein_distance.py
@@ -5,25 +5,97 @@
     Copyright (C) 2019 Inria
 
     Modification(s):
+      - 2020/07 Théo Lacombe: Added tests about handling essential parts in diagrams.
       - YYYY/MM Author: Description of the modification
 """
 
-from gudhi.wasserstein.wasserstein import _proj_on_diag
+from gudhi.wasserstein.wasserstein import _proj_on_diag, _finite_part, _handle_essential_parts, _get_essential_parts
+from gudhi.wasserstein.wasserstein import _warn_infty
 from gudhi.wasserstein import wasserstein_distance as pot
 from gudhi.hera import wasserstein_distance as hera
 import numpy as np
 import pytest
 
+
 __author__ = "Theo Lacombe"
 __copyright__ = "Copyright (C) 2019 Inria"
 __license__ = "MIT"
 
+
 def test_proj_on_diag():
     dgm = np.array([[1., 1.], [1., 2.], [3., 5.]])
     assert np.array_equal(_proj_on_diag(dgm), [[1., 1.], [1.5, 1.5], [4., 4.]])
     empty = np.empty((0, 2))
     assert np.array_equal(_proj_on_diag(empty), empty)
 
+
+def test_finite_part():
+    diag = np.array([[0, 1], [3, 5], [2, np.inf], [3, np.inf], [-np.inf, 8], [-np.inf, 12], [-np.inf, -np.inf],
+                     [np.inf, np.inf], [-np.inf, np.inf], [-np.inf, np.inf]])
+    assert np.array_equal(_finite_part(diag), [[0, 1], [3, 5]])
+
+
+def test_handle_essential_parts():
+    diag1 = np.array([[0, 1], [3, 5],
+                      [2, np.inf], [3, np.inf],
+                      [-np.inf, 8], [-np.inf, 12],
+                      [-np.inf, -np.inf],
+                      [np.inf, np.inf],
+                      [-np.inf, np.inf], [-np.inf, np.inf]])
+
+    diag2 = np.array([[0, 2], [3, 5],
+                      [2, np.inf], [4, np.inf],
+                      [-np.inf, 8], [-np.inf, 11],
+                      [-np.inf, -np.inf],
+                      [np.inf, np.inf],
+                      [-np.inf, np.inf], [-np.inf, np.inf]])
+
+    diag3 = np.array([[0, 2], [3, 5],
+                      [2, np.inf], [4, np.inf], [6, np.inf],
+                      [-np.inf, 8], [-np.inf, 11],
+                      [-np.inf, -np.inf],
+                      [np.inf, np.inf],
+                      [-np.inf, np.inf], [-np.inf, np.inf]])
+
+    c, m = _handle_essential_parts(diag1, diag2, order=1)
+    assert c == pytest.approx(2, 0.0001)  # Note: here c is only the cost due to essential part (thus 2, not 3)
+    # Similarly, the matching only corresponds to essential parts.
+    # Note that (-inf,-inf) and (+inf,+inf) coordinates are matched to the diagonal.
+    assert np.array_equal(m, [[4, 4], [5, 5], [2, 2], [3, 3], [8, 8], [9, 9], [6, -1], [7, -1], [-1, 6], [-1, 7]])
+
+    c, m = _handle_essential_parts(diag1, diag3, order=1)
+    assert c == np.inf
+    assert (m is None)
+
+
+def test_get_essential_parts():
+    diag1 = np.array([[0, 1], [3, 5], [2, np.inf], [3, np.inf], [-np.inf, 8], [-np.inf, 12], [-np.inf, -np.inf],
+                     [np.inf, np.inf], [-np.inf, np.inf], [-np.inf, np.inf]])
+
+    diag2 = np.array([[0, 1], [3, 5], [2, np.inf], [3, np.inf]])
+
+    res  = _get_essential_parts(diag1)
+    res2 = _get_essential_parts(diag2)
+    assert np.array_equal(res[0], [4, 5])
+    assert np.array_equal(res[1], [2, 3])
+    assert np.array_equal(res[2], [8, 9])
+    assert np.array_equal(res[3], [6]   )
+    assert np.array_equal(res[4], [7]   )
+
+    assert np.array_equal(res2[0], []    )
+    assert np.array_equal(res2[1], [2, 3])
+    assert np.array_equal(res2[2], []    )
+    assert np.array_equal(res2[3], []    )
+    assert np.array_equal(res2[4], []    )
+
+
+def test_warn_infty():
+    assert _warn_infty(matching=False)==np.inf
+    c, m = _warn_infty(matching=True)
+    assert (c == np.inf)
+    assert (m is None)
+
+
 def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_matching=True):
     diag1 = np.array([[2.7, 3.7], [9.6, 14.0], [34.2, 34.974]])
     diag2 = np.array([[2.8, 4.45], [9.5, 14.1]])
@@ -64,7 +136,7 @@ def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_mat
 
         assert wasserstein_distance(diag4, diag5) == np.inf
         assert wasserstein_distance(diag5, diag6, order=1, internal_p=np.inf) == approx(4.)
-
+        assert wasserstein_distance(diag5, emptydiag) == np.inf
 
     if test_matching:
         match = wasserstein_distance(emptydiag, emptydiag, matching=True, internal_p=1., order=2)[1]
@@ -78,6 +150,31 @@ def _basic_wasserstein(wasserstein_distance, delta, test_infinity=True, test_mat
         match = wasserstein_distance(diag1, diag2, matching=True, internal_p=2., order=2.)[1]
         assert np.array_equal(match, [[0, 0], [1, 1], [2, -1]])
 
+    if test_matching and test_infinity:
+        diag7 = np.array([[0, 3], [4, np.inf], [5, np.inf]])
+        diag8 = np.array([[0,1], [0, np.inf], [-np.inf, -np.inf], [np.inf, np.inf]])
+        diag9 = np.array([[-np.inf, -np.inf], [np.inf, np.inf]])
+        diag10 = np.array([[0,1], [-np.inf, -np.inf], [np.inf, np.inf]])
+
+        match = wasserstein_distance(diag5, diag6, matching=True, internal_p=2., order=2.)[1]
+        assert np.array_equal(match, [[0, -1], [-1,0], [-1, 1], [1, 2]])
+        match = wasserstein_distance(diag5, diag7, matching=True, internal_p=2., order=2.)[1]
+        assert (match is None)
+        cost, match = wasserstein_distance(diag7, emptydiag, matching=True, internal_p=2., order=2.3)
+        assert (cost == np.inf)
+        assert (match is None)
+        cost, match = wasserstein_distance(emptydiag, diag7, matching=True, internal_p=2.42, order=2.)
+        assert (cost == np.inf)
+        assert (match is None)
+        cost, match = wasserstein_distance(diag8, diag9, matching=True, internal_p=2., order=2.)
+        assert (cost == np.inf)
+        assert (match is None)
+        cost, match = wasserstein_distance(diag9, diag10, matching=True, internal_p=1., order=1.)
+        assert (cost == 1)
+        assert (match == [[0, -1],[1, -1],[-1, 0], [-1, 1], [-1, 2]]) # type 4 and 5 are match to the diag anyway.
+        cost, match = wasserstein_distance(diag9, emptydiag, matching=True, internal_p=2., order=2.)
+        assert (cost == 0.)
+        assert (match == [[0, -1], [1, -1]])
 
 
 def hera_wrap(**extra):
@@ -85,15 +182,19 @@ def hera_wrap(**extra):
         return hera(*kargs,**kwargs,**extra)
     return fun
 
+
 def pot_wrap(**extra):
     def fun(*kargs,**kwargs):
         return pot(*kargs,**kwargs,**extra)
     return fun
 
+
 def test_wasserstein_distance_pot():
-    _basic_wasserstein(pot, 1e-15, test_infinity=False, test_matching=True)
-    _basic_wasserstein(pot_wrap(enable_autodiff=True), 1e-15, test_infinity=False, test_matching=False)
+    _basic_wasserstein(pot, 1e-15, test_infinity=False, test_matching=True)  # pot with its standard args
+    _basic_wasserstein(pot_wrap(enable_autodiff=True, keep_essential_parts=False), 1e-15, test_infinity=False, test_matching=False)
+
 
 def test_wasserstein_distance_hera():
     _basic_wasserstein(hera_wrap(delta=1e-12), 1e-12, test_matching=False)
     _basic_wasserstein(hera_wrap(delta=.1), .1, test_matching=False)
+