适当的 numpy/scipy 函数来插入在单纯形（非规则网格）上定义的函数

Posted 2023-02-18

【中文标题】适当的 numpy/scipy 函数来插入在单纯形（非规则网格）上定义的函数

【英文标题】：Appropriate numpy/scipy function to interpolate function defined on simplex (non regular grid)

【发布时间】：2019-01-28 03:29:24

【问题描述】：

我在 3 维单纯形上定义了一个函数。即点 x,y,z 的集合，每个点都在 0 和 1 之间，使得 x + y + z = 1.0

例如，如果我为每个 x、y 和 z 考虑 4 个点，那么我将得到一个如下所示的 (10, 3) numpy 数组（每行总和正好为 1）：

points = array([[0.        , 0.        , 1.        ],
       [0.        , 0.33333333, 0.66666667],
       [0.        , 0.66666667, 0.33333333],
       [0.        , 1.        , 0.        ],
       [0.33333333, 0.        , 0.66666667],
       [0.33333333, 0.33333333, 0.33333333],
       [0.33333333, 0.66666667, 0.        ],
       [0.66666667, 0.        , 0.33333333],
       [0.66666667, 0.33333333, 0.        ],
       [1.        , 0.        , 0.        ]])

我添加了生成单纯形的便利函数：

def generate_simplex_3dims(n_per_dim):
    xlist = np.linspace(0.0, 1.0, n_per_dim)
    ylist = np.linspace(0.0, 1.0, n_per_dim)
    zlist = np.linspace(0.0, 1.0, n_per_dim)
    return np.array([[x, y, z] for x in xlist for y in ylist for z in zlist
                     if np.allclose(x+y+z, 1.0)])

我还将为这些点设置值。例如，让我们生成如下值：

def approx_this_f(x, y, z):
    return 2*x - y + 5*z

values = np.empty(len(points))
for i, point in enumerate(points):
    values[i] = approx_this_f(point[0], point[1],
                         point[2])

我的目标是获得一个interpolated_f，我可以用它来评估像interpolated_f([0.3, 0.5, 0.2]) 或interpolated_f(0.3, 0.5, 0.2) 这样的单纯形中的任意点。

我查看了文档，但不明白这里合适的插值器是什么，因为我的网格点是在单纯形上定义的并且我想要返回一个插值函数。

我试过scipy.interpolate.griddata，它只适用于method='nearest'，这个返回一个值数组，但我需要一个插值函数。我在 scipy 上看到了其他函数，它们返回一个插值函数，但似乎只适用于常规网格。

谢谢！

---- 以griddata 为例，以防万一 ------

from scipy.interpolate import griddata
xi = generate_simplex_3dims(n_per_dim=20) #Generates lots of points
interpolated_grid = griddata(points, values, xi,
         method='linear') #this fails
interpolated_grid = griddata(points, values, xi,
         method='nearest') #this works, but returns a grid, not a function

method=linear 抛出错误，但是，更多

【问题讨论】：

regularGridInterpolator 在这里有用吗？docs.scipy.org/doc/scipy-0.16.1/reference/generated/…

@GhasemNaddaf 我想我不能使用它。函数的域是一个单纯形，我认为它不能写成规则网格（但如果可以的话，那就太好了！）

【参考方案1】：

您需要一种方法，该方法 (a) 接受非结构化 N 维数据点（N > 2），并且 (b) 返回一个可调用对象。阅读docs 我看到两个选项

LinearNDInterpolator - 仅限分段线性插值。 Rbf - 使用径向基函数；平滑但不会像分段线性插值器那样尊重数据的单调性或最大/最小值。

尝试两者，然后决定哪种方式更适合您的目的。

【讨论】：

谢谢，这很有用！我赞成你的回答，但我发现，因为域是单纯形，我应该删除一个维度。我会写一个答案来解释

【参考方案2】：

感谢@user6655984 的回答，我知道该怎么做（谢谢！）

我想了很多，我很确定（我很乐意得到纠正）：

单纯形域意味着它不能有规则的网格（空间的某些部分没有值） 我需要删除一个维度才能使插值起作用。由于在 3D 单纯形中，元素需要加起来为 1，我将有 [0.4, 0.3, 0.3]，但我不能有 [0.4, 0.3, 0.5]，所以实际上该函数将只有值 on 3D 的 2D 子空间！

让我们与问题中的设置相同

def approx_this_f(x, y, z):
    return 2*x - y + 5*z

#Uses the function defined in the question
simplex_points = generate_simplex_3dims(10)
values = np.empty(len(simplex_points))
for i, lambda_0 in enumerate(simplex_points):
    values[i] = approx_this_f(lambda_0[0], lambda_0[1],
                         lambda_0[2])

由于第 2 条评论，以下内容不起作用：

from scipy.interpolate import LinearNDInterpolator

interpolated = LinearNDInterpolator(simplex_points, values)

并抛出此错误

QhullError: qhull precision warning: 
The initial hull is narrow (cosine of min. angle is 0.9999999999999999).
Is the input lower dimensional (e.g., on a plane in 3-d)?  Qhull may
produce a wide facet. 

所以我需要传递少一维的点（即只有第 1 列和第 2 列，而不是第三列）：

interpolated = LinearNDInterpolator(simplex_points[:, 0:2], values)

现在我们可以在其他点上评估它

#Silly code to make the original function take a matrix
def approx_this_f_vec(array):
    res = np.empty(len(array))
    for row in range(len(array)):
         res[row] = approx_this_f(*array[row])
     return res

points_test = src.generate_simplex_3dims(50) #1275 new points
interpolated_vals = interpolated_f(points_test[:, 0:2])
real_values = approx_this_f_vec(points_test)

print((interpolated_vals - real_values).max())

给出1.77e-15，这意味着插值做得很好！

【参考方案3】：

这是一个使用 simpllex_grid 和线性插值的简单解决方案：

import dit
import numpy as np
from scipy.interpolate import LinearNDInterpolator
import matplotlib.pyplot as plt


dim = 3 #dimensions
res= 10  #resolution
grid = list(dit.simplex_grid(dim, res, using=np.array))
grid = list(map(list, zip(*grid)))  #transpose
x = grid[0]
y = grid[1]
z = np.hypot(x, y) #Hypothenuse function (as a function of only 2 variables of the simplex)
X = np.linspace(min(x), max(x)) #Creates a mesh with may points (around 50)
Y = np.linspace(min(y), max(y))
X, Y = np.meshgrid(X, Y)  # 2D grid for interpolation
interp = LinearNDInterpolator(list(zip(x, y)), z)
Z = interp(X, Y) #This is the function! :D (returns nah if outside of sample)
plt.pcolormesh(X, Y, Z, shading='auto')
plt.plot(x, y, ".",color="black", label="input point")
plt.legend()
plt.colorbar()
plt.axis("equal")
plt.show()

图像输出可以在这里找到： [1]：https://i.stack.imgur.com/8K01u.png