Python Scipy для 2D экстраполированной функции сплайна?

Question

Я хотел бы написать экстраполированную сплайн-функцию для 2D-матрицы. Теперь у меня есть экстраполированная функция сплайна для 1D-массивов, как показано ниже. scipy.interpolate.InterpolatedUnivariateSpline() используется.

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 [0], y = y0 [0]; и когда x> x0 [-1], y = y0 [-1]. Здесь, предполагая, что x0 находится в порядке возрастания.

Я хочу иметь подобную функцию, экстраполированный сплайн для работы с 2D матриц с использованием np.select() как в extrapolated_spline_1D. Я думал, что scipy.interpolate.RectBivariateSpline() может помочь, но я не уверен, как это сделать.

Для справки, моя текущая версия extrapolated_spline_2D является очень медленно. Основная идея заключается в:

(1) во-первых, с учетом 1D массивы x0, y0 и 2D массив z2d0 в качестве входных данных, что делает nx0 extrapolated_spline_1D функции, y0_spls, каждый из которых выступает за z2d0 слоя, определенной на у0;

(2) второй, для точки (x, y) не на сетке, вычисляя значения nx0, каждый равен y0_spls [i] (y);

(3) третий, фитинг (х0, y0_spls [I] (Y)) с extrapolated_spline_1D к x_spl и возвращение x_spl (х) в качестве конечного результата.

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

Я думаю, что я пришел с ответом сам, который использует scipy.interpolate.RectBivariateSpline() и более чем в 10 раз быстрее, чем мой старый друг.

Вот функция extrapolated_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 

Я не знаю, что это значит.