'Query the kd-tree for nearest neighbors
Parameters
x : array_like, last dimension self.m
An array of points to query.
k : int, optional
The number of nearest neighbors to return.
eps : nonnegative float, optional
Return approximate nearest neighbors; the kth returned value
is guaranteed to be no further than (1+eps) times the
distance to the real kth nearest neighbor.
p : float, 1<=p<=infinity, optional
Which Minkowski p-norm to use.
1 is the sum-of-absolute-values "Manhattan" distance
2 is the usual Euclidean distance
infinity is the maximum-coordinate-difference distance
distance_upper_bound : nonnegative float, optional
Return only neighbors within this distance. This is used to prune
tree searches, so if you are doing a series of nearest-neighbor
queries, it may help to supply the distance to the nearest neighbor
of the most recent point.
Returns
d : float or array of floats
The distances to the nearest neighbors.
If x has shape tuple+(self.m,), then d has shape tuple if
k is one, or tuple+(k,) if k is larger than one. Missing
neighbors (e.g. when k > n or distance_upper_bound is
given) are indicated with infinite distances.  If k is None,
then d is an object array of shape tuple, containing lists
of distances. In either case the hits are sorted by distance
(nearest first).
i : integer or array of integers
The locations of the neighbors in self.data. i is the same
shape as d.
Examples
>>> from scipy import spatial
>>> x, y = np.mgrid[0:5, 2:8]
>>> tree = spatial.KDTree(list(zip(x.ravel(), y.ravel())))
>>> tree.data
array([[0, 2],
[0, 3],
[0, 4],
[0, 5],
[0, 6],
[0, 7],
[1, 2],
[1, 3],
[1, 4],
[1, 5],
[1, 6],
[1, 7],
[2, 2],
[2, 3],
[2, 4],
[2, 5],
[2, 6],
[2, 7],
[3, 2],
[3, 3],
[3, 4],
[3, 5],
[3, 6],
[3, 7],
[4, 2],
[4, 3],
[4, 4],
[4, 5],
[4, 6],
[4, 7]])
>>> pts = np.array([[0, 0], [2.1, 2.9]])
>>> tree.query(pts)
(array([ 2.        ,  0.14142136]), array([ 0, 13]))
>>> tree.query(pts[0])
(2.0, 0)'

def query(self, x, k=1, eps=0, p=2, distance_upper_bound=np.inf):

    x = np.asarray(x)
    if (np.shape(x)[(-1)] != self.m):
        raise ValueError(('x    must    consist    of    vectors    of    length    %d    but    has    shape    %s' % (self.m, np.shape(x))))
    if (p < 1):
        raise ValueError('Only    p-norms    with    1<=p<=infinity    permitted')
    retshape = np.shape(x)[:(-1)]
    if (retshape != ()):
        if (k is None):
            dd = np.empty(retshape, dtype=object)
            ii = np.empty(retshape, dtype=object)
        elif (k > 1):
            dd = np.empty((retshape + (k,)), dtype=float)
            dd.fill(np.inf)
            ii = np.empty((retshape + (k,)), dtype=int)
            ii.fill(self.n)
        elif (k == 1):
            dd = np.empty(retshape, dtype=float)
            dd.fill(np.inf)
            ii = np.empty(retshape, dtype=int)
            ii.fill(self.n)
        else:
            raise ValueError('Requested    %s    nearest    neighbors;    acceptable    numbers    are    integers    greater    than    or    equal    to    one,    or    None')
        for c in np.ndindex(retshape):
            hits = self.__query(x[c], k=k, eps=eps, p=p, distance_upper_bound=distance_upper_bound)
            if (k is None):
                dd[c] = [d for (d, i) in hits]
                ii[c] = [i for (d, i) in hits]
            elif (k > 1):
                for j in range(len(hits)):
                    (dd[(c + (j,))], ii[(c + (j,))]) = hits[j]
            elif (k == 1):
                if (len(hits) > 0):
                    (dd[c], ii[c]) = hits[0]
                else:
                    dd[c] = np.inf
                    ii[c] = self.n
        return (dd, ii)
    else:
        hits = self.__query(x, k=k, eps=eps, p=p, distance_upper_bound=distance_upper_bound)
        if (k is None):
            return ([d for (d, i) in hits], [i for (d, i) in hits])
        elif (k == 1):
            if (len(hits) > 0):
                return hits[0]
            else:
                return (np.inf, self.n)
        elif (k > 1):
            dd = np.empty(k, dtype=float)
            dd.fill(np.inf)
            ii = np.empty(k, dtype=int)
            ii.fill(self.n)
            for j in range(len(hits)):
                (dd[j], ii[j]) = hits[j]
            return (dd, ii)
        else:
            raise ValueError('Requested    %s    nearest    neighbors;    acceptable    numbers    are    integers    greater    than    or    equal    to    one,    or    None')
