对不同大小的顺序数组切片进行操作

Question

假设我有以下signal 数组，其中每个值对应于time 数组中的一个时间：

np.random.seed(123)
time = np.array([0,2,4,7,10,11,12,17,21,25,29,30,31,40])  # in seconds
signal = np.random.randint(1,5,len(time))

我想要做的是将signal 数组分割成更小的数组，这样每个切片的时间跨度至少 10 秒。然后对每个切片求和 signal。视觉上：

         |-----sum------||-----sum------||---sum----||--X--
time   = 0,  2,  4,  7, 10, 11, 12, 17, 21, 25, 29, 31, 40
signal = 3,  2,  3,  3,  1,  3,  3,  2,  4,  3,  4,  3,  2

我想要的输出是一个列表，其中包含每 10 秒切片的 signal 的总和：

[12,  # 3+2+3+3+1
 13,  # 1+3+3+2+4
 14]  # 4+3+4+3

请注意，最后 2 个 signal 元素不会被求和，因为那里的时间差小于 10 秒

我写了以下函数：

def count(x, time, epoch=60):
    # calculate time diff
    time = time - time[0]

    # get indices at time boundaries
    num_bins = int(max(time) / epoch)
    inds = [0]

    for i in range(num_bins):
        upper_ind = np.argmax(time >= time[inds[-1]] + epoch)

        if time[upper_ind] - time[inds[-1]] >= epoch:
            inds.append(upper_ind)

    # calculate sums between each boundary
    counts = []
    for i in range(len(inds) - 1):
        lower = inds[i]
        upper = inds[i+1] + 1

        cur_signal = x[lower:upper]

        counts.append(sum(cur_signal))

    return counts

调用者：

counts = count(signal, time, epoch=10)

它可以工作，但是对于大型数组来说它很慢并且很hacky。有没有更有效的方法来做到这一点，也许有一些 numpy 魔法，我不必通过来确定边界，然后再通过来得到总和？

如果有一种方法可以在 2 个时间点之间进行线性插值（即，如果前一个时间点略短于 10 秒，而下一个时间点略大于 10 秒），则可以通过在精确 10 秒处估计 signal 的值来获得奖励积分间隔

编辑：

只是为了从 cmets 那里提供一些额外的信息...

至少 10 秒长意味着切片不能短于 10 秒，但可以更长。我将取大于 10s 的第一个时间点。请参阅上面示例中的第二个切片

边界处的信号值将被计算两次。也就是说，一个切片的结束值就是下一个切片的起始值

Answer 1

编辑

在考虑了一会儿之后，我意识到最好的选择可能不是优雅的 numpy 代码，尤其是在您关心性能的情况下。即使是 @PaulPanzer 的代码，虽然很漂亮，但也依赖于在循环中调用 searchsorted（基于相对昂贵的二进制搜索）。

相反，您可以在一次循环中执行整个算法，无需搜索：

signal = np.array([3,  2,  3,  3,  1,  3,  3,  2,  4,  3,  4,  3,  2])
time = np.array([0,  2,  4,  7, 10, 11, 12, 17, 21, 25, 29, 31, 40])

def count(signal, time, epoch=10):
    counts = []
    total = 0
    timestart = times[0]

    for x,t in zip(signal, time):
        total += x

        if t - timestart >= epoch:
            counts.append(total)
            total = x
            timestart = t

    return counts

count(signal, time)

输出：

[12, 13, 14]

时间

看起来简单循环确实比 numpy/searchsorted/where 方法快了很多。

我的代码：

%%timeit        
count(signal, time)

5.88 µs ± 165 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

@PaulPanzer 的代码：

%%timeit
idx = np.fromiter(iter(accumulate(chain((0,), repeat(10)), lambda now, delta: time.searchsorted(time[now] + delta)).__next__, len(time)), int)
np.add.reduceat(signal[:idx[-1]], idx[:-1]) + signal[idx[1:]]

9.63 µs ± 182 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

@Brenlla 的代码：

%%timeit
out=[]
prev=0
# need to reinitialize the time array since the loop eats it
time = np.array([0,  2,  4,  7, 10, 11, 12, 17, 21, 25, 29, 31, 40])
while True:
    try:
        idx10 = np.where(time >=10)[0][0]
        time-=time[idx10]
        out.append(np.sum(signal[prev:idx10+1]))
        prev=idx10
    except:
        break

30.1 µs ± 502 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Answer 2

这是一种结合使用itertools 和numpy 的方法：

>>> time   = 0,  2,  4,  7, 10, 11, 12, 17, 21, 25, 29, 31, 40
>>> signal = 3,  2,  3,  3,  1,  3,  3,  2,  4,  3,  4,  3,  2
>>> time, signal = map(np.array, (time, signal))
>>> 
>>> idx = np.fromiter(iter(accumulate(chain((0,), repeat(10)), lambda now, delta: time.searchsorted(time[now] + delta)).__next__, len(time)), int)
>>> np.add.reduceat(signal[:idx[-1]], idx[:-1]) + signal[idx[1:]]
array([12, 13, 14])

Answer 3

这是一种矢量化和 Numpythonic 方法：

# time is array([ 2,  2,  4,  7, 10, 11, 12, 17, 21, 25, 29, 30, 31, 40])

# Using broadcasting you can get a 2d array of the difference of all items
# from other items within your array
In [115]: arr = time[:, None] - time
# Then find indices where the difference is less and equal to -10
In [116]: x, y = np.where(arr <= -10)
# find the first occurrences of where for each item the difference is less and equal to -10 
In [117]: first_acc = np.concatenate(([0], np.where(np.diff(x) != 0)[0]  + 1, [x.size]))

# use a recursive generator function to retrieve all the expected indices.
In [118]: def get_ind_rec(ind=0):
     ...:     try:
     ...:         ind = y[first_acc[ind]]
     ...:         yield ind
     ...:         yield from get_ind_rec(ind)
     ...:     except: # IndexError
     ...:         pass
     ...:     
     ...:     

In [119]: list(get_ind_rec())
Out[119]: [6, 9, 13]

现在您可以简单地使用np.split() 根据这些索引拆分signal，并使用map 将sum 应用于所有切片。

Answer 4

也有点老套，但我觉得很容易理解。可能需要用更健壮/优雅的东西替换try ... except

time   = 0,  2,  4,  7, 10, 11, 12, 17, 21, 25, 29, 31, 40
signal = 3,  2,  3,  3,  1,  3,  3,  2,  4,  3,  4,  3,  2
time, signal = map(np.array, (time, signal))

out=[]
prev=0
while True:
    try:
        idx10 = np.where(time >=10)[0][0]
        time-=time[idx10]
        out.append(sum(signal[prev:idx10+1]))
        prev=idx10
    except:
        break

Answer 5

如果编译器是一个选项

如果使用 numpy 和广播解决问题并非易事，这应该始终是一种替代方法。即使矢量化很容易解决问题，您通常也可以获得显着的加速。

随心所欲地循环

import numpy as np
import numba as nb

@nb.njit(fastmath=True)
def count(x, time, epoch=10):
  max_bins=int((time[-1]-time[0]))//epoch
  sum_arr=np.zeros((max_bins),dtype=x.dtype)

  start_time=time[0]
  ii=0
  for i in range(x.shape[0]):
    if (time[i]-start_time) < epoch:
      sum_arr[ii]+=x[i]
    else:
      sum_arr[ii]+=x[i]
      ii+=1
      sum_arr[ii]+=x[i]
      start_time=time[i]

  return sum_arr[0:ii]

编译

在此示例中，我使用 numba，因为它很简单。只需要导入和函数装饰器即可获得一些数量级的加速。

衡量绩效

#create some data
t=np.arange(0,1e6,2)
signal = np.random.randint(1,5,len(t))
sum_arr=count(signal, t, epoch=10)

t1=time.time()
sum_arr_1=your_count(signal, t, epoch=10)
print(time.time()-t1)

#The first call gets about 0.2s compilation overhead
sum_arr_2=count(signal, t, epoch=10)

t1=time.time()
for i in range(1000):
  sum_arr_2=count(signal, t, epoch=10)

print((time.time()-t1)/1000)
np.allclose(sum_arr_1,sum_arr_2)

结果

your_version:13.6s
compiled_version: 0.6ms
np.allclose: True

总之，加速了 20200 倍。