用于2D外推样条函数的Python Scipy

Question

我想为2D矩阵编写一个外推样条函数。我现在有一个1D阵列的外推样条函数，如下所示。使用scipy.interpolate.Interpolated单变量样条线（）。

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，x0是函数的定义位置，y0是相应的值。当x<x0[0]，y=y0[0]；当x>x0[-1]时，y=y0[-1]。这里，假设x0按升序排列。

我想要一个类似的外推样条函数，用于使用np.select（）处理2D矩阵，如extractd_spline_1D中所述。我认为scipy.interpole.RectBivariateSpline（）可能会有所帮助，但我不确定如何做到。

作为参考，我当前版本的外推d_spline_2D非常慢。基本思想是：

（1） 首先，给定1D阵列x0、y0和2D阵列z2d0作为输入，使nx0外推d_spline_1D函数y0_spls，每个函数代表在y0上定义的层z2d0；

（2） 第二，对于不在网格上的点（x，y），计算nx0个值，每个值等于y0_spls[i]（y）；

（3） 第三，拟合（x0，y0_spls[i]（y）），外推spline_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

Answer 1

我在这里做了一项名为GlobalSpline2D的类似工作，它在线性、三次或五次样条曲线下都能完美工作。

它基本上继承了interp2d，并通过插值单变量样条线促进了对二维外推的使用。两者都是病态的内部功能。

其用法应参考文档以及interp2d的调用方法。

Answer 2

我想我自己已经想出了一个答案，它使用了scipy.interpole.RectBivariateSpline（），比我的旧方法快10多倍。

这是函数extracted_spline_2D_new。

def extrapolated_spline_2D_new(x0,y0,z2d0):
    '''    
    x0,y0 : array_like，1-D arrays of coordinates in strictly ascending order. 
    z2d0  : array_like，2-D array of data with shape (x.size,y.size).
    ''' 
    assert z2d0.shape == (x0.shape[0],y0.shape[0])
    spline = scipy.interpolate.RectBivariateSpline(x0,y0,z2d0,kx=3,ky=3)
    '''
    scipy.interpolate.RectBivariateSpline
    x,y : array_like, 1-D arrays of coordinates in strictly ascending order.
    z   : array_like, 2-D array of data with shape (x.size,y.size).
    '''  
    def f(x,y,spline=spline):
        '''
        x and y have the same shape with the output. 
        '''
        x = np.array(x,dtype='f4')
        y = np.array(y,dtype='f4') 
        assert x.shape == y.shape 
        ndim = x.ndim   
        # We want the output to have the same dimension as the input, 
        # and when ndim == 0 or 1, spline(x,y) is always 2D. 
        if   ndim == 0: result = spline(x,y)[0][0]
        elif ndim == 1: 
            result = np.array([spline(x[i],y[i])[0][0] for i in range(len(x))])
        else:           
            result = np.array([spline(x.flatten()[i],y.flatten()[i])[0][0] for i in range(len(x.flatten()))]).reshape(x.shape)         
        return result
    return f

注：在上面的版本中，我逐个计算值，而不是使用下面的代码。

def f(x,y,spline=spline):
    '''
    x and y have the same shape with the output. 
    '''
    x = np.array(x,dtype='f4')
    y = np.array(y,dtype='f4') 
    assert x.shape == y.shape 
    ndim = x.ndim
    if   ndim == 0: result = spline(x,y)[0][0]
    elif ndim == 1: 
         result = spline(x,y).diagonal()
    else:           
         result = spline(x.flatten(),y.flatten()).diagonal().reshape(x.shape)       
    return result

因为当我试图用下面的代码进行计算时，它有时会给出错误信息：

<ipython-input-65-33285fd2319d> in f(x, y, spline)
 29         if   ndim == 0: result = spline(x,y)[0][0]
 30         elif ndim == 1:
---> 31             result = spline(x,y).diagonal()
 32         else:
 33             result = spline(x.flatten(),y.flatten()).diagonal().reshape(x.shape)
/usr/local/lib/python2.7/site-packages/scipy/interpolate/fitpack2.pyc in __call__(self, x, y, mth, dx, dy, grid)
826                 z,ier = dfitpack.bispev(tx,ty,c,kx,ky,x,y)
827                 if not ier == 0:
--> 828                     raise ValueError("Error code returned by bispev: %s" % ier)
829         else:
830             # standard Numpy broadcasting
ValueError: Error code returned by bispev: 10

我不知道这意味着什么。