Python/numpy 棘手的切片问题答案

【问题标题】：Python/numpy tricky slicing problemPython/numpy 棘手的切片问题
【发布时间】：2009-10-19 07:28:08
【问题描述】：

我对一些 numpy 的东西有疑问。我需要一个 numpy 数组以不寻常的方式运行，方法是将切片作为我切片的数据的视图返回，而不是副本。所以这是我想做的一个例子：

假设我们有一个像这样的简单数组：

a = array([1, 0, 0, 0])

我想用数组中的前一个条目更新数组中的连续条目（从左到右移动），使用如下语法：

a[1:] = a[0:3]

这将得到以下结果：

a = array([1, 1, 1, 1])

或者是这样的：

a[1:] = 2*a[:3]
# a = [1,2,4,8]

为了进一步说明，我想要以下行为：

for i in range(len(a)):
    if i == 0 or i+1 == len(a): continue
    a[i+1] = a[i]

除了我想要numpy的速度。

numpy 的默认行为是获取切片的副本，所以我实际得到的是这样的：

a = array([1, 1, 0, 0])

我已经将此数组作为 ndarray 的子类，因此如果需要，我可以对其进行进一步更改，我只需要不断更新右侧的切片，因为它会更新左侧的切片边。

我是在做梦还是这魔法可能？

更新：这都是因为我或多或少地尝试使用 Gauss-Seidel 迭代来解决线性代数问题。这是一个涉及调和函数的特殊情况，我试图避免进入这个，因为它真的没有必要并且可能会进一步混淆事情，但是这里是。

算法是这样的：

while not converged:
    for i in range(len(u[:,0])):
        for j in range(len(u[0,:])):
            # skip over boundary entries, i,j == 0 or len(u)
            u[i,j] = 0.25*(u[i-1,j] + u[i+1,j] + u[i, j-1] + u[i,j+1])

对吗？但是您可以通过两种方式做到这一点，Jacobi 涉及使用其邻居更新每个元素，而不考虑在 while 循环循环之前您已经进行的更新，要在循环中执行此操作，您将复制数组，然后从复制的数组中更新一个数组。然而，Gauss-Seidel 使用您已经为每个 i-1 和 j-1 条目更新的信息，因此不需要副本，循环应该基本上“知道”，因为在每个单个元素更新后重新评估了数组.也就是说，每次调用 u[i-1,j] 或 u[i,j-1] 这样的条目时，前面循环计算的信息都会出现。

我想用一个使用 numpy 切片的简洁代码行来替换这种缓慢而丑陋的嵌套循环情况：

u[1:-1,1:-1] = 0.25(u[:-2,1:-1] + u[2:,1:-1] + u[1:-1,:-2] + u[1:-1,2:])

但结果是 Jacobi 迭代，因为当您获取切片时： u[:,-2,1:-1] 您复制数据，因此切片不知道所做的任何更新。现在 numpy 仍然循环对吗？它不是并行的，它只是一种更快的循环方式，看起来像 python 中的并行操作。我想通过破解 numpy 来利用这种行为，以便在我获取切片时返回指针而不是副本。对？然后每次 numpy 循环时，该切片都会“更新”或实际上只是复制更新中发生的任何事情。为此，我需要将数组两侧的切片作为指针。

无论如何，如果那里有一个非常聪明的人那么棒，但我几乎已经让自己相信唯一的答案是在 C 中循环。

【问题讨论】：

很抱歉，我不太明白你的问题...无论如何，你用 a.copy() 试过了吗？
您的要求在很大程度上是毫无意义的。为什么你可以说 a[0:3] 并且只有 a[0] 的意思？那是无法理解的。
只是对您正在尝试做的事情进行类比，因为大多数人似乎并不理解它：这有点像将摄像机指向显示其自身输出的电视屏幕。所以你想要的是某种递归赋值——但我不相信有任何保证这会稳定为一个常数值。确定在您的情况下确实如此，但通常不是-例如：a[:] = 2*a[:] 将永远循环。所以不，如果没有显式循环和比较直到值的传播完成，你想要的在 numpy 中是不可能的。
@daver：你能用循环发布一个合适的例子吗？现在循环与 a[1:]=a[1] 相同。
@daver，我不认为您了解 numpy 的工作原理。当您说 a[1:] = 2*a[:3] 时，有两个循环。第一个是 2*a[:3] ，它是一个临时数组。然后，第二个循环执行分配。 a[1:] 不是从自身分配的，而是从临时分配的。问题不在于 a[:3] 不是 a 的视图，因为它是，而是 2*a[:3] 是一个完全不同的数组。 Numpy 为您提供 C 的速度，但代价是更多的临时性。一旦你掌握了这一点，你就会明白为什么你不能在 NumPy 中轻松地做你想做的事而不深入研究底层的东西。

标签： python numpy slice

【解决方案1】：

迟到的答案，但这出现在 Google 上，所以我可能指向 OP 想要的文档。您的问题很清楚：使用 NumPy 切片时，会创建临时文件。快速调用 weave.blitz 将您的代码包装起来，以摆脱临时性并获得您想要的行为。

阅读PerformancePython tutorial 的 weave.blitz 部分了解完整详情。

【讨论】：

【解决方案2】：

accumulate 旨在做你想做的事；也就是说，沿数组分配操作。这是一个例子：

from numpy import *

a = array([1,0,0,0])
a[1:] = add.accumulate(a[0:3])
# a = [1, 1, 1, 1]

b = array([1,1,1,1])
b[1:] = multiply.accumulate(2*b[0:3])
# b = [1 2 4 8]

另一种方法是将结果数组显式指定为输入数组。这是一个例子：

c = array([2,0,0,0])
multiply(c[:3], c[:3], c[1:])
# c = [  2   4  16 256]

【讨论】：

【解决方案3】：

只需使用循环。我无法立即想到任何方法来使切片运算符按照您所说的方式运行，除了可能通过子类化 numpy 的 array 并用某种覆盖适当的方法Python voodoo ......但更重要的是，a[1:] = a[0:3] 应该将a 的第一个值复制到接下来的三个插槽中的想法对我来说似乎完全没有意义。我想这很容易让其他查看您的代码的人感到困惑（至少在前几次）。

【讨论】：

不不不，它不会将第一个值复制到接下来的三个值中，而是使用前一个条目中的数据更新每个连续条目。但是，每次更新时，我都希望它知道以前的更新。是的，我可以循环播放，但就我所想的目的而言，它非常缓慢且笨拙。但这是循环中的样子：for i in a:
*该死，对不起，我点击了标签，它更新了我的评论......代码：for i in a: a[i+1] = a[i] 关键是这是二维的数组，是我需要实现的特定数值算法。
gnibbler：这完全忽略了算法的重点，这是一个简单的例子，我正在做的是 Gauss-Seidel 迭代，它通过使用已经存在的数据来推断矩阵中的位置信息在之前的条目中被推断出来。据我了解，在 numpy 的机器深处，它执行这个循环。但是，它会循环切片中原始数据的副本。我希望它遍历切片并随时更新切片。也许我可以这样说得更清楚：a[1:] = 2*a[0:3]，预期结果是：a = [1,2,4,8]。
请修正问题。不要在 cmets 的答案中更正您的问题。请使用修改后的问题描述更新问题。
numpy 的主要优势之一是能够避免昂贵的 python 迭代。（又名广播）。这个问题是完全有效的，并且是期望 numpy 做的合理的事情。

【解决方案4】：

这不是正确的逻辑。我会尽量用字母来解释。

图像array = abcd 以a、b、c、d 作为元素。
现在，array[1:] 表示从位置为1（从0 开始）的元素开始。
在这种情况下：bcd 和array[0:3] 表示从位置0 到第三个字符（位置3-1）在这种情况下：'abc'。

这样写：
array[1:] = array[0:3]

意思是：用abc替换bcd

要获得你想要的输出，现在在 python 中，你应该使用类似的东西：

a[1:] = a[0]

【讨论】：

OP 正在尝试做一些比这更复杂的事情。目标是根据当前操作中更新的元素值更新切片的元素。很难理解，但是比传统的python循环快很多。
我刚刚阅读了他上面的评论。这个问题很不清楚。我帖子的第一部分解释了为什么它不像 OP 认为的那样工作。我正在等待问题改进以尝试解决他的问题。目前还不清楚，因为他想更新它。
是的，最初的问题含糊不清且过于简单化。在评论中找到澄清。我想知道为什么人们害怕编辑原始帖子？
我确实编辑过，现在应该很清楚了。抱歉，这是我的第一篇文章，我想知道为什么这个问题得到了 3 票反对？
我对这个问题发表了评论。我没有对此投反对票，但我认为由于说明不清楚，它被否决了。尝试提供一种替代工作方式来解决问题，或者至少详细说明您希望使用的算法。

【解决方案5】：

它必须与分配切片有关。但是，正如您可能已经知道的那样，运算符确实遵循您的预期行为：

>>> a = numpy.array([1,0,0,0])
>>> a[1:]+=a[:3]
>>> a
array([1, 1, 1, 1])

如果您在示例中的实际问题中已经有零，那么这可以解决它。否则，以额外的成本，通过乘以零或分配为零将它们设置为零，（以更快的为准）

编辑：我有另一个想法。你可能更喜欢这个：

numpy.put(a,[1,2,3],a[:3])

【讨论】：

【解决方案6】：

在执行 setkey 调用时，Numpy 必须检查目标数组是否与输入数组相同。幸运的是，有一些方法可以解决它。首先，我尝试改用numpy.put

In [46]: a = numpy.array([1,0,0,0])

In [47]: numpy.put(a,[1,2,3],a[0:3])

In [48]: a
Out[48]: array([1, 1, 1, 1])

然后从文档中，我尝试使用 flatiters (a.flat)

In [49]: a = numpy.array([1,0,0,0])

In [50]: a.flat[1:] = a[0:3]

In [51]: a
Out[51]: array([1, 1, 1, 1])

但这并不能解决您想到的问题

In [55]: a = np.array([1,0,0,0])

In [56]: a.flat[1:] = 2*a[0:3]

In [57]: a
Out[57]: array([1, 2, 0, 0])

这失败了，因为乘法是在赋值之前完成的，而不是你想要的并行。

Numpy 设计用于在数组中并行重复应用完全相同的操作。要做一些更复杂的事情，除非你能找到像 numpy.cumsum 和 numpy.cumprod 这样的函数来分解它，你必须求助于 scipy.weave 之类的东西或用 C 编写函数。（参见 PerfomancePython页面了解更多详细信息。）（另外，我从未使用过编织，所以我不能保证它会做你想要的。）

【讨论】：

【解决方案7】：

你可以看看 np.lib.stride_tricks。

这些优秀的幻灯片中有一些信息： http://mentat.za.net/numpy/numpy_advanced_slides/

stride_tricks 从幻灯片 29 开始。

虽然我对这个问题并不完全清楚，所以不能提出更具体的建议 - 尽管我可能会在 cython 或 fortran 中使用 f2py 或 weave 来做。我现在更喜欢 fortran，因为当你在 cython 中添加所有必需的类型注释时，我认为它最终看起来不如 fortran 清晰。

这里有这些方法的比较：

万维网。 scipy。 org/PerformancePython

（由于我是新用户，无法发布更多链接）举一个与您的案例相似的示例。

【讨论】：

【解决方案8】：

最后我遇到了和你一样的问题。我不得不求助于 Jacobi 迭代和编织器：

 while (iter_n < max_time_steps):
        expr = "field[1:-1, 1:-1] = (field[2:, 1:-1] "\
                                                      "+ field[:-2, 1:-1]+"\
                                                      "field[1:-1, 2:] +"\
                                                      "field[1:-1, :-2] )/4."                                       

        weave.blitz(expr, check_size=0)

         #Toroidal conditions
        field[:,0] = field[:,self.flow.n_x - 2]
        field[:,self.flow.n_x -1] = field[:,1]

        iter_n = iter_n + 1

它有效且速度快，但不是 Gauss-Seidel，因此收敛可能有点棘手。将 Gauss-Seidel 作为带索引的传统循环的唯一选择。

【讨论】：

【解决方案9】：

我建议使用 cython 而不是在 c 中循环。可能有一些花哨的 numpy 方法可以让你的示例使用很多中间步骤来工作......但是既然你已经知道如何用 c 编写它，只需将那个快速的一点点写成 cython功能，让 cython 的魔力让您轻松完成剩下的工作。

【讨论】：