代码之家 › 专栏 › 技术社区 › jkr

用体素洗牌四维时间序列

indexing numpy python

jkr · 技术社区 · 7 年前

我有一个四维数组,它是一个三维数组的时间序列。我想沿着时间轴洗牌三维数组中的每个点。下面是我使用嵌套 for

import numpy as np

timepoints = 2
x = 4
y = 4
z = 3

vol_1 = np.zeros((x, y, z))
vol_2 = np.ones((x, y, z))
timeseries = np.array((vol_1, vol_2))

timeseries.shape  # (2, 4, 4, 3)

# One voxel over time.
timeseries[:, 0, 0, 0]

for xx in range(x):
    for yy in range(y):
        for zz in range(z):
            np.random.shuffle(timeseries[:, xx, yy, zz])

2 回复 | 直到 7 年前

Divakar 7 年前

我们可以沿着第一个轴生成所有无序索引,然后简单地使用 advanced-indexing 获得随机版本。现在,为了得到所有的乱序索引,我们可以生成一个与输入数组形状相同的随机数组,并沿着第一个轴得到argsort索引。这一点以前已经探讨过,因为 here

因此,我们将有一个这样的矢量化实现-

m,n,r,p = a.shape # a is the input array
idx = np.random.rand(*a.shape).argsort(0)
out = a[idx, np.arange(n)[:,None,None], np.arange(r)[:,None], np.arange(p)]

1) 输入4D阵列:

In [711]: a
Out[711]: 
array([[[[60, 22, 34],
         [29, 18, 79]],

        [[11, 69, 41],
         [75, 30, 30]]],


       [[[63, 61, 42],
         [70, 56, 57]],

        [[70, 98, 71],
         [29, 93, 96]]]])

2) 使用拟议方法沿第一轴索引生成的随机索引:

In [712]: idx
Out[712]: 
array([[[[1, 0, 1],
         [0, 1, 1]],

        [[0, 0, 1],
         [1, 0, 1]]],


       [[[0, 1, 0],
         [1, 0, 0]],

        [[1, 1, 0],
         [0, 1, 0]]]])

In [713]: out
Out[713]: 
array([[[[63, 22, 42],
         [29, 56, 57]],

        [[11, 69, 71],
         [29, 30, 96]]],


       [[[60, 61, 34],
         [70, 18, 79]],

        [[70, 98, 41],
         [75, 93, 30]]]])

仔细观察,我们会发现 63 在 a[0,0,0,0] 60 在 a[1,0,0,0] 由于 idx 1 和 0 分别位于 idx公司 . 接下来, 22 和 61 留在他们的地方,因为值为和等等

运行时测试

In [726]: timeseries = np.random.rand(10,10,10,10)

In [727]: %timeit org_app(timeseries)
100 loops, best of 3: 5.24 ms per loop

In [728]: %timeit proposed_app(timeseries)
1000 loops, best of 3: 289 Âµs per loop

In [729]: timeseries = np.random.rand(50,50,50,50)

In [730]: %timeit org_app(timeseries)
1 loop, best of 3: 720 ms per loop

In [731]: %timeit proposed_app(timeseries)
1 loop, best of 3: 426 ms per loop

在大尺寸情况下,创建随机阵列的成本被证明是该方法的瓶颈,但仍显示出比原始循环版本有很好的加速。

norok2 7 年前

我添加这个作为一个答案,因为它不适合评论,因为它只是@Divakar的优秀答案之上的一个小添加:

def divakar(a):
    m,n,r,p = a.shape # a is the input array
    idx = np.random.rand(*a.shape).argsort(0)
    return a[idx, np.arange(n)[:,None,None], np.arange(r)[:,None], np.arange(p)]

a = np.random.rand(50,50,50,50)
%timeit divakar(a)
560 ms Â± 2.62 ms per loop (mean Â± std. dev. of 7 runs, 1 loop each)

def norok2(a):
    shape = a.shape
    idx = np.random.rand(*a.shape).argsort(0).reshape(shape[0], -1)
    return a.reshape(shape[0], -1)[idx, np.arange(shape[1] * shape[2] * shape[3])].reshape(shape)

a = np.random.rand(50,50,50,50)
%timeit norok2(a)
495 ms Â± 1.61 ms per loop (mean Â± std. dev. of 7 runs, 1 loop each)

与OP的提案相比:

def jakub(a):
    t, x, y, z = a.shape
    for xx in range(x):
        for yy in range(y):
            for zz in range(z):
                np.random.shuffle(a[:, xx, yy, zz])


%timeit jakub(a)
2 s Â± 30.8 ms per loop (mean Â± std. dev. of 7 runs, 1 loop each)

顺便提一下,我提出的修改更容易扩展到n维阵列和任意洗牌轴,例如:

import numpy as np
import functools

def shuffle_axis(arr, axis=0):
    arr = np.swapaxes(arr, 0, axis)
    shape = arr.shape
    i = np.random.rand(*shape).argsort(0).reshape(shape[0], -1)
    return arr.reshape(shape[0], -1)[i, np.arange(functools.reduce(lambda x, y: x * y, shape[1:]))].reshape(shape).swapaxes(axis, 0)

a = np.random.rand(50,50,50,50)
%timeit shuffle_axis(a)
499 ms Â± 2.1 ms per loop (mean Â± std. dev. of 7 runs, 1 loop each)

编辑再次访问

...时间安排并不比随机化更糟糕:

a = np.random.rand(50,50,50,50)
%timeit np.random.shuffle(a.ravel())
310 ms Â± 1.84 ms per loop (mean Â± std. dev. of 7 runs, 1 loop each)

这应该是该问题任何解决方案性能的某种下限(但它不能解决OP问题)。