我喜欢最小二乘匹配具有许多已知信号形状的数据(浮点数的数组)。我的代码可以运行,但是对于我计划执行的许多运行来说太慢了:
import numpy
import time
samples = 50000
width_signal = 100
data = numpy.random.normal(0, 1, samples)
signal = numpy.random.normal(0, 1, width_signal) # Placeholder
t0 = time.clock()
for i in range(samples - width_signal):
data_chunk = data[i:i + width_signal]
residuals = data_chunk - signal
squared_residuals = residuals**2
summed_residuals = numpy.sum(squared_residuals)
t1 = time.clock()
print('Time elapsed (sec)', t1-t0)
编辑:更正了一个错误:首先平方残差,然后求和。
在我的计算机上运行大约需要0.2秒。由于我有许多数据集和信号形状,所以这太慢了。我的特定问题不允许使用典型的MCMC方法,因为信号形状差异太大。它必须是蛮力的。
典型的数据量为50,000浮点数,信号为100。这些可以相差几倍。
我的测试表明:
numpy.sum(residuals)
的总和占90%的时间。我尝试了Python的sum(residuals)
,对于较小的数组(〜<50个元素),速度更快,而对于较大的数组,它的速度较慢。我应该插入if
条件吗? numpy.roll()
而不是直接获取数据,因此.roll()
速度较慢。 问题:
参考方案
基于 Compute mean squared, absolute deviation and custom similarity measure - Python/NumPy
中提出的各种方法,我们希望在这里解决问题。
方法#1
我们可以利用基于 np.lib.stride_tricks.as_strided
的 scikit-image's view_as_windows
来获取滑动窗口,从而在这里有了我们的第一个解决方案,就像这样-
from skimage.util import view_as_windows
d = view_as_windows(data,(width_signal))-signal # diffs
out = np.einsum('ij,ij->i',d,d)
More info on use of as_strided
based view_as_windows
.
方法2
再次根据该答案文章中的矩阵乘法技巧,我们可以改善性能,如下所示:
def MSD_strided(data, signal):
w = view_as_windows(data,(width_signal))
return (w**2).sum(1) + (signal**2).sum(0) - 2*w.dot(signal)
方法#3
我们将通过统一过滤和卷积来改进方法2-
from scipy.ndimage.filters import uniform_filter
def MSD_uniffilt_conv(data, signal):
hW = width_signal//2
l = len(data)-len(signal)+1
parte1 = uniform_filter(data**2,width_signal)[hW:hW+l]*width_signal
parte3 = np.convolve(data, signal[::-1],'valid')
return parte1 + (signal**2).sum(0) - 2*parte3
基准化
发布样本的时间-
In [117]: %%timeit
...: for i in range(samples - width_signal + 1):
...: data_chunk = data[i:i + width_signal]
...: residuals = data_chunk - signal
...: squared_residuals = residuals**2
...: summed_residuals = numpy.sum(squared_residuals)
1 loop, best of 3: 239 ms per loop
In [118]: %%timeit
...: d = view_as_windows(data,(width_signal))-signal
...: np.einsum('ij,ij->i',d,d)
100 loops, best of 3: 11.1 ms per loop
In [209]: %timeit MSD_strided(data, signal)
10 loops, best of 3: 18.4 ms per loop
In [210]: %timeit MSD_uniffilt_conv(data, signal)
1000 loops, best of 3: 1.71 ms per loop
~140x
加快了第三个速度!
我知道这似乎是一个荒谬的问题,但是我必须在与部门中其他人共享的计算服务器上定期运行作业,当我开始10个作业时,我真的希望它只占用10个核心而不是更多;我不在乎每次运行一个内核所需的时间是否更长:我只是不想让它侵犯其他人的领土,这将需要我放弃工作等等。我只想拥有10个核心,仅此而已。更具体地说,我在基于Python 2.7.3和numpy 1.6.1的Redh…
numpy.savetxt“元组索引超出范围”? - python我试图在文本文件中写几行,这是我使用的代码:import numpy as np # Generate some test data data = np.arange(0.0,1000.0,50.0) with file('test.txt', 'w') as outfile: outfile.write('…
Python / Scipy过滤器离散化 - python我目前正在尝试从Matlab转向Python,并在多个方面取得了成功。但是,我经常使用的Matlab信号处理工具箱中的一个函数是impinvar函数,用于从其模拟版本计算数字滤波器。在Scipy.signal中,我仅发现bilinear函数可以执行类似的操作。但是,与Matlab bilinear function相比,它不需要可选参数来对频率进行一些预变形…
Python GPU资源利用 - python我有一个Python脚本在某些深度学习模型上运行推理。有什么办法可以找出GPU资源的利用率水平?例如,使用着色器,float16乘法器等。我似乎在网上找不到太多有关这些GPU资源的文档。谢谢! 参考方案 您可以尝试在像Renderdoc这样的GPU分析器中运行pyxthon应用程序。它将分析您的跑步情况。您将能够获得有关已使用资源,已用缓冲区,不同渲染状态上…
Python:图像处理可产生皱纹纸效果 - python也许很难描述我的问题。我正在寻找Python中的算法,以在带有某些文本的白色图像上创建皱纹纸效果。我的第一个尝试是在带有文字的图像上添加一些真实的皱纹纸图像(具有透明度)。看起来不错,但副作用是文本没有真正起皱。所以我正在寻找更好的解决方案,有什么想法吗?谢谢 参考方案 除了使用透明性之外,假设您有两张相同尺寸的图像,一张在皱纹纸上明亮,一张在白色背景上有深…