Returns indices that perform nearest interpolation given a
    set of points. Similar to scipy.interpolate.griddata with
    method set to "nearest". In fact, it's construction is based
    on that function and related types.

    Parameters:
    -----------

    xs -- ndarray of floats, shape (n,D)
          Data point coordinates. Can either be an array of
          shape (n,D), or a tuple of `ndim` arrays.
    XS -- ndarray of floats, shape (M, D)
          Points at which to interpolate/sample data.

def nearest_indices(xs, XS):
    '''
    Returns indices that perform nearest interpolation given a
    set of points. Similar to scipy.interpolate.griddata with
    method set to "nearest". In fact, it's construction is based
    on that function and related types.

    Parameters:
    -----------

    xs -- ndarray of floats, shape (n,D)
          Data point coordinates. Can either be an array of
          shape (n,D), or a tuple of `ndim` arrays.
    XS -- ndarray of floats, shape (M, D)
          Points at which to interpolate/sample data.
    '''
    xs=_ndim_coords_from_arrays(xs);
    XS = _ndim_coords_from_arrays(XS, ndim=xs.shape[1])
    if XS.shape[-1] != xs.shape[1]:
        raise ValueError("dimensions of the points don't the sample points")
    #the magical kd-tree
    _,i = cKDTree(xs).query(XS)
    return i;