沿轴的 2D 数组的 Numpy 总和 = 1,浮动范围
Numpy sum of 2D array along axis=1, floating range
我想在第二个轴上对二维数组求和,但范围是可变的。没有向量化它是:`
import numpy as np
nx = 3
ny = 5
a = np.ones((nx, ny))
left_bnd = np.array([0, 1, 0])
right_bnd = np.array([2, 2, 4])
b = np.zeros(nx)
for jx in range(nx):
b[jx] = np.sum(a[jx, left_bnd[jx]: right_bnd[jx]])
print(b)
输出 b 是 [2. 1. 4.]
我喜欢矢量化循环,有点
b = np.sum(a[:, left_bnd[:]: right_bnd[:], axis=1)
为了加快计算速度,因为我的“n”一般是1e6几个。不幸的是我找不到合适的工作语法。
在 for 循环中手动求和的 jitted numba
实现大约快 100 倍。在 numba
函数内使用带有切片的 np.sum
只有一半的速度。此解决方案假定所有切片都在有效范围内。
为基准测试生成足够大的样本数据
import numpy as np
import numba as nb
np.random.seed(42) # just for reproducibility
n, m = 5000, 100
a = np.random.rand(n,m)
bnd_l, bnd_r = np.sort(np.random.randint(m+1, size=(n,2))).T
与 numba
合影。请确保通过 运行 函数对编译的热代码进行基准测试至少两次。
@nb.njit
def slice_sum(a, bnd_l, bnd_r):
b = np.zeros(a.shape[0])
for j in range(a.shape[0]):
for i in range(bnd_l[j], bnd_r[j]):
b[j] += a[j,i]
return b
slice_sum(a, bnd_l, bnd_r)
输出
# %timeit 1000 loops, best of 5: 297 µs per loop
array([ 4.31060848, 35.90684722, 38.03820523, ..., 37.9578962 ,
3.61011028, 6.53631388])
在 python 循环中使用 numpy
(这是一个很好、简单的实现)
b = np.zeros(n)
for j in range(n):
b[j] = np.sum(a[ j, bnd_l[j] : bnd_r[j] ])
b
输出
# %timeit 10 loops, best of 5: 29.2 ms per loop
array([ 4.31060848, 35.90684722, 38.03820523, ..., 37.9578962 ,
3.61011028, 6.53631388])
验证结果是否相等
np.testing.assert_allclose(slice_sum(a, bnd_l, bnd_r), b)
这是一个纯 numpy 解决方案,其速度接近已发布的 numba 解决方案。它利用 reduceat
但设置非常复杂。
def slice_sum_np(a, left_bnd, right_bnd):
nx, ny = a.shape
linear_indices = np.c_[left_bnd, right_bnd] + ny * np.arange(nx)[:,None]
sums = np.add.reduceat(a.ravel(), linear_indices.ravel())[::2]
# account for reduceat special case
sums[left_bnd >= right_bnd] = 0
return sums
我想在第二个轴上对二维数组求和,但范围是可变的。没有向量化它是:`
import numpy as np
nx = 3
ny = 5
a = np.ones((nx, ny))
left_bnd = np.array([0, 1, 0])
right_bnd = np.array([2, 2, 4])
b = np.zeros(nx)
for jx in range(nx):
b[jx] = np.sum(a[jx, left_bnd[jx]: right_bnd[jx]])
print(b)
输出 b 是 [2. 1. 4.] 我喜欢矢量化循环,有点
b = np.sum(a[:, left_bnd[:]: right_bnd[:], axis=1)
为了加快计算速度,因为我的“n”一般是1e6几个。不幸的是我找不到合适的工作语法。
在 for 循环中手动求和的 jitted numba
实现大约快 100 倍。在 numba
函数内使用带有切片的 np.sum
只有一半的速度。此解决方案假定所有切片都在有效范围内。
为基准测试生成足够大的样本数据
import numpy as np
import numba as nb
np.random.seed(42) # just for reproducibility
n, m = 5000, 100
a = np.random.rand(n,m)
bnd_l, bnd_r = np.sort(np.random.randint(m+1, size=(n,2))).T
与 numba
合影。请确保通过 运行 函数对编译的热代码进行基准测试至少两次。
@nb.njit
def slice_sum(a, bnd_l, bnd_r):
b = np.zeros(a.shape[0])
for j in range(a.shape[0]):
for i in range(bnd_l[j], bnd_r[j]):
b[j] += a[j,i]
return b
slice_sum(a, bnd_l, bnd_r)
输出
# %timeit 1000 loops, best of 5: 297 µs per loop
array([ 4.31060848, 35.90684722, 38.03820523, ..., 37.9578962 ,
3.61011028, 6.53631388])
在 python 循环中使用 numpy
(这是一个很好、简单的实现)
b = np.zeros(n)
for j in range(n):
b[j] = np.sum(a[ j, bnd_l[j] : bnd_r[j] ])
b
输出
# %timeit 10 loops, best of 5: 29.2 ms per loop
array([ 4.31060848, 35.90684722, 38.03820523, ..., 37.9578962 ,
3.61011028, 6.53631388])
验证结果是否相等
np.testing.assert_allclose(slice_sum(a, bnd_l, bnd_r), b)
这是一个纯 numpy 解决方案,其速度接近已发布的 numba 解决方案。它利用 reduceat
但设置非常复杂。
def slice_sum_np(a, left_bnd, right_bnd):
nx, ny = a.shape
linear_indices = np.c_[left_bnd, right_bnd] + ny * np.arange(nx)[:,None]
sums = np.add.reduceat(a.ravel(), linear_indices.ravel())[::2]
# account for reduceat special case
sums[left_bnd >= right_bnd] = 0
return sums