scipy中重采样数组的索引答案

【问题标题】：Indices of resampled array in scipyscipy中重采样数组的索引
【发布时间】：2013-01-30 22:47:21
【问题描述】：

我有两个长度相同的一维数组，包含一个时间序列和值序列，例如

t = linspace(0, 5, 5) # [0, 1.25, 2.5, 3.75, 5]
x = array(range(10, 25)) # [10, 11, 12, 13, 14]

例如，我必须使用不同的采样时间点重新采样 x 数组（具有相同的起点和终点，但可以有任意数量的元素）

r = linspace(0, 5, 4) # [ 0, 1.667, 3.333, 5]
x2 = resample(t, x, r) # [10, 11, 12, 14]

即r的每一个时间点都放在t的两个时间点之间，我想求出这两个时间点在t中的下一点的索引。从索引数组中，可以得到x的相对点。

我想要一个基于向量的解决方案，没有循环，可能使用 scipy 的运算符。如果使用 scipy 的函数会更好。

编辑：这是我需要的代码，但更短、更快且基于矢量的解决方案会更好。我找不到一个（直到尝试）。

def resample(t, r):
    i, j, k = 0, 1, 0
    s = []
    while j < len(t):
        if t[i] <= r[k] < t[j]:
            s.append(i)
            k += 1
        else:
            i += 1
            j += 1
    s.append(len(t) - 1)
    return array(s)

【问题讨论】：

标签： python numpy python-3.x scipy

【解决方案1】：

以下两个小函数中的第二个可以完成你想要的：

def resample_up(t, x, r) :
    return x[np.argmax(r[:, None] <= t, axis=1)]

def resample_down(t, x, r) :
    return x[::-1][np.argmax(r[:, None] >= t[::-1], axis=1)]

>>> resample_up(t, x, r)
array([10, 12, 13, 14])
>>> resample_down(t, x, r)
array([10, 11, 12, 14])

如果您发现很难弄清楚发生了什么，以下可能会有所帮助：

>>> r[:, None] <= t
array([[ True,  True,  True,  True,  True],
       [False, False,  True,  True,  True],
       [False, False, False,  True,  True],
       [False, False, False, False,  True]], dtype=bool)
>>> r[:, None] >= t[::-1]
array([[False, False, False, False,  True],
       [False, False, False,  True,  True],
       [False, False,  True,  True,  True],
       [ True,  True,  True,  True,  True]], dtype=bool)

然后np.argmax 返回每行中第一次出现True 的索引。

编辑很难让它比一行代码短，但是对于大型数组，性能会受到影响，因为索引查找永远不会在循环的早期中断。因此，对于非常大的数组，使用 python 循环扫描数组可能会更快。对于较小的则没有：

In [2]: %timeit resample_up(t, x, r)
100000 loops, best of 3: 7.32 us per loop

In [3]: %timeit resample_down(t, x, r)
100000 loops, best of 3: 8.44 us per loop

In [4]: %timeit resample(t, x, r) # modified version of the OP's taking also x
100000 loops, best of 3: 13.7 us per loop

【讨论】：

很好，它使用我不知道的语法来索引数组！

【解决方案2】：

您可以尝试在scipy.interpolate 中使用interp1d 函数，将kind 参数指定为zero。使用你的数组：

>>> from scipy.interpolate import interp1d
>>> f = interp1d(t,x,kind="zero")
>>> f(r)
array((10, 11, 12, 13))

请注意，“重新采样”数组中的最后一个元素是 13，而不是您在问题中要求的 14，而是f(5.001) = 14 (*)。只要“重采样”数组与原始数组中的一个点匹配，插值函数就是不连续的。

(*) 如果要在t 范围之外重新采样，则需要在interp1d 调用中指定关键字参数bounds_error=False。

【讨论】：

很好，感谢您添加的详细信息，但最后一点必须一致。当然，可以人为地引入另一个，但我相信这可能会导致抽样中的（非常小的）错误。

【解决方案3】：

numpy.interp 是一个快速简单的分段线性插值器：

from __future__ import division
import numpy as np

t = np.linspace(0, 5, 5)  # [0, 1.25, 2.5, 3.75, 5]
x = np.array(range(10, 15))  # [10, 11, 12, 13, 14]
r = np.linspace(0, 5, 4)  # [ 0, 1.667, 3.333, 5]

print "np.interp:", np.interp( r, t, x )
    # [ 10.    11.33  12.67  14.  ]
xint = np.arange( len(t) )
print "r to int:", np.interp( r, t, xint ).astype(int)
    # [0 1 2 4]

【讨论】：