将 numpy 数组的组名映射到索引的最快方法是什么？
What is the fastest way to map group names of numpy array to indices?
我正在使用激光雷达的 3D 点云。这些点由 numpy 数组给出，如下所示：
points = np.array([[61651921, 416326074, 39805], [61605255, 416360555, 41124], [61664810, 416313743, 39900], [61664837, 416313749, 39910], [61674456, 416316663, 39503], [61651933, 416326074, 39802], [61679969, 416318049, 39500], [61674494, 416316677, 39508], [61651908, 416326079, 39800], [61651908, 416326087, 39802], [61664845, 416313738, 39913], [61674480, 416316668, 39503], [61679996, 416318047, 39510], [61605290, 416360572, 41118], [61605270, 416360565, 41122], [61683939, 416313004, 41052], [61683936, 416313033, 41060], [61679976, 416318044, 39509], [61605279, 416360555, 41109], [61664837, 416313739, 39915], [61674487, 416316666, 39505], [61679961, 416318035, 39503], [61683943, 416313004, 41054], [61683930, 416313042, 41059]])
我想将我的数据分组为大小为 50*50*50 的多维数据集，以便每个多维数据集都保留我的 points 它包含的一些可哈希索引和 numpy 索引。为了进行拆分，我将 cubes = points \ 50 分配给输出到：
cubes = np.array([[1233038, 8326521, 796], [1232105, 8327211, 822], [1233296, 8326274, 798], [1233296, 8326274, 798], [1233489, 8326333, 790], [1233038, 8326521, 796], [1233599, 8326360, 790], [1233489, 8326333, 790], [1233038, 8326521, 796], [1233038, 8326521, 796], [1233296, 8326274, 798], [1233489, 8326333, 790], [1233599, 8326360, 790], [1232105, 8327211, 822], [1232105, 8327211, 822], [1233678, 8326260, 821], [1233678, 8326260, 821], [1233599, 8326360, 790], [1232105, 8327211, 822], [1233296, 8326274, 798], [1233489, 8326333, 790], [1233599, 8326360, 790], [1233678, 8326260, 821], [1233678, 8326260, 821]])
我想要的输出如下所示：
{(1232105, 8327211, 822): [1, 13, 14, 18]),
(1233038, 8326521, 796): [0, 5, 8, 9],
(1233296, 8326274, 798): [2, 3, 10, 19],
(1233489, 8326333, 790): [4, 7, 11, 20],
(1233599, 8326360, 790): [6, 12, 17, 21],
(1233678, 8326260, 821): [15, 16, 22, 23]}
我的真实点云包含多达数亿个 3D 点。进行这种分组最快的方法是什么？
我已经尝试了大多数不同的解决方案。下面是时间消耗的比较，假设点的大小约为 2000 万，不同立方体的大小约为 100 万：
Pandas [tuple(elem) -> np.array(dtype=int64)]
import pandas as pd
print(pd.DataFrame(cubes).groupby([0,1,2]).indices)
#takes 9sec
默认 [elem.tobytes() 或元组 -> 列表]
#thanks @abc:
result = defaultdict(list)
for idx, elem in enumerate(cubes):
result[elem.tobytes()].append(idx) # takes 20.5sec
# result[elem[0], elem[1], elem[2]].append(idx) #takes 27sec
# result[tuple(elem)].append(idx) # takes 50sec
numpy_indexed [整数 -> np.array]
# thanks @Eelco Hoogendoorn for his library
values = npi.group_by(cubes).split(np.arange(len(cubes)))
result = dict(enumerate(values))
# takes 9.8sec
Pandas + 降维[int -> np.array(dtype=int64)]
# thanks @Divakar for showing numexpr library:
import numexpr as ne
def dimensionality_reduction(cubes):
#cubes = cubes - np.min(cubes, axis=0) #in case some coords are negative
cubes = cubes.astype(np.int64)
s0, s1 = cubes[:,0].max()+1, cubes[:,1].max()+1
d = {'s0':s0,'s1':s1,'c0':cubes[:,0],'c1':cubes[:,1],'c2':cubes[:,2]}
c1D = ne.evaluate('c0+c1*s0+c2*s0*s1',d)
return c1D
cubes = dimensionality_reduction(cubes)
result = pd.DataFrame(cubes).groupby([0]).indices
# takes 2.5 seconds
可以下载 cubes.npz 文件 here 并使用命令
cubes = np.load('cubes.npz')['array']
查看表演时间。
您可能只是迭代并将每个元素的索引添加到相应的列表中。
from collections import defaultdict
res = defaultdict(list)
for idx, elem in enumerate(cubes):
#res[tuple(elem)].append(idx)
res[elem.tobytes()].append(idx)
可以通过使用 tobytes() 而不是将键转换为元组来进一步改进运行时间。
每组索引数量不变
方法 #1
我们可以执行 dimensionality-reduction 将 cubes 缩减为一维数组。这是基于将给定立方体数据映射到 n 维网格以计算线性索引等价物，详细讨论 。然后，基于这些线性索引的唯一性，我们可以分离唯一组及其对应的索引。因此，按照这些策略，我们将有一个解决方案，就像这样 -
N = 4 # number of indices per group
c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)
sidx = c1D.argsort()
indices = sidx.reshape(-1,N)
unq_groups = cubes[indices[:,0]]
# If you need in a zipped dictionary format
out = dict(zip(map(tuple,unq_groups), indices))
备选方案 #1 : 如果 cubes 中的整数值太大，我们可能需要 dimensionality-reduction 使得维度更短选择范围作为主轴。因此，对于这些情况，我们可以修改缩减步骤以获得 c1D，就像这样 -
s1,s2 = cubes[:,:2].max(0)+1
s = np.r_[s2,1,s1*s2]
c1D = cubes.dot(s)
方法 #2
接下来，我们可以使用 Cython-powered kd-tree for quick nearest-neighbor lookup 来获取最近的相邻索引，从而像这样解决我们的问题 -
from scipy.spatial import cKDTree
idx = cKDTree(cubes).query(cubes, k=N)[1] # N = 4 as discussed earlier
I = idx[:,0].argsort().reshape(-1,N)[:,0]
unq_groups,indices = cubes[I],idx[I]
一般情况：每组索引数量可变
我们将通过一些拆分来扩展基于 argsort 的方法以获得我们想要的输出，就像这样 -
c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)
sidx = c1D.argsort()
c1Ds = c1D[sidx]
split_idx = np.flatnonzero(np.r_[True,c1Ds[:-1]!=c1Ds[1:],True])
grps = cubes[sidx[split_idx[:-1]]]
indices = [sidx[i:j] for (i,j) in zip(split_idx[:-1],split_idx[1:])]
# If needed as dict o/p
out = dict(zip(map(tuple,grps), indices))
使用 cubes 组的一维版本作为键
我们将扩展前面列出的方法，使用 cubes 的组作为键，以简化创建字典的过程，并使其高效，就像这样 -
def numpy1(cubes):
c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)
sidx = c1D.argsort()
c1Ds = c1D[sidx]
mask = np.r_[True,c1Ds[:-1]!=c1Ds[1:],True]
split_idx = np.flatnonzero(mask)
indices = [sidx[i:j] for (i,j) in zip(split_idx[:-1],split_idx[1:])]
out = dict(zip(c1Ds[mask[:-1]],indices))
return out
接下来，我们将使用 numba 包进行迭代并获得最终的可哈希字典输出。随之而来的是，将有两种解决方案 - 一种使用 numba 分别获取键和值，主调用将压缩并转换为 dict，而另一种将创建 numba-supported dict 类型和因此主调用函数不需要额外的工作。
因此，我们将有第一个 numba 解决方案：
from numba import njit
@njit
def _numba1(sidx, c1D):
out = []
n = len(sidx)
start = 0
grpID = []
for i in range(1,n):
if c1D[sidx[i]]!=c1D[sidx[i-1]]:
out.append(sidx[start:i])
grpID.append(c1D[sidx[start]])
start = i
out.append(sidx[start:])
grpID.append(c1D[sidx[start]])
return grpID,out
def numba1(cubes):
c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)
sidx = c1D.argsort()
out = dict(zip(*_numba1(sidx, c1D)))
return out
第二个 numba 解决方案为：
from numba import types
from numba.typed import Dict
int_array = types.int64[:]
@njit
def _numba2(sidx, c1D):
n = len(sidx)
start = 0
outt = Dict.empty(
key_type=types.int64,
value_type=int_array,
)
for i in range(1,n):
if c1D[sidx[i]]!=c1D[sidx[i-1]]:
outt[c1D[sidx[start]]] = sidx[start:i]
start = i
outt[c1D[sidx[start]]] = sidx[start:]
return outt
def numba2(cubes):
c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)
sidx = c1D.argsort()
out = _numba2(sidx, c1D)
return out
时间 cubes.npz 数据 -
In [4]: cubes = np.load('cubes.npz')['array']
In [5]: %timeit numpy1(cubes)
...: %timeit numba1(cubes)
...: %timeit numba2(cubes)
2.38 s ± 14.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
2.13 s ± 25.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
1.8 s ± 5.95 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
备选方案 #1： 我们可以使用 numexpr 实现进一步加速，以便大型数组计算 c1D，就像这样 -
import numexpr as ne
s0,s1 = cubes[:,0].max()+1,cubes[:,1].max()+1
d = {'s0':s0,'s1':s1,'c0':cubes[:,0],'c1':cubes[:,1],'c2':cubes[:,2]}
c1D = ne.evaluate('c0+c1*s0+c2*s0*s1',d)
这将适用于所有需要 c1D.
的地方
您可以使用 Cython：
%%cython -c-O3 -c-march=native -a
#cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True, infer_types=True
import math
import cython as cy
cimport numpy as cnp
cpdef groupby_index_dict_cy(cnp.int32_t[:, :] arr):
cdef cy.size_t size = len(arr)
result = {}
for i in range(size):
key = arr[i, 0], arr[i, 1], arr[i, 2]
if key in result:
result[key].append(i)
else:
result[key] = [i]
return result
但它不会让你比 Pandas 更快，尽管它是之后最快的（也许是基于 numpy_index 的解决方案），并且不会带来内存损失它。
到目前为止已提出的集合是 here.
在 OP 的机器中，执行时间应该接近 ~12 秒。
我正在使用激光雷达的 3D 点云。这些点由 numpy 数组给出，如下所示：
points = np.array([[61651921, 416326074, 39805], [61605255, 416360555, 41124], [61664810, 416313743, 39900], [61664837, 416313749, 39910], [61674456, 416316663, 39503], [61651933, 416326074, 39802], [61679969, 416318049, 39500], [61674494, 416316677, 39508], [61651908, 416326079, 39800], [61651908, 416326087, 39802], [61664845, 416313738, 39913], [61674480, 416316668, 39503], [61679996, 416318047, 39510], [61605290, 416360572, 41118], [61605270, 416360565, 41122], [61683939, 416313004, 41052], [61683936, 416313033, 41060], [61679976, 416318044, 39509], [61605279, 416360555, 41109], [61664837, 416313739, 39915], [61674487, 416316666, 39505], [61679961, 416318035, 39503], [61683943, 416313004, 41054], [61683930, 416313042, 41059]])
我想将我的数据分组为大小为 50*50*50 的多维数据集，以便每个多维数据集都保留我的 points 它包含的一些可哈希索引和 numpy 索引。为了进行拆分，我将 cubes = points \ 50 分配给输出到：
cubes = np.array([[1233038, 8326521, 796], [1232105, 8327211, 822], [1233296, 8326274, 798], [1233296, 8326274, 798], [1233489, 8326333, 790], [1233038, 8326521, 796], [1233599, 8326360, 790], [1233489, 8326333, 790], [1233038, 8326521, 796], [1233038, 8326521, 796], [1233296, 8326274, 798], [1233489, 8326333, 790], [1233599, 8326360, 790], [1232105, 8327211, 822], [1232105, 8327211, 822], [1233678, 8326260, 821], [1233678, 8326260, 821], [1233599, 8326360, 790], [1232105, 8327211, 822], [1233296, 8326274, 798], [1233489, 8326333, 790], [1233599, 8326360, 790], [1233678, 8326260, 821], [1233678, 8326260, 821]])
我想要的输出如下所示：
{(1232105, 8327211, 822): [1, 13, 14, 18]),
(1233038, 8326521, 796): [0, 5, 8, 9],
(1233296, 8326274, 798): [2, 3, 10, 19],
(1233489, 8326333, 790): [4, 7, 11, 20],
(1233599, 8326360, 790): [6, 12, 17, 21],
(1233678, 8326260, 821): [15, 16, 22, 23]}
我的真实点云包含多达数亿个 3D 点。进行这种分组最快的方法是什么？
我已经尝试了大多数不同的解决方案。下面是时间消耗的比较，假设点的大小约为 2000 万，不同立方体的大小约为 100 万：
Pandas [tuple(elem) -> np.array(dtype=int64)]
import pandas as pd
print(pd.DataFrame(cubes).groupby([0,1,2]).indices)
#takes 9sec
默认 [elem.tobytes() 或元组 -> 列表]
#thanks @abc:
result = defaultdict(list)
for idx, elem in enumerate(cubes):
result[elem.tobytes()].append(idx) # takes 20.5sec
# result[elem[0], elem[1], elem[2]].append(idx) #takes 27sec
# result[tuple(elem)].append(idx) # takes 50sec
numpy_indexed [整数 -> np.array]
# thanks @Eelco Hoogendoorn for his library
values = npi.group_by(cubes).split(np.arange(len(cubes)))
result = dict(enumerate(values))
# takes 9.8sec
Pandas + 降维[int -> np.array(dtype=int64)]
# thanks @Divakar for showing numexpr library:
import numexpr as ne
def dimensionality_reduction(cubes):
#cubes = cubes - np.min(cubes, axis=0) #in case some coords are negative
cubes = cubes.astype(np.int64)
s0, s1 = cubes[:,0].max()+1, cubes[:,1].max()+1
d = {'s0':s0,'s1':s1,'c0':cubes[:,0],'c1':cubes[:,1],'c2':cubes[:,2]}
c1D = ne.evaluate('c0+c1*s0+c2*s0*s1',d)
return c1D
cubes = dimensionality_reduction(cubes)
result = pd.DataFrame(cubes).groupby([0]).indices
# takes 2.5 seconds
可以下载 cubes.npz 文件 here 并使用命令
cubes = np.load('cubes.npz')['array']
查看表演时间。
您可能只是迭代并将每个元素的索引添加到相应的列表中。
from collections import defaultdict
res = defaultdict(list)
for idx, elem in enumerate(cubes):
#res[tuple(elem)].append(idx)
res[elem.tobytes()].append(idx)
可以通过使用 tobytes() 而不是将键转换为元组来进一步改进运行时间。
每组索引数量不变
方法 #1
我们可以执行 dimensionality-reduction 将 cubes 缩减为一维数组。这是基于将给定立方体数据映射到 n 维网格以计算线性索引等价物，详细讨论
N = 4 # number of indices per group
c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)
sidx = c1D.argsort()
indices = sidx.reshape(-1,N)
unq_groups = cubes[indices[:,0]]
# If you need in a zipped dictionary format
out = dict(zip(map(tuple,unq_groups), indices))
备选方案 #1 : 如果 cubes 中的整数值太大，我们可能需要 dimensionality-reduction 使得维度更短选择范围作为主轴。因此，对于这些情况，我们可以修改缩减步骤以获得 c1D，就像这样 -
s1,s2 = cubes[:,:2].max(0)+1
s = np.r_[s2,1,s1*s2]
c1D = cubes.dot(s)
方法 #2
接下来，我们可以使用 Cython-powered kd-tree for quick nearest-neighbor lookup 来获取最近的相邻索引，从而像这样解决我们的问题 -
from scipy.spatial import cKDTree
idx = cKDTree(cubes).query(cubes, k=N)[1] # N = 4 as discussed earlier
I = idx[:,0].argsort().reshape(-1,N)[:,0]
unq_groups,indices = cubes[I],idx[I]
一般情况：每组索引数量可变
我们将通过一些拆分来扩展基于 argsort 的方法以获得我们想要的输出，就像这样 -
c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)
sidx = c1D.argsort()
c1Ds = c1D[sidx]
split_idx = np.flatnonzero(np.r_[True,c1Ds[:-1]!=c1Ds[1:],True])
grps = cubes[sidx[split_idx[:-1]]]
indices = [sidx[i:j] for (i,j) in zip(split_idx[:-1],split_idx[1:])]
# If needed as dict o/p
out = dict(zip(map(tuple,grps), indices))
使用 cubes 组的一维版本作为键
我们将扩展前面列出的方法，使用 cubes 的组作为键，以简化创建字典的过程，并使其高效，就像这样 -
def numpy1(cubes):
c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)
sidx = c1D.argsort()
c1Ds = c1D[sidx]
mask = np.r_[True,c1Ds[:-1]!=c1Ds[1:],True]
split_idx = np.flatnonzero(mask)
indices = [sidx[i:j] for (i,j) in zip(split_idx[:-1],split_idx[1:])]
out = dict(zip(c1Ds[mask[:-1]],indices))
return out
接下来，我们将使用 numba 包进行迭代并获得最终的可哈希字典输出。随之而来的是，将有两种解决方案 - 一种使用 numba 分别获取键和值，主调用将压缩并转换为 dict，而另一种将创建 numba-supported dict 类型和因此主调用函数不需要额外的工作。
因此，我们将有第一个 numba 解决方案：
from numba import njit
@njit
def _numba1(sidx, c1D):
out = []
n = len(sidx)
start = 0
grpID = []
for i in range(1,n):
if c1D[sidx[i]]!=c1D[sidx[i-1]]:
out.append(sidx[start:i])
grpID.append(c1D[sidx[start]])
start = i
out.append(sidx[start:])
grpID.append(c1D[sidx[start]])
return grpID,out
def numba1(cubes):
c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)
sidx = c1D.argsort()
out = dict(zip(*_numba1(sidx, c1D)))
return out
第二个 numba 解决方案为：
from numba import types
from numba.typed import Dict
int_array = types.int64[:]
@njit
def _numba2(sidx, c1D):
n = len(sidx)
start = 0
outt = Dict.empty(
key_type=types.int64,
value_type=int_array,
)
for i in range(1,n):
if c1D[sidx[i]]!=c1D[sidx[i-1]]:
outt[c1D[sidx[start]]] = sidx[start:i]
start = i
outt[c1D[sidx[start]]] = sidx[start:]
return outt
def numba2(cubes):
c1D = np.ravel_multi_index(cubes.T, cubes.max(0)+1)
sidx = c1D.argsort()
out = _numba2(sidx, c1D)
return out
时间 cubes.npz 数据 -
In [4]: cubes = np.load('cubes.npz')['array']
In [5]: %timeit numpy1(cubes)
...: %timeit numba1(cubes)
...: %timeit numba2(cubes)
2.38 s ± 14.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
2.13 s ± 25.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
1.8 s ± 5.95 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
备选方案 #1： 我们可以使用 numexpr 实现进一步加速，以便大型数组计算 c1D，就像这样 -
import numexpr as ne
s0,s1 = cubes[:,0].max()+1,cubes[:,1].max()+1
d = {'s0':s0,'s1':s1,'c0':cubes[:,0],'c1':cubes[:,1],'c2':cubes[:,2]}
c1D = ne.evaluate('c0+c1*s0+c2*s0*s1',d)
这将适用于所有需要 c1D.
您可以使用 Cython：
%%cython -c-O3 -c-march=native -a
#cython: language_level=3, boundscheck=False, wraparound=False, initializedcheck=False, cdivision=True, infer_types=True
import math
import cython as cy
cimport numpy as cnp
cpdef groupby_index_dict_cy(cnp.int32_t[:, :] arr):
cdef cy.size_t size = len(arr)
result = {}
for i in range(size):
key = arr[i, 0], arr[i, 1], arr[i, 2]
if key in result:
result[key].append(i)
else:
result[key] = [i]
return result
但它不会让你比 Pandas 更快，尽管它是之后最快的（也许是基于 numpy_index 的解决方案），并且不会带来内存损失它。
到目前为止已提出的集合是 here.
在 OP 的机器中，执行时间应该接近 ~12 秒。