Merge master and resolve commits

14 files changed, 715 insertions, 121 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_bottleneck_distance.py b/src/python/test/test_bottleneck_distance.py
index 6915bea8..07fcc9cc 100755
--- a/src/python/test/test_bottleneck_distance.py
+++ b/src/python/test/test_bottleneck_distance.py
@@ -25,3 +25,15 @@ def test_basic_bottleneck():
     assert gudhi.bottleneck_distance(diag1, diag2, 0.1) == pytest.approx(0.75, abs=0.1)
     assert gudhi.hera.bottleneck_distance(diag1, diag2, 0) == 0.75
     assert gudhi.hera.bottleneck_distance(diag1, diag2, 0.1) == pytest.approx(0.75, rel=0.1)
+
+    import numpy as np
+
+    # Translating both diagrams along the diagonal should not affect the distance,
+    # checks that negative numbers are not an issue
+    diag1 = np.array(diag1) - 100
+    diag2 = np.array(diag2) - 100
+
+    assert gudhi.bottleneck_distance(diag1, diag2) == 0.75
+    assert gudhi.bottleneck_distance(diag1, diag2, 0.1) == pytest.approx(0.75, abs=0.1)
+    assert gudhi.hera.bottleneck_distance(diag1, diag2, 0) == 0.75
+    assert gudhi.hera.bottleneck_distance(diag1, diag2, 0.1) == pytest.approx(0.75, rel=0.1)
diff --git a/src/python/test/test_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 e5c211a0..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
@@ -39,11 +53,37 @@ def test_multiple():
     d2 = BottleneckDistance(epsilon=0.00001).fit_transform(l1)
     d3 = pairwise_persistence_diagram_distances(l1, l1b, e=0.00001, n_jobs=4)
     assert d1 == pytest.approx(d2)
-    assert d3 == pytest.approx(d2, abs=1e-5) # Because of 0 entries (on the diagonal)
+    assert d3 == pytest.approx(d2, abs=1e-5)  # Because of 0 entries (on the diagonal)
     d1 = pairwise_persistence_diagram_distances(l1, l2, metric="wasserstein", order=2, internal_p=2)
     d2 = WassersteinDistance(order=2, internal_p=2, n_jobs=4).fit(l2).transform(l1)
     print(d1.shape, d2.shape)
-    assert d1 == pytest.approx(d2, rel=.02)
+    assert d1 == pytest.approx(d2, rel=0.02)
+
+
+# Test sorted values as points order can be inverted, and sorted test is not documentation-friendly
+# Note the test below must be up to date with the Atol class documentation
+def test_atol_doc():
+    a = np.array([[1, 2, 4], [1, 4, 0], [1, 0, 4]])
+    b = np.array([[4, 2, 0], [4, 4, 0], [4, 0, 2]])
+    c = np.array([[3, 2, -1], [1, 2, -1]])
+
+    atol_vectoriser = Atol(quantiser=KMeans(n_clusters=2, random_state=202006))
+    # Atol will do
+    # X = np.concatenate([a,b,c])
+    # kmeans = KMeans(n_clusters=2, random_state=202006).fit(X) 
+    # kmeans.labels_ will be : array([1, 0, 1, 0, 0, 1, 0, 0])
+    first_cluster = np.asarray([a[0], a[2], b[2]])
+    second_cluster = np.asarray([a[1], b[0], b[2], c[0], c[1]])
+
+    # Check the center of the first_cluster and second_cluster are in Atol centers
+    centers = atol_vectoriser.fit(X=[a, b, c]).centers
+    np.isclose(centers, first_cluster.mean(axis=0)).all(1).any() 
+    np.isclose(centers, second_cluster.mean(axis=0)).all(1).any() 
+
+    vectorization = atol_vectoriser.transform(X=[a, b, c])
+    assert np.allclose(vectorization[0], atol_vectoriser(a))
+    assert np.allclose(vectorization[1], atol_vectoriser(b))
+    assert np.allclose(vectorization[2], atol_vectoriser(c))
 
 
 def test_dummy_atol():
@@ -53,8 +93,78 @@ def test_dummy_atol():
 
     for weighting_method in ["cloud", "iidproba"]:
         for contrast in ["gaussian", "laplacian", "indicator"]:
-            atol_vectoriser = Atol(quantiser=KMeans(n_clusters=1, random_state=202006), weighting_method=weighting_method, contrast=contrast)
+            atol_vectoriser = Atol(
+                quantiser=KMeans(n_clusters=1, random_state=202006),
+                weighting_method=weighting_method,
+                contrast=contrast,
+            )
             atol_vectoriser.fit([a, b, c])
             atol_vectoriser(a)
             atol_vectoriser.transform(X=[a, b, c])
 
+
+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_rips_complex.py b/src/python/test/test_rips_complex.py
index b86e7498..a2f43a1b 100755
--- a/src/python/test/test_rips_complex.py
+++ b/src/python/test/test_rips_complex.py
@@ -133,3 +133,24 @@ def test_filtered_rips_from_distance_matrix():
 
     assert simplex_tree.num_simplices() == 8
     assert simplex_tree.num_vertices() == 4
+
+
+def test_sparse_with_multiplicity():
+    points = [
+        [3, 4],
+        [0.1, 2],
+        [0.1, 2],
+        [0.1, 2],
+        [0.1, 2],
+        [0.1, 2],
+        [0.1, 2],
+        [0.1, 2],
+        [0.1, 2],
+        [0.1, 2],
+        [0.1, 2],
+        [3, 4.1],
+    ]
+    rips = RipsComplex(points=points, sparse=0.01)
+    simplex_tree = rips.create_simplex_tree(max_dimension=2)
+    assert simplex_tree.num_simplices() == 7
+    diag = simplex_tree.persistence()
diff --git a/src/python/test/test_simplex_tree.py b/src/python/test/test_simplex_tree.py
index 33b0ac99..d8173b52 100755
--- a/src/python/test/test_simplex_tree.py
+++ b/src/python/test/test_simplex_tree.py
@@ -8,7 +8,8 @@
       - YYYY/MM Author: Description of the modification
 """
 
-from gudhi import SimplexTree
+from gudhi import SimplexTree, __GUDHI_USE_EIGEN3
+import numpy as np
 import pytest
 
 __author__ = "Vincent Rouvreau"
@@ -353,11 +354,167 @@ def test_collapse_edges():
 
     assert st.num_simplices() == 10
 
-    st.collapse_edges()
-    assert st.num_simplices() == 9
-    assert st.find([1, 3]) == False
-    for simplex in st.get_skeleton(0): 
-        assert simplex[1] == 1. 
+    if __GUDHI_USE_EIGEN3:
+        st.collapse_edges()
+        assert st.num_simplices() == 9
+        assert st.find([1, 3]) == False
+        for simplex in st.get_skeleton(0):
+            assert simplex[1] == 1.
+    else:
+        # If no Eigen3, collapse_edges throws an exception
+        with pytest.raises(RuntimeError):
+            st.collapse_edges()
+
+def test_reset_filtration():
+    st = SimplexTree()
+    
+    assert st.insert([0, 1, 2], 3.) == True
+    assert st.insert([0, 3], 2.) == True
+    assert st.insert([3, 4, 5], 3.) == True
+    assert st.insert([0, 1, 6, 7], 4.) == True
+
+    # Guaranteed by construction
+    for simplex in st.get_simplices():
+        assert st.filtration(simplex[0]) >= 2.
+    
+    # dimension until 5 even if simplex tree is of dimension 3 to test the limits
+    for dimension in range(5, -1, -1):
+        st.reset_filtration(0., dimension)
+        for simplex in st.get_skeleton(3):
+            print(simplex)
+            if len(simplex[0]) < (dimension) + 1:
+                assert st.filtration(simplex[0]) >= 2.
+            else:
+                assert st.filtration(simplex[0]) == 0.
+
+def test_boundaries_iterator():
+    st = SimplexTree()
+
+    assert st.insert([0, 1, 2, 3], filtration=1.0) == True
+    assert st.insert([1, 2, 3, 4], filtration=2.0) == True
+
+    assert list(st.get_boundaries([1, 2, 3])) == [([1, 2], 1.0), ([1, 3], 1.0), ([2, 3], 1.0)]
+    assert list(st.get_boundaries([2, 3, 4])) == [([2, 3], 1.0), ([2, 4], 2.0), ([3, 4], 2.0)]
+    assert list(st.get_boundaries([2])) == []
+
+    with pytest.raises(RuntimeError):
+        list(st.get_boundaries([]))
+
+    with pytest.raises(RuntimeError):
+        list(st.get_boundaries([0, 4])) # (0, 4) does not exist
+
+    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)
+
+def test_equality_operator():
+    st1 = SimplexTree()
+    st2 = SimplexTree()
+
+    assert st1 == st2
+
+    st1.insert([1,2,3], 4.)
+    assert st1 != st2
+
+    st2.insert([1,2,3], 4.)
+    assert st1 == st2
+
+def test_simplex_tree_deep_copy():
+    st = SimplexTree()
+    st.insert([1, 2, 3], 0.)
+    # compute persistence only on the original
+    st.compute_persistence()
+
+    st_copy = st.copy()
+    assert st_copy == st
+    st_filt_list = list(st.get_filtration())
+
+    # check persistence is not copied
+    assert st.__is_persistence_defined() == True
+    assert st_copy.__is_persistence_defined() == False
+
+    # remove something in the copy and check the copy is included in the original
+    st_copy.remove_maximal_simplex([1, 2, 3])
+    a_filt_list = list(st_copy.get_filtration())
+    assert len(a_filt_list) < len(st_filt_list)
+
+    for a_splx in a_filt_list:
+        assert a_splx in st_filt_list
+    
+    # test double free
+    del st
+    del st_copy
+
+def test_simplex_tree_deep_copy_constructor():
+    st = SimplexTree()
+    st.insert([1, 2, 3], 0.)
+    # compute persistence only on the original
+    st.compute_persistence()
+
+    st_copy = SimplexTree(st)
+    assert st_copy == st
+    st_filt_list = list(st.get_filtration())
+
+    # check persistence is not copied
+    assert st.__is_persistence_defined() == True
+    assert st_copy.__is_persistence_defined() == False
+
+    # remove something in the copy and check the copy is included in the original
+    st_copy.remove_maximal_simplex([1, 2, 3])
+    a_filt_list = list(st_copy.get_filtration())
+    assert len(a_filt_list) < len(st_filt_list)
+
+    for a_splx in a_filt_list:
+        assert a_splx in st_filt_list
+    
+    # test double free
+    del st
+    del st_copy
+
+def test_simplex_tree_constructor_exception():
+    with pytest.raises(TypeError):
+        st = SimplexTree(other = "Construction from a string shall raise an exception")
 
 def test_expansion_with_blocker():
     st=SimplexTree()
diff --git a/src/python/test/test_subsampling.py b/src/python/test/test_subsampling.py
index 31f64e32..4019852e 100755
--- a/src/python/test/test_subsampling.py
+++ b/src/python/test/test_subsampling.py
@@ -141,12 +141,16 @@ def test_simple_sparsify_points():
     # assert gudhi.sparsify_point_set(points = [], min_squared_dist = 0.0) == []
     # assert gudhi.sparsify_point_set(points = [], min_squared_dist = 10.0) == []
     assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=0.0) == point_set
-    assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=1.0) == point_set
-    assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=2.0) == [
+    assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=0.999) == point_set
+    assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=1.001) == [
         [0, 1],
         [1, 0],
     ]
-    assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=2.01) == [[0, 1]]
+    assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=1.999) == [
+        [0, 1],
+        [1, 0],
+    ]
+    assert gudhi.sparsify_point_set(points=point_set, min_squared_dist=2.001) == [[0, 1]]
 
     assert (
         len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=0.0))
@@ -157,11 +161,11 @@ def test_simple_sparsify_points():
         == 5
     )
     assert (
-        len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=40.0))
+        len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=40.1))
         == 4
     )
     assert (
-        len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=90.0))
+        len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=89.9))
         == 3
     )
     assert (
@@ -169,7 +173,7 @@ def test_simple_sparsify_points():
         == 2
     )
     assert (
-        len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=325.0))
+        len(gudhi.sparsify_point_set(off_file="subsample.off", min_squared_dist=324.9))
         == 2
     )
     assert (
diff --git a/src/python/test/test_tomato.py b/src/python/test/test_tomato.py
index ecab03c4..c571f799 100755
--- a/src/python/test/test_tomato.py
+++ b/src/python/test/test_tomato.py
@@ -37,7 +37,7 @@ def test_tomato_1():
     t = Tomato(metric="euclidean", graph_type="radius", r=4.7, k=4)
     t.fit(a)
     assert t.max_weight_per_cc_.size == 2
-    assert np.array_equal(t.neighbors_, [[0, 1, 2], [0, 1, 2], [0, 1, 2], [3, 4, 5, 6], [3, 4, 5], [3, 4, 5], [3, 6]])
+    assert t.neighbors_ == [[0, 1, 2], [0, 1, 2], [0, 1, 2], [3, 4, 5, 6], [3, 4, 5], [3, 4, 5], [3, 6]]
     t.plot_diagram()
 
     t = Tomato(graph_type="radius", r=4.7, k=4, symmetrize_graph=True)
diff --git a/src/python/test/test_wasserstein_distance.py b/src/python/test/test_wasserstein_distance.py
index 90d26809..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,39 +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)
 
-def test_wasserstein_distance_grad():
-    import torch
-
-    diag1 = torch.tensor([[2.7, 3.7], [9.6, 14.0], [34.2, 34.974]], requires_grad=True)
-    diag2 = torch.tensor([[2.8, 4.45], [9.5, 14.1]], requires_grad=True)
-    diag3 = torch.tensor([[2.8, 4.45], [9.5, 14.1]], requires_grad=True)
-    assert diag1.grad is None and diag2.grad is None and diag3.grad is None
-    dist12 = pot(diag1, diag2, internal_p=2, order=2, enable_autodiff=True)
-    dist30 = pot(diag3, torch.tensor([]), internal_p=2, order=2, enable_autodiff=True)
-    dist12.backward()
-    dist30.backward()
-    assert not torch.isnan(diag1.grad).any() and not torch.isnan(diag2.grad).any() and not torch.isnan(diag3.grad).any()
-    diag4 = torch.tensor([[0., 10.]], requires_grad=True)
-    diag5 = torch.tensor([[1., 11.], [3., 4.]], requires_grad=True)
-    dist45 = pot(diag4, diag5, internal_p=1, order=1, enable_autodiff=True)
-    assert dist45 == 3.
-    dist45.backward()
-    assert np.array_equal(diag4.grad, [[-1., -1.]])
-    assert np.array_equal(diag5.grad, [[1., 1.], [-1., 1.]])
-    diag6 = torch.tensor([[5., 10.]], requires_grad=True)
-    pot(diag6, diag6, internal_p=2, order=2, enable_autodiff=True).backward()
-    # https://github.com/jonasrauber/eagerpy/issues/6
-    # assert np.array_equal(diag6.grad, [[0., 0.]])
diff --git a/src/python/test/test_wasserstein_with_tensors.py b/src/python/test/test_wasserstein_with_tensors.py
new file mode 100755
index 00000000..e3f1411a
--- /dev/null
+++ b/src/python/test/test_wasserstein_with_tensors.py
@@ -0,0 +1,47 @@
+""" This file is part of the Gudhi Library - https://gudhi.inria.fr/ - which is released under MIT.
+    See file LICENSE or go to https://gudhi.inria.fr/licensing/ for full license details.
+    Author(s):       Mathieu Carriere
+
+    Copyright (C) 2020 Inria
+
+    Modification(s):
+      - YYYY/MM Author: Description of the modification
+"""
+
+from gudhi.wasserstein import wasserstein_distance as pot
+import numpy as np
+import torch
+import tensorflow as tf
+
+def test_wasserstein_distance_grad():
+    diag1 = torch.tensor([[2.7, 3.7], [9.6, 14.0], [34.2, 34.974]], requires_grad=True)
+    diag2 = torch.tensor([[2.8, 4.45], [9.5, 14.1]], requires_grad=True)
+    diag3 = torch.tensor([[2.8, 4.45], [9.5, 14.1]], requires_grad=True)
+    assert diag1.grad is None and diag2.grad is None and diag3.grad is None
+    dist12 = pot(diag1, diag2, internal_p=2, order=2, enable_autodiff=True)
+    dist30 = pot(diag3, torch.tensor([]), internal_p=2, order=2, enable_autodiff=True)
+    dist12.backward()
+    dist30.backward()
+    assert not torch.isnan(diag1.grad).any() and not torch.isnan(diag2.grad).any() and not torch.isnan(diag3.grad).any()
+    diag4 = torch.tensor([[0., 10.]], requires_grad=True)
+    diag5 = torch.tensor([[1., 11.], [3., 4.]], requires_grad=True)
+    dist45 = pot(diag4, diag5, internal_p=1, order=1, enable_autodiff=True)
+    assert dist45 == 3.
+    dist45.backward()
+    assert np.array_equal(diag4.grad, [[-1., -1.]])
+    assert np.array_equal(diag5.grad, [[1., 1.], [-1., 1.]])
+    diag6 = torch.tensor([[5., 10.]], requires_grad=True)
+    pot(diag6, diag6, internal_p=2, order=2, enable_autodiff=True).backward()
+    # https://github.com/jonasrauber/eagerpy/issues/6
+    # assert np.array_equal(diag6.grad, [[0., 0.]])
+
+def test_wasserstein_distance_grad_tensorflow():
+    with tf.GradientTape() as tape:
+        diag4 = tf.convert_to_tensor(tf.Variable(initial_value=np.array([[0., 10.]]),           trainable=True))
+        diag5 = tf.convert_to_tensor(tf.Variable(initial_value=np.array([[1., 11.], [3., 4.]]), trainable=True))
+        dist45 = pot(diag4, diag5, internal_p=1, order=1, enable_autodiff=True)
+        assert dist45 == 3.
+
+    grads = tape.gradient(dist45, [diag4, diag5])
+    assert np.array_equal(grads[0].values, [[-1., -1.]])
+    assert np.array_equal(grads[1].values, [[1., 1.], [-1., 1.]])
\ No newline at end of file