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_ diff --git a/test/test_circular.py b/test/test_circular.py index 4d205b7..7e0141a 100644 --- a/test/test_circular.py +++ b/test/test_circular.py @@ -97,6 +97,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