0

data 是一个包含 2500 个时间序列的矩阵。我需要随着时间的推移对每个时间序列进行平均，丢弃围绕峰值记录的数据点（在间隔 tspike-dt*10...tspike+10*dt 中）。每个神经元的尖峰时间数量是可变的，并存储在一个包含 2500 个条目的字典中。我当前的代码迭代神经元和尖峰时间并将掩码值设置为 NaN。然后调用bottleneck.nanmean()。但是，此代码在当前版本中速度较慢，我想知道是否有更快的解决方案。谢谢！

import bottleneck
import numpy as np
from numpy.random import rand, randint

t = 1
dt = 1e-4
N = 2500
dtbin = 10*dt

data = np.float32(ones((N, t/dt)))
times = np.arange(0,t,dt)
spiketimes = dict.fromkeys(np.arange(N))
for key in spiketimes:
  spiketimes[key] = rand(randint(100))

means = np.empty(N)

for i in range(N):        
  spike_times = spiketimes[i]
  datarow = data[i]
  if len(spike_times) > 0:
    for spike_time in spike_times:                        
      start=max(spike_time-dtbin,0)
      end=min(spike_time+dtbin,t)
      idx = np.all([times>=start,times<=end],0)
      datarow[idx] = np.NaN
  means[i] = bottleneck.nanmean(datarow)
4

2 回答 2

0

而不是使用nanmean你可以只索引你需要和使用的值mean。

means[i] = data[ (times<start) | (times>end) ].mean()

如果我误解了并且您确实需要索引，您可以尝试

means[i] = data[numpy.logical_not( np.all([times>=start,times<=end],0) )].mean()

同样在您可能不想使用的代码中if len(spike_times) > 0（我假设您在每次迭代中删除尖峰时间，否则该语句将始终为真并且您将有一个无限循环），仅使用for spike_time in spike_times.

于 2012-08-03T19:28:21.020 回答
0

代码中的绝大多数处理时间都来自这一行：

idx = np.all([times>=start,times<=end],0)

这是因为对于每个尖峰，您都在将每个值与开始和结束时间进行比较。由于您在此示例中具有统一的时间步长（我认为这在您的数据中也是如此），因此简单地计算开始和结束索引要快得多：

# This replaces the last loop in your example:
for i in range(N):        
    spike_times = spiketimes[i]
    datarow = data[i]
    if len(spike_times) > 0:
        for spike_time in spike_times:
            start=max(spike_time-dtbin,0)
            end=min(spike_time+dtbin,t)
            #idx = np.all([times>=start,times<=end],0)
            #datarow[idx] = np.NaN
            datarow[int(start/dt):int(end/dt)] = np.NaN
    ## replaced this with equivalent for testing
    means[i] = datarow[~np.isnan(datarow)].mean()

这将我的运行时间从 ~100s 减少到 ~1.5s。您还可以通过在spike_times 上对循环进行矢量化来节省更多时间。这样做的效果将取决于您的数据的特征（对于高峰值率应该是最有效的）：

kernel = np.ones(20, dtype=bool)
for i in range(N):        
    spike_times = spiketimes[i]
    datarow = data[i]
    mask = np.zeros(len(datarow), dtype=bool)
    indexes = (spike_times / dt).astype(int)
    mask[indexes] = True  
    mask = np.convolve(mask, kernel)[10:-9]

    means[i] = datarow[~mask].mean()
于 2012-08-04T20:52:35.537 回答