Repository: scipy/scipy
Issue category: Performance Issue

Issue:
Massive performance regression in cKDTree.query with L_inf distance in scipy 0.17
Hi,

While running some of my own code, I noticed a major slowdown in scipy 0.17 which I tracked to the cKDtree.query call with L_inf distance.
I did a git bisect and tracked the commit that caused the issue -- the commit is :
58a10c2688398ab440de8d6d7b5fba29671df182
Here is the example program that finishes in ~ 10 seconds before this commit and pretty much never finishes after it

```
import scipy.spatial.ckdtree
import numpy as np
arr1 = np.random.uniform(size=(6000000, 3))
tree2 = scipy.spatial.ckdtree.cKDTree(arr1)
distances, ind = tree2.query(arr1, 1, p=np.inf)
```

I also checked that if I change the distance from L_inf to L_2, then program finishes pretty quickly. I also tried to stop the program with gdb a few times to see where the code is mostly spending the time on and it seems to be here scipy/spatial/ckdtree/src/query.cxx:311. I also tried to dig into the problematic commit but I don't know enough about ckdtree internals to suggest a patch. 

Thanks, 
 Sergey

If relevant, I'm seeing this issue on a both systems that I tried  64bit RHEL6.7 with gcc 4.4 and Debian 7.9 with gcc-4.7. In all the cases it is python 2.7



Hints:

Patch:
diff --git a/benchmarks/benchmarks/spatial.py b/benchmarks/benchmarks/spatial.py
index f3005f6d70bb..0825bd5dcfa8 100644
--- a/benchmarks/benchmarks/spatial.py
+++ b/benchmarks/benchmarks/spatial.py
@@ -106,27 +106,74 @@ class Neighbors(Benchmark):
         [1, 2, np.inf],
         [0.2, 0.5],
         BOX_SIZES, LEAF_SIZES,
+        ['cKDTree', 'cKDTree_weighted'],
     ]
-    param_names = ['(m, n1, n2)', 'p', 'probe radius', 'boxsize', 'leafsize']
+    param_names = ['(m, n1, n2)', 'p', 'probe radius', 'boxsize', 'leafsize', 'cls']
 
-    def setup(self, mn1n2, p, probe_radius, boxsize, leafsize):
+    def setup(self, mn1n2, p, probe_radius, boxsize, leafsize, cls):
         m, n1, n2 = mn1n2
 
         self.data1 = np.random.uniform(size=(n1, m))
         self.data2 = np.random.uniform(size=(n2, m))
 
+        self.w1 = np.ones(n1)
+        self.w2 = np.ones(n2)
+
         self.T1 = cKDTree(self.data1, boxsize=boxsize, leafsize=leafsize)
         self.T2 = cKDTree(self.data2, boxsize=boxsize, leafsize=leafsize)
 
-    def time_sparse_distance_matrix(self, mn1n2, p, probe_radius, boxsize, leafsize):
+    def time_sparse_distance_matrix(self, mn1n2, p, probe_radius, boxsize, leafsize, cls):
         self.T1.sparse_distance_matrix(self.T2, probe_radius, p=p)
 
-    def time_count_neighbors(self, mn1n2, p, probe_radius, boxsize, leafsize):
+    def time_count_neighbors(self, mn1n2, p, probe_radius, boxsize, leafsize, cls):
         """
         Count neighbors kd-tree
-        dim | # points T1 | # points T2 | p | probe radius |  BoxSize | LeafSize
+        dim | # points T1 | # points T2 | p | probe radius |  BoxSize | LeafSize | cls
+        """
+
+        if cls != 'cKDTree_weighted':
+            self.T1.count_neighbors(self.T2, probe_radius, p=p)
+        else:
+            self.T1.count_neighbors(self.T2, probe_radius, self_weights=self.w1, other_weights=self.w2, p=p)
+
+class CNeighbors(Benchmark):
+    params = [
+        [
+          (2,1000,1000),
+          (8,1000,1000),
+          (16,1000,1000)
+        ],
+        [2, 10, 100, 400, 1000],
+    ]
+    param_names = ['(m, n1, n2)', 'Nr']
+
+    def setup(self, mn1n2, Nr):
+        m, n1, n2 = mn1n2
+
+        data1 = np.random.uniform(size=(n1, m))
+        data2 = np.random.uniform(size=(n2, m))
+        self.w1 = np.ones(len(data1))
+        self.w2 = np.ones(len(data2))
+ 
+        self.T1d = cKDTree(data1, leafsize=1)
+        self.T2d = cKDTree(data2, leafsize=1)
+        self.T1s = cKDTree(data1, leafsize=8)
+        self.T2s = cKDTree(data2, leafsize=8)
+        self.r = np.linspace(0, 0.5, Nr)
+
+    def time_count_neighbors_deep(self, mn1n2, Nr):
+        """
+        Count neighbors for a very deep kd-tree
+        dim | # points T1 | # points T2 | Nr
+        """
+        self.T1d.count_neighbors(self.T2d, self.r)
+
+    def time_count_neighbors_shallow(self, mn1n2, Nr):
+        """
+        Count neighbors for a shallow kd-tree
+        dim | # points T1 | # points T2 | Nr
         """
-        self.T1.count_neighbors(self.T2, probe_radius, p=p)
+        self.T1s.count_neighbors(self.T2s, self.r)
 
 def generate_spherical_points(num_points):
         # generate uniform points on sphere (see:
diff --git a/scipy/spatial/ckdtree.pyx b/scipy/spatial/ckdtree.pyx
index 9096c580a849..6cf244e6c46a 100644
--- a/scipy/spatial/ckdtree.pyx
+++ b/scipy/spatial/ckdtree.pyx
@@ -346,6 +346,10 @@ cdef extern from "ckdtree_methods.h":
                          np.float64_t *mins, 
                          int _median, 
                          int _compact)
+
+    object build_weights(ckdtree *self, 
+                         np.float64_t *node_weights,
+                         np.float64_t *weights)
        
     object query_knn(const ckdtree *self, 
                      np.float64_t *dd, 
@@ -364,15 +368,27 @@ cdef extern from "ckdtree_methods.h":
                        const np.float64_t p, 
                        const np.float64_t eps,
                        vector[ordered_pair] *results)
-                       
-    object count_neighbors(const ckdtree *self,
+
+    object count_neighbors_unweighted(const ckdtree *self,
                            const ckdtree *other,
                            np.intp_t     n_queries,
                            np.float64_t  *real_r,
                            np.intp_t     *results,
-                           np.intp_t     *idx,
-                           const np.float64_t p)
-                           
+                           const np.float64_t p,
+                           int cumulative)
+
+    object count_neighbors_weighted(const ckdtree *self,
+                           const ckdtree *other,
+                           np.float64_t  *self_weights,
+                           np.float64_t  *other_weights,
+                           np.float64_t  *self_node_weights,
+                           np.float64_t  *other_node_weights,
+                           np.intp_t     n_queries,
+                           np.float64_t  *real_r,
+                           np.float64_t     *results,
+                           const np.float64_t p,
+                           int cumulative)
+
     object query_ball_point(const ckdtree *self,
                             const np.float64_t *x,
                             const np.float64_t r,
@@ -438,8 +454,9 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
         brute-force. Default: 16.
     compact_nodes : bool, optional    
         If True, the kd-tree is built to shrink the hyperrectangles to
-        the actual data range. This usually gives a more compact tree and 
-        faster queries at the expense of longer build time. Default: True.
+        the actual data range. This usually gives a more compact tree that 
+        is robust against degenerated input data and gives faster queries 
+        at the expense of longer build time. Default: True.
     copy_data : bool, optional
         If True the data is always copied to protect the kd-tree against 
         data corruption. Default: False.
@@ -454,6 +471,11 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
         into :math:`[0, L_i)`. A ValueError is raised if any of the data is
         outside of this bound.
 
+    Attributes
+    ----------
+    size : int
+        The number of nodes in the tree.
+
     See Also
     --------
     KDTree : Implementation of `cKDTree` in pure Python
@@ -477,6 +499,7 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
         readonly object          boxsize
         np.ndarray               boxsize_data
         np.float64_t             *raw_boxsize_data
+        readonly np.intp_t       size
 
     def __cinit__(cKDTree self):
         self.tree_buffer = NULL        
@@ -499,7 +522,6 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
 
         if boxsize is None:
             self.boxsize = None
-            self.raw_boxsize_data = NULL
             self.boxsize_data = None
         else:
             boxsize_arr = np.empty(2 * self.m, dtype=np.float64)
@@ -507,7 +529,6 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
             boxsize_arr[self.m:] = 0.5 * boxsize_arr[:self.m]
             # FIXME: how to use a matching del if new is used?
             self.boxsize_data = boxsize_arr
-            self.raw_boxsize_data = <np.float64_t*> np.PyArray_DATA(boxsize_arr)
             self.boxsize = boxsize_arr[:self.m].copy()
             if (self.data >= self.boxsize[None, :]).any():
                 raise ValueError("Some input data are greater than the size of the periodic box.")
@@ -518,10 +539,7 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
         self.mins = np.ascontiguousarray(np.amin(self.data,axis=0), dtype=np.float64)
         self.indices = np.ascontiguousarray(np.arange(self.n,dtype=np.intp))
 
-        self.raw_data = <np.float64_t*> np.PyArray_DATA(self.data)
-        self.raw_maxes = <np.float64_t*> np.PyArray_DATA(self.maxes)
-        self.raw_mins = <np.float64_t*> np.PyArray_DATA(self.mins)
-        self.raw_indices = <np.intp_t*> np.PyArray_DATA(self.indices)
+        self._pre_init()
 
         _compact = 1 if compact_nodes else 0
         _median = 1 if balanced_tree else 0
@@ -539,12 +557,11 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
                 _median, _compact)
         finally:
             PyMem_Free(tmp)
-                
+
         self._median_workspace = None
         
         # set up the tree structure pointers
-        self.ctree = tree_buffer_root(self.tree_buffer)
-        self._post_init(self.ctree)
+        self._post_init()
         
         # make the tree viewable from Python
         self.tree = cKDTreeNode()
@@ -553,9 +570,32 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
         self.tree._indices = self.indices
         self.tree.level = 0
         self.tree._setup()
+
+    cdef int _pre_init(cKDTree self) except -1:
+
+        # finalize the pointers from array attributes
+
+        self.raw_data = <np.float64_t*> np.PyArray_DATA(self.data)
+        self.raw_maxes = <np.float64_t*> np.PyArray_DATA(self.maxes)
+        self.raw_mins = <np.float64_t*> np.PyArray_DATA(self.mins)
+        self.raw_indices = <np.intp_t*> np.PyArray_DATA(self.indices)
+
+        if self.boxsize_data is not None:
+            self.raw_boxsize_data = <np.float64_t*>np.PyArray_DATA(self.boxsize_data)
+
+        return 0        
+
+    cdef int _post_init(cKDTree self) except -1:
+        # finalize the tree points, this calls _post_init_traverse
         
-                
-    cdef int _post_init(cKDTree self, ckdtreenode *node) except -1:
+        self.ctree = tree_buffer_root(self.tree_buffer)
+
+        # set the size attribute after tree_buffer is built
+        self.size = self.tree_buffer.size()
+
+        return self._post_init_traverse(self.ctree)
+         
+    cdef int _post_init_traverse(cKDTree self, ckdtreenode *node) except -1:
         # recurse the tree and re-initialize
         # "less" and "greater" fields
         if node.split_dim == -1:
@@ -565,8 +605,9 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
         else:
             node.less = self.ctree + node._less
             node.greater = self.ctree + node._greater
-            self._post_init(node.less)
-            self._post_init(node.greater)
+            self._post_init_traverse(node.less)
+            self._post_init_traverse(node.greater)
+                
         return 0
         
 
@@ -617,18 +658,18 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
         -------
         d : array of floats
             The distances to the nearest neighbors. 
-            If x has shape tuple+(self.m,), then d has shape tuple+(len(k),).
+            If ``x`` has shape ``tuple+(self.m,)``, then ``d`` has shape ``tuple+(k,)``.
             When k == 1, the last dimension of the output is squeezed.
             Missing neighbors are indicated with infinite distances.
         i : ndarray of ints
-            The locations of the neighbors in self.data.
-            If `x` has shape tuple+(self.m,), then `i` has shape tuple+(len(k),).
+            The locations of the neighbors in ``self.data``.
+            If ``x`` has shape ``tuple+(self.m,)``, then ``i`` has shape ``tuple+(k,)``.
             When k == 1, the last dimension of the output is squeezed.
-            Missing neighbors are indicated with self.n.
+            Missing neighbors are indicated with ``self.n``.
 
         Notes
         -----
-        If the KD-Tree is periodic, the position :py:code:`x` is wrapped into the
+        If the KD-Tree is periodic, the position ``x`` is wrapped into the
         box.
         
         When the input k is a list, a query for arange(max(k)) is performed, but
@@ -1045,7 +1086,7 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
         r : positive float
             The maximum distance.
         p : float, optional
-            Which Minkowski norm to use.  `p` has to meet the condition
+            Which Minkowski norm to use.  ``p`` has to meet the condition
             ``1 <= p <= infinity``.
         eps : float, optional
             Approximate search.  Branches of the tree are not explored
@@ -1076,46 +1117,183 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
         else:
             raise ValueError("Invalid output type") 
 
+    def _build_weights(cKDTree self, object weights):
+        """
+        _build_weights(weights)
+
+        Compute weights of nodes from weights of data points. This will sum
+        up the total weight per node. This function is used internally.
+
+        Parameters
+        ----------
+        weights : array_like
+            weights of data points; must be the same length as the data points.
+            currently only scalar weights are supported. Therefore the weights
+            array must be 1 dimensional.
+
+        Returns
+        -------
+        node_weights : array_like
+            total weight for each KD-Tree node.
+
+        """
+        cdef np.intp_t num_of_nodes
+        cdef np.ndarray[np.float64_t, ndim=1, mode="c"] node_weights
+        cdef np.ndarray[np.float64_t, ndim=1, mode="c"] proper_weights
+
+        num_of_nodes = self.tree_buffer.size();
+        node_weights = np.empty(num_of_nodes, dtype=np.float64)
+
+        # FIXME: use templates to avoid the type conversion 
+        proper_weights = np.ascontiguousarray(weights, dtype=np.float64)
+
+        if len(proper_weights) != self.n:
+            raise ValueError('Number of weights differ from the number of data points')
+
+        build_weights(<ckdtree*> self, <np.float64_t*>np.PyArray_DATA(node_weights),
+                            <np.float64_t*> np.PyArray_DATA(proper_weights))
+
+        return node_weights
 
     # ---------------
     # count_neighbors
     # ---------------
 
     @cython.boundscheck(False)
-    def count_neighbors(cKDTree self, cKDTree other, object r, np.float64_t p=2.):
+    def count_neighbors(cKDTree self, cKDTree other, object r, np.float64_t p=2., 
+                        object weights=None, int cumulative=True):
         """
-        count_neighbors(self, other, r, p=2.)
+        count_neighbors(self, other, r, p=2., weights=None, cumulative=True)
 
-        Count how many nearby pairs can be formed.
+        Count how many nearby pairs can be formed. (pair-counting)
 
         Count the number of pairs (x1,x2) can be formed, with x1 drawn
-        from self and x2 drawn from `other`, and where
+        from self and x2 drawn from ``other``, and where
         ``distance(x1, x2, p) <= r``.
-        This is the "two-point correlation" described in Gray and Moore 2000,
-        "N-body problems in statistical learning", and the code here is based
-        on their algorithm.
+
+        Data points on self and other are optionally weighted by the ``weights``
+        argument. (See below)
+
+        The algorithm we implement here is based on [1]_. See notes for further discussion.
 
         Parameters
         ----------
         other : cKDTree instance
-            The other tree to draw points from.
+            The other tree to draw points from, can be the same tree as self.
         r : float or one-dimensional array of floats
             The radius to produce a count for. Multiple radii are searched with
-            a single tree traversal.
+            a single tree traversal. 
+            If the count is non-cumulative(``cumulative=False``), ``r`` defines 
+            the edges of the bins, and must be non-decreasing.
         p : float, optional 
-            1<=p<=infinity, default 2.0
-            Which Minkowski p-norm to use
+            1<=p<=infinity. 
+            Which Minkowski p-norm to use.
+            Default 2.0.
+        weights : tuple, array_like, or None, optional
+            If None, the pair-counting is unweighted.
+            If given as a tuple, weights[0] is the weights of points in ``self``, and
+            weights[1] is the weights of points in ``other``; either can be None to 
+            indicate the points are unweighted.
+            If given as an array_like, weights is the weights of points in ``self``
+            and ``other``. For this to make sense, ``self`` and ``other`` must be the
+            same tree. If ``self`` and ``other`` are two different trees, a ``ValueError``
+            is raised.
+            Default: None
+        cumulative : bool, optional
+            Whether the returned counts are cumulative. When cumulative is set to ``False``
+            the algorithm is optimized to work with a large number of bins (>10) specified
+            by ``r``. When ``cumulative`` is set to True, the algorithm is optimized to work
+            with a small number of ``r``. Default: True
 
         Returns
         -------
-        result : int or 1-D array of ints
-            The number of pairs.
+        result : scalar or 1-D array
+            The number of pairs. For unweighted counts, the result is integer.
+            For weighted counts, the result is float.
+            If cumulative is False, ``result[i]`` contains the counts with
+            ``(-inf if i == 0 else r[i-1]) < R <= r[i]``
+
+        Notes
+        -----
+        Pair-counting is the basic operation used to calculate the two point
+        correlation functions from a data set composed of position of objects.
+
+        Two point correlation function measures the clustering of objects and
+        is widely used in cosmology to quantify the large scale structure
+        in our Universe, but it may be useful for data analysis in other fields
+        where self-similar assembly of objects also occur.
+
+        The Landy-Szalay estimator for the two point correlation function of
+        ``D`` measures the clustering signal in ``D``. [2]_
+
+        For example, given the position of two sets of objects,
+
+        - objects ``D`` (data) contains the clustering signal, and
+
+        - objects ``R`` (random) that contains no signal,
+
+        .. math::
+
+             \\xi(r) = \\frac{<D, D> - 2 f <D, R> + f^2<R, R>}{f^2<R, R>},
+
+        where the brackets represents counting pairs between two data sets
+        in a finite bin around ``r`` (distance), corresponding to setting
+        `cumulative=False`, and ``f = float(len(D)) / float(len(R))`` is the
+        ratio between number of objects from data and random.
+
+        The algorithm implemented here is loosely based on the dual-tree
+        algorithm described in [1]_. We switch between two different
+        pair-cumulation scheme depending on the setting of ``cumulative``.
+        The computing time of the method we use when for
+        ``cumulative == False`` does not scale with the total number of bins.
+        The algorithm for ``cumulative == True`` scales linearly with the
+        number of bins, though it is slightly faster when only
+        1 or 2 bins are used. [5]_.
+
+        As an extension to the naive pair-counting,
+        weighted pair-counting counts the product of weights instead
+        of number of pairs.
+        Weighted pair-counting is used to estimate marked correlation functions
+        ([3]_, section 2.2),
+        or to properly calculate the average of data per distance bin
+        (e.g. [4]_, section 2.1 on redshift).
+
+        .. [1] Gray and Moore,
+               "N-body problems in statistical learning",
+               Mining the sky, 2000,
+               https://arxiv.org/abs/astro-ph/0012333
+
+        .. [2] Landy and Szalay,
+               "Bias and variance of angular correlation functions",
+               The Astrophysical Journal, 1993,
+               http://adsabs.harvard.edu/abs/1993ApJ...412...64L
+
+        .. [3] Sheth, Connolly and Skibba,
+               "Marked correlations in galaxy formation models",
+               Arxiv e-print, 2005,
+               https://arxiv.org/abs/astro-ph/0511773
+
+        .. [4] Hawkins, et al.,
+               "The 2dF Galaxy Redshift Survey: correlation functions,
+               peculiar velocities and the matter density of the Universe",
+               Monthly Notices of the Royal Astronomical Society, 2002,
+               http://adsabs.harvard.edu/abs/2003MNRAS.346...78H
+
+        .. [5] https://github.com/scipy/scipy/pull/5647#issuecomment-168474926
+
         """
         cdef: 
             int r_ndim
             np.intp_t n_queries, i
             np.ndarray[np.float64_t, ndim=1, mode="c"] real_r
-            np.ndarray[np.intp_t, ndim=1, mode="c"] results, idx
+            np.ndarray[np.float64_t, ndim=1, mode="c"] fresults
+            np.ndarray[np.intp_t, ndim=1, mode="c"] iresults
+            np.ndarray[np.float64_t, ndim=1, mode="c"] w1, w1n
+            np.ndarray[np.float64_t, ndim=1, mode="c"] w2, w2n
+            np.float64_t *w1p
+            np.float64_t *w1np
+            np.float64_t *w2p
+            np.float64_t *w2np
 
         # Make sure trees are compatible
         if self.m != other.m:
@@ -1129,6 +1307,10 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
             raise ValueError("r must be either a single value or a "
                              "one-dimensional array of values")
         real_r = np.array(r, ndmin=1, dtype=np.float64, copy=True)
+        if not cumulative:
+            if (real_r[:-1] > real_r[1:]).any():
+                raise ValueError("r must be non-decreasing for non-cumulative counting.");
+        real_r, uind, inverse = np.unique(real_r, return_inverse=True, return_index=True)
         n_queries = real_r.shape[0]
 
         # Internally, we represent all distances as distance ** p
@@ -1137,21 +1319,69 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
                 if not ckdtree_isinf(real_r[i]):
                     real_r[i] = real_r[i] ** p
 
-        results = np.zeros(n_queries, dtype=np.intp)
-        idx = np.arange(n_queries, dtype=np.intp)
-        
-        count_neighbors(<ckdtree*> self, <ckdtree*> other, n_queries,
-                        &real_r[0], &results[0], &idx[0], p)
-        
+        if weights is None:
+            self_weights = other_weights = None
+        elif isinstance(weights, tuple):
+            self_weights, other_weights = weights
+        else:
+            self_weights = other_weights = weights
+            if other is not self:
+                raise ValueError("Two different trees are used. Specify weights for both in a tuple.")
+
+        if self_weights is None and other_weights is None:
+            int_result = True
+            # unweighted, use the integer arithmetics
+            results = np.zeros(n_queries + 1, dtype=np.intp)
+
+            iresults = results
+            count_neighbors_unweighted(<ckdtree*> self, <ckdtree*> other, n_queries,
+                            &real_r[0], &iresults[0], p, cumulative)
+
+        else:
+            int_result = False
+            # weighted / half weighted, use the floating point arithmetics
+            if self_weights is not None:
+                w1 = np.ascontiguousarray(self_weights, dtype=np.float64)
+                w1n = self._build_weights(w1)
+                w1p = <np.float64_t*> np.PyArray_DATA(w1)
+                w1np = <np.float64_t*> np.PyArray_DATA(w1n)
+            else:
+                w1p = NULL
+                w1np = NULL
+            if other_weights is not None:
+                w2 = np.ascontiguousarray(other_weights, dtype=np.float64)
+                w2n = other._build_weights(w2)
+                w2p = <np.float64_t*> np.PyArray_DATA(w2)
+                w2np = <np.float64_t*> np.PyArray_DATA(w2n)
+            else:
+                w2p = NULL
+                w2np = NULL
+
+            results = np.zeros(n_queries + 1, dtype=np.float64)
+            fresults = results
+            count_neighbors_weighted(<ckdtree*> self, <ckdtree*> other,
+                                    w1p, w2p, w1np, w2np,
+                                    n_queries,
+                                    &real_r[0], &fresults[0], p, cumulative)
+
+        results2 = np.zeros(inverse.shape, results.dtype)
+        if cumulative:
+            # copy out the results (taking care of duplication and sorting)
+            results2[...] = results[inverse]
+        else:
+            # keep the identical ones zero
+            # this could have been done in a more readable way.
+            results2[uind] = results[inverse][uind]
+        results = results2
+
         if r_ndim == 0:
-            if results[0] <= <np.intp_t> LONG_MAX:
+            if int_result and results[0] <= <np.intp_t> LONG_MAX:
                 return int(results[0])
             else:
                 return results[0]
         else:
             return results
-
-
+    
     # ----------------------
     # sparse_distance_matrix
     # ----------------------
@@ -1236,21 +1466,15 @@ cdef public class cKDTree [object ckdtree, type ckdtree_type]:
         # unpack the state
         (tree, self.data, self.n, self.m, self.leafsize, 
             self.maxes, self.mins, self.indices, self.boxsize, self.boxsize_data) = state
+
+        # set raw pointers
+        self._pre_init()
         
         # copy kd-tree buffer 
         unpickle_tree_buffer(self.tree_buffer, tree)    
         
-        # set raw pointers
-        self.raw_data = <np.float64_t*>np.PyArray_DATA(self.data)
-        self.raw_maxes = <np.float64_t*>np.PyArray_DATA(self.maxes)
-        self.raw_mins = <np.float64_t*>np.PyArray_DATA(self.mins)
-        self.raw_indices = <np.intp_t*>np.PyArray_DATA(self.indices)
-        if self.boxsize_data is not None:
-            self.raw_boxsize_data = <np.float64_t*>np.PyArray_DATA(self.boxsize_data)
- 
         # set up the tree structure pointers
-        self.ctree = tree_buffer_root(self.tree_buffer)
-        self._post_init(self.ctree)
+        self._post_init()
         
         # make the tree viewable from Python
         self.tree = cKDTreeNode()
diff --git a/scipy/spatial/ckdtree/src/build.cxx b/scipy/spatial/ckdtree/src/build.cxx
index 0e3236c7d6f9..1770fcc80722 100644
--- a/scipy/spatial/ckdtree/src/build.cxx
+++ b/scipy/spatial/ckdtree/src/build.cxx
@@ -18,18 +18,11 @@
 #include "ordered_pair.h"
 #include "ckdtree_methods.h"
 #include "cpp_exc.h"
+#include "cpp_utils.h"
 #include "partial_sort.h"
 
 
 
-inline ckdtreenode*
-tree_buffer_root(std::vector<ckdtreenode> *buf)
-{
-    std::vector<ckdtreenode> &tmp = *buf;
-    return &tmp[0];
-}
-
-
 static npy_intp 
 build(ckdtree *self, npy_intp start_idx, npy_intp end_idx,
       npy_float64 *maxes, npy_float64 *mins, 
@@ -50,13 +43,15 @@ build(ckdtree *self, npy_intp start_idx, npy_intp end_idx,
     node_index = self->tree_buffer->size() - 1;
     root = tree_buffer_root(self->tree_buffer);
     n = root + node_index;
-                
+    memset(n, 0, sizeof(n[0]));
+
+    n->start_idx = start_idx;
+    n->end_idx = end_idx;
+    n->children = end_idx - start_idx;
+    
     if (end_idx-start_idx <= self->leafsize) {
         /* below brute force limit, return leafnode */
         n->split_dim = -1;
-        n->children = end_idx - start_idx;
-        n->start_idx = start_idx;
-        n->end_idx = end_idx;
         return node_index;
     }
     else {
@@ -99,9 +94,6 @@ build(ckdtree *self, npy_intp start_idx, npy_intp end_idx,
              * return leafnode
              */
             n->split_dim = -1;
-            n->children = end_idx - start_idx;
-            n->start_idx = start_idx;
-            n->end_idx = end_idx;
             return node_index;
         }   
         
@@ -199,7 +191,6 @@ build(ckdtree *self, npy_intp start_idx, npy_intp end_idx,
         n->_greater = _greater;
         n->less = root + _less;
         n->greater = root + _greater;
-        n->children = n->less->children + n->greater->children;      
         n->split_dim = d;
         n->split = split;
         
@@ -236,3 +227,67 @@ build_ckdtree(ckdtree *self, npy_intp start_idx, npy_intp end_idx,
         Py_RETURN_NONE;
     }
 }
+
+static npy_float64
+add_weights(ckdtree *self, 
+           npy_float64 *node_weights, 
+           npy_intp node_index, 
+           npy_float64 *weights)
+{
+      
+    npy_intp *indices = (npy_intp *)(self->raw_indices);
+
+    ckdtreenode *n, *root;
+
+    root = tree_buffer_root(self->tree_buffer); 
+
+    n = root + node_index;
+    
+    npy_float64 sum = 0;
+
+    if (n->split_dim != -1) {
+        /* internal nodes; recursively calculate the total weight */
+        npy_float64 left, right;
+        left = add_weights(self, node_weights, n->_less, weights);
+        right = add_weights(self, node_weights, n->_greater, weights);
+        sum = left + right;
+    } else {
+        npy_intp i;
+
+        /* Leaf nodes */
+        for (i = n->start_idx; i < n->end_idx; ++i) {
+            sum += weights[indices[i]];
+        }
+    }
+
+    node_weights[node_index] = sum;
+    return sum;
+}
+        
+        
+extern "C" PyObject*
+build_weights (ckdtree *self, npy_float64 *node_weights, npy_float64 *weights)
+{
+    
+    /* release the GIL */
+    NPY_BEGIN_ALLOW_THREADS
+    {
+        try {
+            add_weights(self, node_weights, 0, weights);        
+        } 
+        catch(...) {
+            translate_cpp_exception_with_gil();
+        }
+    }
+    /* reacquire the GIL */
+    NPY_END_ALLOW_THREADS
+
+    if (PyErr_Occurred()) 
+        /* true if a C++ exception was translated */
+        return NULL;
+    else {
+        /* return None if there were no errors */
+        Py_RETURN_NONE;
+    }
+}
+
diff --git a/scipy/spatial/ckdtree/src/ckdtree_decl.h b/scipy/spatial/ckdtree/src/ckdtree_decl.h
index 72df1924f734..74aa79f66860 100644
--- a/scipy/spatial/ckdtree/src/ckdtree_decl.h
+++ b/scipy/spatial/ckdtree/src/ckdtree_decl.h
@@ -26,7 +26,6 @@ struct ckdtreenode {
 };
 
 #ifdef CKDTREE_METHODS_IMPL
-
 struct ckdtree {
     PyObject_HEAD
     // vtab pointer is present as long as cKDTree has cdef methods
@@ -51,7 +50,7 @@ struct ckdtree {
     const PyObject      *boxsize;
     const PyArrayObject *boxsize_data;
     const npy_float64   *raw_boxsize_data;
+    const npy_intp size;
 };
-
 #endif
 #endif
diff --git a/scipy/spatial/ckdtree/src/ckdtree_methods.h b/scipy/spatial/ckdtree/src/ckdtree_methods.h
index f07723abec50..78190861ed51 100644
--- a/scipy/spatial/ckdtree/src/ckdtree_methods.h
+++ b/scipy/spatial/ckdtree/src/ckdtree_methods.h
@@ -232,6 +232,8 @@ CKDTREE_EXTERN PyObject*
 build_ckdtree(ckdtree *self, npy_intp start_idx, npy_intp end_idx,
               npy_float64 *maxes, npy_float64 *mins, int _median, int _compact);
 
+extern "C" PyObject*
+build_weights (ckdtree *self, npy_float64 *node_weights, npy_float64 *weights);
 
 /* Query methods in C++ for better speed and GIL release */
 
@@ -256,13 +258,26 @@ query_pairs(const ckdtree *self,
             std::vector<ordered_pair> *results);
             
 CKDTREE_EXTERN PyObject*
-count_neighbors(const ckdtree *self,
+count_neighbors_unweighted(const ckdtree *self,
                 const ckdtree *other,
                 npy_intp n_queries,
                 npy_float64 *real_r,
                 npy_intp *results,
-                npy_intp *idx, 
-                const npy_float64 p);
+                const npy_float64 p,
+                int cumulative);
+
+CKDTREE_EXTERN PyObject*
+count_neighbors_weighted(const ckdtree *self,
+                const ckdtree *other,
+                npy_float64 *self_weights, 
+                npy_float64 *other_weights, 
+                npy_float64 *self_node_weights, 
+                npy_float64 *other_node_weights, 
+                npy_intp n_queries,
+                npy_float64 *real_r,
+                npy_float64 *results,
+                const npy_float64 p,
+                int cumulative);
                                
 CKDTREE_EXTERN PyObject*
 query_ball_point(const ckdtree *self,
diff --git a/scipy/spatial/ckdtree/src/count_neighbors.cxx b/scipy/spatial/ckdtree/src/count_neighbors.cxx
index ce7616180d7c..a07e0721d2ef 100644
--- a/scipy/spatial/ckdtree/src/count_neighbors.cxx
+++ b/scipy/spatial/ckdtree/src/count_neighbors.cxx
@@ -6,6 +6,7 @@
 #include <cstring>
 
 #include <vector>
+#include <algorithm>
 #include <string>
 #include <sstream>
 #include <new>
@@ -19,184 +20,249 @@
 #include "cpp_exc.h"
 #include "rectangle.h"
 
-template <typename MinMaxDist> static void
-traverse(const ckdtree *self, const ckdtree *other,
-         npy_intp n_queries, npy_float64 *r, npy_intp *results, npy_intp *idx,
-         const ckdtreenode *node1, const ckdtreenode *node2,
-         RectRectDistanceTracker<MinMaxDist> *tracker)
+struct WeightedTree {
+    const ckdtree *tree;
+    npy_float64 *weights; 
+    npy_float64 *node_weights; 
+};
+
+struct CNBParams
 {
+    npy_float64 *r; 
+    void * results; /* will be casted inside */
+    WeightedTree self, other;
+    int cumulative;
+};
 
-    const ckdtreenode *lnode1;
-    const ckdtreenode *lnode2;
-    npy_float64 d;
-    npy_intp *old_idx;
-    npy_intp old_n_queries, l, i, j;
-    
-    /* 
+template <typename MinMaxDist, typename WeightType, typename ResultType> static void
+traverse(
+    RectRectDistanceTracker<MinMaxDist> *tracker,
+    const CNBParams *params, 
+    npy_float64 *start, npy_float64 *end, 
+    const ckdtreenode *node1,
+    const ckdtreenode *node2)
+{
+    static void (* const next)(RectRectDistanceTracker<MinMaxDist> *tracker,
+            const CNBParams *params, 
+            npy_float64 *start, npy_float64 *end, 
+            const ckdtreenode *node1,
+            const ckdtreenode *node2) = traverse<MinMaxDist, WeightType, ResultType>;
+
+    ResultType *results = (ResultType*) params->results;
+
+    /*
      * Speed through pairs of nodes all of whose children are close
      * and see if any work remains to be done
      */
-    
-    old_idx = idx;
-    
-    std::vector<npy_intp> inner_idx(n_queries);
-    idx = &inner_idx[0];
 
-    old_n_queries = n_queries;
-    n_queries = 0;
-    
-    for (i=0; i<old_n_queries; ++i) {
-        if (tracker->max_distance < r[old_idx[i]])
-            results[old_idx[i]] += node1->children * node2->children;
-        else if (tracker->min_distance <= r[old_idx[i]]) {
-            idx[n_queries] = old_idx[i];
-            ++n_queries;
+    npy_float64 * new_start = std::lower_bound(start, end, tracker->min_distance);
+    npy_float64 * new_end = std::lower_bound(start, end, tracker->max_distance);
+
+
+    /* since max_distance >= min_distance, end < start never happens */
+    if (params->cumulative) {
+        npy_float64 * i;
+        if (new_end != end) {
+            ResultType nn = WeightType::get_weight(&params->self, node1)
+                          * WeightType::get_weight(&params->other, node2);
+
+            for (i = new_end; i < end; ++i) {
+                results[i - params->r] += nn;
+            }
+        }
+        /* any bins larger than end have been correctly counted, thus
+         * thus we can truncate the queries in future of this branch of the traversal*/
+        start = new_start;
+        end = new_end;
+    } else {
+        start = new_start;
+        end = new_end;
+
+        if (end == start) {
+            ResultType nn = WeightType::get_weight(&params->self, node1)
+                          * WeightType::get_weight(&params->other, node2);
+            results[start - params->r] += nn;
         }
     }
-    
-    if (n_queries > 0) {
-        /* OK, need to probe a bit deeper */
-        if (node1->split_dim == -1) {  /* 1 is leaf node */
-            lnode1 = node1;
-            if (node2->split_dim == -1) {  /* 1 & 2 are leaves */
-                lnode2 = node2;
-                const npy_float64 p = tracker->p;
-                const npy_float64 tmd = tracker->max_distance;                
-                const npy_float64 *sdata = self->raw_data;
-                const npy_intp *sindices = self->raw_indices;
-                const npy_float64 *odata = other->raw_data;
-                const npy_intp *oindices = other->raw_indices;
-                const npy_intp m = self->m;
-                const npy_intp start1 = lnode1->start_idx;
-                const npy_intp start2 = lnode2->start_idx;
-                const npy_intp end1 = lnode1->end_idx;
-                const npy_intp end2 = lnode2->end_idx;
-                
-                prefetch_datapoint(sdata + sindices[start1] * m, m);
-                
-                if (start1 < end1)
-                    prefetch_datapoint(sdata + sindices[start1+1] * m, m);
-                                        
-                /* brute-force */
-                for (i = start1; i < end1; ++i) {
-                    
-                    if (i < end1-2)
-                        prefetch_datapoint(sdata + sindices[i+2] * m, m);
-                                      
-                    prefetch_datapoint(odata + oindices[start2] * m, m);
-                        
-                    if (start2 < end2)
-                        prefetch_datapoint(odata + oindices[start2+1] * m, m);
-                  
-                    for (j = start2; j < end2; ++j) {
-                     
-                        if (j < end2-2)
-                            prefetch_datapoint(odata + oindices[j+2] * m, m);
-                 
-                        d = MinMaxDist::distance_p(self,
-                                sdata + sindices[i] * m,
-                                odata + oindices[j] * m,
-                                p, m, tmd);
+
+    if (end == start) {
+        /* this pair falls into exactly one bin, no need to probe deeper. */
+        return;
+    }
+
+    /* OK, need to probe a bit deeper */
+    if (node1->split_dim == -1) {  /* 1 is leaf node */
+        if (node2->split_dim == -1) {  /* 1 & 2 are leaves */
+            npy_intp i, j;
+            const npy_float64 p = tracker->p;
+            const npy_float64 tmd = tracker->max_distance;                
+            const npy_float64 *sdata = params->self.tree->raw_data;
+            const npy_intp *sindices = params->self.tree->raw_indices;
+            const npy_float64 *odata = params->other.tree->raw_data;
+            const npy_intp *oindices = params->other.tree->raw_indices;
+            const npy_intp m = params->self.tree->m;
+            const npy_intp start1 = node1->start_idx;
+            const npy_intp start2 = node2->start_idx;
+            const npy_intp end1 = node1->end_idx;
+            const npy_intp end2 = node2->end_idx;
+
+            prefetch_datapoint(sdata + sindices[start1] * m, m);
+
+            if (start1 < end1 - 1)
+                prefetch_datapoint(sdata + sindices[start1+1] * m, m);
+
+            /* brute-force */
+            for (i = start1; i < end1; ++i) {
+
+                if (i < end1 - 2)
+                    prefetch_datapoint(sdata + sindices[i+2] * m, m);
+
+                prefetch_datapoint(odata + oindices[start2] * m, m);
+
+                if (start2 < end2 - 1)
+                    prefetch_datapoint(odata + oindices[start2+1] * m, m);
+
+                for (j = start2; j < end2; ++j) {
+
+                    if (j < end2 - 2)
+                        prefetch_datapoint(odata + oindices[j+2] * m, m);
+
+                    npy_float64 d = MinMaxDist::distance_p(params->self.tree,
+                            sdata + sindices[i] * m,
+                            odata + oindices[j] * m,
+                            p, m, tmd);
+
+                    if (params->cumulative) {
                         /*
                          * I think it's usually cheaper to test d against all 
                          * r's than to generate a distance array, sort it, then
                          * search for all r's via binary search
                          */
-                        for (l=0; l<n_queries; ++l) {
-                            if (d <= r[idx[l]]) results[idx[l]] += 1;
+                        npy_float64 * l;
+                        for (l = start; l < end; ++l) {
+                            if (d <= *l) {
+                                results[l - params->r] += WeightType::get_weight(&params->self, sindices[i])
+                                                        * WeightType::get_weight(&params->other, sindices[j]);
+                            }
                         }
+                    } else {
+                        const npy_float64 *l = std::lower_bound(start, end, d);
+                        results[l - params->r] += WeightType::get_weight(&params->self, sindices[i])
+                                                * WeightType::get_weight(&params->other, sindices[j]);
                     }
                 }
             }
-            else {  /* 1 is a leaf node, 2 is inner node */
-                tracker->push_less_of(2, node2);
-                traverse(self, other, n_queries, r, results, idx,
-                    node1, node2->less, tracker);
-                tracker->pop();
-
-                tracker->push_greater_of(2, node2);
-                traverse(self, other, n_queries, r, results, idx,
-                    node1, node2->greater, tracker);
-                tracker->pop();
-            }
         }
-        else { /* 1 is an inner node */
-            if (node2->split_dim == -1) {
-                /* 1 is an inner node, 2 is a leaf node */
-                tracker->push_less_of(1, node1);
-                traverse(self, other, n_queries, r, results, idx,
-                    node1->less, node2, tracker);
-                tracker->pop();
+        else {  /* 1 is a leaf node, 2 is inner node */
+            tracker->push_less_of(2, node2);
+            next(tracker, params, start, end, node1, node2->less);
+            tracker->pop();
+
+            tracker->push_greater_of(2, node2);
+            next(tracker, params, start, end, node1, node2->greater);
+            tracker->pop();
+        }
+    }
+    else { /* 1 is an inner node */
+        if (node2->split_dim == -1) {
+            /* 1 is an inner node, 2 is a leaf node */
+            tracker->push_less_of(1, node1);
+            next(tracker, params, start, end, node1->less, node2);
+            tracker->pop();
+            
+            tracker->push_greater_of(1, node1);
+            next(tracker, params, start, end, node1->greater, node2);
+            tracker->pop();
+        }
+        else { /* 1 and 2 are inner nodes */
+            tracker->push_less_of(1, node1);
+            tracker->push_less_of(2, node2);
+            next(tracker, params, start, end, node1->less, node2->less);
+            tracker->pop();
                 
-                tracker->push_greater_of(1, node1);
-                traverse(self, other, n_queries, r, results, idx,
-                    node1->greater, node2, tracker);
-                tracker->pop();
-            }
-            else { /* 1 and 2 are inner nodes */
-                tracker->push_less_of(1, node1);
-                tracker->push_less_of(2, node2);
-                traverse(self, other, n_queries, r, results, idx,
-                    node1->less, node2->less, tracker);
-                tracker->pop();
-                    
-                tracker->push_greater_of(2, node2);
-                traverse(self, other, n_queries, r, results, idx,
-                    node1->less, node2->greater, tracker);
-                tracker->pop();
-                tracker->pop();
-                    
-                tracker->push_greater_of(1, node1);
-                tracker->push_less_of(2, node2);
-                traverse(self, other, n_queries, r, results, idx,
-                    node1->greater, node2->less, tracker);
-                tracker->pop();
-                    
-                tracker->push_greater_of(2, node2);
-                traverse(self, other, n_queries, r, results, idx,
-                    node1->greater, node2->greater, tracker);
-                tracker->pop();
-                tracker->pop();
-            }
+            tracker->push_greater_of(2, node2);
+            next(tracker, params, start, end, node1->less, node2->greater);
+            tracker->pop();
+            tracker->pop();
+                
+            tracker->push_greater_of(1, node1);
+            tracker->push_less_of(2, node2);
+            next(tracker, params, start, end, node1->greater, node2->less);
+            tracker->pop();
+                
+            tracker->push_greater_of(2, node2);
+            next(tracker, params, start, end, node1->greater, node2->greater);
+            tracker->pop();
+            tracker->pop();
         }
     }
 }
 
-
-extern "C" PyObject*
-count_neighbors(const ckdtree *self, const ckdtree *other,
-                npy_intp n_queries, npy_float64 *real_r, npy_intp *results,
-                npy_intp *idx, const npy_float64 p)
+template <typename WeightType, typename ResultType> void
+count_neighbors(struct CNBParams *params,
+                npy_intp n_queries, const npy_float64 p)
 {
 
+    const ckdtree *self = params->self.tree;
+    const ckdtree *other = params->other.tree;
+
 #define HANDLE(cond, kls) \
-    if(cond) { \
+    if (cond) { \
         RectRectDistanceTracker<kls> tracker(self, r1, r2, p, 0.0, 0.0);\
-        traverse(self, other, n_queries, real_r, results, idx, \
-                 self->ctree, other->ctree, &tracker); \
+        traverse<kls, WeightType, ResultType>(&tracker, params, params->r, params->r+n_queries, \
+                 self->ctree, other->ctree); \
     } else
 
+    Rectangle r1(self->m, self->raw_mins, self->raw_maxes);
+    Rectangle r2(other->m, other->raw_mins, other->raw_maxes);
+    
+    if (NPY_LIKELY(self->raw_boxsize_data == NULL)) {
+        HANDLE(NPY_LIKELY(p == 2), MinkowskiDistP2)
+        HANDLE(p == 1, MinkowskiDistP1)
+        HANDLE(ckdtree_isinf(p), MinkowskiDistPinf)
+        HANDLE(1, MinkowskiDistPp) 
+        {}
+    } else {
+        HANDLE(NPY_LIKELY(p == 2), BoxMinkowskiDistP2)
+        HANDLE(p == 1, BoxMinkowskiDistP1)
+        HANDLE(ckdtree_isinf(p), BoxMinkowskiDistPinf)
+        HANDLE(1, BoxMinkowskiDistPp) 
+        {}
+    }
+}
+
+struct Unweighted {
+    /* the interface for accessing weights of unweighted data. */
+    static inline npy_intp
+    get_weight(const WeightedTree *wt, const ckdtreenode * node)
+    {
+        return node->children;
+    }
+    static inline npy_intp
+    get_weight(const WeightedTree *wt, const npy_intp i)
+    {
+        return 1;
+    }
+};
+
+extern "C" PyObject*
+count_neighbors_unweighted(const ckdtree *self, const ckdtree *other,
+                npy_intp n_queries, npy_float64 *real_r, npy_intp *results,
+                const npy_float64 p, int cumulative) {
+
+    CNBParams params = {0};
+
+    params.r = real_r;
+    params.results = (void*) results;
+    params.self.tree = self;
+    params.other.tree = other;
+    params.cumulative = cumulative;
+
     /* release the GIL */
     NPY_BEGIN_ALLOW_THREADS   
     {
         try {
-            
-            Rectangle r1(self->m, self->raw_mins, self->raw_maxes);
-            Rectangle r2(other->m, other->raw_mins, other->raw_maxes);
-            
-            if(NPY_LIKELY(self->raw_boxsize_data == NULL)) {
-                HANDLE(NPY_LIKELY(p == 2), MinkowskiDistP2)
-                HANDLE(p == 1, MinkowskiDistP1)
-                HANDLE(ckdtree_isinf(p), MinkowskiDistPinf)
-                HANDLE(1, MinkowskiDistPp) 
-                {}
-            } else {
-                HANDLE(NPY_LIKELY(p == 2), BoxMinkowskiDistP2)
-                HANDLE(p == 1, BoxMinkowskiDistP1)
-                HANDLE(ckdtree_isinf(p), BoxMinkowskiDistPinf)
-                HANDLE(1, BoxMinkowskiDistPp) 
-                {}
-            }
+            count_neighbors<Unweighted, npy_intp>(&params, n_queries, p);
         } 
         catch(...) {
             translate_cpp_exception_with_gil();
@@ -214,4 +280,66 @@ count_neighbors(const ckdtree *self, const ckdtree *other,
     }
 }
 
+struct Weighted {
+    /* the interface for accessing weights of weighted data. */
+    static inline npy_float64
+    get_weight(const WeightedTree *wt, const ckdtreenode * node)
+    {
+        return (wt->weights != NULL)
+           ? wt->node_weights[node - wt->tree->ctree]
+           : node->children;
+    }
+    static inline npy_float64
+    get_weight(const WeightedTree *wt, const npy_intp i)
+    {
+        return (wt->weights != NULL)?wt->weights[i]:1;
+    }
+};
+
+
+extern "C" PyObject*
+count_neighbors_weighted(const ckdtree *self, const ckdtree *other,
+                npy_float64 *self_weights, npy_float64 *other_weights, 
+                npy_float64 *self_node_weights, npy_float64 *other_node_weights, 
+                npy_intp n_queries, npy_float64 *real_r, npy_float64 *results,
+                const npy_float64 p, int cumulative)
+{
+
+    CNBParams params = {0};
+
+    params.r = real_r;
+    params.results = (void*) results;
+    params.cumulative = cumulative;
+
+    params.self.tree = self;
+    params.other.tree = other;
+    if (self_weights) {
+        params.self.weights = self_weights;
+        params.self.node_weights = self_node_weights;
+    }
+    if (other_weights) {
+        params.other.weights = other_weights;
+        params.other.node_weights = other_node_weights;
+    }
+    /* release the GIL */
+    NPY_BEGIN_ALLOW_THREADS   
+    {
+        try {
+            count_neighbors<Weighted, npy_float64>(&params, n_queries, p);
+        } 
+        catch(...) {
+            translate_cpp_exception_with_gil();
+        }
+    }  
+    /* reacquire the GIL */
+    NPY_END_ALLOW_THREADS
+
+    if (PyErr_Occurred()) 
+        /* true if a C++ exception was translated */
+        return NULL;
+    else {
+        /* return None if there were no errors */
+        Py_RETURN_NONE;
+    }
+}
 
diff --git a/scipy/spatial/ckdtree/src/query.cxx b/scipy/spatial/ckdtree/src/query.cxx
index d67140e80753..00131665eb68 100644
--- a/scipy/spatial/ckdtree/src/query.cxx
+++ b/scipy/spatial/ckdtree/src/query.cxx
@@ -302,12 +302,12 @@ query_single_point(const ckdtree *self,
                 const npy_intp *indices = self->raw_indices;
                 
                 prefetch_datapoint(data+indices[start_idx]*m, m);
-                if (start_idx < end_idx)
+                if (start_idx < end_idx - 1)
                     prefetch_datapoint(data+indices[start_idx+1]*m, m);
                 
                 for (i=start_idx; i<end_idx; ++i) {
                        
-                    if (i < end_idx-2)
+                    if (i < end_idx - 2)
                         prefetch_datapoint(data+indices[i+2]*m, m);
                 
                     d = MinMaxDist::distance_p(self, data+indices[i]*m, x, p, m, distance_upper_bound);
diff --git a/scipy/spatial/ckdtree/src/query_ball_point.cxx b/scipy/spatial/ckdtree/src/query_ball_point.cxx
index a1d081d4b440..7a170d3883ae 100644
--- a/scipy/spatial/ckdtree/src/query_ball_point.cxx
+++ b/scipy/spatial/ckdtree/src/query_ball_point.cxx
@@ -71,12 +71,12 @@ traverse_checking(const ckdtree *self,
         const npy_intp end = lnode->end_idx;
             
         prefetch_datapoint(data + indices[start] * m, m);
-        if (start < end)
+        if (start < end - 1)
             prefetch_datapoint(data + indices[start+1] * m, m);
                 
         for (i = start; i < end; ++i) {
             
-            if (i < end-2)
+            if (i < end -2 )
                 prefetch_datapoint(data + indices[i+2] * m, m);
            
             d = MinMaxDist::distance_p(self, data + indices[i] * m, tpt, p, m, tub);
diff --git a/scipy/spatial/ckdtree/src/query_ball_tree.cxx b/scipy/spatial/ckdtree/src/query_ball_tree.cxx
index 990c653d9308..b25a0e520a0a 100644
--- a/scipy/spatial/ckdtree/src/query_ball_tree.cxx
+++ b/scipy/spatial/ckdtree/src/query_ball_tree.cxx
@@ -99,24 +99,24 @@ traverse_checking(const ckdtree *self, const ckdtree *other,
         
             prefetch_datapoint(sdata + sindices[start1] * m, m);
             
-            if (start1 < end1)
+            if (start1 < end1 - 1)
                 prefetch_datapoint(sdata + sindices[start1+1] * m, m);                                    
             
             for (i = start1; i < end1; ++i) {
             
-                if (i < end1-2)
+                if (i < end1 - 2)
                     prefetch_datapoint(sdata + sindices[i+2] * m, m);
                                   
                 prefetch_datapoint(odata + oindices[start2] * m, m);
                     
-                if (start2 < end2)
+                if (start2 < end2 - 1)
                     prefetch_datapoint(odata + oindices[start2+1] * m, m);
                         
                 results_i = results[sindices[i]];        
                                                                 
                 for (j = start2; j < end2; ++j) {
                 
-                    if (j < end2-2)
+                    if (j < end2 - 2)
                         prefetch_datapoint(odata + oindices[j+2] * m, m);
                 
                     d = MinMaxDist::distance_p(
diff --git a/scipy/spatial/ckdtree/src/query_pairs.cxx b/scipy/spatial/ckdtree/src/query_pairs.cxx
index cc39814f0b3d..166208cb81e4 100644
--- a/scipy/spatial/ckdtree/src/query_pairs.cxx
+++ b/scipy/spatial/ckdtree/src/query_pairs.cxx
@@ -112,12 +112,12 @@ traverse_checking(const ckdtree *self,
             const npy_intp end2 = lnode2->end_idx;
             
             prefetch_datapoint(data+indices[start1]*m, m);
-            if (start1 < end1)
+            if (start1 < end1 - 1)
                prefetch_datapoint(data+indices[start1+1]*m, m);
             
             for(i = start1; i < end1; ++i) {
             
-                if (i < end1-2)
+                if (i < end1 - 2)
                      prefetch_datapoint(data+indices[i+2]*m, m);
                                
                 /* Special care here to avoid duplicate pairs */
@@ -126,13 +126,14 @@ traverse_checking(const ckdtree *self,
                 else
                     min_j = start2;
                                
-                prefetch_datapoint(data+indices[min_j]*m, m);
                 if (min_j < end2)
+                    prefetch_datapoint(data+indices[min_j]*m, m);
+                if (min_j < end2 - 1)
                     prefetch_datapoint(data+indices[min_j+1]*m, m);
                             
                 for (j = min_j; j < end2; ++j) {
                                         
-                    if (j < end2-2)
+                    if (j < end2 - 2)
                         prefetch_datapoint(data+indices[j+2]*m, m);
                                         
                     d = MinMaxDist::distance_p(
diff --git a/scipy/spatial/ckdtree/src/sparse_distances.cxx b/scipy/spatial/ckdtree/src/sparse_distances.cxx
index 939ca56e22be..766b43af9850 100644
--- a/scipy/spatial/ckdtree/src/sparse_distances.cxx
+++ b/scipy/spatial/ckdtree/src/sparse_distances.cxx
@@ -47,21 +47,21 @@ traverse(const ckdtree *self, const ckdtree *other,
             const npy_intp end2 = node2->end_idx;
                         
             prefetch_datapoint(sdata + sindices[start1] * m, m);
-            if (start1 < end1)
+            if (start1 < end1 - 1)
                prefetch_datapoint(sdata + sindices[start1+1] * m, m);                         
                         
             for (npy_intp i = start1; i < end1; ++i) {
             
-                if (i < end1-2)
+                if (i < end1 - 2)
                      prefetch_datapoint(sdata + sindices[i+2] * m, m);
             
                 prefetch_datapoint(odata + oindices[start2] * m, m);
-                if (start2 < end2)
+                if (start2 < end2 - 1)
                     prefetch_datapoint(sdata + oindices[start2+1] * m, m);
                     
                 for (npy_intp j = start2; j < end2; ++j) {
                 
-                    if (j < end2-2)
+                    if (j < end2 - 2)
                         prefetch_datapoint(odata + oindices[j+2] * m, m);
                 
                     npy_float64 d = MinMaxDist::distance_p(


Test patch:
diff --git a/scipy/spatial/tests/test_kdtree.py b/scipy/spatial/tests/test_kdtree.py
index e5f85440fa86..6c1a51772e20 100644
--- a/scipy/spatial/tests/test_kdtree.py
+++ b/scipy/spatial/tests/test_kdtree.py
@@ -4,13 +4,17 @@
 from __future__ import division, print_function, absolute_import
 
 from numpy.testing import (assert_equal, assert_array_equal,
-    assert_almost_equal, assert_array_almost_equal, assert_, run_module_suite)
+    assert_almost_equal, assert_array_almost_equal, assert_,
+    assert_raises,
+    run_module_suite)
 
 import numpy as np
 from scipy.spatial import KDTree, Rectangle, distance_matrix, cKDTree
 from scipy.spatial.ckdtree import cKDTreeNode
 from scipy.spatial import minkowski_distance
 
+import itertools
+
 def distance_box(a, b, p, boxsize):
     diff = a - b
     diff[diff > 0.5 * boxsize] -= boxsize
@@ -1065,19 +1069,11 @@ def test_ckdtree_box():
 
 def test_ckdtree_box_upper_bounds():
     data = np.linspace(0, 2, 10).reshape(-1, 1)
-    try:
-        cKDTree(data, leafsize=1, boxsize=1.0)
-    except ValueError:
-        return
-    raise AssertionError("ValueError is not raised")
+    assert_raises(ValueError, cKDTree, data, leafsize=1, boxsize=1.0)
 
 def test_ckdtree_box_lower_bounds():
     data = np.linspace(-1, 1, 10)
-    try:
-        cKDTree(data, leafsize=1, boxsize=1.0)
-    except ValueError:
-        return
-    raise AssertionError("ValueError is not raised")
+    assert_raises(ValueError, cKDTree, data, leafsize=1, boxsize=1.0)
 
 def simulate_periodic_box(kdtree, data, k, boxsize):
     dd = []
@@ -1105,6 +1101,10 @@ def simulate_periodic_box(kdtree, data, k, boxsize):
     
 def test_ckdtree_memuse():
     # unit test adaptation of gh-5630
+
+    # NOTE: this will fail when run via valgrind,
+    # because rss is no longer a reliable memory usage indicator.
+
     try:
         import resource
     except ImportError:
@@ -1139,6 +1139,71 @@ def test_ckdtree_memuse():
     # ideally zero leaks, but errors might accidentally happen
     # outside cKDTree
     assert_(num_leaks < 10)
+
+def test_ckdtree_weights():
+
+    data = np.linspace(0, 1, 4).reshape(-1, 1)
+    tree1 = cKDTree(data, leafsize=1)
+    weights = np.ones(len(data), dtype='f4')
+
+    nw = tree1._build_weights(weights)
+    assert_array_equal(nw, [4, 2, 1, 1, 2, 1, 1])
+
+    assert_raises(ValueError, tree1._build_weights, weights[:-1])
+
+    for i in range(10):
+        # since weights are uniform, these shall agree:
+        c1 = tree1.count_neighbors(tree1, np.linspace(0, 10, i))
+        c2 = tree1.count_neighbors(tree1, np.linspace(0, 10, i), 
+                weights=(weights, weights))
+        c3 = tree1.count_neighbors(tree1, np.linspace(0, 10, i), 
+                weights=(weights, None))
+        c4 = tree1.count_neighbors(tree1, np.linspace(0, 10, i), 
+                weights=(None, weights))
+        c5 = tree1.count_neighbors(tree1, np.linspace(0, 10, i), 
+                weights=weights)
+
+        assert_array_equal(c1, c2)
+        assert_array_equal(c1, c3)
+        assert_array_equal(c1, c4)
+
+    for i in range(len(data)):
+        # this tests removal of one data point by setting weight to 0
+        w1 = weights.copy()
+        w1[i] = 0
+        data2 = data[w1 != 0]
+        w2 = weights[w1 != 0]
+        tree2 = cKDTree(data2)
+
+        c1 = tree1.count_neighbors(tree1, np.linspace(0, 10, 100),
+                weights=(w1, w1))
+        # "c2 is correct"
+        c2 = tree2.count_neighbors(tree2, np.linspace(0, 10, 100))
+
+        assert_array_equal(c1, c2)
+
+        #this asserts for two different trees, singular weights
+        # crashes
+        assert_raises(ValueError, tree1.count_neighbors,
+            tree2, np.linspace(0, 10, 100), weights=w1)
+
+def test_ckdtree_count_neighbous_multiple_r():
+    n = 2000
+    m = 2
+    np.random.seed(1234)
+    data = np.random.normal(size=(n, m))
+    kdtree = cKDTree(data, leafsize=1)
+    r0 = [0, 0.01, 0.01, 0.02, 0.05]
+    i0 = np.arange(len(r0))
+    n0 = kdtree.count_neighbors(kdtree, r0)
+    nnc = kdtree.count_neighbors(kdtree, r0, cumulative=False)
+    assert_equal(n0, nnc.cumsum())
+
+    for i, r in zip(itertools.permutations(i0), 
+                    itertools.permutations(r0)):
+        # permute n0 by i and it shall agree 
+        n = kdtree.count_neighbors(kdtree, r)
+        assert_array_equal(n, n0[list(i)])
     
 def test_len0_arrays():
     # make sure len-0 arrays are handled correctly
@@ -1186,6 +1251,38 @@ def test_len0_arrays():
     z = np.empty(shape=(0,2), dtype=np.intp)
     assert_array_equal(y, z)
 
+def test_ckdtree_duplicated_inputs():
+    # check ckdtree with duplicated inputs
+    n = 1024 
+    for m in range(1, 8):
+        data = np.concatenate([
+            np.ones((n // 2, m)) * 1,
+            np.ones((n // 2, m)) * 2], axis=0)
+
+        # it shall not divide more than 3 nodes.
+        # root left (1), and right (2)
+        kdtree = cKDTree(data, leafsize=1)
+        assert_equal(kdtree.size, 3)
+
+        kdtree = cKDTree(data)
+        assert_equal(kdtree.size, 3)
+
+        # if compact_nodes are disabled, the number
+        # of nodes is n (per leaf) + (m - 1)* 2 (splits per dimension) + 1
+        # and the root
+        kdtree = cKDTree(data, compact_nodes=False, leafsize=1)
+        assert_equal(kdtree.size, n + m * 2 - 1)
+
+def test_ckdtree_noncumulative_nondecreasing():
+    # check ckdtree with duplicated inputs
+
+    # it shall not divide more than 3 nodes.
+    # root left (1), and right (2)
+    kdtree = cKDTree([[0]], leafsize=1)
+
+    assert_raises(ValueError, kdtree.count_neighbors,
+        kdtree, [0.1, 0], cumulative=False)
+
 def test_short_knn():
 
     # The test case is based on github: #6425 by @SteveDoyle2