Coordinates on new points

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,

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)

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_

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