我想为二维矩阵编写一个外推样条函数。我现在拥有的是一维数组的外推样条函数，如下所示。scipy.interpolate.InterpolatedUnivariateSpline() http://docs.scipy.org/doc/scipy-0.16.0/reference/generated/scipy.interpolate.InterpolatedUnivariateSpline.html用来。
import numpy as np
import scipy as sp
def extrapolated_spline_1D(x0,y0):
x0 = np.array(x0)
y0 = np.array(y0)
assert x0.shape == y.shape
spline = sp.interpolate.InterpolatedUnivariateSpline(x0,y0)
def f(x, spline=spline):
return np.select(
[(x<x0[0]), (x>x0[-1]), np.ones_like(x,dtype='bool')],
[np.zeros_like(x)+y0[0], np.zeros_like(x)+y0[-1], spline(x)])
return f
它采用 x0（定义函数的位置）和 y0（相应的值）。当x x0[-1] 时，y = y0[-1]。这里，假设x0按升序排列。
我想要一个类似的外推样条函数来处理二维矩阵np.select() http://docs.scipy.org/doc/numpy-1.10.1/reference/generated/numpy.select.html as in 外推样条_1D。我想scipy.interpolate.RectBivariateSpline() http://docs.scipy.org/doc/scipy/reference/generated/scipy.interpolate.RectBivariateSpline.html#scipy.interpolate.RectBivariateSpline可能有帮助，但我不知道该怎么做。
作为参考，我当前的版本外推样条_2D is 非常慢。
基本思想是：
(1) 首先，给定一维数组 x0、y0 和二维数组 z2d0 作为输入，使得 nx0外推样条_1D函数y0_spls，每个函数代表y0上定义的一层z2d0；
(2) 其次，对于不在网格上的点(x,y)，计算nx0个值，每个值等于y0_spls[i](y)；
(3) 第三，拟合 (x0, y0_spls[i](y)) 与外推样条_1D到 x_spl 并返回 x_spl(x) 作为最终结果。
def extrapolated_spline_2D(x0,y0,z2d0):
'''
x0,y0 : array_like, 1-D arrays of coordinates in strictly monotonic order.
z2d0 : array_like, 2-D array of data with shape (x.size,y.size).
'''
nx0 = x0.shape[0]
ny0 = y0.shape[0]
assert z2d0.shape == (nx0,ny0)
# make nx0 splines, each of which stands for a layer of z2d0 on y0
y0_spls = [extrapolated_spline_1D(y0,z2d0[i,:]) for i in range(nx0)]
def f(x, y):
'''
f takes 2 arguments at the same time --> x, y have the same dimention
Return: a numpy ndarray object with the same shape of x and y
'''
x = np.array(x,dtype='f4')
y = np.array(y,dtype='f4')
assert x.shape == y.shape
ndim = x.ndim
if ndim == 0:
'''
Given a point on the xy-plane.
Make ny = 1 splines, each of which stands for a layer of new_xs on x0
'''
new_xs = np.array([y0_spls[i](y) for i in range(nx0)])
x_spl = extrapolated_spline_1D(x0,new_xs)
result = x_spl(x)
elif ndim == 1:
'''
Given a 1-D array of points on the xy-plane.
'''
ny = len(y)
new_xs = np.array([y0_spls[i](y) for i in range(nx0)]) # new_xs.shape = (nx0,ny)
x_spls = [extrapolated_spline_1D(x0,new_xs[:,i]) for i in range(ny)]
result = np.array([x_spls[i](x[i]) for i in range(ny)])
else:
'''
Given a multiple dimensional array of points on the xy-plane.
'''
x_flatten = x.flatten()
y_flatten = y.flatten()
ny = len(y_flatten)
new_xs = np.array([y0_spls[i](y_flatten) for i in range(nx0)])
x_spls = [extrapolated_spline_1D(x0,new_xs[:,i]) for i in range(ny)]
result = np.array([x_spls[i](x_flatten[i]) for i in range(ny)]).reshape(y.shape)
return result
return f
