Scipy Fast 1-D插值,无任何循环

Question

我有两个2D数组,x(ni,nj)和y(ni,nj),我需要在一个轴上进行插值.我想为每个ni插入最后一个轴.

我写

import numpy as np
from scipy.interpolate import interp1d

z = np.asarray([200,300,400,500,600])
out = []
for i in range(ni):
    f = interp1d(x[i,:], y[i,:], kind='linear')
    out.append(f(z))
out = np.asarray(out)

但是,我认为如果数组大小太大,这种方法由于循环而效率低且速度慢.像这样插入多维数组的最快方法是什么？有没有办法在没有循环的情况下执行线性和三次插值？谢谢.

Answer 1

你提出的方法确实有一个python循环,所以对于大的值,ni它会变慢.也就是说,除非你有大片,否则你ni不必担心.

我使用以下代码创建了示例输入数据:

def sample_data(n_i, n_j, z_shape) :
    x = np.random.rand(n_i, n_j) * 1000
    x.sort()
    x[:,0] = 0
    x[:, -1] = 1000
    y = np.random.rand(n_i, n_j)
    z = np.random.rand(*z_shape) * 1000
    return x, y, z

并用这两个版本的线性插值测试了它们:

def interp_1(x, y, z) :
    rows, cols = x.shape
    out = np.empty((rows,) + z.shape, dtype=y.dtype)
    for j in xrange(rows) :
        out[j] =interp1d(x[j], y[j], kind='linear', copy=False)(z)
    return out

def interp_2(x, y, z) :
    rows, cols = x.shape
    row_idx = np.arange(rows).reshape((rows,) + (1,) * z.ndim)
    col_idx = np.argmax(x.reshape(x.shape + (1,) * z.ndim) > z, axis=1) - 1
    ret = y[row_idx, col_idx + 1] - y[row_idx, col_idx]
    ret /= x[row_idx, col_idx + 1] - x[row_idx, col_idx]
    ret *= z - x[row_idx, col_idx]
    ret += y[row_idx, col_idx]
    return ret

interp_1是Dave的回答,是代码的优化版本.interp_2是一个线性插值的矢量化实现,它可以避免任何python循环.编写类似这样的东西需要对numpy中的广播和索引有充分的理解,而且有些东西的优化程度要低于那些interp1d.一个主要的例子是找到插入值的bin:interp1d一旦找到bin,肯定会提前摆脱循环,上面的函数是将值与所有bin进行比较.

因此,结果将是非常依赖于什么n_i以及n_j是,甚至多久你的数组z值进行插值是.如果n_j它很小而且n_i很大,那么你应该期待它的优势interp_2,而且interp_1如果它是相反的话.较小的z应该是一个优势interp_2,更长的interp_1.

其实我已经超时都用各种方法n_i和n_j用于z形状(5,)和(50,),这里是图:

因此,对于z形状来说,(5,)你应该interp_2随时随地n_j < 1000和interp_1其他地方一起去.毫不奇怪,z形状的门槛是不同的(50,),现在已经存在n_j < 100.看来你应该坚持使用你的代码似乎很有诱惑力n_j * len(z) > 5000,但interp_2如果不是这样的话就改变它,但是在那个陈述中有大量的推断!如果你想进一步尝试自己,这里是我用来生成图形的代码.

n_s = np.logspace(1, 3.3, 25)
int_1 = np.empty((len(n_s),) * 2)
int_2 = np.empty((len(n_s),) * 2)
z_shape = (5,)

for i, n_i in enumerate(n_s) :
    print int(n_i)
    for j, n_j in enumerate(n_s) :
        x, y, z = sample_data(int(n_i), int(n_j), z_shape)
        int_1[i, j] = min(timeit.repeat('interp_1(x, y, z)',
                                        'from __main__ import interp_1, x, y, z',
                                        repeat=10, number=1))
        int_2[i, j] = min(timeit.repeat('interp_2(x, y, z)',
                                        'from __main__ import interp_2, x, y, z',
                                        repeat=10, number=1))

cs = plt.contour(n_s, n_s, np.transpose(int_1-int_2))
plt.clabel(cs, inline=1, fontsize=10)
plt.xlabel('n_i')
plt.ylabel('n_j')
plt.title('timeit(interp_2) - timeit(interp_1), z.shape=' + str(z_shape))
plt.show()