根据单元格索引填充 numpy 数组答案

【问题标题】：Populate numpy array based on cell indices根据单元格索引填充 numpy 数组
【发布时间】：2018-11-30 16:08:59
【问题描述】：

我正在尝试创建一个 3d 数组，其单元格条目将根据单元格索引进行计算。具体来说，我想要那个单元格(i,j,k) = sqrt(i+j+k)。

使用以下 for 循环很容易做到这一点：

N=10
A=np.zeros((N,N,N))

for i in range(N):
    for j in range(N):
        for k in range(N):
            A[i][j][k] = np.sqrt(i+j+k)

我想知道 numpy 是否具有使这些嵌套 for 循环变得多余的内置函数。

【问题讨论】：

标签： python numpy

【解决方案1】：

最简单和高效的方法是使用np.ogrid 使用开放网格，然后执行相关操作 -

I,J,K = np.ogrid[:N,:N,:N]
A = np.sqrt(I+J+K)

或者np.sum 用于单线开放网格的广播总和 -

A = np.sqrt(np.sum(np.ogrid[:N,:N,:N]))

【讨论】：

太棒了！我会在几个小时内接受你的回答......也许其他人有更蟒蛇式的方式。我对此表示怀疑
我不知道np.sum 可以做到这一点。（将不强制到单个数组的东西相加）。不过，似乎比内置的 sum 慢一些（通过微小的恒定开销）。
@DouglasJamesBock 我不认为有人可以提供更好的答案，至少我不能。 PS：我从 Divakar 学了numpy

【解决方案2】：

您可以使用np.arange，然后使用np.newaxis 创建不同的维度。用一个简单的sum 和np.sqrt 来完成这项工作：

arr = np.arange(N)
A = np.sqrt(arr + arr[:,np.newaxis]+ arr[:,np.newaxis,np.newaxis])

你得到相同的结果：

N = 10
arr = np.arange(N)
A = np.sqrt(arr + arr[:,np.newaxis]+ arr[:,np.newaxis,np.newaxis])
B = np.sqrt(np.sum(np.ogrid[:N,:N,:N]))
print ((A==B).all())
#True

这种方法比使用np.ogrid快一点：

N = 10
%timeit arr = np.arange(N); A = np.sqrt(arr + arr[:,np.newaxis]+ arr[:,np.newaxis,np.newaxis])
#18.6 µs ± 3.29 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
%timeit A = np.sqrt(np.sum(np.ogrid[:N,:N,:N]))
#58.5 µs ± 8.01 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

【讨论】：

这太厚脸皮了——利用所有迭代都在N 的事实 :) 通过避免为所有三个迭代器创建范围，您可以获得性能。
@Divakar 谢谢 :) 实际上，即使你不创建 arr 而是每次都直接创建 np.arange(N)，它仍然比 ogrid 快。所以我想如果你在每个维度上都有不同的长度，它仍然会更快。也许你增加N就不会这样了，我没有尝试更高N
Ben.T, @Divakar 许多numpy 便利功能都有可怕的开销...

【解决方案3】：

这对于大型N 更快，但可能被视为作弊 ;-)

充分利用高度规律性和重复性的模式，节省了大量的平方根计算。

def cheat(N):
    values = np.sqrt(np.arange(3*N-2))
    result = np.lib.stride_tricks.as_strided(values, (N, N, N), 3*values.strides)
    return np.ascontiguousarray(result)

如果您可以使用非连续的只读（通过所有实际方式）视图，这可以大大加快：

def cheat_nc_view(N):
    values = np.sqrt(np.arange(3*N-2))
    return np.lib.stride_tricks.as_strided(values, (N, N, N), 3*values.strides)

供参考：

def cheek(N):
    arr = np.arange(N)
    return np.sqrt(arr + arr[:,np.newaxis] + arr[:,np.newaxis,np.newaxis])

>>> np.all(cheek(20) == cheat(20))
True
>>> np.all(cheek(200) == cheat_nc_view(200))
True

时间安排：

>>> timeit(lambda: cheek(20), number=1000)
0.05387042500660755
>>> timeit(lambda: cheat(20), number=1000)
0.020798540994292125
>>> timeit(lambda: cheat_nc_view(20), number=1000)
0.010791150998556986

>>> timeit(lambda: cheek(200), number=100)
6.823299437994137
>>> timeit(lambda: cheat(200), number=100)
2.0583883369981777
>>> timeit(lambda: cheat_nc_view(200), number=100)
0.0014881940151099116

【讨论】：