python中的多个1d插值

我正在模拟CCD阵列中的陷阱。目前我正在使用NumPy和Scipy，并且我已经能够矢量化大部分使我加快速度的呼叫。目前我的代码中的瓶颈是我必须从我的代码的内部循环中从大量不同的插值中检索一个数字。这一特定步骤占用计算时间的约97％。python中的多个1d插值

我做了我的问题在这里的一个简单的例子：

import numpy as np 
from scipy.interpolate import interp1d 

# the CCD array containing values from 0-100 
array = np.random.random(200)*100 

# a number of traps at different positions in the CCD array 
n_traps = 100 
trap_positions = np.random.randint(0,200,n_traps) 

# xvalues for the interpolations 
xval = [0,10,100] 
# each trap has y values corresponding to the x values 
trap_yvals = [np.random.random(3)*100 for _ in range(n_traps)] 
# The xval-to-yval interpolation is made for each trap 
yval_interps = [interp1d(xval,yval) for yval in trap_yvals] 

# moving the trap positions down over the array 
for i in range(len(array)): 
    # calculating new trap position 
    new_trap_pos = trap_positions+i 
    # omitting traps that are outside array 
    trap_inside_array = new_trap_pos < len(array) 
    # finding the array_vals (corresponding to the xvalues in the interpolations) 
    array_vals = array[new_trap_pos[trap_inside_array]] 

    # retrieving the interpolated y-values (this is the bottleneck) 
    yvals = np.array([yval_interps[trap_inside_array[t]](array_vals[t]) 
         for t in range(len(array_vals))]) 

    # some more operations using yvals

是否有此可以优化，也许用用Cython或类似的方法吗？

来源

2016-02-19 Skottfelt

请参阅[this]（http://stackoverflow.com/questions/23909266/interpolation-and-extrapolation-for-large-arrays/34289911#34289911），而不是interp1d ，使用InterpolatedUnivariateSpline。改进性能达几个数量级 –

另一个重大改进是将数组作为参数传递给interp1d/InterpolatedUnivariateSpline，而不是在可能的情况下循环使用单个值 –

@ M.T：感谢关于使用InterpolatedUnivariateSpline加速的提示。我循环单个值的原因是每个值都需要从不同的插值中提取出来，而且我找不到解决方法。 – Skottfelt

我已经仔细考虑了这一点，我想我已经找到了一个我想要分享的很好的解决方案，虽然这意味着我会回答我自己的问题。

首先，我明白了，而不是使用其中一个scipy.interpolation函数，我可以找到两个值之间的插值。这可以用这个小功能

from bisect import bisect_left 

def two_value_interpolation(x,y,val): 
    index = bisect_left(x,val) 
    _xrange = x[index] - x[index-1] 
    xdiff = val - x[index-1] 
    modolo = xdiff/_xrange 
    ydiff = y[index] - y[index-1] 
    return y[index-1] + modolo*ydiff

这给了我一些加速进行，但我想看看我是否能做得更好，所以我移植的功能，用Cython，并且增加了遍历所有的陷阱等等我没有这样做，在Python代码：

# cython: boundscheck=False 
# cython: wraparound=False 
# cython: cdivision=True 

import numpy as np 
cimport numpy as np 

def two_value_interpolation_c(np.ndarray[np.float64_t] x, 
           np.ndarray[np.float64_t, ndim=2] y, 
           np.ndarray[np.float64_t] val_array): 
    cdef unsigned int index, trap 
    cdef unsigned int ntraps=val_array.size 
    cdef long double val, _xrange, xdiff, modolo, ydiff 
    cdef np.ndarray y_interp = np.zeros(ntraps, dtype=np.float64) 

    for trap in range(ntraps): 
     index = 0 
     val = val_array[trap] 
     while x[index] <= val: 
      index += 1 

     _xrange = x[index] - x[index-1] 
     xdiff = val - x[index-1] 
     modolo = xdiff/_xrange 
     ydiff = y[trap,index] - y[trap,index-1] 
     y_interp[trap] = y[trap,index-1] + modolo*ydiff 
    return y_interp

我跑了不同方法的一些时序（使用一些较大的阵列和更多的陷阱比原来的问题表示）：

使用原来的方法，即interp1d :(最好的3）15.1秒

使用InterpolatedUnivariateSpline（K = 1），而不是如interp1d通过@MT建议：（最好的3）7.25秒

使用two_value_interpolation功能：（最好3）1.34秒

使用用Cython实现two_value_interpolation_c :(最好的3）0.113秒

来源

2016-02-24 19:19:37 Skottfelt

python中的多个1d插值

回答

相关问题