flann returns the squared distance when called with 'euclidean' distance -> fix

fix radius nn graph for spdist

Author: naspert

pygsp/graphs/nngraphs/nngraph.py

@@ -49,6 +49,7 @@ def _knn_flann(X, num_neighbors, dist_type, order):
     # do not allow it
     if dist_type == 'max_dist':
         raise ValueError('FLANN and max_dist is not supported')
+    
     pfl = _import_pfl()
     pfl.set_distance_type(dist_type, order=order)
     flann = pfl.FLANN()
@@ -58,6 +59,8 @@ def _knn_flann(X, num_neighbors, dist_type, order):
     # seems to work best).
     NN, D = flann.nn(X, X, num_neighbors=(num_neighbors + 1), 
                      algorithm='kdtree')
+    if dist_type == 'euclidean': # flann returns squared distances
+        return NN, np.sqrt(D)
     return NN, D
 
 def _radius_sp_kdtree(X, epsilon, dist_type, order=0):
@@ -86,8 +89,9 @@ def _radius_sp_pdist(X, epsilon, dist_type, order):
     NN = []
     for k in range(N):
         v = pd[k, pdf[k, :]]
+        d = pd[k, :].argsort()
         # use the same conventions as in scipy.distance.kdtree
-        NN.append(v.argsort())
+        NN.append(d[0:len(v)])
         D.append(np.sort(v))
 
     return NN, D
@@ -98,21 +102,32 @@ def _radius_flann(X, epsilon, dist_type, order=0):
     # do not allow it
     if dist_type == 'max_dist':
         raise ValueError('FLANN and max_dist is not supported')
-        
     pfl = _import_pfl()
+    
     pfl.set_distance_type(dist_type, order=order)
     flann = pfl.FLANN()
     flann.build_index(X)
 
     D = []
     NN = []
     for k in range(N):
-        nn, d = flann.nn_radius(X[k, :], epsilon)
+        nn, d = flann.nn_radius(X[k, :], epsilon*epsilon)
         D.append(d)
         NN.append(nn)
     flann.delete_index()
+    if dist_type == 'euclidean': # flann returns squared distances
+        return NN, np.sqrt(D)
     return NN, D
 
+def center_input(X, N):
+    return X - np.kron(np.ones((N, 1)), np.mean(X, axis=0))
+
+def rescale_input(X, N, d):
+     bounding_radius = 0.5 * np.linalg.norm(np.amax(X, axis=0) -
+                                                   np.amin(X, axis=0), 2)
+     scale = np.power(N, 1. / float(min(d, 3))) / 10.
+     return X * scale / bounding_radius
+
 class NNGraph(Graph):
     r"""Nearest-neighbor graph from given point cloud.
 
@@ -207,14 +222,10 @@ def __init__(self, Xin, NNtype='knn', backend='scipy-kdtree', center=True,
         Xout = self.Xin
 
         if self.center:
-            Xout = self.Xin - np.kron(np.ones((N, 1)),
-                                      np.mean(self.Xin, axis=0))
+            Xout = center_input(Xout, N)
 
         if self.rescale:
-            bounding_radius = 0.5 * np.linalg.norm(np.amax(Xout, axis=0) -
-                                                   np.amin(Xout, axis=0), 2)
-            scale = np.power(N, 1. / float(min(d, 3))) / 10.
-            Xout *= scale / bounding_radius
+            Xout = rescale_input(Xout, N, d)
 
 
 

pygsp/tests/test_graphs.py

@@ -187,6 +187,7 @@ def test_nngraph(self):
 
         for cur_backend in backends:
             for dist_type in dist_types:
+                #print("backend={} dist={}".format(cur_backend, dist_type))
                 if cur_backend == 'flann' and dist_type == 'max_dist':
                     self.assertRaises(ValueError, graphs.NNGraph, Xin, 
                                       NNtype='knn', backend=cur_backend, 
@@ -199,6 +200,20 @@ def test_nngraph(self):
                                    backend=cur_backend, 
                                    dist_type=dist_type, order=order)
 
+    def test_nngraph_consistency(self):
+        #Xin = np.arange(180).reshape(60, 3)
+        Xin = np.random.uniform(-5, 5, (60, 3))
+        dist_types = ['euclidean', 'manhattan', 'max_dist', 'minkowski']
+        backends = ['scipy-kdtree', 'flann']
+        num_neighbors=5
+        
+        G = graphs.NNGraph(Xin, NNtype='knn', 
+                           backend='scipy-pdist', k=num_neighbors)
+         for cur_backend in backends:
+            for dist_type in dist_types:
+                Gt = graphs.NNGraph(Xin, NNtype='knn', 
+                                    backend=cur_backend, k=num_neighbors)
+        
     def test_bunny(self):
         graphs.Bunny()