From 3160df09d8bab660bb2c443aefad90df0af1dabc Mon Sep 17 00:00:00 2001 From: wreise Date: Tue, 17 Oct 2023 18:41:49 +0200 Subject: [PATCH 1/4] Toroidal and circular coordinates on new points --- dreimac/circularcoords.py | 6 +++++- dreimac/emcoords.py | 19 +++++++++++++------ dreimac/toroidalcoords.py | 11 ++++++++--- 3 files changed, 26 insertions(+), 10 deletions(-) diff --git a/dreimac/circularcoords.py b/dreimac/circularcoords.py index 2992663..02945bb 100644 --- a/dreimac/circularcoords.py +++ b/dreimac/circularcoords.py @@ -40,6 +40,7 @@ def __init__( def get_coordinates( self, + X_query=None, perc=0.5, cocycle_idx=0, partunity_fn=PartUnity.linear, @@ -52,6 +53,8 @@ def get_coordinates( Parameters ---------- + X_query: ndarray(M, d) or None + A point cloud to compute the circular coordinates on. If None, uses self.X. perc : float Percent coverage. Must be between 0 and 1. cocycle_idx : integer @@ -69,12 +72,13 @@ def get_coordinates( Returns ------- - thetas : ndarray(N) + thetas : ndarray(M) Circular coordinates """ return ToroidalCoords.get_coordinates( self, + X_query, perc, [cocycle_idx], partunity_fn, diff --git a/dreimac/emcoords.py b/dreimac/emcoords.py index f71d985..3ed9c79 100644 --- a/dreimac/emcoords.py +++ b/dreimac/emcoords.py @@ -3,6 +3,7 @@ """ import numpy as np from scipy.sparse.linalg import lsqr +from scipy.spatial.distance import cdist import time from .utils import CohomologyUtils from ripser import ripser @@ -43,10 +44,11 @@ def __init__(self, X, n_landmarks, distance_matrix, prime, maxdim, verbose): if verbose: tic = time.time() print("Doing TDA...") + self.distance_matrix = distance_matrix if distance_matrix is False: - ripser_metric_input = X + ripser_metric_input = X elif X.shape[0] == X.shape[1]: - ripser_metric_input = X + ripser_metric_input = X else: ripser_metric_input = X[:,:X.shape[0]] res = ripser( @@ -149,7 +151,7 @@ def get_cover_radius(self, perc, cohomdeath_rips, cohombirth_rips, standard_rang return self.r_cover_, self.rips_threshold_ - def get_covering_partition(self, r_cover, partunity_fn): + def get_covering_partition(self, r_cover, partunity_fn, X_query=None): """ Create the open covering U = {U_1,..., U_{s+1}} and partition of unity @@ -159,15 +161,20 @@ def get_covering_partition(self, r_cover, partunity_fn): Covering radius partunity_fn: (dist_land_data, r_cover) -> phi A function from the distances of each landmark to a bump function + X_query: ndarray(M, d) + A point cloud to compute the circular coordinates on. If None, uses self.X. Returns ------- - varphi: ndarray(n_data, dtype=float) + varphi: ndarray(M, dtype=float) varphi_j(b) = phi_j(b)/(phi_1(b) + ... + phi_{n_landmarks}(b)), - ball_indx: ndarray(n_data, dtype=int) + ball_indx: ndarray(M, dtype=int) The index of the first open set each data point belongs to """ - dist_land_data = self.dist_land_data_ + if X_query is None: + dist_land_data = self.dist_land_data_ + else: # calculate the distance between the landmarks and the query + dist_land_data = cdist(self.X_[self.idx_land_], X_query) U = dist_land_data < r_cover phi = np.zeros_like(dist_land_data) phi[U] = partunity_fn(dist_land_data[U], r_cover) diff --git a/dreimac/toroidalcoords.py b/dreimac/toroidalcoords.py index 9dc4196..1f4fb75 100644 --- a/dreimac/toroidalcoords.py +++ b/dreimac/toroidalcoords.py @@ -44,6 +44,7 @@ def __init__( def get_coordinates( self, + X_query=None, perc=0.5, cocycle_idxs=[0], partunity_fn=PartUnity.linear, @@ -55,6 +56,8 @@ def get_coordinates( Parameters ---------- + X_query: ndarray(M, d) or None + A point cloud to compute the toroidal coordinates on. If None, uses self.X. perc : float Percent coverage. Must be between 0 and 1. cocycle_idx : integer @@ -72,11 +75,13 @@ def get_coordinates( Returns ------- - thetas : ndarray(n, N) + thetas : ndarray(n, M) List of circular coordinates, with n the length of cocycle_idxs """ - + if (X_query is not None) and self.distance_matrix: + raise NotImplementedError("To compute coordinates on points other than the input, " + "initialize ToroidalCoords with points instead of a distance matrix.") # get representative cocycles and the intersection of their supports homological_dimension = 1 cohomdeaths, cohombirths, cocycles = zip( @@ -107,7 +112,7 @@ def get_coordinates( ) # compute partition of unity and choose a cover element for each data point - varphi, ball_indx = EMCoords.get_covering_partition(self, r_cover, partunity_fn) + varphi, ball_indx = EMCoords.get_covering_partition(self, r_cover, partunity_fn, X_query) # compute boundary matrix dist_land_land = self.dist_land_land_ From 2c32134027774f108ee22c308fc80b02bfcb7682 Mon Sep 17 00:00:00 2001 From: wreise Date: Wed, 25 Oct 2023 14:23:46 +0200 Subject: [PATCH 2/4] Tests for circular --- test/test_circular.py | 32 +++++++++++++++++++++++++++++++- 1 file changed, 31 insertions(+), 1 deletion(-) diff --git a/test/test_circular.py b/test/test_circular.py index 4d205b7..f934bea 100644 --- a/test/test_circular.py +++ b/test/test_circular.py @@ -1,4 +1,5 @@ import numpy as np +from pytest import approx from scipy.spatial import KDTree from dreimac import CircularCoords, ToroidalCoords, GeometryExamples @@ -97,6 +98,35 @@ def test_toroidal_coordinates_less_energy(self): np.linalg.norm(tc.gram_matrix_), np.linalg.norm(tc.original_gram_matrix_) ) + def test_toroidal_consistent_on_query(self): + """ + Check that the coordinates computed on the query point-cloud are the same as + those computed on the initial points, when the two are equal up to permutation. + """ + X = GeometryExamples.bullseye() + n_landmarks = 300 + tc = ToroidalCoords(X, n_landmarks=n_landmarks) + perc = 0.1; cohomology_classes = [1] + t_coords_ref = tc.get_coordinates(perc=perc, cocycle_idxs=cohomology_classes)[0] + + indices = np.random.randint(low=0, high=X.shape[0], size=(20,)).astype(int) + X_query = X[indices] + t_coords_new = tc.get_coordinates(X_query=X_query, perc=perc, cocycle_idxs=cohomology_classes)[0] + assert _less_than_or_equal_with_tolerance(0., np.array([_circle_distance(x,y) for x,y in zip(t_coords_ref[indices], t_coords_new)])) + + def test_toroidal_coordinates_continuity(self): + X = GeometryExamples.bullseye() + n_landmarks = 300 + tc = ToroidalCoords(X, n_landmarks=n_landmarks) + perc = 0.1; cohomology_classes = [1] + indices = np.random.randint(low=0, high=X.shape[0], size=(20,)).astype(int) + + X_subsample = X[indices] + t_coords_ref = tc.get_coordinates(X_query=X_subsample, perc=perc, cocycle_idxs=cohomology_classes)[0] + + X_noisy = X_subsample + 0.01*(np.random.rand(*X_subsample.shape)-0.5) + t_coords_new = tc.get_coordinates(X_query=X_noisy, perc=perc, cocycle_idxs=cohomology_classes)[0] + assert all(_less_than_or_equal_with_tolerance(0., np.array([_circle_distance(x,y) for x,y in zip(t_coords_ref, t_coords_new)]))) def test_trefoil(self): """Check that circular coordinates returns a continuous map, even when the lifted @@ -110,7 +140,7 @@ def test_trefoil(self): cc = CircularCoords(X, 300, prime=prime) coords = cc.get_coordinates( perc=large_perc, - cocycle_idx=0, + cocycle_idx=[0], check_cocycle_condition=True, ) assert len(coords) == len(X) From a08695bf44f39991e9dcff4f376e7ea8be6b3ca1 Mon Sep 17 00:00:00 2001 From: wreise Date: Wed, 25 Oct 2023 14:31:40 +0200 Subject: [PATCH 3/4] Remove un-necessary import --- test/test_circular.py | 1 - 1 file changed, 1 deletion(-) diff --git a/test/test_circular.py b/test/test_circular.py index f934bea..04597f7 100644 --- a/test/test_circular.py +++ b/test/test_circular.py @@ -1,5 +1,4 @@ import numpy as np -from pytest import approx from scipy.spatial import KDTree from dreimac import CircularCoords, ToroidalCoords, GeometryExamples From 470e450a61dfc9f1fbad558dd5f7f453e6a27fbc Mon Sep 17 00:00:00 2001 From: wreise Date: Mon, 30 Oct 2023 10:34:29 +0100 Subject: [PATCH 4/4] Revert erroneous changes --- test/test_circular.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/test/test_circular.py b/test/test_circular.py index 04597f7..7e0141a 100644 --- a/test/test_circular.py +++ b/test/test_circular.py @@ -139,7 +139,7 @@ def test_trefoil(self): cc = CircularCoords(X, 300, prime=prime) coords = cc.get_coordinates( perc=large_perc, - cocycle_idx=[0], + cocycle_idx=0, check_cocycle_condition=True, ) assert len(coords) == len(X)