理解 numpy 中奇怪的布尔二维数组索引行为答案

【问题标题】：Understanding weird boolean 2d-array indexing behavior in numpy理解 numpy 中奇怪的布尔二维数组索引行为
【发布时间】：2011-10-19 11:46:28
【问题描述】：

为什么会这样：

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

【问题讨论】：

标签： python numpy

【解决方案1】：

Numpy 参考文档的 page on indexing 包含答案，但需要仔细阅读。

这里的答案是，使用布尔值索引等同于使用通过首先使用np.nonzero 转换布尔数组获得的整数数组进行索引。因此，对于布尔数组m1、m2

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

【讨论】：

这真的没有意义。我会在邮件列表中询问为什么会这样。
scipy.org/Cookbook/Indexing p。 14 on Multidimenional Boolean Indexing 说“看看 numpy 的掩码数组工具......显而易见的方法没有给出正确的答案。” （该文档写得很好，需要更新。）
@denis，大约在 2013 年，该文档确实很好地解释了它。但是，如果你用 google numpy 逻辑索引，出现的文档是docs.scipy.org/doc/numpy/reference/arrays.indexing.html，几乎没有解释。
如果m1.nonzero()[0].shape == m2.nonzero()[0].shape 成功，至少在当前版本中。
而且不等同于[a[i,i] ...]，而是等同于[a[i,j] for i,j in zip(m1.nonzero()[0], m2.nonzero()[0])]

【解决方案2】：

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 或硬件版本而变化。

【讨论】：