题 理解numpy中奇怪的布尔2d数组索引行为


为什么这样做:

a=np.random.rand(10,20)
x_range=np.arange(10)
y_range=np.arange(20)

a_tmp=a[x_range<5,:]
b=a_tmp[:,np.in1d(y_range,[3,4,8])]

这不是:

a=np.random.rand(10,20)
x_range=np.arange(10)
y_range=np.arange(20)    

b=a[x_range<5,np.in1d(y_range,[3,4,8])]

15
2017-10-19 11:46


起源




答案:


Numpy参考文档 关于索引的页面 包含答案,但需要仔细阅读。

这里的答案是使用布尔值进行索引等效于使用首先使用布尔数组转换得到的整数数组进行索引 np.nonzero。因此,使用布尔数组 m1m2

a[m1, m2] == a[m1.nonzero(), m2.nonzero()]

哪个(当它成功时,即, m1.nonzero().shape == m2.nonzero().shape)相当于:

[a[i, i] for i in range(a.shape[0]) if m1[i] and m2[i]]

我不确定为什么它被设计成这样工作 - 通常,这是  你想要什么。

为了获得更直观的结果,您可以改为

a[np.ix_(m1, m2)]

产生相当于的结果

[[a[i,j] for j in range(a.shape[1]) if m2[j]] for i in range(a.shape[0]) if m1[i]]

19
2017-10-19 12:14



这真的没有意义。我会在maillist中询问为什么会这样。 - tillsten
scipy.org/Cookbook/Indexing 页。 14关于多维布尔索引说“看看numpy的蒙面数组工具......显而易见的方法没有给出正确答案。” (该文件写得很好,需要更新。) - denis
@denis,大约在2013年该文件确实解释得相当好。但是,如果你谷歌numpy逻辑索引,出现的文件是 docs.scipy.org/doc/numpy/reference/arrays.indexing.html 而且几乎没有解释。 - John
如果成功的话 m1.nonzero()[0].shape == m2.nonzero()[0].shape,至少在当前版本中。 - kasal
它并不等同于 [a[i,i] ...] 但是 [a[i,j] for i,j in zip(m1.nonzero()[0], m2.nonzero()[0])] - kasal


替代 np.ix_ 是将布尔数组转换为整数数组(使用 np.nonzero()),然后使用 np.newaxis 创建正确形状的数组以利用广播。

import numpy as np

a=np.random.rand(10,20)
x_range=np.arange(10)
y_range=np.arange(20)

a_tmp=a[x_range<5,:]
b_correct=a_tmp[:,np.in1d(y_range,[3,4,8])]

m1=(x_range<5).nonzero()[0]
m2=np.in1d(y_range,[3,4,8]).nonzero()
b=a[m1[:,np.newaxis], m2]
assert np.allclose(b,b_correct)

b2=a[np.ix_(x_range<5,np.in1d(y_range,[3,4,8]))]
assert np.allclose(b2,b_correct)

np.ix_ 往往比双索引慢。 长格式解决方案似乎更快一点:

长表

In [83]: %timeit a[(x_range<5).nonzero()[0][:,np.newaxis], (np.in1d(y_range,[3,4,8])).nonzero()[0]]
10000 loops, best of 3: 131 us per loop

双索引

In [85]: %timeit a[x_range<5,:][:,np.in1d(y_range,[3,4,8])]
10000 loops, best of 3: 144 us per loop

使用np.ix_

In [84]: %timeit a[np.ix_(x_range<5,np.in1d(y_range,[3,4,8]))]
10000 loops, best of 3: 160 us per loop

注意:在您的计算机上测试这些计时是个好主意,因为排名可能会根据您的Python,numpy或硬件版本而改变。


4
2017-10-19 12:59