我们从Python开源项目中，提取了以下50个代码示例，用于说明如何使用scipy.spatial.cKDTree()。
def eliminate_overlapping_locations(f, separation): """ Makes sure that no position is within `separation` from each other, by deleting one of the that are to close to each other. """ separation = validate_tuple(separation, f.shape[1]) assert np.greater(separation, 0).all() # Rescale positions, so that pairs are identified below a distance of 1. f = f / separation while True: duplicates = cKDTree(f, 30).query_pairs(1) if len(duplicates) == 0: break to_drop = [] for pair in duplicates: to_drop.append(pair[1]) f = np.delete(f, to_drop, 0) return f * separation
def darts(n, xx, yy, rr, dst): """ get at most n random, uniformly distributed, points in a circle. centered at (xx,yy), with radius rr. points are no closer to each other than dst. """ ## remove new nodes that are too close to other ## new nodes visited = set() dartsxy = random_points_in_circle(n, xx, yy, rr) tree = kdt(dartsxy) near = tree.query_ball_point(dartsxy, dst) jj = [] for j,n in enumerate(near): if len(visited.intersection(n))<1: jj.append(j) visited.add(j) res = dartsxy[jj,:] return res
def darts_rect(n, xx, yy, w=1, h=1, dst=0): ## remove new nodes that are too close to other ## new nodes visited = set() dartsxy = random_points_in_rectangle(n, xx, yy, w, h) tree = kdt(dartsxy) near = tree.query_ball_point(dartsxy, dst) jj = [] for j,n in enumerate(near): if len(visited.intersection(n))<1: jj.append(j) visited.add(j) res = dartsxy[jj,:] return res
def get_neighbors(self, pos, delta): def _pos2coor(pos): a, b, c = np.array(self.m) x, y, z = pos coor = a*x + b*y + c*z # an array return tuple(coor) def p_gen(): for z in range(self.depth): for x in range(self.width): for y in range(self.length): yield(x, y, z) point = _pos2coor(pos) # w = self.width # l = self.length coor_map = {p: _pos2coor(p) for p in p_gen()} del coor_map[pos] points = list(coor_map.values()) points_tree = cKDTree(points) ind = points_tree.query_ball_point(point, delta) neighbors = itemgetter(*ind)(list(coor_map.keys())) return set(neighbors)
def get_neighbors(self, pos, delta): def _pos2coor(pos): a, b = np.array(self.m) x, y = pos coor = a*x + b*y # an array return tuple(coor) def p_gen(): for x in range(self.width): for y in range(self.length): yield(x, y) point = _pos2coor(pos) # w = self.width # l = self.length coor_map = {p : _pos2coor(p) for p in p_gen()} del coor_map[pos] points = list(coor_map.values()) points_tree = cKDTree(points) ind = points_tree.query_ball_point(point, delta) neighbors = itemgetter(*ind)(list(coor_map.keys())) return set(neighbors)
def _pquery(scheduler, data, ndata, ndim, leafsize, x, nx, d, i, k, eps, p, dub, ierr): try: _data = shmem_as_nparray(data).reshape((ndata, ndim)) _x = shmem_as_nparray(x).reshape((nx, ndim)) _d = shmem_as_nparray(d).reshape((nx, k)) _i = shmem_as_nparray(i).reshape((nx, k)) kdtree = cKDTree(_data, leafsize=leafsize) for s in scheduler: d_out, i_out = kdtree.query(_x[s, :], k=k, eps=eps, p=p, distance_upper_bound=dub) m_d = d_out.shape[0] m_i = i_out.shape[0] _d[s, :], _i[s, :] = d_out.reshape(m_d, 1), i_out.reshape(m_i, 1) except: ierr.value += 1
def __init__(self, mode=1, precision_mode=2): self.mode = mode if precision_mode == 0: loc_path = RG_FILE_1000 elif precision_mode == 1: loc_path = RG_FILE_5000 else: loc_path = RG_FILE_15000 if not os.path.exists(loc_path): loc_path = rel_path(loc_path) # coordinates, locations = self.extract(path) coordinates, self.locations = self.extract(loc_path) if mode == 1: # Single-process self.tree = KDTree(coordinates) else: # Multi-process self.tree = KDTree_MP.cKDTree_MP(coordinates)
def mi(x, y, k=3, base=2): """ Mutual information of x and y x, y should be a list of vectors, e.g. x = [[1.3], [3.7], [5.1], [2.4]] if x is a one-dimensional scalar and we have four samples """ assert len(x) == len(y), "Lists should have same length" assert k <= len(x) - 1, "Set k smaller than num. samples - 1" intens = 1e-10 # small noise to break degeneracy, see doc. x = [list(p + intens * nr.rand(len(x[0]))) for p in x] y = [list(p + intens * nr.rand(len(y[0]))) for p in y] points = zip2(x, y) # Find nearest neighbors in joint space, p=inf means max-norm tree = ss.cKDTree(points) dvec = [tree.query(point, k + 1, p=float('inf'))[0][k] for point in points] a, b, c, d = avgdigamma(x, dvec), avgdigamma(y, dvec), digamma(k), digamma(len(x)) return (-a - b + c + d) / log(base)
def cmi(x, y, z, k=3, base=2): """ Mutual information of x and y, conditioned on z x, y, z should be a list of vectors, e.g. x = [[1.3], [3.7], [5.1], [2.4]] if x is a one-dimensional scalar and we have four samples """ assert len(x) == len(y), "Lists should have same length" assert k <= len(x) - 1, "Set k smaller than num. samples - 1" intens = 1e-10 # small noise to break degeneracy, see doc. x = [list(p + intens * nr.rand(len(x[0]))) for p in x] y = [list(p + intens * nr.rand(len(y[0]))) for p in y] z = [list(p + intens * nr.rand(len(z[0]))) for p in z] points = zip2(x, y, z) # Find nearest neighbors in joint space, p=inf means max-norm tree = ss.cKDTree(points) dvec = [tree.query(point, k + 1, p=float('inf'))[0][k] for point in points] a, b, c, d = avgdigamma(zip2(x, z), dvec), avgdigamma(zip2(y, z), dvec), avgdigamma(z, dvec), digamma(k) return (-a - b + c + d) / log(base)
def kldiv(x, xp, k=3, base=2): """ KL Divergence between p and q for x~p(x), xp~q(x) x, xp should be a list of vectors, e.g. x = [[1.3], [3.7], [5.1], [2.4]] if x is a one-dimensional scalar and we have four samples """ assert k <= len(x) - 1, "Set k smaller than num. samples - 1" assert k <= len(xp) - 1, "Set k smaller than num. samples - 1" assert len(x[0]) == len(xp[0]), "Two distributions must have same dim." d = len(x[0]) n = len(x) m = len(xp) const = log(m) - log(n - 1) tree = ss.cKDTree(x) treep = ss.cKDTree(xp) nn = [tree.query(point, k + 1, p=float('inf'))[0][k] for point in x] nnp = [treep.query(point, k, p=float('inf'))[0][k - 1] for point in x] return (const + d * np.mean(map(log, nnp)) - d * np.mean(map(log, nn))) / log(base) # DISCRETE ESTIMATORS
def load_data(fname='data/gps_data/gps_points.csv'): """ Given a file that contains gps points, this method creates different data structures :param fname: the name of the input file, as generated by QMIC :return: data_points (list of gps positions with their metadata), raw_points (coordinates only), points_tree is the KDTree structure to enable searching the points space """ data_points = list() raw_points = list() with open(fname, 'r') as f: f.readline() for line in f: if len(line) < 10: continue vehicule_id, timestamp, lat, lon, speed, angle = line.split(',') pt = GpsPoint(vehicule_id=vehicule_id, timestamp=timestamp, lat=lat, lon=lon, speed=speed,angle=angle) data_points.append(pt) raw_points.append(pt.get_coordinates()) points_tree = cKDTree(raw_points) return np.array(data_points), np.array(raw_points), points_tree
def connect_graphs(self, sets_orig, edges_orig): ''' Returns the edges needed to connect unconnected graphs (sets of nodes) given a set of sets of nodes, select the master_graph (the biggest) one and search the shortest edges to connect the other sets of nodes ''' master_graph = max(sets_orig, key=len) sets = sets_orig.copy() edges = np.array([], dtype=object) sets.remove(master_graph) master_tree = cKDTree(self.nodes[list(master_graph)]) for s in sets: x = np.array(list(s)) nearests = np.array([master_tree.query(self.nodes[v]) for v in x]) tails = nearests[ nearests[:, 0].argsort()][:, 1][:self.max_neighbours] heads = x[nearests[:, 0].argsort()][:self.max_neighbours] for head, tail in zip(heads, tails): edges = np.append(edges, None) edges[-1] = (head, tail) edges = np.append(edges, None) edges[-1] = (tail, head) return edges
def tsp(self, start): res = [] nodes = np.array([i for i in range(self.n)]) visited = np.empty(self.n, dtype=np.bool) visited.fill(False) visited[start] = True node = start while False in visited: tree = cKDTree(self.nodes[~visited]) nearest = tree.query(self.nodes[node], k=1)[1] t = nodes[~visited][nearest] res.append([node, t]) visited[t] = True node = t return res + [node, start]
def KL_entropy(x,k=5): ''' Estimate the entropy H(X) from samples {x_i}_{i=1}^N Using Kozachenko-Leonenko (KL) estimator Input: x: 2D list of size N*d_x k: k-nearest neighbor parameter Output: one number of H(X) ''' assert k <= len(x)-1, "Set k smaller than num. samples - 1" N = len(x) d = len(x[0]) tree = ss.cKDTree(x) knn_dis = [tree.query(point,k+1,p=2)[0][k] for point in x] ans = -log(k) + log(N) + vd(d,2) return ans + d*np.mean(map(log,knn_dis))
def _3LNN_1_KSG_mi(data,split,k=5,tr=30): ''' Estimate the mutual information I(X;Y) from samples {x_i,y_i}_{i=1}^N Using I(X;Y) = H_{LNN}(X) + H_{LNN}(Y) - H_{LNN}(X;Y) with "KSG trick" where H_{LNN} is the LNN entropy estimator with order 1 Input: data: 2D list of size N*(d_x + d_y) split: should be d_x, splitting the data into two parts, X and Y k: k-nearest neighbor parameter tr: number of sample used for computation Output: one number of I(X;Y) ''' assert split >=1, "x must have at least one dimension" assert split <= len(data[0]) - 1, "y must have at least one dimension" x = data[:,:split] y = data[:,split:] tree_xy = ss.cKDTree(data) knn_dis = [tree_xy.query(point,k+1,p=2)[0][k] for point in data] return LNN_1_entropy(x,k,tr,bw=knn_dis) + LNN_1_entropy(y,k,tr,bw=knn_dis) - LNN_1_entropy(data,k,tr,bw=knn_dis)
def mi(x, y, k=3, base=2): """Mutual information of x and y. x,y should be a list of vectors, e.g. x = [[1.3], [3.7], [5.1], [2.4]] if x is a one-dimensional scalar and we have four samples. """ assert len(x)==len(y), 'Lists should have same length.' assert k <= len(x) - 1, 'Set k smaller than num samples - 1.' intens = 1e-10 # Small noise to break degeneracy, see doc. x = [list(p + intens*nr.rand(len(x[0]))) for p in x] y = [list(p + intens*nr.rand(len(y[0]))) for p in y] points = zip2(x,y) # Find nearest neighbors in joint space, p=inf means max-norm. tree = ss.cKDTree(points) dvec = [tree.query(point, k+1, p=float('inf'))[0][k] for point in points] a = avgdigamma(x,dvec) b = avgdigamma(y,dvec) c = digamma(k) d = digamma(len(x)) return (-a-b+c+d) / log(base)
def cmi(x, y, z, k=3, base=2): """Mutual information of x and y, conditioned on z x,y,z should be a list of vectors, e.g. x = [[1.3], [3.7], [5.1], [2.4]] if x is a one-dimensional scalar and we have four samples """ assert len(x)==len(y), 'Lists should have same length.' assert k <= len(x) - 1, 'Set k smaller than num samples - 1.' intens = 1e-10 # Small noise to break degeneracy, see doc. x = [list(p + intens*nr.rand(len(x[0]))) for p in x] y = [list(p + intens*nr.rand(len(y[0]))) for p in y] z = [list(p + intens*nr.rand(len(z[0]))) for p in z] points = zip2(x,y,z) # Find nearest neighbors in joint space, p=inf means max-norm. tree = ss.cKDTree(points) dvec = [tree.query(point, k+1, p=float('inf'))[0][k] for point in points] a = avgdigamma(zip2(x,z), dvec) b = avgdigamma(zip2(y,z), dvec) c = avgdigamma(z,dvec) d = digamma(k) return (-a-b+c+d) / log(base)
def kldiv(x, xp, k=3, base=2): """KL Divergence between p and q for x~p(x),xp~q(x). x, xp should be a list of vectors, e.g. x = [[1.3],[3.7],[5.1],[2.4]] if x is a one-dimensional scalar and we have four samples """ assert k <= len(x) - 1, 'Set k smaller than num samples - 1.' assert k <= len(xp) - 1, 'Set k smaller than num samples - 1.' assert len(x[0]) == len(xp[0]), 'Two distributions must have same dim.' d = len(x[0]) n = len(x) m = len(xp) const = log(m) - log(n-1) tree = ss.cKDTree(x) treep = ss.cKDTree(xp) nn = [tree.query(point, k+1, p=float('inf'))[0][k] for point in x] nnp = [treep.query(point, k, p=float('inf'))[0][k-1] for point in x] return (const + d * np.mean(map(log, nnp)) \ - d * np.mean(map(log, nn))) / log(base) # DISCRETE ESTIMATORS
def calc_vbar(self): # Calculates the average velocities from neighboring birds depending on max_dist and max_num. my_tree = spatial.cKDTree(self.current_points) dist, indexes = my_tree.query(self.current_points, k=self._max_num) ri = np.zeros((len(self.current_points), self._max_num), dtype=int) rk = np.zeros_like(ri, dtype=int) good_inds = (dist < self._max_dist) ri[good_inds] = indexes[good_inds] rk[good_inds] = 1 # I should get the angle and average ori = np.arctan2(self.current_v[:, 1], self.current_v[:, 0]) mean_n = [] for i in range(len(ri)): nei = ri[i][np.where(rk[i] == 1)[0]] mm = (np.arctan2(np.sum(np.sin(ori[nei])), np.sum(np.cos(ori[nei])))) % (2 * np.pi) mean_n.append(mm) vbar = np.array(zip(np.cos(mean_n), np.sin(mean_n))) return vbar, [np.array(ri), np.array(rk)]
def genDSM_building(BD_footprint,PC_data,h,w,n,ground_height,affine_par,upper_bd_par): [A,D,B,E,C,F] = np.loadtxt(affine_par) BD_footprint = cv2.imread(BD_footprint,0) point_cloud = pd.read_csv(PC_data,names=['X', 'Y', 'Z', 'R','G','B'],delim_whitespace=True) X = point_cloud.X.copy() Y = point_cloud.Y.copy() # Inverse the Y-axis Z = point_cloud.Z.copy() del point_cloud Dsm_arr = np.zeros((h, w)) # cKDtree tree = spatial.cKDTree(list(zip(Y, X))) #Affine Trans yv, xv = np.where(BD_footprint>127) XA = A*xv+B*yv+C YA = D*xv+E*yv+F pts = np.dstack((YA,XA)) dis, loc = tree.query(pts, k=n ,distance_upper_bound=upper_bd_par) #building use max Dsm_arr[yv,xv] = Z[loc[:,:,:].ravel()].values.reshape(-1, n).max(axis=1) Dsm_arr=np.float32(Dsm_arr) return Dsm_arr
def voronoi_tessellation(x, y, xnode, ynode, scale): """ Computes (Weighted) Voronoi Tessellation of the pixels grid """ if scale[0] == 1: # non-weighted VT tree = cKDTree(np.column_stack([xnode, ynode])) classe = tree.query(np.column_stack([x, y]))[1] else: if x.size < 1e4: classe = np.argmin(((x[:, None] - xnode)**2 + (y[:, None] - ynode)**2)/scale**2, axis=1) else: # use for loop to reduce memory usage classe = np.zeros(x.size, dtype=int) for j, (xj, yj) in enumerate(zip(x, y)): classe[j] = np.argmin(((xj - xnode)**2 + (yj - ynode)**2)/scale**2) return classe #----------------------------------------------------------------------
def remove_close(points, face_index, radius): ''' Given an (n, m) set of points where n=(2|3) return a list of points where no point is closer than radius ''' from scipy.spatial import cKDTree as KDTree tree = KDTree(points) consumed = np.zeros(len(points), dtype=np.bool) unique = np.zeros(len(points), dtype=np.bool) for i in range(len(points)): if consumed[i]: continue neighbors = tree.query_ball_point(points[i], r=radius) consumed[neighbors] = True unique[i] = True return points[unique], face_index[unique]
def remove_close(points, radius): ''' Given an (n, m) set of points where n=(2|3) return a list of points where no point is closer than radius ''' from scipy.spatial import cKDTree as KDTree tree = KDTree(points) consumed = np.zeros(len(points), dtype=np.bool) unique = np.zeros(len(points), dtype=np.bool) for i in range(len(points)): if consumed[i]: continue neighbors = tree.query_ball_point(points[i], r=radius) consumed[neighbors] = True unique[i] = True return points[unique]
def remove_close_set(points_fixed, points_reduce, radius): ''' Given two sets of points and a radius, return a set of points that is the subset of points_reduce where no point is within radius of any point in points_fixed ''' from scipy.spatial import cKDTree as KDTree tree_fixed = KDTree(points_fixed) tree_reduce = KDTree(points_reduce) reduce_duplicates = tree_fixed.query_ball_tree(tree_reduce, r = radius) reduce_duplicates = np.unique(np.hstack(reduce_duplicates).astype(int)) reduce_mask = np.ones(len(points_reduce), dtype=np.bool) reduce_mask[reduce_duplicates] = False points_clean = points_reduce[reduce_mask] return points_clean
def _get_neighbours_ix(obs_coords, raw_coords, nnear): """ Returns <nnear> neighbour indices per <obs_coords> coordinate pair Parameters ---------- obs_coords : array of float of shape (num_points,ndim) in the neighbourhood of these coordinate pairs we look for neighbours raw_coords : array of float of shape (num_points,ndim) from these coordinate pairs the neighbours are selected nnear : integer number of neighbours to be selected per coordinate pair of obs_coords """ # plant a tree tree = cKDTree(raw_coords) # return nearest neighbour indices return tree.query(obs_coords, k=nnear)[1]
def _get_neighbours(obs_coords, raw_coords, raw, nnear): """ Returns <nnear> neighbour values per <obs_coords> coordinate pair Parameters ---------- obs_coords : array of float of shape (num_points,ndim) in the neighbourhood of these coordinate pairs we look for neighbours raw_coords : array of float of shape (num_points,ndim) from these coordinate pairs the neighbours are selected raw : array of float of shape (num_points,...) this is the data corresponding to the coordinate pairs raw_coords nnear : integer number of neighbours to be selected per coordinate pair of obs_coords """ # plant a tree tree = cKDTree(raw_coords) # retrieve nearest neighbour indices ix = tree.query(obs_coords, k=nnear)[1] # return the values of the nearest neighbours return raw[ix]
def test_point_correctness(): import itertools stencil = [-1, 0, 1] ndim = 3 n = 2000 stencil = itertools.product(*[stencil]*ndim) stencil = np.array(list(stencil)).astype(np.int32) points = (np.random.rand(n, ndim) * [1, 2, 3]).astype(np.float32) scale = 0.1 spec = GridSpec(points, float(scale)) offsets = spec.stencil(stencil).astype(np.int32) grid = PointGrid(spec, points, offsets) pairs = grid.pairs() from scipy.spatial import cKDTree tree = cKDTree(points) tree_pairs = tree.query_pairs(scale, output_type='ndarray') print(tree_pairs) print(pairs) assert np.alltrue(npi.unique(tree_pairs) == npi.unique(np.sort(pairs, axis=1)))
def hp_kd_tree(nside=set_default_nside(), leafsize=100): """ Generate a KD-tree of healpixel locations Parameters ---------- nside : int A valid healpix nside leafsize : int (100) Leafsize of the kdtree Returns ------- tree : scipy kdtree """ hpid = np.arange(hp.nside2npix(nside)) ra, dec = _hpid2RaDec(nside, hpid) x, y, z = treexyz(ra, dec) tree = kdtree(list(zip(x, y, z)), leafsize=leafsize, balanced_tree=False, compact_nodes=False) return tree
def _set_altaz_fields(self, leafsize=10): """ Have a fixed grid of alt,az pointings to use """ tmp = read_fields() names = ['alt', 'az'] types = [float, float] self.fields = np.zeros(tmp.size, dtype=list(zip(names, types))) self.fields['alt'] = tmp['dec'] self.fields['az'] = tmp['RA'] x, y, z = treexyz(self.fields['az'], self.fields['alt']) self.field_tree = kdtree(list(zip(x, y, z)), leafsize=leafsize, balanced_tree=False, compact_nodes=False) hpids = np.arange(hp.nside2npix(self.nside)) self.reward_ra, self.reward_dec = _hpid2RaDec(self.nside, hpids)
def _spin_fields(self, lon=None, lat=None): """Spin the field tesselation """ if lon is None: lon = np.random.rand()*np.pi*2 if lat is None: lat = np.random.rand()*np.pi*2 # rotate longitude ra = (self.fields['RA'] + lon) % (2.*np.pi) dec = self.fields['dec'] + 0 # Now to rotate ra and dec about the x-axis x, y, z = thetaphi2xyz(self.fields['RA'], self.fields['dec']+np.pi/2.) xp, yp, zp = rotx(lat, x, y, z) theta, phi = xyz2thetaphi(xp, yp, zp) dec = phi - np.pi/2 ra = theta + np.pi self.fields['RA'] = ra self.fields['dec'] = dec # Rebuild the kdtree with the new positions # XXX-may be doing some ra,dec to conversions xyz more than needed. self._hp2fieldsetup(ra, dec)
def fit(self, X, w): if len(X) == 0: raise NotEnoughParticles("Fitting not possible.") self.X_arr = X.as_matrix() ctree = cKDTree(X) _, indices = ctree.query(X, k=min(self.k + 1, X.shape[0])) covs, inv_covs, dets = list(zip(*[self._cov_and_inv(n, indices) for n in range(X.shape[0])])) self.covs = sp.array(covs) self.inv_covs = sp.array(inv_covs) self.determinants = sp.array(dets) self.normalization = sp.sqrt( (2 * sp.pi) ** self.X_arr.shape[1] * self.determinants) if not sp.isreal(self.normalization).all(): raise Exception("Normalization not real") self.normalization = sp.real(self.normalization)
def __init__(self, init_num, size, stp): self.num = init_num self.size = size self.stp = stp self.one = 1.0/size self.rad = 0.04 edge = 0.2 self.xy = edge + random((init_num, 2))*(1.0-2*edge) self.v = zeros((init_num, 2), 'float') self.dx = zeros((init_num, 2), 'float') self.random_velocity(self.stp) self.tree = kdtree(self.xy) self.slowdown = 0.99999
def get_cKDTree(): try: from scipy.spatial import cKDTree except ImportError: cKDTree = None return cKDTree # getheader gives a 87a header and a color palette (two elements in a list). # getdata()[0] gives the Image Descriptor up to (including) "LZW min code size". # getdatas()[1:] is the image data itself in chuncks of 256 bytes (well # technically the first byte says how many bytes follow, after which that # amount (max 255) follows).
def quantize_with_scipy(self, image): w,h = image.size px = np.asarray(image).copy() px2 = px[:,:,:3].reshape((w*h,3)) cKDTree = get_cKDTree() kdtree = cKDTree(self.colormap[:,:3],leafsize=10) result = kdtree.query(px2) colorindex = result[1] print("Distance: %1.2f" % (result[0].sum()/(w*h)) ) px2[:] = self.colormap[colorindex,:3] return Image.fromarray(px).convert("RGB").quantize(palette=self.paletteImage())
def compute_index(codebook): return cKDTree(codebook)
def quantize_with_scipy(self, image): w, h = image.size px = np.asarray(image).copy() px2 = px[:, :, :3].reshape((w * h, 3)) cKDTree = get_cKDTree() kdtree = cKDTree(self.colormap[:, :3], leafsize=10) result = kdtree.query(px2) colorindex = result[1] print("Distance: %1.2f" % (result[0].sum() / (w * h))) px2[:] = self.colormap[colorindex, :3] return Image.fromarray(px).convert("RGB").quantize(palette=self.paletteImage())
def invalidate_distances(self): self._tree = cKDTree(self.params.X) self._distances, self._neighbors = None, None
def where_close(pos, separation, intensity=None): """ Returns indices of features that are closer than separation from other features. When intensity is given, the one with the lowest intensity is returned: else the most topleft is returned (to avoid randomness) To be implemented in trackpy v0.4""" if len(pos) == 0: return [] separation = validate_tuple(separation, pos.shape[1]) if any([s == 0 for s in separation]): return [] # Rescale positions, so that pairs are identified below a distance # of 1. pos_rescaled = pos / separation duplicates = cKDTree(pos_rescaled, 30).query_pairs(1 - 1e-7) if len(duplicates) == 0: return [] index_0 = np.fromiter((x[0] for x in duplicates), dtype=int) index_1 = np.fromiter((x[1] for x in duplicates), dtype=int) if intensity is None: to_drop = np.where(np.sum(pos_rescaled[index_0], 1) > np.sum(pos_rescaled[index_1], 1), index_1, index_0) else: intensity_0 = intensity[index_0] intensity_1 = intensity[index_1] to_drop = np.where(intensity_0 > intensity_1, index_1, index_0) edge_cases = intensity_0 == intensity_1 if np.any(edge_cases): index_0 = index_0[edge_cases] index_1 = index_1[edge_cases] to_drop[edge_cases] = np.where(np.sum(pos_rescaled[index_0], 1) > np.sum(pos_rescaled[index_1], 1), index_1, index_0) return np.unique(to_drop)
def sort_positions(actual, expected): assert_equal(len(actual), len(expected)) tree = cKDTree(actual) devs, argsort = tree.query([expected]) return devs, actual[argsort][0]
def __init__(self, X=None, z=None, leafsize=10): if not X is None: self.tree = cKDTree(X, leafsize=leafsize ) if not z is None: self.z = z
def _make_kdtree(self, x): self.kdtree = cKDTree(mpiops.random_full_points(x, Napprox=self.nodes)) if not np.isfinite(self.kdtree.query(x, k=self.k)[0]).all(): log.warning('Kdtree computation encountered problem. ' 'Not enough neighbors available to compute ' 'kdtree. Printing kdtree for debugging purpose') raise ValueError('Computed kdtree is not fully populated.' 'Not enough valid neighbours available.')
def __init__(self, train, weights=None, truncation=3.0, nmin=4, factor=1.0): """ Parameters ---------- train : np.ndarray The training data set. Should be a 1D array of samples or a 2D array of shape (n_samples, n_dim). weights : np.ndarray, optional An array of weights. If not specified, equal weights are assumed. truncation : float, optional The maximum deviation (in sigma) to use points in the KDE nmin : int, optional The minimum number of points required to estimate the density factor : float, optional Send bandwidth to this factor of the data estimate """ self.truncation = truncation self.nmin = nmin self.train = train if len(train.shape) == 1: train = np.atleast_2d(train).T self.num_points, self.num_dim = train.shape if weights is None: weights = np.ones(self.num_points) self.weights = weights self.mean = np.average(train, weights=weights, axis=0) dx = train - self.mean cov = np.atleast_2d(np.cov(dx.T, aweights=weights)) self.A = np.linalg.cholesky(np.linalg.inv(cov)) # The sphere-ifying transform self.d = np.dot(dx, self.A) # Sphere-ified data self.tree = spatial.cKDTree(self.d) # kD tree of data self.sigma = 2.0 * factor * np.power(self.num_points, -1. / (4 + self.num_dim)) # Starting sigma (bw) of Gauss self.sigma_fact = -0.5 / (self.sigma * self.sigma) # Cant get normed probs to work atm, turning off for now as I don't need normed pdfs for contours # self.norm = np.product(np.diagonal(self.A)) * (2 * np.pi) ** (-0.5 * self.num_dim) # prob norm # self.scaling = np.power(self.norm * self.sigma, -self.num_dim)
def entropy(x, k=3, base=2): """ The classic K-L k-nearest neighbor continuous entropy estimator x should be a list of vectors, e.g. x = [[1.3], [3.7], [5.1], [2.4]] if x is a one-dimensional scalar and we have four samples """ assert k <= len(x) - 1, "Set k smaller than num. samples - 1" d = len(x[0]) N = len(x) intens = 1e-10 # small noise to break degeneracy, see doc. x = [list(p + intens * nr.rand(len(x[0]))) for p in x] tree = ss.cKDTree(x) nn = [tree.query(point, k + 1, p=float('inf'))[0][k] for point in x] const = digamma(N) - digamma(k) + d * log(2) return (const + d * np.mean(map(log, nn))) / log(base)
def avgdigamma(points, dvec): # This part finds number of neighbors in some radius in the marginal space # returns expectation value of <psi(nx)> N = len(points) tree = ss.cKDTree(points) avg = 0. for i in range(N): dist = dvec[i] # subtlety, we don't include the boundary point, # but we are implicitly adding 1 to kraskov def bc center point is included num_points = len(tree.query_ball_point(points[i], dist - 1e-15, p=float('inf'))) avg += digamma(num_points) / N return avg
def maskMeshByDistance(x,y,d,grid): """Filters all (x,y) coordinates that are more than d in meshgrid given some actual coordinates (x,y). Args: x (numpy.ndarray): x-coordinates. y (numpy.ndarray): y-coordinates. d (float): Maximum distance. grid (numpy.ndarray): Numpy meshgrid. Returns: idxs (list): List of booleans. """ #Convert x/y into useful array xy=np.vstack((x,y)) #Compute distances to nearest neighbors tree = cKDTree(xy.T) xi = _ndim_coords_from_arrays(tuple(grid)) dists, indexes = tree.query(xi) #Get indices of distances that are further apart than d idxs = (dists > d) return idxs,
