使用scipy.interpolate.interpn插入N維數組

Question

假設我有依賴4個變量的數據：a，b，c和d。我想插入返回一個二維數組，它對應於a和b的單個值，以及c和d的值數組。但是，數組大小不必相同。具體來說，我的數據來自晶體管模擬。電流取決於這裏的4個變量。我想繪製一個參數變化。參數上的點數比橫軸的點數少很多。

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) 

以下工作：

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

等做到這一點：

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

但不是這樣的：

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

正如你可以看到，在過去的情況下， ，最後兩個數組的大小在xi是不同的。 scipy不支持這種功能嗎？或者我錯誤地使用了interpn？我需要這個創建一個情節，其中一個xi是一個參數，而另一個是橫軸。

Answer 1

我會嘗試在2D中向您解釋這一點，以便您更好地瞭解發生了什麼。首先，我們來創建一個線性數組來測試。

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) 

這給我們：

然後，讓我們說你要沿直線插補，即一個點沿着第一維度，而是沿着第二個維度的所有點。這些點顯然不在原始數組(x1, x2)中。假設我們想插入到位於x1軸兩點之間的點x1 = 3.5。

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

這使您可以根據需要的結果：注意，黑點正確地趴在飛機上，在3.5的X1值。

注意，大部分的「魔力」，並回答你的問題，就在於這兩條線：

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

我已經解釋過這個elsewhere的工作。簡而言之，它所做的就是創建一個大小爲10x2的數組，其中包含要插入arr處的10個點的座標。 （該職位，這一次之間的唯一區別是，我已經寫了解釋np.mgrid，這是一個捷徑，爲一堆arange在寫np.meshgrid）

爲了您的4x4x4x4情況下，你可能需要的東西像這樣：

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) 

希望有所幫助！