我对 NumPy/SciPy 插值方法感到困惑。我使用 LinearNDInterpolator 实现了 3D 线性插值，但我发现它非常慢。然后我用纯 Python 写了一个强力三重 for 循环方法，令人惊讶的是它给了我 1000 倍的加速。我也试了一下 Numba 包，但结果并没有更快。

根据我在互联网上找到的任何来源，与 NumPy/SciPy 和 Numba 相比，Python 循环应该超慢。但这不是我所看到的。

我发布了我运行的整个源代码。我在我的机器上得到了这些时间:

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]

最佳答案

这两个函数在内部循环非常不同，numpyInterp 遍历广播数组的每个元素，而您的 pythonInterp 假设数据在网格上并且只遍历每个元素方面。所以实际上发生的是一个循环是 O(N^3)，而另一个是 O(3N)，这解释了您看到的加速。

您可以使用 scipy.ndimage 中的插值方法，因为您的数据位于规则网格上，这应该更快。

关于python - For 循环似乎比 NumPy/SciPy 3D 插值更快，我们在Stack Overflow上找到一个类似的问题： https://stackoverflow.com/questions/22593277/