NumPy: Как рассчитать кусочно-линейный интерполант по нескольким осям

Question

С учетом следующего ndarray t -

In [26]: t.shape                                                                                     
Out[26]: (3, 3, 2)

In [27]: t                                                                                           
Out[27]: 
array([[[ 0,  1],
        [ 2,  3],
        [ 4,  5]],

       [[ 6,  7],
        [ 8,  9],
        [10, 11]],

       [[12, 13],
        [14, 15],
        [16, 17]]])

этот кусочно-линейный интерполятор для точек t[:, 0, 0] можно оценить для [0 , 0.66666667, 1.33333333, 2.] следующим образом, используя numpy .interp -

In [38]: x = np.linspace(0, t.shape[0]-1, 4)                                                         

In [39]: x                                                                                           
Out[39]: array([0.        , 0.66666667, 1.33333333, 2.        ])

In [30]: xp = np.arange(t.shape[0])                                                                  

In [31]: xp                                                                                          
Out[31]: array([0, 1, 2])

In [32]: fp = t[:,0,0]                                                                               

In [33]: fp                                                                                          
Out[33]: array([ 0,  6, 12])

In [40]: np.interp(x, xp, fp)                                                                        
Out[40]: array([ 0.,  4.,  8., 12.])

Каким образом можно эффективно рассчитать и вернуть все интерполанты для всех значений fp -

array([[[ 0,  1],
        [ 2,  3],
        [ 4,  5]],

       [[ 4,  5],
        [ 6,  7],
        [ 8,  9]],

       [[ 8,  9],
        [10, 11],
        [12, 13]],

       [[12, 13],
        [14, 15],
        [16, 17]]])

Answer 1

Поскольку интерполяция равна 1d с изменением значений y, она должна выполняться для каждого 1d среза t. Вероятно, это быстрее, чем l oop явно, но точнее, чем l oop, используя np.apply_along_axis

import numpy as np

t = np.arange( 18 ).reshape(3,3,2)

x = np.linspace( 0, t.shape[0]-1, 4)
xp = np.arange(t.shape[0])


def interfunc( arr ):
    """ Function interpolates a 1d array. """
    return np.interp( x, xp, arr )

np.apply_along_axis( interfunc, 0, t ) # apply function along axis 0

"""  Result
array([[[ 0.,  1.],
        [ 2.,  3.],
        [ 4.,  5.]],

       [[ 4.,  5.],
        [ 6.,  7.],
        [ 8.,  9.]],

       [[ 8.,  9.],
        [10., 11.],
        [12., 13.]],

       [[12., 13.],
        [14., 15.],
        [16., 17.]]]) """

С явными циклами

result = np.zeros((4,3,2))
for c in range(t.shape[1]):
    for p in range(t.shape[2]):
       result[:,c,p] = np.interp( x, xp, t[:,c,p])

На моей машине второй вариант работает пополам время.

Изменить для использования np.nditer

Поскольку результат и параметр имеют разные формы, мне кажется, что мне нужно создать два объекта np.nditer, один для параметр и один для результата. Это моя первая попытка использовать nditer для чего-либо, чтобы это могло быть слишком сложным.

def test( t ):                                                              
    ts = t.shape              
    result = np.zeros((ts[0]+1,ts[1],ts[2]))
    param = np.nditer( [t], ['external_loop'], ['readonly'],  order = 'F')
    with np.nditer( [result], ['external_loop'], ['writeonly'],  order = 'F') as res:       
        for p, r in zip( param, res ):
            r[:] = interfunc(p)
    return result

Это немного медленнее, чем явные циклы, и менее легко следовать, чем любое из других решений.

Answer 2

По запросу @Tis Chris, здесь есть решение, использующее np.nditer с флагом multi_index, но я предпочитаю явный метод вложенных циклов for выше, потому что он на 10% быстрее

In [29]: t = np.arange( 18 ).reshape(3,3,2)                                                          

In [30]: ax0old = np.arange(t.shape[0])                                                              

In [31]: ax0new = np.linspace(0, t.shape[0]-1, 4)                                                    

In [32]: tnew = np.zeros((len(ax0new), t.shape[1], t.shape[2]))                                      

In [33]: it = np.nditer(t[0], flags=['multi_index'])                                                 

In [34]: for _ in it: 
    ...:     tnew[:, it.multi_index[0], it.multi_index[1]] = np.interp(ax0new, ax0old, t[:, it.multi_
    ...: index[0], it.multi_index[1]]) 
    ...:                                                                                             

In [35]: tnew                                                                                        
Out[35]: 
array([[[ 0.,  1.],
        [ 2.,  3.],
        [ 4.,  5.]],

       [[ 4.,  5.],
        [ 6.,  7.],
        [ 8.,  9.]],

       [[ 8.,  9.],
        [10., 11.],
        [12., 13.]],

       [[12., 13.],
        [14., 15.],
        [16., 17.]]])

Answer 3

Вы можете попробовать scipy.interpolate.interp1d:

from scipy.interpolate import interp1d
import numpy as np

 t = np.array([[[ 0,  1],
                [ 2,  3],
                [ 4,  5]],

               [[ 6,  7],
                [ 8,  9],
                [10, 11]],

               [[12, 13],
                [14, 15],
                [16, 17]]])

# for the first slice
f = interp1d(np.arange(t.shape[0]), t[..., 0], axis=0)
# returns a function which you call with values within range np.arange(t.shape[0])

# data used for interpolation
t[..., 0]
>>> array([[ 0,  2,  4],
           [ 6,  8, 10],
           [12, 14, 16]])

f(1)
>>> array([ 6.,  8., 10.])

f(1.5)
>>> array([ 9., 11., 13.])