More tweaks to ripser test

src/python/test/test_sklearn_rips_persistence.py

@@ -12,9 +12,9 @@
 from gudhi.sklearn.rips_persistence import RipsPersistence
 from gudhi import RipsComplex, SimplexTree
 from gudhi._ripser import _lower, _full, _sparse, _lower_to_coo, _lower_cone_radius
+from gudhi import bottleneck_distance
 from scipy.sparse import coo_matrix
 from scipy.spatial.distance import cdist
-from scipy.spatial import cKDTree
 import numpy as np
 import random
 import pytest
@@ -94,13 +94,11 @@ def test_big():
     RipsPersistence(range(25), threshold=10).fit_transform([X])
 
 
-def test_ripser_interfaces():
-    primes = [2, 3, 5, 7, 11, 13, 17, 19, 23, 29]
-    field = primes[random.randint(0, 9)]
+def cmp_rips(point_cloud):
+    primes = [2, 3, 11, 17, 29]  # Small list so 2 is often selected
+    field = random.choice(primes)
     print(f"random prime = {field}")
 
-    point_cloud = points.sphere(n_samples=random.randint(100, 150), ambient_dim=2)
-    # point_cloud = np.random.rand(random.randint(100, 150), random.randint(2, 3))
     print(f"nb points = {len(point_cloud)}, dim = {point_cloud.shape[1]}")
     dists = cdist(point_cloud, point_cloud)
 
@@ -109,26 +107,34 @@ def test_ripser_interfaces():
     assert cr == _lower_cone_radius(dists)
     assert cr < 2.0
 
-    ## As there is no easy way to force the use of SimplexTree, let's build it
+    ## Compute with the SimplexTree
     stree = RipsComplex(distance_matrix=dists).create_simplex_tree(max_dimension=2)
     stree.compute_persistence(homology_coeff_field=field, persistence_dim_max=True)
     dgm0 = stree.persistence_intervals_in_dimension(0)
     dgm1 = stree.persistence_intervals_in_dimension(1)
 
+    # Compute with Ripser
     dgm = _full(dists, max_dimension=1, max_edge_length=float("inf"), homology_coeff_field=field)
-    np.testing.assert_almost_equal(dgm0, dgm[0])
-    np.testing.assert_almost_equal(dgm1, dgm[1])
+    # The order of the intervals may differ, so we cannot compare the arrays with np.testing.assert_almost_equal
+    assert bottleneck_distance(dgm0, dgm[0]) < 1e-8
+    assert bottleneck_distance(dgm1, dgm[1]) < 1e-8
 
     dgm = _lower(dists, max_dimension=1, max_edge_length=float("inf"), homology_coeff_field=field)
-    np.testing.assert_almost_equal(dgm0, dgm[0])
-    np.testing.assert_almost_equal(dgm1, dgm[1])
+    assert bottleneck_distance(dgm0, dgm[0]) < 1e-8
+    assert bottleneck_distance(dgm1, dgm[1]) < 1e-8
 
-    # From a coo matrix
+    # Convert to coo matrix
     n = len(point_cloud)
     dists_copy = np.array(dists, copy=True)
     dists_copy[np.triu_indices_from(dists_copy)] = 0  # Keep only the lower entries
     dists_sparse = coo_matrix(dists_copy)
-    ## As there is no easy way to force the use of SimplexTree, let's build it
+    s1 = _lower_to_coo(dists, float("inf"))
+    s2 = _lower_to_coo(dists_copy, float("inf"))
+    s1 = coo_matrix((s1[2], (s1[0], s1[1])), shape=(n, n))
+    s2 = coo_matrix((s2[2], (s2[0], s2[1])), shape=(n, n))
+    assert np.array_equal(s1.toarray(), s2.toarray())
+    assert np.array_equal(dists_sparse.toarray(), s1.toarray())
+    ## Compute with the SimplexTree
     stree = SimplexTree()
     stree.insert_batch(np.arange(n).reshape(1, -1), np.zeros(n))
     stree.insert_edges_from_coo_matrix(dists_sparse)
@@ -137,6 +143,7 @@ def test_ripser_interfaces():
     sp_dgm0 = stree.persistence_intervals_in_dimension(0)
     sp_dgm1 = stree.persistence_intervals_in_dimension(1)
 
+    # Compute with Ripser
     sp_dgm = _sparse(
         dists_sparse.row,
         dists_sparse.col,
@@ -146,7 +153,12 @@ def test_ripser_interfaces():
         max_edge_length=float("inf"),
         homology_coeff_field=field,
     )
-    np.testing.assert_almost_equal(sp_dgm0, sp_dgm[0])
-    np.testing.assert_almost_equal(sp_dgm1, sp_dgm[1])
-    np.testing.assert_almost_equal(sp_dgm0, dgm0)
-    np.testing.assert_almost_equal(sp_dgm1, dgm1)
+    assert bottleneck_distance(sp_dgm0, sp_dgm[0]) < 1e-8
+    assert bottleneck_distance(sp_dgm1, sp_dgm[1]) < 1e-8
+    assert bottleneck_distance(sp_dgm0, dgm0) < 1e-8
+    assert bottleneck_distance(sp_dgm1, dgm1) < 1e-8
+
+
+def test_ripser_interfaces():
+    cmp_rips(points.sphere(n_samples=random.randint(100, 150), ambient_dim=2))
+    cmp_rips(np.random.rand(random.randint(100, 150), random.randint(2, 3)))