循环似乎比NumPy/SciPy 3D插值快
我对NumPy和SciPy的插值方法感到有些困惑。我用LinearNDInterpolator实现了3D线性插值，但发现速度非常慢。于是我用纯Python写了一个暴力的三重循环方法，结果竟然快了1000倍。我也尝试了Numba这个包，但速度并没有变快。
根据我在网上找到的资料，Python的循环应该比NumPy、SciPy和Numba慢得多。但我看到的情况并不是这样。
我把我运行的完整源代码贴出来了。在我的机器上，我得到了这些时间：
Numpy ready: 3.94499993324 s, result[0]= 0.480961746817
Python for loop...
Python ready: 0.0299999713898 s, result[0]= 0.480961746817
Numba for loop...
Numba 0 ready: 0.223000049591 s, result[0]= 0.480961746817
Numba for loop...
Numba 1 ready: 0.0360000133514 s, result[0]= 0.480961746817
我使用的是Anaconda Python 2.7。我这里缺少了什么呢？
import numpy
import scipy.interpolate
import time
from numba import jit
# x: a (40,) numpy array of ordered ints
# y: a (30,) numpy array of ordered ints
# z: a (10,) numpy array of ordered ints
# values: a (10,30,40) numpy array of floats
# targetxs: a (NP,) numpy array of random floats
# targetys: a (NP,) numpy array of random floats
# targetzs: a (NP,) numpy array of random floats
NP=1000
def numpyInterp(x,y,z,values,targetxs,targetys,targetzs):
start=time.time()
zz, yy, xx = numpy.broadcast_arrays(z,y[:,numpy.newaxis],x[:,numpy.newaxis,numpy.newaxis])
grid=numpy.reshape(numpy.array([zz,yy,xx]).swapaxes(1,3),(3,-1)).T
values3D=numpy.reshape(values,-1)
print 'Reshape matrix: ',time.time()-start
start=time.time()
f=scipy.interpolate.LinearNDInterpolator(grid,values3D)
print 'Interpolation: ',time.time()-start
#start=time.time()
#result1=[f(targetzs[i],targetys[i],targetxs[i]) for i in range(len(targetzs))]
#print 'Evaluation (list comprehension): ',time.time()-start
# I found that map is slightly (not much) faster on my machine than list comprehension
start=time.time()
result=numpy.squeeze(map(f,targetzs,targetys,targetxs))
print 'Evaluation (map): ',time.time()-start
return result
def pythonInterp(x,y,z,values,targetxs,targetys,targetzs):
nx=len(x)
ny=len(y)
nz=len(z)
ntarget=targetxs.shape[0]
result=numpy.zeros((ntarget,))
for targ in range(ntarget):
westix=len(x)-2
eastix=len(x)-1
for ix in range(1,nx):
if targetxs[targ] <= x[ix]:
westix=ix-1
eastix=ix
break
southiy=len(y)-2
northiy=len(y)-1
for iy in range(1,ny):
if targetys[targ] <= y[iy]:
southiy=iy-1
northiy=iy
break
upiz=len(z)-1
downiz=len(z)-2
for iz in range(1,nz):
if targetzs[targ] <= z[iz]:
downiz=iz-1
upiz=iz
break
xratio=(targetxs[targ]-x[westix])/(x[eastix]-x[westix])
yratio=(targetys[targ]-y[southiy])/(y[northiy]-y[southiy])
lowerresult=values[downiz,southiy,westix]+(values[downiz,southiy,eastix]-values[downiz,southiy,westix])*xratio+(values[downiz,northiy,westix]-values[downiz,southiy,westix])*yratio+(values[downiz,northiy,eastix]-values[downiz,northiy,westix]-values[downiz,southiy,eastix]+values[downiz,southiy,westix])*xratio*yratio
upperresult=values[upiz,southiy,westix]+(values[upiz,southiy,eastix]-values[upiz,southiy,westix])*xratio+(values[upiz,northiy,westix]-values[upiz,southiy,westix])*yratio+(values[upiz,northiy,eastix]-values[upiz,northiy,westix]-values[upiz,southiy,eastix]+values[upiz,southiy,westix])*xratio*yratio
result[targ]=lowerresult+(upperresult-lowerresult)*(targetzs[targ]-z[downiz])/(z[upiz]-z[downiz])
return result
@jit
def numbaInterp(x,y,z,values,targetxs,targetys,targetzs):
nx=len(x)
ny=len(y)
nz=len(z)
ntarget=targetxs.shape[0]
result=numpy.zeros((ntarget,))
for targ in range(ntarget):
westix=len(x)-2
eastix=len(x)-1
for ix in range(1,nx):
if targetxs[targ] <= x[ix]:
westix=ix-1
eastix=ix
break
southiy=len(y)-2
northiy=len(y)-1
for iy in range(1,ny):
if targetys[targ] <= y[iy]:
southiy=iy-1
northiy=iy
break
upiz=len(z)-1
downiz=len(z)-2
for iz in range(1,nz):
if targetzs[targ] <= z[iz]:
downiz=iz-1
upiz=iz
break
xratio=(targetxs[targ]-x[westix])/(x[eastix]-x[westix])
yratio=(targetys[targ]-y[southiy])/(y[northiy]-y[southiy])
lowerresult=values[downiz,southiy,westix]+(values[downiz,southiy,eastix]-values[downiz,southiy,westix])*xratio+(values[downiz,northiy,westix]-values[downiz,southiy,westix])*yratio+(values[downiz,northiy,eastix]-values[downiz,northiy,westix]-values[downiz,southiy,eastix]+values[downiz,southiy,westix])*xratio*yratio
upperresult=values[upiz,southiy,westix]+(values[upiz,southiy,eastix]-values[upiz,southiy,westix])*xratio+(values[upiz,northiy,westix]-values[upiz,southiy,westix])*yratio+(values[upiz,northiy,eastix]-values[upiz,northiy,westix]-values[upiz,southiy,eastix]+values[upiz,southiy,westix])*xratio*yratio
result[targ]=lowerresult+(upperresult-lowerresult)*(targetzs[targ]-z[downiz])/(z[upiz]-z[downiz])
return result
# Declare input data grid coordinates
z=numpy.arange(10000,100001,10000) # 10
y=numpy.arange(30,60) # 30
x=numpy.arange(0,40) # 40
# Initialize values (pointwise sin)
zz, yy, xx = numpy.broadcast_arrays(z,y[:,numpy.newaxis],x[:,numpy.newaxis,numpy.newaxis])
grid=numpy.array([zz,yy,xx]).swapaxes(1,3)[0,:,:,:]
values=numpy.sin(grid)
# Initialize points for interpolation
targetxs=numpy.random.random((NP,))*40
targetys=numpy.random.random((NP,))*30+30
targetzs=numpy.random.random((NP,))*90000+10000
# Running functions
start=time.time()
print 'Numpy...'
a=numpyInterp(x,y,z,values,targetxs,targetys,targetzs)
print 'Numpy ready: ',time.time()-start,' s, result[0]= ',a[0]
start=time.time()
print 'Python for loop...'
a=pythonInterp(x,y,z,values,targetxs,targetys,targetzs)
print 'Python ready: ',time.time()-start,' s, result[0]= ',a[0]
for i in range(5):
start=time.time()
print 'Numba for loop...'
a=numbaInterp(x,y,z,values,targetxs,targetys,targetzs)
print 'Numba ',i,' ready: ',time.time()-start,' s, result[0]= ',a[0]
1 个回答
1
这两个函数的内部循环方式差别很大，numpyInterp 会遍历广播数组中的每一个元素，而你的 pythonInterp 则假设数据是在一个网格上，只遍历每个维度。所以实际上，一个循环的复杂度是 O(N^3)，而另一个是 O(3N)，这就解释了你看到的速度提升。
你可以使用来自 scipy.ndimage 的插值方法，因为你的数据是在一个规则的网格上，这样应该会更快。