Evaluar la salida de la interpolación 2D a lo largo de una curva

Tengo datos z muestreados desde una función 2D f en los puntos de cuadrícula x, y , como en z = f(x, y) .

Es fácil interpolar f con scipy.interp2d través de f = interp2d(x, y, z ).

Sin embargo, evaluar f(x, y) devuelve una cuadrícula 2D completa como si lo hubiera hecho

 xx, yy = np.meshgrid(x, y) f(xx, yy) 

El comportamiento que quiero es simplemente devolver los valores [f(x[i], y[i]) for i in range(len(x))] , que creo que es el comportamiento de casi cualquier otro método en numpy.

La razón por la que quiero este comportamiento es que estoy buscando el camino trazado a lo largo de la superficie de f durante el “tiempo” por el par (t, u(t)) .

También es sorprendente que np.diag(f(t, u(t))) difiera de np.array([f(ti, u(ti)) for ti in t]) , por lo que no me queda claro cómo llegar a la ruta f(t, u(t)) partir de lo que se devuelve a través de interp2d .

EDIT: Acerca de diag , pensé que parecía que deberíamos tener np.diag(f(t, u(t))) == np.array([f(ti, u(ti)) for ti in t]) , pero ese no es el caso

Ejemplo completo:

 def f(t, u): return (t**2) * np.exp((u**2) / (1 + u**2)) x = np.linspace(0, 1, 250) xx, yy = np.meshgrid(x, x) z = f(xx, yy) f = scipy.interpolate.interp2d(x, y, z) print(f(x, y)) print(np.array([f(xi, yi)[0] for xi, yi in zip(x, y)])) 

Me gustaría que la salida de ambas declaraciones print sea ​​la misma.

El método interp2d devuelve un objeto cuyo método de llamada espera que los vectores x, y sean coordenadas de una cuadrícula rectangular. Y la razón por la que no obtiene los valores deseados de la diagonal de la matriz devuelta es que primero ordena x, y.

Pero hay una solución, que también usé en la consulta de múltiples puntos en una spline bivariada en la base de B-spline . Despues de ejecutar

 import scipy.interpolate as si f = si.interp2d(x, y, z) 

evalúe f no llamándolo, sino pasando sus propiedades tck , seguidas de sus coordenadas x, y, al método interno de bispeu . Me gusta esto:

 print(si.dfitpack.bispeu(f.tck[0], f.tck[1], f.tck[2], f.tck[3], f.tck[4], x, y)[0]) 

Lo anterior devuelve lo mismo que el ciclo lento.

 print(np.array([f(xi, yi)[0] for xi, yi in zip(x, y)])) 

Explicación

El objeto f es secretamente una B-spline de orden 1. Los parámetros de spline (nudos, coeficientes, orden) están contenidos en su propiedad tck y pueden ser utilizados directamente por rutinas de orden inferior para el efecto deseado.

(Idealmente, el método de llamada de f tendría una grid parámetros booleanos que estableceríamos en Falso para hacerle saber que no queremos una evaluación de cuadrícula. Lamentablemente, no está implementado).

Parece que interp2d() comporta como lo hace porque así fue concebida la función Fortran correspondiente. La única solución a esto (en la que puedo pensar) es llamar a f en pares de coordenadas:

 [f(*p)[0] for p in zip(x, y)]