Python scipy.spatial 模块，KDTree() 实例源码

def computeDistMatrices(emb1, emb2, gold):
    correct = 0
    MRR = float(0.0)
    vob2 = []
    matrix2 = []
    for i in emb2.keys():
        vob2.append(i)
        matrix2.append(emb2[i])
    matrix2 = np.array(matrix2)

    kdtree = KDTree(matrix2, leafsize=100)
    for gold_en, gold_trans in gold.items():
        d, index = kdtree.query(emb1[gold_en], k=10)
        for i in index:
            if vob2[i] in gold_trans:
                correct += 1
                MRR += float(1/(i+1))
                break

    print('\nfinished!\n')
    print('{}/{} % age {}'.format(correct, len(gold.keys()),
          float(correct/len(gold.keys()))))
    print('{}/{} MRR {}'.format(MRR, len(gold.keys()),
          float(MRR/len(gold.keys()))))
    print("dist matrix done!")

def __init__(self, csv_file = './data/king_county_data_geocoded.csv', data = None, values = None):
    if (data is None and csv_file is not None):
      df = pd.read_csv(csv_file)
      self.values = df['AppraisedValue']
      df = df.drop('AppraisedValue', 1)
      df = (df - df.mean()) / (df.max() - df.min())
      self.df = df
      self.df = self.df[['lat', 'long', 'SqFtLot']]

    elif (data is not None and values is not None):
      self.df = data
      self.values = values
    else:
      raise ValueError("Must have either csv_file or data set")

    self.n = len(self.df)

    self.kdtree = KDTree(self.df)
    self.metric = np.mean

    # TODO: set k to a number, try a few numbers out
    # self.k = None

def __init__(self, r, az, sitecoords, proj, x, y, nnear=9):
        self.nnear = nnear
        self.az = az
        self.r = r
        self.x = x
        self.y = y
        # compute the centroid coordinates in lat/lon
        bin_lon, bin_lat = georef.polar2centroids(r, az, sitecoords)
        # reproject the centroids to cartesian map coordinates
        binx, biny = georef.reproject(bin_lon, bin_lat,
                                      projection_target=proj)
        self.binx, self.biny = binx.ravel(), biny.ravel()
        # compute the KDTree
        tree = KDTree(list(zip(self.binx, self.biny)))
        # query the tree for nearest neighbours
        self.dist, self.ix = tree.query(list(zip(x, y)), k=nnear)

def pairwise_angles(lines):

    # line_starts = np.asarray([list(x[0]) for x in lines])
    line_ends = np.asarray([list(x[1]) for x in lines])
    tree = KDTree(line_ends)

    angles = []
    for i, line_A in enumerate(lines):

        # Find nearest line
        dists, j = tree.query(line_A[0])

        if i == j or j > len(lines):
            continue

        line_B = lines[j]

        # Calculate angle between vectors
        a0, a1 = np.asarray(line_A)
        b0, b1 = np.asarray(line_B)
        angles.append(angle_between(a1 - a0, b1 - b0))

    return angles

def within(self, x, ctrs, kdtree=None):
        """Checks which cubes `x` falls within. Uses a K-D Tree to
        accelerate the search if provided."""

        if kdtree is None:
            # If no KDTree is provided, execute a brute-force search
            # over all cubes.
            nctrs = len(ctrs)
            within = np.array([max(abs(ctrs[i] - x)) <= self.hside
                               for i in range(nctrs)], dtype='bool')
            idxs = np.arange(nctrs)[within]
        else:
            # If a KDTree is provided, find all points within r (`hside`).
            idxs = kdtree.query_ball_point(x, self.hside, p=np.inf, eps=0)

        return idxs

def _friends_leaveoneout_radius(points, ftype):
    """Internal method used to compute the radius (half-side-length) for each
    ball (cube) used in :class:`RadFriends` (:class:`SupFriends`) using
    leave-one-out (LOO) cross-validation."""

    # Construct KDTree to enable quick nearest-neighbor lookup for
    # our resampled objects.
    kdtree = spatial.KDTree(points)

    if ftype == 'balls':
        # Compute radius to two nearest neighbors (self + neighbor).
        dists, ids = kdtree.query(points, k=2, eps=0, p=2)
    elif ftype == 'cubes':
        # Compute half-side-length to two nearest neighbors (self + neighbor).
        dists, ids = kdtree.query(points, k=2, eps=0, p=np.inf)

    dist = dists[:, 1]  # distances to LOO nearest neighbor

    return dist

def update(self, pointvol):
        """Update the N-sphere radii using the current set of live points."""

        # Initialize a K-D Tree to assist nearest neighbor searches.
        kdtree = spatial.KDTree(self.live_u)

        # Check if we should use the provided pool for updating.
        if self.use_pool_update:
            pool = self.pool
        else:
            pool = None

        # Update the N-spheres.
        self.radfriends.update(self.live_u, pointvol=pointvol,
                               rstate=self.rstate, bootstrap=self.bootstrap,
                               pool=pool, kdtree=kdtree)
        if self.enlarge != 1.:
            self.radfriends.scale_to_vol(self.radfriends.vol_ball *
                                         self.enlarge)

        return copy.deepcopy(self.radfriends)

def propose_unif(self):
        """Propose a new live point by sampling *uniformly* within
        the union of N-spheres defined by our live points."""

        # Initialize a K-D Tree to assist nearest neighbor searches.
        kdtree = spatial.KDTree(self.live_u)

        while True:
            # Sample a point `u` from the union of N-spheres along with the
            # number of overlapping spheres `q` at point `u`.
            u, q = self.radfriends.sample(self.live_u, rstate=self.rstate,
                                          return_q=True, kdtree=kdtree)

            # Check if our sample is within the unit cube.
            if self._check_unit_cube(u):
                # Accept the point with probability 1/q to account for
                # overlapping balls.
                if q == 1 or self.rstate.rand() < 1.0 / q:
                    break  # if successful, we're done!

        # Define the axes of the N-sphere.
        ax = np.identity(self.npdim) * self.radfriends.radius

        return u, ax

def propose_unif(self):
        """Propose a new live point by sampling *uniformly* within
        the collection of N-cubes defined by our live points."""

        # Initialize a K-D Tree to assist nearest neighbor searches.
        kdtree = spatial.KDTree(self.live_u)

        while True:
            # Sample a point `u` from the union of N-cubes along with the
            # number of overlapping cubes `q` at point `u`.
            u, q = self.supfriends.sample(self.live_u, rstate=self.rstate,
                                          return_q=True, kdtree=kdtree)

            # Check if our point is within the unit cube.
            if self._check_unit_cube(u):
                # Accept the point with probability 1/q to account for
                # overlapping cubes.
                if q == 1 or self.rstate.rand() < 1.0 / q:
                    break  # if successful, we're done!

        # Define the axes of our N-cube.
        ax = np.identity(self.npdim) * self.supfriends.hside

        return u, ax

def buildNNDataStructure(self):
        """Builds a nearest neighbor data structure.  User doesn't need to
        call this unless the self.problems attribute was changed manually."""
        if len(self.problemFeatures)==0 or len(self.featureNames)==0:
            return
        try:
            from sklearn.neighbors import NearestNeighbors,BallTree
            from scipy.spatial import KDTree
            with self.lock:
                try:
                    farray = self.problemFeatures.array
                except AttributeError:
                    farray = np.array(self.problemFeatures.items)
                if self.metricTransform is not None:
                    farray = np.dot(farray,self.metricTransform)
                #self.nn = NearestNeighbors(n_neighbors=1,algorithm="auto").fit(farray)
                self.nn = BallTree(farray)
                #self.nn = KDTree(farray)
                self.nnBuildSize = len(self.problemFeatures)
        except ImportError:
            print "IKDatabase: Warning, scikit-learn is not installed, queries will be much slower"
            with self.lock:
                self.nn = None
                self.nnBuildSize = 0
        return

def warpCloud( xyc, sourceGridPoints, targetGridPoints, warpQuality=9 ):

    sourceTree = KDTree(sourceGridPoints, leafsize=10)
    warpedXYC = []  
    for c in xyc:
        nearestEdge = sourceTree.query(c,k=warpQuality)
        nx = 0.0
        ny = 0.0
        ws = 0.0
        for i in range(warpQuality):
            p = targetGridPoints[nearestEdge[1][i]]
            w = nearestEdge[0][i]
            if w == 0.0:
                nx = p[0]
                ny = p[1]
                ws = 1.0
                break
            else:
                w = 1.0 / w
                nx += w * p[0]
                ny += w * p[1]
                ws += w

        warpedXYC.append([nx/ws,ny/ws])

    warpedXYC = np.array(warpedXYC)
    return warpedXYC

def match_points(current_points, prior_points, distance_cutoff):
    """
    Takes in an nxd input vector of d-dimensional Euclidean coordinates representing the current dataset
    and an mxd input vector of d-dimensional Euclidean cooridnates representing the prior dataset. 
    Output gives, for each row in _current_points, an mx2 vector that gives
      - the index number from _prior_points in the first column, 
      - and distance matched in second.
    If the matched distance is greater than the cutoff, consider the pair unmatched.
    Unmatched points get -1 for the "matched" index number and the cutoff value for the distance (infinity).
    """

    # Initialize matched indices to -1 and distances to the cutoff value.
    matches = np.ones((current_points.shape[0], 2))*-1
    matches[:,1] = distance_cutoff

    # Initialize index numbers for current points
    current_idx = np.asarray(range(current_points.shape[0]))
    prior_idx = np.asarray(range(prior_points.shape[0]))

    # Generate kd trees
    curr_kd_tree = spatial.KDTree(current_points)
    prior_kd_tree = spatial.KDTree(prior_points)

    # Compute closest keypoint from current->prior and from prior->current
    matches_a = prior_kd_tree.query(current_points)
    matches_b = curr_kd_tree.query(prior_points)

    # Mutual matches are the positive matches within the distance cutoff. All others unmatched.
    potential_matches = matches_b[1][matches_a[1]]
    matched_indices = np.equal(potential_matches, current_idx)

    # Filter out matches that are more than the distance cutoff away.
    in_bounds = (matches_a[0] <= distance_cutoff)
    matched_indices = np.multiply(matched_indices, in_bounds)

    # Add the matching data to the output
    matches[current_idx[matched_indices],0] = prior_idx[matches_a[1]][matched_indices].astype(np.int)
    matches[current_idx[matched_indices],1] = matches_a[0][matched_indices]

    return matches

def init_framefield(self):
        """Initialize the frame field based on surface curvature and normals."""
        boundary_frames = []
        boundary_ids = {}
        # The frame field is initialized at the boundary,
        # based on the curvature cross-field and normals.
        for fi, face in enumerate(self.surf_mesh.faces):
            # Retrieve the tet this face belongs to.
            ti = self.tet_mesh.adjacent_elements[self.surf_mesh.face_map.inv[fi]][0]
            tet = self.tet_mesh.elements[ti]
            # Ignore faces which have 0 curvature.
            if np.isclose(self.surf_mesh.k1[face[0]], 0) and np.isclose(self.surf_mesh.k2[face[0]], 0):
                continue
            # @TODO(aidan) Find actual face values, not vertex values.
            uvw = np.hstack((np.vstack(self.surf_mesh.pdir1[face[0]]),
                             np.vstack(self.surf_mesh.pdir2[face[0]]),
                             np.vstack(self.surf_mesh.vertex_normals[face[0]])))
            boundary_frames.append(Frame(uvw, tet_centroid(self.tet_mesh, ti)))
            boundary_ids[ti] = len(boundary_frames) - 1

        # Prepare a KDTree of boundary frame coords for quick spatial queries.
        tree = spatial.KDTree(np.vstack([frame.location for frame in boundary_frames]))

        # Now propagate the boundary frames throughout the tet mesh.
        # Each tet frame takes the value of its nearest boundary tet.
        for ti, tet in enumerate(self.tet_mesh.elements):
            location = tet_centroid(self.tet_mesh, ti)
            if ti in boundary_ids:
                self.frames.append(Frame(boundary_frames[boundary_ids[ti]].uvw, location, True))
            else:
                nearest_ti = tree.query(location)[1] # Find closest boundary frame
                self.frames.append(Frame(boundary_frames[nearest_ti].uvw, location, False))

def _create_histogram(self, image, hist_size, codebook):
        histogram = np.zeros(hist_size)
        descriptors = self.sift.detectAndCompute(image)
        tree = spatial.KDTree(codebook)

        for i in xrange(len(descriptors)):
            histogram[tree.query(descriptors[i])[1]] += 1

        return normalize(histogram[:, np.newaxis], axis=0).ravel()

def create_histogram(image, hist_size=2048, codebook=[], detectAndCompute=SIFT_create().detectAndCompute):
    histogram = np.zeros(hist_size)
    descriptors = detectAndCompute(image, window_size=None)
    tree = spatial.KDTree(codebook)

    for i in xrange(len(descriptors)):
        histogram[tree.query(descriptors[i])[1]] += 1

    return normalize(histogram[:, np.newaxis], axis=0).ravel()

def calculate_rms_with_rotran(atoms1: np.ndarray, atoms2: np.ndarray, rotran: tuple) -> Tuple[float, float, Tuple]:
    # apply rotran on atoms2
    atoms2 = np.dot(atoms2, rotran[0]) + rotran[1]

    # fix atoms1 in place, use KDTree for distances
    kd_tree_a = KDTree(atoms1)
    kd_tree_b = KDTree(atoms2)
    # get nearest neighbours for atoms2 in atoms1 using KDTree, and also the other way around
    distances, indexes_a = kd_tree_b.query(atoms1)
    _, indexes_b = kd_tree_a.query(atoms2)

    # pick only the pairs that are both closest to each other
    closest_indexes = np.where(indexes_b[indexes_a] == list(range(len(indexes_a))))
    smoothing_rotran = calculate_rotran(atoms1[closest_indexes], atoms2[indexes_a[closest_indexes]])
    atoms2 = np.dot(atoms2, smoothing_rotran[0]) + smoothing_rotran[1]

    kd_tree_b = KDTree(atoms2)
    distances, indexes_a = kd_tree_b.query(atoms1)
    _, indexes_b = kd_tree_a.query(atoms2)

    # pick only the pairs that are both closest to each other
    closest_indexes = np.where(indexes_b[indexes_a] == list(range(len(indexes_a))))
    distances = distances[closest_indexes]

    # return RMSD ( sqrt(1/N * SUM(distance(xi, yi)^2)) )
    post_rms = math.sqrt(np.sum(np.power(distances, 2)) / float(len(distances)))

    psi = 100.0 * np.sum(distances <= 4.0) / float(min(len(atoms1), len(atoms2)))

    return post_rms, psi, smoothing_rotran

def generate_neighbourhoods(self):
        coords = self.root.subtree_coords
        kd_tree = KDTree(coords)
        self.root.set_neighbours(kd_tree)

def build_kdtree(self, grid='node'):
        """Builds the kdtree for the specified grid"""
        from scipy.spatial import KDTree
        if not hasattr(self, '_kd_trees'):
            self._kd_trees = {'node': None,
                              'edge1': None,
                              'edge2': None,
                              'center': None}
        lon, lat = self._get_grid_vars(grid)
        if lon is None or lat is None:
            raise ValueError("{0}_lon and {0}_lat must be defined in order to "
                             "create and use KDTree for this grid".format(grid))
        lin_points = np.column_stack((lon.ravel(), lat.ravel()))
        self._kd_trees[grid] = KDTree(lin_points, leafsize=4)

def _generate_invariants(sources):
    """Return an array of (unique) invariants derived from the array `sources`.
Return an array of the indices of `sources` that correspond to each invariant,
arranged as described in _arrangetriplet.
"""
    from scipy.spatial import KDTree
    from itertools import combinations
    from functools import partial
    arrange = partial(_arrangetriplet, sources=sources)

    inv = []
    triang_vrtx = []
    coordtree = KDTree(sources)
    for asrc in sources:
        __, indx = coordtree.query(asrc, NUM_NEAREST_NEIGHBORS)

        # Generate all possible triangles with the 5 indx provided, and store
        # them with the order (a, b, c) defined in _arrangetriplet
        all_asterism_triang = [arrange(vertex_indices=list(cmb))
                               for cmb in combinations(indx, 3)]
        triang_vrtx.extend(all_asterism_triang)

        inv.extend([_invariantfeatures(*sources[triplet])
                    for triplet in all_asterism_triang])

    # Remove here all possible duplicate triangles
    uniq_ind = [pos for (pos, elem) in enumerate(inv)
                if elem not in inv[pos + 1:]]
    inv_uniq = _np.array(inv)[uniq_ind]
    triang_vrtx_uniq = _np.array(triang_vrtx)[uniq_ind]

    return inv_uniq, triang_vrtx_uniq

def test_find_sources(self):
        srcs = aa._find_sources(self.image_ref)

        from scipy.spatial import KDTree
        ref_coordtree = KDTree(self.star_ref_pos)

        # Compare here srcs list with self.star_ref_pos
        num_sources = 0
        for asrc in srcs:
            found_source = ref_coordtree.query_ball_point(asrc, 3)
            if found_source:
                num_sources += 1
        fraction_found = float(num_sources) / float(len(srcs))
        self.assertGreater(fraction_found, 0.85)

def findNearestCell(hf,AreaName,gageCoord):
    # Finds the nearest cell to the gageCoordinate from the HDF file

    CellPts = hf['Geometry']['2D Flow Areas'][AreaName]['Cells Center Coordinate'][:] # All Coordinates which are in 2D Flow Area with name AreaName
    fPtInd = hf['Geometry']['2D Flow Areas'][AreaName]['Cells FacePoint Indexes'] # All FacePoint Indices which correspond to each 'cell center point'
    # the reason this next step is necessary is that CellPts returns coordinates for face points as well as center points, and we need the index of the cell point.
    distances,indices = spatial.KDTree(CellPts).query(gageCoord, k=7) # finds the closest points (cell center or face)
    # Now to eliminate FacePoints and only have cell centers left
    i = 0
    for cell in range(0, len(indices)):
        if fPtInd[indices[i]][2] == -1:
            distances=np.delete(distances,i)
            indices=np.delete(indices,i)
        else:
            i+=1
            continue
    if indices.size: # If the indices list is not empty
        CellInd = indices[np.where(distances==min(distances))] # Choose the index that's left with the shortest distance
        CellDist = distances[np.where(distances==min(distances))]  # Displays the distance left
    else:
        CellInd = [0] # Be careful of this
        CellDist = [9999] # Seeing this value in the Distances array will notify you that none of the cells in the 2D flow area were associated with a cell point. 
                        # This is likely because the gage coordinates are too far from a 2D area to be considered "close" to a cell point 
                        # and so face points are all that are being rendered as "close"
    return CellInd[0], CellDist[0]  # The index of the cell center which is closest to the gage coordinates

# Finds the absolute closest cell to the gage coordinates given. Returns the cell index to be used to access
# output data and the distance of that cell's center to the gage coordinates.

def _friends_bootstrap_radius(args):
    """Internal method used to compute the radius (half-side-length) for each
    ball (cube) used in :class:`RadFriends` (:class:`SupFriends`) using
    bootstrapping."""

    # Unzipping.
    points, ftype = args
    rstate = np.random

    # Resampling.
    npoints, ndim = points.shape
    idxs = rstate.randint(npoints, size=npoints)  # resample
    idx_in = np.unique(idxs)  # selected objects
    sel = np.ones(npoints, dtype='bool')
    sel[idx_in] = False
    idx_out = np.arange(npoints)[sel]  # "missing" objects
    if len(idx_out) < 2:  # edge case
        idx_out = np.append(idx_out, [0, 1])
    points_in, points_out = points[idx_in], points[idx_out]

    # Construct KDTree to enable quick nearest-neighbor lookup for
    # our resampled objects.
    kdtree = spatial.KDTree(points_in)

    if ftype == 'balls':
        # Compute distances from our "missing" points its closest neighbor
        # among the resampled points using the Euclidean norm
        # (i.e. "radius" of n-sphere).
        dists, ids = kdtree.query(points_out, k=1, eps=0, p=2)
    elif ftype == 'cubes':
        # Compute distances from our "missing" points its closest neighbor
        # among the resampled points using the Euclidean norm
        # (i.e. "half-side-length" of n-cube).
        dists, ids = kdtree.query(points_out, k=1, eps=0, p=np.inf)

    # Conservative upper-bound on radius.
    dist = max(dists)

    return dist

def compute_nn(regions, frameEnd, F=15, L=4, redirect=False):
    """
    Compute transition matrix using nearest neighbors
    Input:
        regions: (k,d): k regions with d-dim feature
        frameEnd: (n,): indices of regions where frame ends: 0-indexed, included
        F: temporal radius: nn to be searched in (2F+1) frames around curr frame
        L: nearest neighbors to be found per frame on an average
    Output: transM: (k,k)
    """
    sTime = time.time()
    M = L * (2 * F + 1)
    k, _ = regions.shape
    n = frameEnd.shape[0]
    transM = np.zeros((k, k), dtype=np.float)

    # Build 0-1 nn graph based on L2 distance using KDTree
    for i in range(n):
        # build KDTree with 2F+1 frames around frame i
        startF = max(0, i - F)
        startF = 1 + frameEnd[startF - 1] if startF > 0 else 0
        endF = frameEnd[min(n - 1, i + F)]
        tree = KDTree(regions[startF:1 + endF], leafsize=100)

        # find nn for regions in frame i
        currStartF = 1 + frameEnd[i - 1] if i > 0 else 0
        currEndF = frameEnd[i]
        distNN, nnInd = tree.query(regions[currStartF:1 + currEndF], M)
        nnInd += startF
        currInd = np.mgrid[currStartF:1 + currEndF, 0:M][0]
        transM[currInd, nnInd] = distNN
        if not redirect:
            sys.stdout.write('NearestNeighbor computation: [% 5.1f%%]\r' %
                                (100.0 * float((i + 1) / n)))
            sys.stdout.flush()

    eTime = time.time()
    print('NearestNeighbor computation finished: %.2f s' % (eTime - sTime))

    return transM

def simulate(self):
    tdoa = []   

    x, y = np.mgrid[-self.range:self.range:self.step, -self.range:self.range:self.step]
    self.points = zip(x.ravel(), y.ravel())
    #ds = SupervisedDataSet(6,2)

    for point in self.points:
        #if (point[0] >
        tdoa.append(self.tdoa(point))
        #ds.addSample(self.tdoa(point), point)

    self.tree = spatial.KDTree(tdoa)

def _prune_blobs(blobs_array, overlap):
    """Eliminated blobs with area overlap.

    Parameters
    ----------
    blobs_array : ndarray
        A 2d array with each row representing 3 (or 4) values,
        ``(row, col, sigma)`` or ``(pln, row, col, sigma)`` in 3D,
        where ``(row, col)`` (``(pln, row, col)``) are coordinates of the blob
        and ``sigma`` is the standard deviation of the Gaussian kernel which
        detected the blob.
    overlap : float
        A value between 0 and 1. If the fraction of area overlapping for 2
        blobs is greater than `overlap` the smaller blob is eliminated.

    Returns
    -------
    A : ndarray
        `array` with overlapping blobs removed.
    """

    # iterating again might eliminate more blobs, but one iteration suffices
    # for most cases
    if len(blobs_array) == 0:
        return np.array([])
    sigma = blobs_array[:, -1].max()
    distance = 2 * sigma * sqrt(blobs_array.shape[1] - 1)
    try:
        tree = spatial.cKDTree(blobs_array[:, :-1])
        pairs = np.array(list(tree.query_pairs(distance)))
    except AttributeError:  # scipy 0.9, min requirements
        tree = spatial.KDTree(blobs_array[:, :-1])
        pairs = np.array(list(tree.query_pairs(distance)))
    if len(pairs) == 0:
        return blobs_array
    else:
        for (i, j) in pairs:
            blob1, blob2 = blobs_array[i], blobs_array[j]
            if _blob_overlap(blob1, blob2) > overlap:
                if blob1[-1] > blob2[-1]:
                    blob2[-1] = 0
                else:
                    blob1[-1] = 0

    return np.array([b for b in blobs_array if b[-1] > 0])

def estimate_i_alpha(y, co):
    """ Estimate i_alpha = \int p^{\alpha}(y)dy.

    The Renyi and Tsallis entropies are simple functions of this quantity. 

    Parameters
    ----------
    y : (number of samples, dimension)-ndarray
        One row of y corresponds to one sample.
    co : cost object; details below.
    co.knn_method : str
                    kNN computation method; 'cKDTree' or 'KDTree'.
    co.k : int, >= 1
           k-nearest neighbors.
    co.eps : float, >= 0
             the k^th returned value is guaranteed to be no further than 
             (1+eps) times the distance to the real kNN.
    co.alpha : float
               alpha in the definition of i_alpha

    Returns
    -------
    i_alpha : float
              Estimated i_alpha value.

    Examples
    --------
    i_alpha = estimate_i_alpha(y,co)

    """

    num_of_samples, dim = y.shape
    distances_yy = knn_distances(y, y, True, co.knn_method, co.k, co.eps,
                                 2)[0]
    v = volume_of_the_unit_ball(dim)

    # Solution-1 (normal k):
    c = (gamma(co.k)/gamma(co.k + 1 - co.alpha))**(1 / (1 - co.alpha))

    # Solution-2 (if k is 'extreme large', say self.k=180 [ =>
    #            gamma(self.k)=inf], then use this alternative form of
    #            'c', after importing gammaln). Note: we used the
    #            'gamma(a) / gamma(b) = exp(gammaln(a) - gammaln(b))'
    #            identity.
    # c = exp(gammaln(co.k) - gammaln(co.k+1-co.alpha))**(1 / (1-co.alpha))

    s = sum(distances_yy[:, co.k-1]**(dim * (1 - co.alpha)))
    i_alpha = \
        (num_of_samples - 1) / num_of_samples * v**(1 - co.alpha) * \
        c**(1 - co.alpha) * s / (num_of_samples - 1)**co.alpha

    return i_alpha

def onclick(self,event):
        ix, iy = event.xdata, event.ydata
        #Make sure the user clicks on the map
        #If they don't, ix and iy will be None. 
        #This will result in an error in the row/col calculations
        if (ix != None and iy != None):
            if (self.plotObj.dataSet.runType != 'IDEAL'):
                lon, lat  = self.plotObj.dataSet.map[self.plotObj.currentGrid-1](
                            self.plotObj.dataSet.glons[self.plotObj.currentGrid-1],
                            self.plotObj.dataSet.glats[self.plotObj.currentGrid-1])
                if (self.plotObj.dataSet.dsetname == 'Sounding'):
                    a = np.vstack((lat,lon)).T
                    pt = [iy, ix]
                    distance, index = spatial.KDTree(a).query(pt)
                    #print(distance)
                    #print(index)
                    #print(a[index])
                    #print(iy,ix)
                    col = index
                    row = index
                    print(self.plotObj.dataSet.stat_id[index])
                else:                    
                    clon,clat = self.plotObj.dataSet.map[self.plotObj.currentGrid-1](
                                self.plotObj.dataSet.lon0[self.plotObj.currentGrid-1],
                                self.plotObj.dataSet.lat0[self.plotObj.currentGrid-1])
                    tmp = np.argsort(np.abs(lat[:,int(self.plotObj.dataSet.nx[self.plotObj.currentGrid -1]/2)] - clat))[0]
                    tmp2 = np.argsort(np.abs(lon[tmp,:] - clon))[0]
                    row = np.argsort(np.abs(lat[:,tmp2] - iy))[0]
                    col  = np.argsort(np.abs(lon[row,:] - ix))[0]
            else:                

                row = int(iy*1000/self.plotObj.dataSet.dy[self.plotObj.currentGrid-1] +
                         (self.plotObj.dataSet.ny[self.plotObj.currentGrid-1]-1)/2)
                col = int(ix*1000/self.plotObj.dataSet.dx[self.plotObj.currentGrid-1] +
                         (self.plotObj.dataSet.nx[self.plotObj.currentGrid-1]-1)/2)

            print( col, row )
            coords = [ix, iy]

            if len(coords) == 2:
                self.plotObj.col = col
                self.plotObj.row = row
                self.plotObj.replot2d = False
                self.plotObj.newvar = True
                self.drawPlot()
            else:
                self.plotObj.col = None
                self.plotObj.row = None
        else:
            self.errorClick()

    #Throw an error for a bad click

def get_closest_station(latitude, longitude, minumum_recent_data=20140000, 
                        match_max=100):
    '''Query function to find the nearest weather station to a particular 
    set of coordinates. Optionally allows for a recent date by which the 
    station is required to be still active at.

    Parameters
    ----------
    latitude : float
        Latitude to search for nearby weather stations at, [degrees]
    longitude : float
        Longitude to search for nearby weather stations at, [degrees]
    minumum_recent_data : int, optional
        Date that the weather station is required to have more recent
        weather data than; format YYYYMMDD; set this to 0 to not restrict data
        by date.
    match_max : int, optional
        The maximum number of results in the KDTree to search for before 
        applying the filtering criteria; an internal parameter which is
        increased automatically if the default value is insufficient [-]

    Returns
    -------
    station : IntegratedSurfaceDatabaseStation
        Instance of IntegratedSurfaceDatabaseStation which was nearest
        to the requested coordinates and with sufficiently recent data
        available [-]

    Notes
    -----
    Searching for 100 stations is a reasonable choice as it takes, ~70 
    microseconds vs 50 microsecond to find only 1 station. The search does get 
    slower as more points are requested. Bad data is returned from a KDTree
    search if more points are requested than are available.

    Examples
    --------
    >>> get_closest_station(51.02532675, -114.049868485806, 20150000)
    <Weather station registered in the Integrated Surface Database, name CALGARY INTL CS, country CA, USAF 713930.0, WBAN None, coords (51.1, -114.0) Weather data from 2004 to 2017>
    '''
    # Both station strings may be important
    # Searching for 100 stations is fine, 70 microseconds vs 50 microsecond for 1
    # but there's little point for more points, it gets slower.
    # bad data is returned if k > station_count
    distances, indexes = kd_tree.query([latitude, longitude], k=min(match_max, station_count)) 
    #
    for i in indexes:
        latlon = _latlongs[i]
        enddate = stations[i].END
        # Iterate for all indexes until one is found whose date is current
        if enddate > minumum_recent_data:
            return stations[i]
    if match_max < station_count:
        return get_closest_station(latitude, longitude, minumum_recent_data=minumum_recent_data, match_max=match_max*10)
    raise Exception('Could not find a station with more recent data than '
                    'specified near the specified coordinates.')


# This should be agressively cached