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

Question

我在 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])

我的objective是为了得到一个interpolated_f，我可以用它来评估单纯形中的任意点，如interpolated_f([0.3, 0.5, 0.2])或interpolated_f(0.3, 0.5, 0.2)。

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

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

谢谢！

---- 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抛出了错误，但是，更多的是

Answer 1

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

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

两者都试一下，然后决定哪种更适合您的目的。

Answer 2

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

我仔细考虑了一下，我很确定（认为我很乐意得到纠正）：

单纯形域意味着它不能有规则的网格（space的某些部分没有值）
我需要删除一维以使插值起作用。因为在 3D 单纯形中，元素需要加起来为 1，所以我会有 [0.4, 0.3, 0.3]，但我不能有 [0.4, 0.3, 0.5]，所以实际上函数只会有值3D 的 2D subspace!

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

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 列，而不是第 3 列）：

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 这意味着插值效果很好！

Answer 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

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

Appropriate numpy/scipy function to interpolate function defined on simplex (non regular grid)

python

interpolation

numpy

scipy

numerical-methods