Usando scipy.interpolate.interpn para interpolar una matriz N-Dimensional

Supongamos que tengo datos que dependen de 4 variables: a, b, c y d. Quiero interpolar para devolver una matriz 2D que corresponde a un solo valor de a y b, y una matriz de valores para cy d. Sin embargo, el tamaño de la matriz no tiene por qué ser el mismo. Para ser específicos, mis datos provienen de una simulación de transistor. La stream depende de 4 variables aquí. Quiero trazar una variación paramétrica. El número de puntos en el parámetro es mucho menor que el número de puntos para el eje horizontal.

import numpy as np from scipy.interpolate import interpn arr = np.random.random((4,4,4,4)) x1 = np.array([0, 1, 2, 3]) x2 = np.array([0, 10, 20, 30]) x3 = np.array([0, 10, 20, 30]) x4 = np.array([0, .1, .2, .30]) points = (x1, x2, x3, x4) 

Los siguientes trabajos:

 xi = (0.1, 9, np.transpose(np.linspace(0, 30, 4)), np.linspace(0, 0.3, 4)) result = interpn(points, arr, xi) 

y también lo hace esto:

 xi = (0.1, 9, 24, np.linspace(0, 0.3, 4)) result = interpn(points, arr, xi) 

pero no esto:

 xi = (0.1, 9, np.transpose(np.linspace(0, 30, 3)), np.linspace(0, 0.3, 4)) result = interpn(points, arr, xi) 

Como puede ver, en el último caso, el tamaño de las dos últimas matrices en xi es diferente. ¿Este tipo de funcionalidad no es compatible con scipy o estoy usando interpn incorrectamente? Necesito esto para crear una gráfica donde uno de los xi ‘s es un parámetro, mientras que el otro es el eje horizontal.

Trataré de explicarte esto en 2D para que tengas una mejor idea de lo que está sucediendo. Primero, vamos a crear una matriz lineal para probar.

 import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from matplotlib import cm # Set up grid and array of values x1 = np.arange(10) x2 = np.arange(10) arr = x1 + x2[:, np.newaxis] # Set up grid for plotting X, Y = np.meshgrid(x1, x2) # Plot the values as a surface plot to depict fig = plt.figure() ax = fig.gca(projection='3d') surf = ax.plot_surface(X, Y, arr, rstride=1, cstride=1, cmap=cm.jet, linewidth=0, alpha=0.8) fig.colorbar(surf, shrink=0.5, aspect=5) 

Esto nos da: gráfico de superficie de valores

Entonces, digamos que quiere interpolar a lo largo de una línea, es decir, un punto a lo largo de la primera dimensión, pero todos los puntos a lo largo de la segunda dimensión. Estos puntos no están en las matrices originales (x1, x2) obviamente. Supongamos que queremos interpolar a un punto x1 = 3.5 , que se encuentra entre dos puntos en el eje x1.

 from scipy.interpolate import interpn interp_x = 3.5 # Only one value on the x1-axis interp_y = np.arange(10) # A range of values on the x2-axis # Note the following two lines that are used to set up the # interpolation points as a 10x2 array! interp_mesh = np.array(np.meshgrid(interp_x, interp_y)) interp_points = np.rollaxis(interp_mesh, 0, 3).reshape((10, 2)) # Perform the interpolation interp_arr = interpn((x1, x2), arr, interp_points) # Plot the result ax.scatter(interp_x * np.ones(interp_y.shape), interp_y, interp_arr, s=20, c='k', depthshade=False) plt.xlabel('x1') plt.ylabel('x2') plt.show() 

Esto le da el resultado deseado: tenga en cuenta que los puntos negros se encuentran correctamente en el plano, con un valor x1 de 3.5 . parcela de superficie de puntos interpolados

Tenga en cuenta que la mayor parte de la “magia”, y la respuesta a su pregunta, se encuentra en estas dos líneas:

 interp_mesh = np.array(np.meshgrid(interp_x, interp_y)) interp_points = np.rollaxis(interp_mesh, 0, 3).reshape((10, 2)) 

He explicado el funcionamiento de este en otro lugar . En resumen, lo que hace es crear una matriz de tamaño 10×2, que contiene las coordenadas de los 10 puntos que desea interpolar. (La única diferencia entre esa publicación y esta es que he escrito esa explicación para np.mgrid , que es un atajo para escribir np.meshgrid para un grupo de arange s.)

Para su estuche 4x4x4x4, probablemente necesitará algo como esto:

 interp_mesh = np.meshgrid([0.1], [9], np.linspace(0, 30, 3), np.linspace(0, 0.3, 4)) interp_points = np.rollaxis(interp_mesh, 0, 5) interp_points = interp_points.reshape((interp_mesh.size // 4, 4)) result = interpn(points, arr, interp_points) 

¡Espero que ayude!