使用 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 = 3.5，它位于 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()

这将为您提供所需的结果：请注意，黑点正确位于平面上，x1 值为3.5。

请注意，大部分“魔法”以及您问题的答案都在于以下两行：

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写了那个解释，这是为一堆aranges写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)

希望有帮助！

Answer 2

在上面的评论中，Ashwith Rego 表示他需要将结果重塑为

np.reshape(result,len(interp_x1),len(interp_x2),len(interp_x3),len(interp_x4))

获得所需形状的interpn 输出。这可能具有所需的形状，但值不在正确的位置。这可能被忽略了，因为在上面的 4D 示例中所有 4 个维度的长度都是相等的。

我做了一个最小的例子来突出我对正确塑造interpn结果的问题的解决方案。这个例子也比上面的一些答案更好地转化为其他维度情况（3D、2D、5D 等......），所以我认为它是对其他回复的有用补充。

import numpy as np
from scipy.interpolate import interpn

# Original data
x1 = np.arange(2)
x2 = np.arange(3)
x3 = np.arange(4)
x4 = np.arange(5)
v = np.reshape(np.arange(len(x1)*len(x2)*len(x3)*len(x4)), (len(x1),len(x2),len(x3),len(x4)))

# New interpolation grid (which for comparison is the same as the origianl)
x1i = x1
x2i = x2
x3i = x3
x4i = x4

# Reshaping trick by stackoverflow user Praveen
interp_mesh = np.array(np.meshgrid(x1i,x2i,x3i,x4i))
ndim = interp_mesh.ndim - 1
interp_points = np.rollaxis(interp_mesh, 0, ndim+1)
interp_points = interp_points.reshape((interp_mesh.size // ndim,ndim))

vi = interpn((x1,x2,x3,x4), v, interp_points) 

# My shaping of the result of interpn
res = np.zeros((len(x1i),len(x2i),len(x3i),len(x4i)))
cnt = 0
for x2i_idx in range(len(x2i)):
    for x1i_idx in range(len(x1i)):
        for x3i_idx in range(len(x3i)):
            for x4i_idx in range(len(x4i)):
                res[x1i_idx,x2i_idx,x3i_idx,x4i_idx] = vi[cnt] 
                cnt += 1

由于本例中的新旧网格是相同的，所以v 和vres 也是如此

print(v-vres)
[[[[0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]]

  [[0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]]

   ...

  [[0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]]

  [[0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]
   [0. 0. 0. 0. 0.]]]]

所有元素都是0，这意味着v和vres相等，即形状和值是正确的。但是，打印最初建议的内容会给出非零值

print(v-np.reshape(vi,len(x1i),len(x2i),len(x3i),len(x4i))
[[[[  0.   0.   0.   0.   0.]
   [  0.   0.   0.   0.   0.]
   [  0.   0.   0.   0.   0.]
   [  0.   0.   0.   0.   0.]]

  [[-40. -40. -40. -40. -40.]
   [-40. -40. -40. -40. -40.]
   [-40. -40. -40. -40. -40.]
   [-40. -40. -40. -40. -40.]]

...

  [[ 40.  40.  40.  40.  40.]
   [ 40.  40.  40.  40.  40.]
   [ 40.  40.  40.  40.  40.]
   [ 40.  40.  40.  40.  40.]]

  [[  0.   0.   0.   0.   0.]
   [  0.   0.   0.   0.   0.]
   [  0.   0.   0.   0.   0.]
   [  0.   0.   0.   0.   0.]]]]

这意味着结果具有正确的形状，但值在错误的位置。

我想不出一个 np.roll 或 np.rolldim 或 np.reshape() 会将 vi 重塑为 (len(x1i),len(x2i),len(x3i),len(x4i)，但是，它确实有效。也许有人可以编写一个比我的四（！）循环（或 N 循环，其中 N 是插值的维数）更好的聪明命令？