在Python中使用数组索引切片2D数组

Question

我正在处理一维numpy数组的切片。为了选择这些切片，我将它们的索引存储在数组中。例如，我有：

mat = np.zeros([xdim,ydim], float)
xmin = np.array([...]) # Array of minimum indices in x
xmax = np.array([...]) # Array of maximum indices in x
ymin = np.array([...]) # Array of minimum indices in y
ymax = np.array([...]) # Array of maximum indices in y
value = np.array([...]) # Values

... 仅表示先前计算的一些整数。所有数组都是明确定义的，并且长度约为 ~265,000。我想要做的事情如下：

mat[xmin:xmax, ymin:ymax] += value

以这样的方式，我将拥有以下第一项：

mat[xmin[0]:xmax[0], ymin[0]:ymax[0]] += value[0]
mat[xmin[1]:xmax[1], ymin[1]:ymax[1]] += value[1]

对于数组的约265000个元素，以此类推。不幸的是，我的写法不起作用，并且抛出错误：IndexError: invalid slice。

我一直在尝试使用np.meshgrid，如此处所建议： NumPy: use 2D index array from argmin in a 3D slice，但它还没有为我工作。另外，我正在寻找一种Pythonic的方法来做到这一点，避免for循环。

非常感谢任何帮助！

谢谢！

- fbecerra

---

np.array() 是无效的。 - zhangxaochen

抱歉造成困惑，我添加了一些文本以使它更清晰明了。 - fbecerra

你的切片是否重叠？mat 多大了？ - Jaime

是的，它们可能会重叠。 mat的大小是由用户定义的，因此它可以非常小或非常大。 - fbecerra

你可能想尝试像Cython或Numba这样的工具，而不是寻找切片技巧。 - user2357112

谢谢建议！我考虑用Cython/Numba编写代码，但目前只需要使用numpy。 - fbecerra

Answer 1

我认为如果不使用Cython或类似的东西，很难找到令人满意的矢量化解决方案。让我先概述一下纯numpy解决方案可能是什么样子，这应该可以说明为什么这可能不是一个非常好的方法。

首先，让我们看一个一维情况。在numpy中，一堆切片并没有太多作用，因此第一个任务是将它们扩展为单个索引。假设你的数组是:

mat = np.zeros((10,))
x_min = np.array([2, 5, 3, 1])
x_max = np.array([5, 9, 8, 7])
value = np.array([0.2, 0.6, 0.1, 0.9])

接下来的代码将切片限制扩展为（可能重复的）索引和值的列表，使用bincount将它们连接起来，并将它们添加到原始的mat中：

x_len = x_max - x_min
x_cum_len = np.cumsum(x_len)
x_idx = np.arange(x_cum_len[-1])
x_idx[x_len[0]:] -= np.repeat(x_cum_len[:-1], x_len[1:])
x_idx += np.repeat(x_min, x_len)
x_val = np.repeat(value, x_len)
x_cumval = np.bincount(x_idx, weights=x_val)
mat[:len(x_cumval)] += x_cumval

>>> mat
array([ 0. ,  0.9,  1.1,  1.2,  1.2,  1.6,  1.6,  0.7,  0.6,  0. ])

虽然这可以扩展到您的二维情况，但这一点都不平凡，而且事情开始变得难以理解：

mat = np.zeros((10, 10))
x_min = np.array([2, 5, 3, 1])
x_max = np.array([5, 9, 8, 7])
y_min = np.array([1, 7, 2, 6])
y_max = np.array([6, 8, 6, 9])
value = np.array([0.2, 0.6, 0.1, 0.9])

x_len = x_max - x_min
y_len = y_max - y_min
total_len = x_len * y_len
x_cum_len = np.cumsum(x_len)
x_idx = np.arange(x_cum_len[-1])
x_idx[x_len[0]:] -= np.repeat(x_cum_len[:-1], x_len[1:])
x_idx += np.repeat(x_min, x_len)
x_val = np.repeat(value, x_len)
y_min_ = np.repeat(y_min, x_len)
y_len_ = np.repeat(y_len, x_len)
y_cum_len = np.cumsum(y_len_)
y_idx = np.arange(y_cum_len[-1])
y_idx[y_len_[0]:] -= np.repeat(y_cum_len[:-1], y_len_[1:])
y_idx += np.repeat(y_min_, y_len_)
x_idx_ = np.repeat(x_idx, y_len_)
xy_val = np.repeat(x_val, y_len_)
xy_idx = np.ravel_multi_index((x_idx_, y_idx), dims=mat.shape)
xy_cumval = np.bincount(xy_idx, weights=xy_val)
mat.ravel()[:len(xy_cumval)] += xy_cumval

这将产生：

>>> mat
array([[ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ],
       [ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0.9,  0.9,  0.9,  0. ],
       [ 0. ,  0.2,  0.2,  0.2,  0.2,  0.2,  0.9,  0.9,  0.9,  0. ],
       [ 0. ,  0.2,  0.3,  0.3,  0.3,  0.3,  0.9,  0.9,  0.9,  0. ],
       [ 0. ,  0.2,  0.3,  0.3,  0.3,  0.3,  0.9,  0.9,  0.9,  0. ],
       [ 0. ,  0. ,  0.1,  0.1,  0.1,  0.1,  0.9,  1.5,  0.9,  0. ],
       [ 0. ,  0. ,  0.1,  0.1,  0.1,  0.1,  0.9,  1.5,  0.9,  0. ],
       [ 0. ,  0. ,  0.1,  0.1,  0.1,  0.1,  0. ,  0.6,  0. ,  0. ],
       [ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0.6,  0. ,  0. ],
       [ 0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ,  0. ]])

但是，如果您有265,000个任意大小的二维切片，则索引数组很快就会达到数百万个项目。必须处理如此多的数据读取和写入可能会抵消使用numpy带来的速度提升。坦率地说，我怀疑这根本不是一个好选项，即使仅仅因为您的代码将变得非常晦涩。

- Jaime