logo标签列表

使用scipy.interpolate.interpn插值N维数组

pythonpython-2.7numpyscipy
11

假设我有依赖于4个变量a、b、c和d的数据。我想要插值返回一个2D数组，它对应于单个值的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 是参数，而另一个是水平轴。
看起来你没有正确使用 interpn... 我假设 points 表示你的网格，即每个维度中的“采样点”，用于表示你所知道的值。arr 保存这些已知值。但是，将它们变成随机数意味着很难检查插值是否有效。尝试将它们全部设置为相等，或者更好的方法是在每个四个维度上进行线性插值。xi 应该是你想要知道 arr 值的点的 _坐标_。xi 应该是一个有 k 行、4 列（4D 数据）的数组，其中 k 表示点的数量。 - Praveen
谢谢Praveen！然而，我的问题是这样的。假设我有依赖于4个变量a、b、c、d的数据。我想要插值返回一个2D数组，该数组对应于单个a和b的值，并且对于c和d的值有一个值数组。但是，数组大小不需要相同。具体来说，我的数据来自晶体管模拟。电流在这里取决于4个变量。我想绘制参数变化。参数上的点数远少于水平轴上的点数。 - Ashwith Rego
下次直接编辑您的问题，而不是在评论中添加大量信息。 - Praveen
顺便提一下，注意np.transpose(np.linspace(0, 30, 3))并不会将其转换为列向量。请阅读有关numpy中1D数组的内容。您需要引入一个新的维度，因此请改用np.linspace(0, 30, 3)[:, np.newaxis]。但在这种情况下，您不需要这样做，因为我在我的答案中已经展示了。 - Praveen
啊，转置是个错误。我在试着调整代码看看哪些会起作用，结果不小心加到了这篇文章里。对此感到抱歉！ - Ashwith Rego
2个回答
15

我会尝试用二维的方式来解释，这样你可以更好地理解发生了什么。首先，让我们创建一个线性数组进行测试。

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)

这给我们带来了：

values的表面图

接下来，假设您想沿一条线进行插值，即在第一个维度上选择一个点，在第二个维度上选择所有点。这些点显然不在原始数组(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))

我在其他地方已经解释了它的工作原理。简而言之，它创建了一个大小为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)

希望这能有所帮助！

我会考虑一段时间。非常好的答案。 - Moritz
Praveen，非常感谢您的帮助。之前没有回复您，我很抱歉。现在才有时间回到这个项目。您所解释的内容是有效的，而且您的解释非常清晰易懂。 - Ashwith Rego
我要补充的是，这之后我还需要一个步骤才能得到我需要的结果 - 我需要重新塑造结果，因为对于我问题中的特定示例，我期望得到一个二维数组，在一般情况下我需要返回。 - Ashwith Rego
我要补充的是，这之后我还需要一个步骤才能得到我需要的结果 - 我需要重新塑造结果，因为对于我问题中的特定示例，我期望得到一个二维数组，在一般情况下，我需要返回 np.reshape(result, len(interp_x1), len(interp_x2), len(interp_x3), len(interp_x4))。 - Ashwith Rego
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是插值的维数)？
1
你在技术上并没有回答他的问题，他的问题是关于4D情况的，而你只提到了3D。 - sarema
1
谢谢您指出这一点。我已经更新了回复和示例，以适用于4D情况，虽然它可以工作，但也凸显了需要更聪明的解决方案，因为我的4D解决方案包含四重循环。 - tyskstil
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接