python中重心坐标的矢量化计算

Question

在科学空间这里有Delaunay函数。文档包括example of how to calculate barycentric coordinates.

下面的代码将使用循环计算重心坐标。你知道吗

points = np.array([(0,0),(0,1),(1,0),(1,1)])
samples = np.array([(0.5,0.5),(0,0),(0.1,0.1)])

dim    = len(points[0])               # determine the dimension of the samples
simp   = Delaunay(points)             # create simplexes for the defined points
s      = simp.find_simplex(samples)   # for each sample, find corresponding simplex for each sample
b0      = np.zeros((len(samples),dim)) # reserve space for each barycentric coordinate
for ii in range(len(samples)):
    b0[ii,:] = simp.transform[s[ii],:dim].dot((samples[ii] - simp.transform[s[ii],dim]).transpose())
coord = np.c_[b0, 1 - b0.sum(axis=1)]

这对于将短样本列表转换为重心坐标是可以的，但是对于非常大的样本列表，性能很差。如何修改它以利用numpy/scipy中的矢量化数学来提高性能？你知道吗

Answer 1

考虑以下修改（对于用numpy方法替换的循环）：

def f_1(points, samples):
    """ original """

    dim = len(points[0])
    simp = ssp.Delaunay(points)
    s = simp.find_simplex(samples)
    b0 = np.zeros((len(samples), dim))

    for ii in range(len(samples)):
        b0[ii, :] = simp.transform[s[ii], :dim].dot(
            (samples[ii] - simp.transform[s[ii], dim]).transpose())
    coord = np.c_[b0, 1 - b0.sum(axis=1)]

    return coord

def f_2(points, samples):
    """ modified """

    simp = ssp.Delaunay(points)
    s = simp.find_simplex(samples)

    b0 = (simp.transform[s, :points.shape[1]].transpose([1, 0, 2]) *
          (samples - simp.transform[s, points.shape[1]])).sum(axis=2).T
    coord = np.c_[b0, 1 - b0.sum(axis=1)]

    return coord

测试用例：

N = 100
points = np.array(list(itertools.product(range(N), repeat=2)))
samples = np.random.rand(100_000, 2) * N

结果：

%timeit f_1(points, samples)
712 ms ± 2.76 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit f_2(points, samples)
422 ms ± 809 µs per loop (mean ± std. dev. of 7 runs, 1 loop each)

对于修改后的版本，行simp.find_simplex(samples)给出了大约95%的运行时间。所以，我想你对矢量化没有别的办法了。为了进一步提高性能，您需要find_simplex方法的另一个实现或解决问题的另一种方法。你知道吗