通过 scipy.interpolate.griddata 无法对 n 维网格进行插值
Interpolation on n-dimensional grid impossible through scipy.interpolate.griddata
我有一个称为值的观察数组，它来自按以下方式定义的函数：
values = np.array([oscillatory(i) for i in range(points.shape[0])])
values 的形状为 (65,1)，points 是一个数组，在该数组中评估振荡函数，它的形状为 (65,7)，7 是一些在 7 维中充当维度的特征 space（7 只是一个任意数字）。
我正在尝试插入一些在此 space 中定义的任意点。
我尝试用以下方法定义这些点：
grid_x = np.random.uniform(0,10, (100,7))
但是没有用。显然是因为网格定义不正确所以我试过了：
grid_x= np.mgrid[-2:2, 5:99 , 4:5, 5:6, 5:7, 6:67, 7:67]
这又不起作用了。
我通过以下方式调用插值函数：
grid_z1 = gdd(points, values, tuple(grid_x))
但我总是遇到一个巨大的错误，我很难理解。
奇怪的是，如果我以随机方式定义点和值，代码将起作用：
points = np.random.uniform(0,10, (65,3))
values = np.random.uniform(0,10,(points.shape[0],1))
grid_x= np.mgrid[0:2, 5:9 , 4:5]
grid_z1 = gdd(points, values, tuple(grid_x))
这里我只是尝试了3而不是7，因为它更快，但原理是一样的。如果我在 7 个维度中定义它，代码也可以工作。
所以我的问题是：
1) 我怎样才能得到 7 个维度的初始代码 运行 ？
2）如果数组具有相同的形状，为什么随机声明与其他声明相比有效？
如有任何帮助，我们将不胜感激。非常感谢。
我收到的错误如下：
Traceback (most recent call last):
File "/usr/lib/python3.6/code.py", line 91, in runcode
exec(code, self.locals)
File "<input>", line 2, in <module>
File "/usr/local/lib/python3.6/dist-packages/scipy/interpolate/ndgriddata.py", line 222, in griddata
rescale=rescale)
File "interpnd.pyx", line 248, in scipy.interpolate.interpnd.LinearNDInterpolator.__init__
File "qhull.pyx", line 1828, in scipy.spatial.qhull.Delaunay.__init__
File "qhull.pyx", line 354, in scipy.spatial.qhull._Qhull.__init__
scipy.spatial.qhull.QhullError: QH6154 Qhull precision error: Initial simplex is flat (facet 1 is coplanar with the interior point)
While executing: | qhull d Qbb Qt Qz Q12 Qc
Options selected for Qhull 2015.2.r 2016/01/18:
run-id 526315618 delaunay Qbbound-last Qtriangulate Qz-infinity-point
Q12-no-wide-dup Qcoplanar-keep _pre-merge _zero-centrum Qinterior-keep
Pgood _max-width 4 Error-roundoff 1.2e-14 _one-merge 1.1e-13
Visible-distance 7.3e-14 U-coplanar-distance 7.3e-14 Width-outside 1.5e-13
_wide-facet 4.4e-13
precision problems (corrected unless 'Q0' or an error)
2 flipped facets
6 degenerate hyperplanes recomputed with gaussian elimination
12 nearly singular or axis-parallel hyperplanes
6 zero divisors during back substitute
126 zero divisors during gaussian elimination
The input to qhull appears to be less than 4 dimensional, or a
computation has overflowed.
Qhull could not construct a clearly convex simplex from points:
- p2(v4): -1.9 5 4 0.012
- p1(v3): -1.9 5 4 0.0054
- p65(v2): 0 5.5 4.5 4
- p64(v1): 2 6 5 3
- p0(v0): -2 5 4 -8.9e-16
The center point is coplanar with a facet, or a vertex is coplanar
with a neighboring facet. The maximum round off error for
computing distances is 1.2e-14. The center point, facets and distances
to the center point are as follows:
center point -0.7625 5.309 4.309 1.407
facet p1 p65 p64 p0 distance= 0
facet p2 p65 p64 p0 distance= -8.9e-16
facet p2 p1 p64 p0 distance= -8.9e-16
facet p2 p1 p65 p0 distance= -8.9e-16
facet p2 p1 p65 p64 distance= -8.9e-16
These points either have a maximum or minimum x-coordinate, or
they maximize the determinant for k coordinates. Trial points
are first selected from points that maximize a coordinate.
The min and max coordinates for each dimension are:
0: -2 2 difference= 4
1: 5 6 difference= 1
2: 4 5 difference= 1
3: -8.882e-16 4 difference= 4
If the input should be full dimensional, you have several options that
may determine an initial simplex:
- use 'QJ' to joggle the input and make it full dimensional
- use 'QbB' to scale the points to the unit cube
- use 'QR0' to randomly rotate the input for different maximum points
- use 'Qs' to search all points for the initial simplex
- use 'En' to specify a maximum roundoff error less than 1.2e-14.
- trace execution with 'T3' to see the determinant for each point.
If the input is lower dimensional:
- use 'QJ' to joggle the input and make it full dimensional
- use 'Qbk:0Bk:0' to delete coordinate k from the input. You should
pick the coordinate with the least range. The hull will have the
correct topology.
- determine the flat containing the points, rotate the points
into a coordinate plane, and delete the other coordinates.
- add one or more points to make the input full dimensional.
我认为问题在于您需要确保维数和这些维的 'domains' 一致。您不会在未采样的地方插值得到好的结果。我认为您遇到的错误与尝试在这些地方进行计算有关。
以下是如何使 scipy.interpolate.griddata 文档中的示例适用于 7 维示例。我正在使用一个更简单的函数，它只是对 points 数据中的 'features' 求和：
import numpy as np
def func(data):
return np.sum(data, axis=1)
grid = np.mgrid[0:1:5j, 0:1:5j, 0:1:5j, 0:1:5j, 0:1:5j, 0:1:5j, 0:1:5j]
points = np.random.rand(100, 7)
values = func(points)
注意网格覆盖了整个坐标范围points。也就是说，由于 points 中的每一列的值都在 0 到 1 范围内，我应该确保在这些相同的坐标上制作网格。
from scipy.interpolate import griddata
grid_z = griddata(points, values, tuple(grid), method='linear')
现在我有了这个：
>>> grid_z[2, 2, 2, 2, 2, :, :]
array([[ nan, nan, nan, nan, nan],
[ nan, 3. , 3.25, 3.5 , nan],
[ nan, 3.25, 3.5 , 3.75, nan],
[ nan, 3.5 , 3.75, 4. , nan],
[ nan, nan, nan, nan, nan]])
注意有很多 NaN。如果你使用nearest作为方法，你总是会得到一个解决方案，但是linear插值当然需要两个东西之间进行插值，所以超立方体的'edges'是无效的（ 7-D space 有很多优势！）。
我有一个称为值的观察数组，它来自按以下方式定义的函数：
values = np.array([oscillatory(i) for i in range(points.shape[0])])
values 的形状为 (65,1)，points 是一个数组，在该数组中评估振荡函数，它的形状为 (65,7)，7 是一些在 7 维中充当维度的特征 space（7 只是一个任意数字）。
我正在尝试插入一些在此 space 中定义的任意点。 我尝试用以下方法定义这些点：
grid_x = np.random.uniform(0,10, (100,7))
但是没有用。显然是因为网格定义不正确所以我试过了：
grid_x= np.mgrid[-2:2, 5:99 , 4:5, 5:6, 5:7, 6:67, 7:67]
这又不起作用了。 我通过以下方式调用插值函数：
grid_z1 = gdd(points, values, tuple(grid_x))
但我总是遇到一个巨大的错误，我很难理解。
奇怪的是，如果我以随机方式定义点和值，代码将起作用：
points = np.random.uniform(0,10, (65,3))
values = np.random.uniform(0,10,(points.shape[0],1))
grid_x= np.mgrid[0:2, 5:9 , 4:5]
grid_z1 = gdd(points, values, tuple(grid_x))
这里我只是尝试了3而不是7，因为它更快，但原理是一样的。如果我在 7 个维度中定义它，代码也可以工作。 所以我的问题是： 1) 我怎样才能得到 7 个维度的初始代码 运行 ？ 2）如果数组具有相同的形状，为什么随机声明与其他声明相比有效？
如有任何帮助，我们将不胜感激。非常感谢。
我收到的错误如下：
Traceback (most recent call last):
File "/usr/lib/python3.6/code.py", line 91, in runcode
exec(code, self.locals)
File "<input>", line 2, in <module>
File "/usr/local/lib/python3.6/dist-packages/scipy/interpolate/ndgriddata.py", line 222, in griddata
rescale=rescale)
File "interpnd.pyx", line 248, in scipy.interpolate.interpnd.LinearNDInterpolator.__init__
File "qhull.pyx", line 1828, in scipy.spatial.qhull.Delaunay.__init__
File "qhull.pyx", line 354, in scipy.spatial.qhull._Qhull.__init__
scipy.spatial.qhull.QhullError: QH6154 Qhull precision error: Initial simplex is flat (facet 1 is coplanar with the interior point)
While executing: | qhull d Qbb Qt Qz Q12 Qc
Options selected for Qhull 2015.2.r 2016/01/18:
run-id 526315618 delaunay Qbbound-last Qtriangulate Qz-infinity-point
Q12-no-wide-dup Qcoplanar-keep _pre-merge _zero-centrum Qinterior-keep
Pgood _max-width 4 Error-roundoff 1.2e-14 _one-merge 1.1e-13
Visible-distance 7.3e-14 U-coplanar-distance 7.3e-14 Width-outside 1.5e-13
_wide-facet 4.4e-13
precision problems (corrected unless 'Q0' or an error)
2 flipped facets
6 degenerate hyperplanes recomputed with gaussian elimination
12 nearly singular or axis-parallel hyperplanes
6 zero divisors during back substitute
126 zero divisors during gaussian elimination
The input to qhull appears to be less than 4 dimensional, or a
computation has overflowed.
Qhull could not construct a clearly convex simplex from points:
- p2(v4): -1.9 5 4 0.012
- p1(v3): -1.9 5 4 0.0054
- p65(v2): 0 5.5 4.5 4
- p64(v1): 2 6 5 3
- p0(v0): -2 5 4 -8.9e-16
The center point is coplanar with a facet, or a vertex is coplanar
with a neighboring facet. The maximum round off error for
computing distances is 1.2e-14. The center point, facets and distances
to the center point are as follows:
center point -0.7625 5.309 4.309 1.407
facet p1 p65 p64 p0 distance= 0
facet p2 p65 p64 p0 distance= -8.9e-16
facet p2 p1 p64 p0 distance= -8.9e-16
facet p2 p1 p65 p0 distance= -8.9e-16
facet p2 p1 p65 p64 distance= -8.9e-16
These points either have a maximum or minimum x-coordinate, or
they maximize the determinant for k coordinates. Trial points
are first selected from points that maximize a coordinate.
The min and max coordinates for each dimension are:
0: -2 2 difference= 4
1: 5 6 difference= 1
2: 4 5 difference= 1
3: -8.882e-16 4 difference= 4
If the input should be full dimensional, you have several options that
may determine an initial simplex:
- use 'QJ' to joggle the input and make it full dimensional
- use 'QbB' to scale the points to the unit cube
- use 'QR0' to randomly rotate the input for different maximum points
- use 'Qs' to search all points for the initial simplex
- use 'En' to specify a maximum roundoff error less than 1.2e-14.
- trace execution with 'T3' to see the determinant for each point.
If the input is lower dimensional:
- use 'QJ' to joggle the input and make it full dimensional
- use 'Qbk:0Bk:0' to delete coordinate k from the input. You should
pick the coordinate with the least range. The hull will have the
correct topology.
- determine the flat containing the points, rotate the points
into a coordinate plane, and delete the other coordinates.
- add one or more points to make the input full dimensional.
我认为问题在于您需要确保维数和这些维的 'domains' 一致。您不会在未采样的地方插值得到好的结果。我认为您遇到的错误与尝试在这些地方进行计算有关。
以下是如何使 scipy.interpolate.griddata 文档中的示例适用于 7 维示例。我正在使用一个更简单的函数，它只是对 points 数据中的 'features' 求和：
import numpy as np
def func(data):
return np.sum(data, axis=1)
grid = np.mgrid[0:1:5j, 0:1:5j, 0:1:5j, 0:1:5j, 0:1:5j, 0:1:5j, 0:1:5j]
points = np.random.rand(100, 7)
values = func(points)
注意网格覆盖了整个坐标范围points。也就是说，由于 points 中的每一列的值都在 0 到 1 范围内，我应该确保在这些相同的坐标上制作网格。
from scipy.interpolate import griddata
grid_z = griddata(points, values, tuple(grid), method='linear')
现在我有了这个：
>>> grid_z[2, 2, 2, 2, 2, :, :]
array([[ nan, nan, nan, nan, nan],
[ nan, 3. , 3.25, 3.5 , nan],
[ nan, 3.25, 3.5 , 3.75, nan],
[ nan, 3.5 , 3.75, 4. , nan],
[ nan, nan, nan, nan, nan]])
注意有很多 NaN。如果你使用nearest作为方法，你总是会得到一个解决方案，但是linear插值当然需要两个东西之间进行插值，所以超立方体的'edges'是无效的（ 7-D space 有很多优势！）。