高效地对大型拆分数组执行 numpy.sum（或 scipy.integrate.simps()）
Perform numpy.sum (or scipy.integrate.simps()) on large splitted array efficiently
让我们考虑一个非常大的 numpy 数组 a (M, N)。
其中 M 通常可以是 1 或 100，N 可以是 10-100,000,000
我们有索引数组，可以沿轴=1 将其拆分为多个 (K = 1,000,000)。
我们希望在每个子数组和 return (M, K) 数组上高效地执行类似沿轴=1 积分（np.sum 采取最简单形式）的操作。
有问题的@Divakar 提出了一个优雅而有效的解决方案[41920367]，但我的理解是它只适用于所有子数组具有相同形状的情况，这允许重塑。
但在我们的例子中，子数组的形状不同，到目前为止，这迫使我在索引上循环……请让我摆脱痛苦……
例子
a = np.random.random((10, 100000000))
ind = np.sort(np.random.randint(10, 9000000, 1000000))
子数组的大小不均匀：
sizes = np.diff(ind)
print(sizes.min(), size.max())
2, 8732
到目前为止，我找到的最好的是：
output = np.concatenate([np.sum(vv, axis=1)[:, None] for vv in np.split(a, ind, axis=1)], axis=1)
numpy 和 scipy 的可能功能请求：
如果循环确实不可避免，至少在 C 语言中的 numpy 和 scipy.integrate.simps（或 romb）函数中完成它可能会加快输出速度。
像
output = np.sum(a, axis=1, split_ind=ind)
output = scipy.integrate.simps(a, x=x, axis=1, split_ind=ind)
output = scipy.integrate.romb(a, x=x, axis=1, split_ind=ind)
非常欢迎！
（其中 x 本身可以拆分，也可以不拆分）
旁注：
在尝试这个例子时，我注意到这些数字几乎总是有一个元素的大小等于 0（sizes.min() 几乎总是零）。
这对我来说很奇怪，因为我们在 10 到 9,000,000 之间选择 10,000 个整数，相同数字出现两次的几率（例如 diff = 0）应该接近 0。它似乎非常接近 1。
这是由于 np.random.randint 背后的算法吗？
你要的是np.add.reduceat
output = np.add.reduceat(a, ind, axis = 1)
output.shape
Out[]: (10, 1000000)
Universal Functions (ufunc) 是 numpy
中非常强大的工具
至于重复索引，那只是Birthday Problem的突然出现。
太棒了！
谢谢！在我的 VM Cent OS 6.9 上，我得到以下结果：
在[71]中：a = np.random.random((10, 10000000))
在 [72] 中：ind = np.unique(np.random.randint(10, 9000000, 100000))
在 [73] 中：ind2 = np.append([0], ind)
In [74]: out = np.concatenate([np.sum(vv, axis=1)[:, None] for vv in np.split(a , ind, 轴=1)], 轴=1)
In [75]: out2 = np.add.reduceat(a, ind2, axis=1)
在 [83] 中：np.allclose(out, out2)
输出[83]：真
In [84]: %timeit out = np.concatenate([np.sum(vv, axis=1)[:, None] for vv in np.split (a, ind, axis=1)], axis=1)
每个循环 2.7 s ± 40.4 ms（7 次运行的平均值 ± 标准偏差，每次 1 个循环）
在 [85] 中：%timeit out2 = np.add.reduceat(a, ind2, axis=1)
每个循环 179 毫秒 ± 15.9 毫秒（7 次运行的平均值 ± 标准偏差，每次 1 个循环）
与列表串联相比，这是一个很好的 93% 的速度增益（或更快的 15 倍）:-)
太棒了！
让我们考虑一个非常大的 numpy 数组 a (M, N)。
其中 M 通常可以是 1 或 100，N 可以是 10-100,000,000
我们有索引数组，可以沿轴=1 将其拆分为多个 (K = 1,000,000)。
我们希望在每个子数组和 return (M, K) 数组上高效地执行类似沿轴=1 积分（np.sum 采取最简单形式）的操作。
有问题的@Divakar 提出了一个优雅而有效的解决方案[41920367]
但在我们的例子中，子数组的形状不同，到目前为止，这迫使我在索引上循环……请让我摆脱痛苦……
例子
a = np.random.random((10, 100000000))
ind = np.sort(np.random.randint(10, 9000000, 1000000))
子数组的大小不均匀：
sizes = np.diff(ind)
print(sizes.min(), size.max())
2, 8732
到目前为止，我找到的最好的是：
output = np.concatenate([np.sum(vv, axis=1)[:, None] for vv in np.split(a, ind, axis=1)], axis=1)
numpy 和 scipy 的可能功能请求：
如果循环确实不可避免，至少在 C 语言中的 numpy 和 scipy.integrate.simps（或 romb）函数中完成它可能会加快输出速度。 像
output = np.sum(a, axis=1, split_ind=ind)
output = scipy.integrate.simps(a, x=x, axis=1, split_ind=ind)
output = scipy.integrate.romb(a, x=x, axis=1, split_ind=ind)
非常欢迎！ （其中 x 本身可以拆分，也可以不拆分）
旁注：
在尝试这个例子时，我注意到这些数字几乎总是有一个元素的大小等于 0（sizes.min() 几乎总是零）。
这对我来说很奇怪，因为我们在 10 到 9,000,000 之间选择 10,000 个整数，相同数字出现两次的几率（例如 diff = 0）应该接近 0。它似乎非常接近 1。
这是由于 np.random.randint 背后的算法吗？
你要的是np.add.reduceat
output = np.add.reduceat(a, ind, axis = 1)
output.shape
Out[]: (10, 1000000)
Universal Functions (ufunc) 是 numpy
至于重复索引，那只是Birthday Problem的突然出现。
太棒了！
谢谢！在我的 VM Cent OS 6.9 上，我得到以下结果：
在[71]中：a = np.random.random((10, 10000000))
在 [72] 中：ind = np.unique(np.random.randint(10, 9000000, 100000))
在 [73] 中：ind2 = np.append([0], ind)
In [74]: out = np.concatenate([np.sum(vv, axis=1)[:, None] for vv in np.split(a , ind, 轴=1)], 轴=1)
In [75]: out2 = np.add.reduceat(a, ind2, axis=1)
在 [83] 中：np.allclose(out, out2)
输出[83]：真
In [84]: %timeit out = np.concatenate([np.sum(vv, axis=1)[:, None] for vv in np.split (a, ind, axis=1)], axis=1)
每个循环 2.7 s ± 40.4 ms（7 次运行的平均值 ± 标准偏差，每次 1 个循环）
在 [85] 中：%timeit out2 = np.add.reduceat(a, ind2, axis=1)
每个循环 179 毫秒 ± 15.9 毫秒（7 次运行的平均值 ± 标准偏差，每次 1 个循环）
与列表串联相比，这是一个很好的 93% 的速度增益（或更快的 15 倍）:-) 太棒了！