具有重复索引的 numpy 数组的矢量化分配 (d[i,j,i,j] = s[i,j])答案

【问题标题】：Vectorized assignment for numpy array with repeated indices (d[i,j,i,j] = s[i,j])具有重复索引的 numpy 数组的矢量化分配 (d[i,j,i,j] = s[i,j])
【发布时间】：2017-10-06 21:29:06
【问题描述】：

如何设置

d[i,j,i,j] = s[i,j]

使用“NumPy”而不使用 for 循环？

我尝试了以下方法：

l1=range(M)
l2=range(N)
d[l1,l2,l1,l2] = s[l1,l2]

【问题讨论】：

d[1,2,3,4] 的值怎么样（所以对于 d[i,j,k,l]，其中 k != i 和 l != j）？
其余值为0
d 的形状是什么？它是用零初始化的吗？如果有，它的形状是什么？
我会对s 的形状感兴趣。是(M, N) 与您用于ranges 的M 和N 吗？

标签： python numpy multidimensional-array indexing

【解决方案1】：

如果您考虑一下，这与创建形状为(m*n, m*n) 的2D 数组并将s 中的值分配到对角线位置相同。要使最终输出为4D，我们只需要在最后进行整形。这基本上是在下面实现-

m,n = s.shape
d = np.zeros((m*n,m*n),dtype=s.dtype)
d.ravel()[::m*n+1] = s.ravel()
d.shape = (m,n,m,n)

运行时测试

方法-

# @MSeifert's solution
def assign_vals_ix(s):    
    d = np.zeros((m, n, m, n), dtype=s.dtype)
    l1 = range(m)
    l2 = range(n)
    d[np.ix_(l1,l2)*2] = s[np.ix_(l1,l2)]
    return d

# Proposed in this post
def assign_vals(s):
    m,n = s.shape
    d = np.zeros((m*n,m*n),dtype=s.dtype)
    d.ravel()[::m*n+1] = s.ravel()
    return d.reshape(m,n,m,n)

# Using a strides based approach
def assign_vals_strides(a):
    m,n = a.shape
    p,q = a.strides

    d = np.zeros((m,n,m,n),dtype=a.dtype)
    out_strides = (q*(n*m*n+n),(m*n+1)*q)
    d_view = np.lib.stride_tricks.as_strided(d, (m,n), out_strides)
    d_view[:] = a
    return d

时间安排 -

In [285]: m,n = 10,10
     ...: s = np.random.rand(m,n)
     ...: d = np.zeros((m,n,m,n))
     ...: 

In [286]: %timeit assign_vals_ix(s)
10000 loops, best of 3: 21.3 µs per loop

In [287]: %timeit assign_vals_strides(s)
100000 loops, best of 3: 9.37 µs per loop

In [288]: %timeit assign_vals(s)
100000 loops, best of 3: 4.13 µs per loop

In [289]: m,n = 20,20
     ...: s = np.random.rand(m,n)
     ...: d = np.zeros((m,n,m,n))


In [290]: %timeit assign_vals_ix(s)
10000 loops, best of 3: 60.2 µs per loop

In [291]: %timeit assign_vals_strides(s)
10000 loops, best of 3: 41.8 µs per loop

In [292]: %timeit assign_vals(s)
10000 loops, best of 3: 35.5 µs per loop

【讨论】：

也可以使用非连续输入 d 和 numpy.lib.stride_tricks 来执行此操作，或者使用 numpy.diagonal 两次。（ravel() 可能会为非连续输入制作副本。）
（不是完全相同的平面、跨步切片分配，但我们可以查看形状像 s 的 d 并对其进行切片分配。）
假设s的形状是(m, n)和d(m*n,m*n)是否明智？我认为这里的问题更普遍，因为在 cmets 中他已经说过"All the rest values are 0"，这意味着形状不一定必须匹配l1 和l2。还是我误解了你的答案？
@user2357112 好主意！实施并计时。不过，它看起来并没有改进扁平视图分配。
@MSeifert：在我看来，该评论（以及它所回复的评论）是在谈论不是i,j,i,j 形式的索引处的元素，例如1,2,3,4。

【解决方案2】：

您可以使用integer array indexing（使用np.ix_ 创建广播索引）：

d[np.ix_(l1,l2)*2] = s[np.ix_(l1,l2)]

第一次必须复制索引时（您需要 [i, j, i, j] 而不仅仅是 [i, j]），这就是我将 np.ix_ 返回的 tuple 乘以 2 的原因。

例如：

>>> d = np.zeros((10, 10, 10, 10), dtype=int)
>>> s = np.arange(100).reshape(10, 10)
>>> l1 = range(3)
>>> l2 = range(5)
>>> d[np.ix_(l1,l2)*2] = s[np.ix_(l1,l2)]

并确保分配了正确的值：

>>> # Assert equality for the given condition
>>> for i in l1:
...     for j in l2:
...         assert d[i, j, i, j] == s[i, j]

>>> # Interactive tests
>>> d[0, 0, 0, 0], s[0, 0]
(0, 0)
>>> d[1, 2, 1, 2], s[1, 2]
(12, 12)
>>> d[2, 0, 2, 0], s[2, 0]
(20, 20)
>>> d[2, 4, 2, 4], s[2, 4]
(24, 24)

【讨论】：