Ndpolator

快速，稀疏网格上的多维线性插值和扩展。

Ndpolator是一个结合插值/外推器，在稀疏（不完全填充）$n$维网格上操作。它在网格定义范围之内和之外估计标量或向量函数值，同时避免需要插补缺失数据或牺牲结构化网格的好处。Ndpolator用C编写以获得速度和可移植性；为了方便，提供了一个使用NumPy数组的Python包装器。

最近更改

ndpolator 1.2.1

• 在网格外顶点中查询点与网格之间的距离计算不正确；已修复。

• find_nearest()函数现在按距离存储最近的点的数组。

ndpolator 1.2

• 现在显式计算并返回网格外查询点与网格之间的距离。

• github工作流程大幅改进。

ndpolator 1.1

• 内存管理改进。

ndpolator 1.0

• 首次发布。

安装

Ndpolator源代码托管在pypi上；您可以通过发出以下命令安装最新版本：

pip install ndpolator

要从GitHub安装ndpolator，请克隆仓库，然后使用pip从本地目录安装

$> git clone https://github.com/aprsa/ndpolator ndpolator
$> cd ndpolator
$> pip install .

安装完成后，您可以通过运行pytest来测试安装

$> cd tests
$> pytest

文档

API参考文档可在gh-pages上找到。

使用示例

为了演示ndpolator的用法，让我们考虑一个具有三个不同顶点大小的轴的3维空间。出于比较的目的，让我们假设我们想要插值和扩展的函数是一个线性标量场

\mathbf a_1 = \{1000, 2000, 3000, 4000, 5000\}, \quad \mathbf a_2 = \{1, 2, 3, 4, 5\}, \quad \mathbf a_3 = \{0.01, 0.02, 0.03, 0.04, 0.05\},

\mathbf F(x, y, z) = \frac{x}{1000} + y + 100 z.

合适的ndpolator实例可以按以下方式初始化和操作

import numpy
import ndpolator

# initialize the axes:
a1 = np.linspace(1000, 5000, 5)
a2 = np.linspace(1, 5, 5)
a3 = np.linspace(0.01, 0.05, 5)

# initialize interpolation space:
ndp = ndpolator.Ndpolator(basic_axes=(a1, a2, a3))

# define a scalar function field and evaluate it across the grid:
def fv(pt):
    return pt[0]/1000 + pt[1] + 100*pt[2]

grid = np.empty((len(ax1), len(ax2), len(ax3), 1))
for i, x in enumerate(ax1):
        for j, y in enumerate(ax2):
            for k, z in enumerate(ax3):
                grid[i, j, k, 0] = fv((x, y, z))

# label the grid ('main') and register it with the ndpolator instance:
ndp.register(table='main', associated_axes=None, grid=grid)

# draw query points randomly within and beyond the definition ranges:
query_pts = np.ascontiguousarray(
    np.vstack((
        np.random.uniform(500, 5500, 1000),
        np.random.uniform(0.5, 5.5, 1000),
        np.random.uniform(0.005, 0.055, 1000))
    ).T
)

# interpolate and extrapolate linearly:
interps = ndp.ndpolate(table='main', query_pts, extrapolation_method='nearest')

目的

多元（$n$维）插值和扩展是数学、统计学和科学中用来估计多维空间中现有数据点之间和之外未知值的技术的。其中，插值涉及估计现有数据点确定的范围内点的函数值，而扩展涉及估计该范围之外的函数值。多元插值有许多稳健的实现，包括k最近邻(Cover & Hart 1967)，自然邻域插值(Sibson 1981)，径向基函数(Hardy 1971)，克里金(Cressie 1990)等。例如，Scipy提供了一个完整的插值模块（scipy.interpolate），它实现了多个多元插值类，包括分段线性、最近邻和径向基函数插值器。不幸的是，Scipy实现的任何方法都不容易用于扩展：最多只能用nan或用户提供的值填充凸包外的值。此外，在规则$n$维网格上操作的内插器不允许有任何缺失数据；这些值要么需要使用非结构化数据内插器进行估计，要么必须完全放弃结构化数据内插器。

Ndpolator旨在填补这一空白：它可以在网格定义范围之内和之外插值和扩展函数值，并且它可以操作不完整的网格。作为额外的好处，ndpolator可以估计标量和向量函数值，并且它可以为位于网格轴上的感兴趣点降低网格维数。它针对速度和可移植性进行了优化（后端是用C语言编写的），并且还提供了一个Python包装器。Ndpolator最初是为PHOEBE食双星建模代码(Prša et al. 2016)开发的，目的是允许插值和扩展恒星大气中的特定强度。鉴于多元插值和扩展领域的空白，ndpolator的开发已从PHOEBE中分离出来，作为独立的软件包提供给社区。

Ndpolator的工作原理

考虑一个标量场或矢量场 $\mathbf{F}$，它在 $N$ 个 $n$ 维点集 $\{\mathbf F (x_1, \dots, x_n)_k\}$ 上采样，其中 $k = 1 \dots N。这些函数值在网格上采样，其中轴线 $\mathbf{a}_k$ 每个网格维度，所以 $\mathsf{X}_{k=1}^{n} \mathbf a_k \equiv \mathbf{a}_1 \times \mathbf{a}_2 \times \dots \times \mathbf{a}_n$ 是一个笛卡尔积，覆盖了整个网格。轴间距不必均匀：顶点可以相隔任何所需的距离，以使网格足够局部线性。如果网格是 完整的，即对于每个网格点都存在一个与之相关的函数值 $\mathbf F(x_1, \dots, x_n)$，则我们有 $N_c = \prod_k l(\mathbf a_k)$ 个函数值样本，其中 $l(\mathbf a)$ 是轴线 $\mathbf a$ 的长度。如果网格是 不完整的，即某些函数值丢失，则 $N < N_c$. Ndpolator 将具有采样值的网格点定义为 节点，将缺失值的网格点定义为 空隙。我们想要估计函数值的点称为 查询点 或 感兴趣点。包围（或相邻于）查询点的最小 $n$ 维子网格称为 超立方体。

单位超立方体变换

ndpolator 的第一个、最基本的原理是，所有插值和外推都是在 单位超立方体 上完成的。在现实世界的应用中，所有轴线都定义在单位区间上通常是很少见的。这可能导致单个轴线上的顶点具有显著不同的数量级。因此，ndpolator 首先将超立方体通过变换归一化为单位超立方体：给定跨越超立方体的连续两个轴线值的集合 $(\mathbf a_{1,p}, \mathbf a_{1,p+1}) \times (\mathbf a_{2,q}, \mathbf a_{2,q+1}) \times \dots \times (\mathbf a_{n,t}, \mathbf a_{n,t+1})$，单位变换将其映射到 $[0, 1] \times [0, 1] \times \dots \times [0, 1] \equiv [0, 1]^n$ 超立方体。所有查询点都受到相同的变换，创建 单位归一化 的查询点。因此，插值（和外推） 总是 在单位超立方体上操作，这在计算上是最经济的，在数值上是最稳定的。这种变换的另一个好处是，外推继承了最近超立方体的网格间距，从而自然地考虑了（可能）不同区域网格中的变量间距。

顺序降维

ndpolator 的第二个操作原理是 顺序降维。考虑一个三维超立方体 \autoref{fig:interpolation}；假设超立方体所有 8 个角上的函数值都已采样，即我们拥有 8 个节点。感兴趣点在左面板中用开放符号表示，以及超立方体面的投影。ndpolator 从最后一个轴线开始，在这个例子中是 $\mathbf a_3$，并沿该轴线将函数值插值到感兴趣点的投影（第二面板）。这些都是 单变量 插值。这个过程产生了 4 个顶点（用开放符号表示），从而将初始维度 $N=3 降低 1，到 $N-1=2。然后重复这个过程（第三面板），这次沿第二个轴线 $\mathbf a_2$，产生 2 个顶点，将维度降低到 1。最后，沿着轴线 $\mathbf a_1 进行最后的插值（右面板），得到单个顶点，即感兴趣点本身。因此，维度降低到 0，并确定函数值。因此，对于一个 $n$ 维超立方体，ndpolator 执行 $\sum_{k=0}^{n-1} 2^k$ 次单变量插值来估计感兴趣点处的 $\mathbf F$，这意味着 $N \log N$ 的时间依赖性。

初始降维

ndpolator的第三个工作原理是初始维度降低。在实际应用中，查询点的某些坐标通常与轴对齐。例如，一个轴可能允许二阶变量的变化，但它的值通常默认为网格上采样的值。当这种情况发生时，每个对齐的轴的初始超立方体维度可以减少1。当查询点与节点重合时，超立方体维度减少到0，不需要进行插值。因此，ndpolator将查询点的每个坐标标记为“在网格上”、“在顶点上”或“超出范围”。当“在顶点上”时，超立方体维度可以立即减少。在这种情况下，时间依赖性减少到$(N-M) \log (N-M)$，其中$M是与轴对齐的坐标数量。

不完整超立方体

ndpolator的第四个工作原理是处理不完整超立方体。如果超立方体的任何角落是空缺的，我们无法进行插值。为此，ndpolator跟踪所有完全定义的$n$维超立方体；当查询点位于不完整超立方体内部时，ndpolator找到最近的完全定义的超立方体，并使用它来估算感兴趣点上的函数值。虽然这在全球范围内仍被视为插值，因为查询点位于网格定义范围内，但估计的函数值严格来说是外推的。请注意，当网格非常稀疏且函数高度非线性时，这可能会导致大量误差累积。在这种情况下，无结构插值技术可能更合适。

外推模式

ndpolator的第五个工作原理是外推。Ndpolator有三种外推方法：none、nearest和linear。当外推方法设置为none时，轴范围外的函数值设置为nan。对于nearest外推方法，ndpolator存储所有节点的列表，并将函数值分配给最近的查询点所在的节点。最后，如果外推方法设置为linear，ndpolator以与处理不完整超立方体等效的方式从最近的完全定义的超立方体进行线性外推。外推方法的选择取决于我们正在估计的多变量函数；如果它是高度非线性的，应避免外推，因此none和nearest可能更合适；如果它是很大程度上线性的或变化缓慢，那么可能需要linear外推方法。Ndpolator是一个线性外推器，因此它不能充分估计非线性多变量函数。

基本轴和相关轴

网格完整性的问题对性能有很大影响；这就是为什么ndpolator的第六个工作原理是区分基本轴和相关轴的原因。其笛卡尔积可能存在空缺的轴被称为基本。对于这些轴，我们需要完整的ndpolator机械来执行插值和外推。另一方面，轴的子集可能在其笛卡尔积中具有所有节点，即它们保证在基本轴采样的所有顶点上被采样；这些被称为相关轴。鉴于它们的采样是确定的，插值和外推可以在没有关于不完整超立方体的担忧的情况下进行——也就是说，只要它们的基本超立方体对应物（由基本轴旋转的超立方体）是完整的。每个相关轴都减少了需要存储的外推查找的超立方体的维度，从而进一步优化性能。

函数值维度

第七个也是最后一个操作原则涉及 函数值维度。大多数插值器假设函数值 $\mathbf F$ 是一个标量；ndpolator 并不做出那个假设。$\mathbf F_r(x_1, \dots, x_n)$ 可以是一个标量或向量，或者任意长度的 $R$（当然，在合理的范围内）。然后，所有节点也必须是 $R$ 维度的。Ndpolator 将分别插值和扩展所有函数值分量，并在感兴趣点提供一个 $R$ 维度的函数值 $\mathbf F 的估计。

超立方体缓存

虽然这并不是 ndpolator 的操作原则的一部分，但 ndpolator 提供了两个辅助函数，import_query_pts() 和 find_hypercubes()，可用于缓存超立方体。这样，调用程序可以分组包含单个超立方体的查询点，并执行批量插值，而无需逐个查询点找到相应的超立方体。虽然索引和超立方体搜索都是二进制的，但在可能的情况下避免查找可以进一步优化运行时间。

API 文档和测试

Ndpolator 采用 GNU 通用公共许可证发布。应用程序编程接口（API）在底层 C 库上可用，请访问 gh-pages。测试套件和自动 API 构建已纳入 github 的持续集成（CI）基础设施。欢迎任何和所有反馈，尤其是问题报告和拉取请求。

致谢

本项目得到美国国家科学基金会资助，资助编号为 #2306996，对此表示衷心感谢。