3

我对一些 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 中循环。

4

9 回答 9

4

Late answer, but this turned up on Google so I probably point to the doc the OP wanted. Your problem is clear: when using NumPy slices, temporaries are created. Wrap your code in a quick call to weave.blitz to get rid of the temporaries and have the behaviour your want.

Read the weave.blitz section of PerformancePython tutorial for full details.

于 2010-04-10T13:59:41.850 回答
2

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

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]
于 2009-10-19T21:13:02.000 回答
1

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

>>> 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])
于 2009-10-19T13:03:14.137 回答
1

这不是正确的逻辑。我会尝试用字母来解释它。

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

写类似：
array[1:] = array[0:3]

意思是：替换bcd为abc

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

a[1:] = a[0]
于 2009-10-19T13:37:10.920 回答
1

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

于 2009-10-19T07:49:35.740 回答
1

最后我想出了和你一样的问题。我不得不求助于 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 作为带索引的传统循环的唯一选择。

于 2013-02-11T02:21:10.013 回答
1

Numpy 必须在执行 setkey 调用时检查目标数组是否与输入数组相同。幸运的是，有一些方法可以解决它。首先，我尝试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.cumsumand这样的函数来分解它，否则你numpy.cumprod将不得不求助于 scipy.weave 之类的东西或用 C 编写函数。（有关更多详细信息，请参见PerfomancePython页面。）（也，我从来没有使用过编织，所以我不能保证它会做你想要的。）

于 2009-10-19T16:04:55.910 回答
1

你可以看看 np.lib.stride_tricks。

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

stride_tricks 从幻灯片 29 开始。

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

这里有这些方法的比较：

万维网。scipy。组织/性能Python

（不能发布更多链接，因为我是新用户）有一个与您的案例相似的示例。

于 2009-10-28T14:41:47.150 回答
0

我建议使用 cython 而不是在 c 中循环。可能有一些奇特的 numpy 方法可以让您的示例使用大量中间步骤来工作......但是既然您已经知道如何用 c 编写它，只需将那个快速的一点点编写为 cython 函数，然后让 cython 的魔法使其余的工作对您来说很容易。

于 2009-10-21T20:54:35.127 回答