如何正确使用numpy中的as_strided(来自np.stride_tricks)?

5
我正在尝试使用numpy.strided_tricks来重塑numpy数组。我正在遵循以下指南:https://dev59.com/pHE95IYBdhLWcg3wCJbk#2487551 我的用例非常相似,不同之处在于我需要步幅为3。
给定这个数组:
a = np.array([1, 2, 3, 4, 5, 6, 7, 8, 9])

我希望能获得:

array([[1, 2, 3],
       [2, 3, 4],
       [3, 4, 5],
       [4, 5, 6],
       [5, 6, 7],
       [6, 7, 8],
       [7, 8, 9]])

以下是我尝试过的内容:

import numpy as np

as_strided = np.lib.stride_tricks.as_strided
a = np.arange(1, 10)

as_strided(a, (len(a) - 2, 3), (3, 3))

array([[                 1,      2199023255552,             131072],
       [     2199023255552,             131072, 216172782113783808],
       [            131072, 216172782113783808,        12884901888],
       [216172782113783808,        12884901888,                768],
       [       12884901888,                768,   1125899906842624],
       [               768,   1125899906842624,           67108864],
       [  1125899906842624,           67108864,                  4]])

我非常确定我完全按照示例操作,但显然并不是。我哪里出了问题?


你为什么认为需要步长为3? - user2357112
@user2357112 我不知道... 根据给出的示例,我认为那就是我需要跨越的方式。 - cs95
看起来这个例子为4字节整数硬编码了步幅为4 - 考虑到它们的输入在不同的操作系统上很容易变成8字节,这不是一个好主意。我要进行编辑。 - user2357112
1
as_strided 允许您访问数组数据缓冲区外的字节。 它不检查步幅和形状是否有效。 使用时请小心。 - hpaulj
3个回答

8

接受的答案(和讨论)非常好,但是为了不想运行自己的测试用例的读者着想,我将尝试说明正在发生的事情:

In [374]: a = np.arange(1,10)
In [375]: as_strided = np.lib.stride_tricks.as_strided

In [376]: a.shape
Out[376]: (9,)
In [377]: a.strides 
Out[377]: (4,)

对于一个连续的一维数组,strides是元素大小,这里为4个字节,即int32。要从一个元素到下一个元素,它向前移动4个字节。

OP尝试的方法:

In [380]: as_strided(a, shape=(7,3), strides=(3,3))
Out[380]: 
array([[        1,       512,    196608],
       [      512,    196608,  67108864],
       [   196608,  67108864,         4],
       [ 67108864,         4,      1280],
       [        4,      1280,    393216],
       [     1280,    393216, 117440512],
       [   393216, 117440512,         7]])

这是按照3个字节步长进行移动,跨越int32边界,并得到大多数不可理解的数字。如果数据类型为bytes或uint8,则可能更容易理解。

相反,使用a.strides*2(元组复制)或(4,4),我们得到所需的数组:

In [381]: as_strided(a, shape=(7,3), strides=(4,4))
Out[381]: 
array([[1, 2, 3],
       [2, 3, 4],
       [3, 4, 5],
       [4, 5, 6],
       [5, 6, 7],
       [6, 7, 8],
       [7, 8, 9]])

列和行都向前移动一个元素,导致产生一个1步长的移动窗口。我们也可以设置shape=(3,7),这将得到3个窗口,每个窗口有7个元素。

In [382]: _.strides
Out[382]: (4, 4)

将步幅更改为(8,4),每个窗口选择2个元素。

In [383]: as_strided(a, shape=(7,3), strides=(8,4))
Out[383]: 
array([[          1,           2,           3],
       [          3,           4,           5],
       [          5,           6,           7],
       [          7,           8,           9],
       [          9,          25, -1316948568],
       [-1316948568,   184787224, -1420192452],
       [-1420192452,           0,           0]])

但是形状不对,向我们展示了原始数据缓冲区末尾的字节。 这可能很危险(我们不知道这些字节是否属于其他对象或数组)。 对于这种大小的数组,我们无法获得完整的2步窗口集。

现在,每行的步骤3元素(3 * 4,4):

In [384]: as_strided(a, shape=(3,3), strides=(12,4))
Out[384]: 
array([[1, 2, 3],
       [4, 5, 6],
       [7, 8, 9]])
In [385]: a.reshape(3,3).strides
Out[385]: (12, 4)

这与3x3变形相同,包括步长。

我们可以设置负步长和0值。 实际上,在具有正步幅的维度上进行负步幅切片会给出负步幅,并且广播通过设置0步幅来实现:

In [399]: np.broadcast_to(a, (2,9))
Out[399]: 
array([[1, 2, 3, 4, 5, 6, 7, 8, 9],
       [1, 2, 3, 4, 5, 6, 7, 8, 9]])
In [400]: _.strides
Out[400]: (0, 4)

In [401]: a.reshape(3,3)[::-1,:]
Out[401]: 
array([[7, 8, 9],
       [4, 5, 6],
       [1, 2, 3]])
In [402]: _.strides
Out[402]: (-12, 4)

然而,负的步幅需要调整原始数组中哪个元素是视图的第一个元素,而as_strided没有此参数。


对于一维数组,步长是元素的大小 - 仅适用于连续的数组。通常情况下,假设连续性是个坏主意,除非你确切知道输入是如何产生的或已经检查了连续性。 - user2357112
太棒了,非常感谢。真的很有启发和帮助。通过所提供的例子,更容易辨别出一个模式。 - cs95
关于负步长,直接通过ndarray构造函数可能更容易,这样您可以提供offset以及shapestrides - user2357112

2

我不知道你为什么认为需要步长为3。实际上,你需要的是一种以字节为单位的步长,它表示一个元素与下一个元素之间的距离。你可以使用a.strides来获取这个步长:

as_strided(a, (len(a) - 2, 3), a.strides*2)

谢谢,这很简单。我其实不知道strides里面放了什么,我以为是3,因为在链接中他们用了4。我看了文档但是理解得不是很多。 - cs95
它表示“在遍历数组时,每个维度要跨越的字节元组。” @user2357112 - 能否添加一些解释? - Vivek Kalyanarangan
出于好奇,如果我想要跨越两个 ([1, 2, 3], [3, 4, 5], ...),我需要 a.strides * 3 吗? - cs95
@cᴏʟᴅsᴘᴇᴇᴅ:我认为你误解了数组的步幅。一个数组的步幅告诉你在任何维度上从一个数组元素移动到下一个元素需要在内存中跨越多少字节。请参阅 https://docs.scipy.org/doc/numpy-1.13.0/reference/generated/numpy.ndarray.strides.html#numpy.ndarray.strides 和 https://docs.scipy.org/doc/numpy-1.13.0/reference/internals.html#internal-organization-of-numpy-arrays。 - user2357112
感谢您的帮助。我已经弄清楚了:as_strided(a, (len(a) - 5, 3), (a.strides[0] *2 , a.strides[0] ))。这个函数让我想起了卷积中的步幅,因为我习惯于使用它,所以我能够很快地掌握它。 - cs95
1
@cᴏʟᴅsᴘᴇᴇᴅ:len(a) - 5 看起来不太对,正确的表达式应该是带有// 2的某个表达式。除此之外,看起来你已经掌握了这个东西。 - user2357112

2
我试图执行类似的操作并遇到了同样的问题。
根据此评论,您的问题是:
  1. 在存储在内存中的元素大小时,您没有考虑到您的元素大小(int32 = 4,可以使用a.dtype.itemsize检查);
  2. 您没有适当地指定要跳过的步幅数量,在您的情况下也为4,因为您只跳过一个元素。
我基于这个答案制作了一个函数,其中我计算给定数组的分段,使用n元素的窗口并指定要重叠的元素数量(由window-number_of_elements_to_skip给出)。
如果有人需要,我在此分享它,因为我花了一些时间才弄清楚stride_tricks的工作原理。
def window_signal(signal, window, overlap):
    """ 
    Windowing function for data segmentation.

    Parameters:
    ------------
    signal: ndarray
            The signal to segment.
    window: int
            Window length, in samples.
    overlap: int
             Number of samples to overlap

    Returns: 
    --------
    nd-array 
            A copy of the signal array with shape (rows, window),
            where row = (N-window)//(window-overlap) + 1
    """
    N = signal.reshape(-1).shape[0] 
    if (window == overlap):
        rows = N//window
        overlap = 0
    else:
        rows = (N-window)//(window-overlap) + 1
        miss = (N-window)%(window-overlap)
        if(miss != 0):
            print('Windowing led to the loss of ', miss, ' samples.')
    item_size = signal.dtype.itemsize 
    strides = (window - overlap) * item_size
    return np.lib.stride_tricks.as_strided(signal, shape=(rows, window),
                                           strides=(strides, item_size))

这种情况的解决方案根据您的代码是: as_strided(a, (len(a) - 2, 3), (4, 4)) 或者,可以使用window_signal函数: window_signal(a, 3, 2) 两者都返回以下数组作为输出:
array([[1, 2, 3],
   [2, 3, 4],
   [3, 4, 5],
   [4, 5, 6],
   [5, 6, 7],
   [6, 7, 8],
   [7, 8, 9]])

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接