Dado el siguiente 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]]])

Este interpolante lineal por partes para los puntos t[:, 0, 0] puede evaluarse para [0 , 0.66666667, 1.33333333, 2.] de la siguiente manera usando 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.])

¿Cómo pueden calcularse y devolverse todos los interpolantes de manera eficiente para todos los valores de 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]]])
0
user2309803 27 abr. 2020 a las 21:56

3 respuestas

La mejor respuesta

Como la interpolación es 1d con el cambio de los valores y, debe ejecutarse para cada 1d segmento de t. Probablemente sea más rápido hacer un bucle explícitamente, pero es más fácil hacerlo usando 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.]]]) """

Con bucles explícitos

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])

En mi máquina, la segunda opción se ejecuta en la mitad del tiempo.

Editar para usar np.nditer

Como el resultado y el parámetro tienen formas diferentes, parece que tengo que crear dos objetos np.nditer, uno para el parámetro y otro para el resultado. Este es mi primer intento de utilizar nditer para cualquier cosa, por lo que podría ser demasiado complicado.

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

Es un poco más lento que los bucles explícitos y menos fácil de seguir que cualquiera de las otras soluciones.

1
Tls Chris 28 abr. 2020 a las 12:48

Puedes probar 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.])
0
Paddy Harrison 27 abr. 2020 a las 20:03

Según lo solicitado por @Tis Chris, aquí hay una solución que utiliza np.nditer con la bandera multi_index pero prefiero el método de bucles anidados explícitos for anterior porque es un 10% más rápido

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.]]])
1
user2309803 29 abr. 2020 a las 05:29