2011-10-19 65 views
14

为什么这个工作:了解怪异布尔二维数组索引行为numpy的

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])] 

回答

19

的numpy的参考文档的page on indexing包含了答案,但需要一点点仔细阅读。

这里答案是索引用布尔等价于通过首先用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]] 
+1

这真的没有意义。我会问在maillist为什么这样。 – tillsten

+1

[scipy.org/Cookbook/Indexing](http://scipy.org/Cookbook/Indexing)p。关于多维布尔索引的14说:“看看numpy的蒙面数组工具......显而易见的方法并不能给出正确的答案。” (该文件写得很好,需要更新。) – denis

+0

@denis,大约在2013年,该文件解释得相当好。但是,如果您使用Google numpy逻辑索引,那么出现的文档是http://docs.scipy.org/doc/numpy/reference/arrays.indexing.html,这一点也没有解释清楚。 – John

4

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_往往比双重索引慢。 长形式解看来是快一点:

长式:使用np.ix_

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

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 [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的,还是硬件的版本中改变一个好主意。