From 896ba1042068951b489ca0bd81bae91f8862e3b3 Mon Sep 17 00:00:00 2001 From: Marc Glisse Date: Fri, 13 Sep 2024 23:32:47 +0200 Subject: [PATCH] More tweaks to ripser test --- .../test/test_sklearn_rips_persistence.py | 46 ++++++++++++------- 1 file changed, 29 insertions(+), 17 deletions(-) diff --git a/src/python/test/test_sklearn_rips_persistence.py b/src/python/test/test_sklearn_rips_persistence.py index 8865f0ff0..928c82c31 100644 --- a/src/python/test/test_sklearn_rips_persistence.py +++ b/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)))