don't include self for knn

Author: scottgigante

graphtools/graphs.py

@@ -85,11 +85,11 @@ def __init__(self, data, knn=5, decay=None,
         if decay is None and bandwidth is not None:
             warnings.warn("`bandwidth` is not used when `decay=None`.",
                           UserWarning)
-        if knn > data.shape[0]:
+        if knn > data.shape[0] - 2:
             warnings.warn("Cannot set knn ({k}) to be greater than "
                           "n_samples ({n}). Setting knn={n}".format(
-                              k=knn, n=data.shape[0]))
-            knn = data.shape[0]
+                              k=knn, n=data.shape[0] - 2))
+            knn = data.shape[0] - 2
         if n_pca is None and data.shape[1] > 500:
             warnings.warn("Building a kNNGraph on data of shape {} is "
                           "expensive. Consider setting n_pca.".format(
@@ -189,7 +189,7 @@ def knn_tree(self):
         except AttributeError:
             try:
                 self._knn_tree = NearestNeighbors(
-                    n_neighbors=self.knn,
+                    n_neighbors=self.knn + 1,
                     algorithm='ball_tree',
                     metric=self.distance,
                     n_jobs=self.n_jobs).fit(self.data_nu)
@@ -201,7 +201,7 @@ def knn_tree(self):
                         self.distance),
                     UserWarning)
                 self._knn_tree = NearestNeighbors(
-                    n_neighbors=self.knn,
+                    n_neighbors=self.knn + 1,
                     algorithm='auto',
                     metric=self.distance,
                     n_jobs=self.n_jobs).fit(self.data_nu)
@@ -219,9 +219,35 @@ def build_kernel(self):
             symmetric matrix with ones down the diagonal
             with no non-negative entries.
         """
-        K = self.build_kernel_to_data(self.data_nu)
+        K = self.build_kernel_to_data(self.data_nu, knn=self.knn + 1)
         return K
 
+    def _check_duplicates(self, distances, indices):
+        if np.any(distances[:, 1] == 0):
+            has_duplicates = distances[:, 1] == 0
+            if np.sum(distances[:, 1:] == 0) < 20:
+                idx = np.argwhere((distances == 0) &
+                                  has_duplicates[:, None])
+                duplicate_ids = np.array(
+                    [[indices[i[0], i[1]], i[0]]
+                     for i in idx if indices[i[0], i[1]] < i[0]])
+                duplicate_ids = duplicate_ids[
+                    np.argsort(duplicate_ids[:, 0])]
+                duplicate_names = ", ".join(["{} and {}".format(i[0], i[1])
+                                             for i in duplicate_ids])
+                warnings.warn(
+                    "Detected zero distance between samples {}. "
+                    "Consider removing duplicates to avoid errors in "
+                    "downstream processing.".format(duplicate_names),
+                    RuntimeWarning)
+            else:
+                warnings.warn(
+                    "Detected zero distance between {} pairs of samples. "
+                    "Consider removing duplicates to avoid errors in "
+                    "downstream processing.".format(
+                        np.sum(np.sum(distances[:, 1:]))),
+                    RuntimeWarning)
+
     def build_kernel_to_data(self, Y, knn=None, bandwidth=None,
                              bandwidth_scale=None):
         """Build a kernel from new input data `Y` to the `self.data`
@@ -281,30 +307,7 @@ def build_kernel_to_data(self, Y, knn=None, bandwidth=None,
             search_knn = min(knn * 20, self.data_nu.shape[0])
             distances, indices = knn_tree.kneighbors(
                 Y, n_neighbors=search_knn)
-            if np.any(distances[:, 1] == 0):
-                has_duplicates = distances[:, 1] == 0
-                if np.sum(distances[:, 1:] == 0) < 20:
-                    idx = np.argwhere((distances == 0) &
-                                      has_duplicates[:, None])
-                    duplicate_ids = np.array(
-                        [[indices[i[0], i[1]], i[0]]
-                         for i in idx if indices[i[0], i[1]] < i[0]])
-                    duplicate_ids = duplicate_ids[
-                        np.argsort(duplicate_ids[:, 0])]
-                    duplicate_names = ", ".join(["{} and {}".format(i[0], i[1])
-                                                 for i in duplicate_ids])
-                    warnings.warn(
-                        "Detected zero distance between samples {}. "
-                        "Consider removing duplicates to avoid errors in "
-                        "downstream processing.".format(duplicate_names),
-                        RuntimeWarning)
-                else:
-                    warnings.warn(
-                        "Detected zero distance between {} pairs of samples. "
-                        "Consider removing duplicates to avoid errors in "
-                        "downstream processing.".format(
-                            np.sum(np.sum(distances[:, 1:]))),
-                        RuntimeWarning)
+            self._check_duplicates(distances, indices)
             tasklogger.log_complete("KNN search")
             tasklogger.log_start("affinities")
             if bandwidth is None:
@@ -338,7 +341,7 @@ def build_kernel_to_data(self, Y, knn=None, bandwidth=None,
                     len(update_idx)))
             if search_knn > self.data_nu.shape[0] / 2:
                 knn_tree = NearestNeighbors(
-                    knn, algorithm='brute',
+                    search_knn, algorithm='brute',
                     n_jobs=self.n_jobs).fit(self.data_nu)
             if len(update_idx) > 0:
                 tasklogger.log_debug(
@@ -771,11 +774,11 @@ def __init__(self, data,
         if knn is None and bandwidth is None:
             raise ValueError(
                 "Either `knn` or `bandwidth` must be provided.")
-        if knn is not None and knn > data.shape[0]:
-            warnings.warn("Cannot set knn ({k}) to be greater than or equal to"
-                          " n_samples ({n}). Setting knn={n}".format(
-                              k=knn, n=data.shape[0] - 1))
-            knn = data.shape[0] - 1
+        if knn is not None and knn > data.shape[0] - 2:
+            warnings.warn("Cannot set knn ({k}) to be greater than "
+                          " n_samples - 2 ({n}). Setting knn={n}".format(
+                              k=knn, n=data.shape[0] - 2))
+            knn = data.shape[0] - 2
         if precomputed is not None:
             if precomputed not in ["distance", "affinity", "adjacency"]:
                 raise ValueError("Precomputed value {} not recognized. "
@@ -918,7 +921,8 @@ def build_kernel(self):
                     "Choose from ['affinity', 'adjacency', 'distance', "
                     "None]".format(self.precomputed))
             if self.bandwidth is None:
-                knn_dist = np.partition(pdx, self.knn, axis=1)[:, :self.knn]
+                knn_dist = np.partition(
+                    pdx, self.knn + 1, axis=1)[:, :self.knn + 1]
                 bandwidth = np.max(knn_dist, axis=1)
             elif callable(self.bandwidth):
                 bandwidth = self.bandwidth(pdx)
@@ -1300,8 +1304,6 @@ def build_kernel_to_data(self, Y, theta=None):
         transformation of the landmarks can be trivially applied to `Y` by
         performing
 
-        TODO: test this.
-
         `transform_Y = transitions.dot(transform)`
 
         Parameters
@@ -1323,52 +1325,6 @@ def build_kernel_to_data(self, Y, theta=None):
             Transition matrix from `Y` to `self.data`
         """
         raise NotImplementedError
-        tasklogger.log_warning("building MNN kernel to theta is experimental")
-        if not isinstance(self.theta, str) and \
-                not isinstance(self.theta, numbers.Number):
-            if theta is None:
-                raise ValueError(
-                    "self.theta is a matrix but theta is not provided.")
-            elif len(theta) != len(self.samples):
-                raise ValueError(
-                    "theta should have one value for every sample")
-
-        Y = self._check_extension_shape(Y)
-        kernel_xy = []
-        kernel_yx = []
-        # don't really need within Y kernel
-        Y_graph = kNNGraph(Y, n_pca=None, knn=0, **(self.knn_args))
-        y_knn = self._weight_knn(sample_size=Y.shape[0])
-        for i, X in enumerate(self.subgraphs):
-            kernel_xy.append(X.build_kernel_to_data(
-                Y, knn=self.weighted_knn[i]))  # kernel X -> Y
-            kernel_yx.append(Y_graph.build_kernel_to_data(
-                X.data_nu, knn=y_knn))  # kernel Y -> X
-        kernel_xy = sparse.hstack(kernel_xy)  # n_cells_y x n_cells_x
-        kernel_yx = sparse.vstack(kernel_yx)  # n_cells_x x n_cells_y
-
-        # symmetrize
-        if theta is not None:
-            # Gamma can be a vector with specific values transitions for
-            # each batch. This allows for technical replicates and
-            # experimental samples to be corrected simultaneously
-            K = np.empty_like(kernel_xy)
-            for i, sample in enumerate(self.samples):
-                sample_idx = self.sample_idx == sample
-                K[:, sample_idx] = theta[i] * \
-                    kernel_xy[:, sample_idx].minimum(
-                        kernel_yx[sample_idx, :].T) + \
-                    (1 - theta[i]) * \
-                    kernel_xy[:, sample_idx].maximum(
-                        kernel_yx[sample_idx, :].T)
-        if self.theta == "+":
-            K = (kernel_xy + kernel_yx.T) / 2
-        elif self.theta == "*":
-            K = kernel_xy.multiply(kernel_yx.T)
-        else:
-            K = self.theta * kernel_xy.minimum(kernel_yx.T) + \
-                (1 - self.theta) * kernel_xy.maximum(kernel_yx.T)
-        return K
 
 
 class kNNLandmarkGraph(kNNGraph, LandmarkGraph):

test/test_exact.py

@@ -103,7 +103,7 @@ def test_k_too_large():
     build_graph(data,
                 n_pca=20,
                 decay=10,
-                knn=len(data) + 1,
+                knn=len(data) - 1,
                 thresh=0)
 
 
@@ -131,7 +131,7 @@ def test_exact_graph():
     np.fill_diagonal(W, 0)
     G = pygsp.graphs.Graph(W)
     G2 = build_graph(data_small, thresh=0, n_pca=n_pca,
-                     decay=a, knn=k, random_state=42,
+                     decay=a, knn=k - 1, random_state=42,
                      bandwidth_scale=bandwidth_scale,
                      use_pygsp=True)
     assert(G.N == G2.N)
@@ -141,7 +141,7 @@ def test_exact_graph():
     assert(isinstance(G2, graphtools.graphs.TraditionalGraph))
     G2 = build_graph(pdx, n_pca=None, precomputed='distance',
                      bandwidth_scale=bandwidth_scale,
-                     decay=a, knn=k, random_state=42, use_pygsp=True)
+                     decay=a, knn=k - 1, random_state=42, use_pygsp=True)
     assert(G.N == G2.N)
     np.testing.assert_equal(G.dw, G2.dw)
     assert((G.W != G2.W).nnz == 0)
@@ -195,7 +195,7 @@ def test_truncated_exact_graph():
     G2 = build_graph(data_small, thresh=thresh,
                      graphtype='exact',
                      n_pca=n_pca,
-                     decay=a, knn=k, random_state=42,
+                     decay=a, knn=k - 1, random_state=42,
                      use_pygsp=True)
     assert(G.N == G2.N)
     np.testing.assert_equal(G.dw, G2.dw)
@@ -204,7 +204,7 @@ def test_truncated_exact_graph():
     assert(isinstance(G2, graphtools.graphs.TraditionalGraph))
     G2 = build_graph(pdx, n_pca=None, precomputed='distance',
                      thresh=thresh,
-                     decay=a, knn=k, random_state=42, use_pygsp=True)
+                     decay=a, knn=k - 1, random_state=42, use_pygsp=True)
     assert(G.N == G2.N)
     np.testing.assert_equal(G.dw, G2.dw)
     assert((G.W != G2.W).nnz == 0)
@@ -252,14 +252,14 @@ def test_truncated_exact_graph_sparse():
     G2 = build_graph(sp.coo_matrix(data_small), thresh=thresh,
                      graphtype='exact',
                      n_pca=n_pca,
-                     decay=a, knn=k, random_state=42,
+                     decay=a, knn=k - 1, random_state=42,
                      use_pygsp=True)
     assert(G.N == G2.N)
     np.testing.assert_allclose(G2.W.toarray(), G.W.toarray())
     assert(isinstance(G2, graphtools.graphs.TraditionalGraph))
     G2 = build_graph(sp.bsr_matrix(pdx), n_pca=None, precomputed='distance',
                      thresh=thresh,
-                     decay=a, knn=k, random_state=42, use_pygsp=True)
+                     decay=a, knn=k - 1, random_state=42, use_pygsp=True)
     assert(G.N == G2.N)
     np.testing.assert_equal(G.dw, G2.dw)
     assert((G.W != G2.W).nnz == 0)
@@ -304,7 +304,7 @@ def test_truncated_exact_graph_no_pca():
     G2 = build_graph(data_small, thresh=thresh,
                      graphtype='exact',
                      n_pca=n_pca,
-                     decay=a, knn=k, random_state=42,
+                     decay=a, knn=k - 1, random_state=42,
                      use_pygsp=True)
     assert(G.N == G2.N)
     np.testing.assert_equal(G.dw, G2.dw)
@@ -314,7 +314,7 @@ def test_truncated_exact_graph_no_pca():
     G2 = build_graph(sp.csr_matrix(data_small), thresh=thresh,
                      graphtype='exact',
                      n_pca=n_pca,
-                     decay=a, knn=k, random_state=42,
+                     decay=a, knn=k - 1, random_state=42,
                      use_pygsp=True)
     assert(G.N == G2.N)
     np.testing.assert_equal(G.dw, G2.dw)
@@ -379,7 +379,7 @@ def test_exact_graph_callable_bandwidth():
     W = np.divide(K, 2)
     np.fill_diagonal(W, 0)
     G = pygsp.graphs.Graph(W)
-    G2 = build_graph(data, n_pca=n_pca, knn=knn,
+    G2 = build_graph(data, n_pca=n_pca, knn=knn - 1,
                      decay=decay, bandwidth=bandwidth,
                      random_state=42,
                      thresh=thresh,
@@ -396,7 +396,7 @@ def test_exact_graph_callable_bandwidth():
     W = np.divide(K, 2)
     np.fill_diagonal(W, 0)
     G = pygsp.graphs.Graph(W)
-    G2 = build_graph(data, n_pca=n_pca, knn=knn,
+    G2 = build_graph(data, n_pca=n_pca, knn=knn - 1,
                      decay=decay, bandwidth=bandwidth,
                      random_state=42,
                      thresh=thresh,
@@ -432,7 +432,7 @@ def test_exact_graph_anisotropy():
     np.fill_diagonal(W, 0)
     G = pygsp.graphs.Graph(W)
     G2 = build_graph(data_small, thresh=0, n_pca=n_pca,
-                     decay=a, knn=k, random_state=42,
+                     decay=a, knn=k - 1, random_state=42,
                      use_pygsp=True, anisotropy=anisotropy)
     assert(isinstance(G2, graphtools.graphs.TraditionalGraph))
     assert(G.N == G2.N)
@@ -441,15 +441,15 @@ def test_exact_graph_anisotropy():
     assert((G.W != G2.W).nnz == 0)
     assert_raises(ValueError, build_graph,
                   data_small, thresh=0, n_pca=n_pca,
-                  decay=a, knn=k, random_state=42,
+                  decay=a, knn=k - 1, random_state=42,
                   use_pygsp=True, anisotropy=-1)
     assert_raises(ValueError, build_graph,
                   data_small, thresh=0, n_pca=n_pca,
-                  decay=a, knn=k, random_state=42,
+                  decay=a, knn=k - 1, random_state=42,
                   use_pygsp=True, anisotropy=2)
     assert_raises(ValueError, build_graph,
                   data_small, thresh=0, n_pca=n_pca,
-                  decay=a, knn=k, random_state=42,
+                  decay=a, knn=k - 1, random_state=42,
                   use_pygsp=True, anisotropy='invalid')
 
 #####################################################
@@ -462,32 +462,32 @@ def test_build_dense_exact_kernel_to_data(**kwargs):
     n = G.data.shape[0]
     K = G.build_kernel_to_data(data[:n // 2, :])
     assert(K.shape == (n // 2, n))
-    K = G.build_kernel_to_data(G.data)
-    assert(np.sum(G.kernel != (K + K.T) / 2) == 0)
-    K = G.build_kernel_to_data(G.data_nu)
-    assert(np.sum(G.kernel != (K + K.T) / 2) == 0)
+    K = G.build_kernel_to_data(G.data, knn=G.knn + 1)
+    np.testing.assert_equal(G.kernel - (K + K.T) / 2, 0)
+    K = G.build_kernel_to_data(G.data_nu, knn=G.knn + 1)
+    np.testing.assert_equal(G.kernel - (K + K.T) / 2, 0)
 
 
 def test_build_dense_exact_callable_bw_kernel_to_data(**kwargs):
     G = build_graph(data, decay=10, thresh=0, bandwidth=lambda x: x.mean(1))
     n = G.data.shape[0]
     K = G.build_kernel_to_data(data[:n // 2, :])
     assert(K.shape == (n // 2, n))
-    K = G.build_kernel_to_data(G.data)
-    assert(np.sum(G.kernel != (K + K.T) / 2) == 0)
-    K = G.build_kernel_to_data(G.data_nu)
-    assert(np.sum(G.kernel != (K + K.T) / 2) == 0)
+    K = G.build_kernel_to_data(G.data, knn=G.knn + 1)
+    np.testing.assert_equal(G.kernel - (K + K.T) / 2, 0)
+    K = G.build_kernel_to_data(G.data_nu, knn=G.knn + 1)
+    np.testing.assert_equal(G.kernel - (K + K.T) / 2, 0)
 
 
 def test_build_sparse_exact_kernel_to_data(**kwargs):
     G = build_graph(data, decay=10, thresh=0, sparse=True)
     n = G.data.shape[0]
     K = G.build_kernel_to_data(data[:n // 2, :])
     assert(K.shape == (n // 2, n))
-    K = G.build_kernel_to_data(G.data)
-    assert(np.sum(G.kernel != (K + K.T) / 2) == 0)
-    K = G.build_kernel_to_data(G.data_nu)
-    assert(np.sum(G.kernel != (K + K.T) / 2) == 0)
+    K = G.build_kernel_to_data(G.data, knn=G.knn + 1)
+    np.testing.assert_equal(G.kernel - (K + K.T) / 2, 0)
+    K = G.build_kernel_to_data(G.data_nu, knn=G.knn + 1)
+    np.testing.assert_equal(G.kernel - (K + K.T) / 2, 0)
 
 
 def test_exact_interpolate():