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
python类cKDTree()的实例源码
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 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, 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)))
