在 Python 中评估极坐标网格上的插值函数的最快方法是什么？

Question

我有两个笛卡尔变量的插值函数（使用 RectBivariateSpline 创建），我正在寻找在极坐标网格上评估该函数的最快方法。

我通过首先在 r（径向坐标）和 t（angular 坐标）中定义空间来解决这个问题，从中创建一个网格，将其转换为meshgrid 到笛卡尔坐标，然后通过循环笛卡尔 meshgrid 来评估每个点的函数。下面的代码演示了这一点。

import numpy as np
import scipy as sp
from scipy.interpolate import RectBivariateSpline

# this shows the type of data/function I'm working with:

n = 1000
x = np.linspace(-10, 10, n)
y = np.linspace(-10, 10, n)
z = np.random.rand(n,n)

fun = RectBivariateSpline(x, y, z)

# this defines the polar grid and converts it to a Cartesian one:

nr = 1000
nt = 360

r = np.linspace(0, 10, nr)
t = np.linspace(0, 2*np.pi, nt)

R, T = np.meshgrid(r, t, indexing = 'ij')

kx = R*np.cos(T)
ky = R*np.sin(T)

# this evaluates the interpolated function over the converted grid:

eval = np.empty((nr, nt))

for i in range(0, nr):
    for j in range(0, nt):
        eval[i][j] = fun(kx[i][j], ky[i][j])

这样，我得到了一个数组，其元素与R和T数组匹配，其中i对应于R，j 到 T。这很重要，因为对于每个半径，我需要对 angular 坐标上的评估函数求和。

这种方法有效，但速度慢得可怕...实际上，我使用的数组比此处的数组大得多。我正在寻找一种方法来加快速度。

我注意到的一件事是，可以将两个一维数组提交给一个 2 变量函数，并返回在输入点的每个可能组合处评估的函数的二维数组。但是，因为我的函数不是极坐标函数，所以我不能只将我的径向和 angular 数组提交给该函数。理想情况下，有一种方法可以将插值函数转换为接受极坐标参数，但我认为这是不可能的。

我应该进一步指出，首先我无法根据径向坐标定义函数：我使用的数据是从离散傅里叶变换输出的，这需要 rectangular单网格数据。

如有任何帮助，我们将不胜感激！

Answer 1

通过检查RectBivariateSplinehere的__call__方法，可以使用grid=False选项来避免这里的慢双循环。

仅此一项就可以将您提供的示例的速度提高一个数量级。我希望在更大的数据集上加速会更好。

此外，方法之间的答案与预期相同。

import numpy as np
import scipy as sp
from scipy.interpolate import RectBivariateSpline

# this shows the type of data/function I'm working with:

n = 1000
x = np.linspace(-10, 10, n)
y = np.linspace(-10, 10, n)
z = np.random.rand(n,n)

fun = RectBivariateSpline(x, y, z)

# this defines the polar grid and converts it to a Cartesian one:

nr = 1000
nt = 360

r = np.linspace(0, 10, nr)
t = np.linspace(0, 2*np.pi, nt)

R, T = np.meshgrid(r, t, indexing = 'ij')

kx = R*np.cos(T)
ky = R*np.sin(T)

# this evaluates the interpolated function over the converted grid:
def evaluate_slow(kx, ky):
    eval = np.empty((nr, nt))

    for i in range(0, nr):
        for j in range(0, nt):
            eval[i][j] = fun(kx[i][j], ky[i][j])
    return eval

def evaluate_fast(kx, ky):
    eval = fun(kx.ravel(), ky.ravel(), grid=False)
    return eval  

%timeit evaluate_slow(kx, ky)
%timeit evaluate_fast(kx, ky)

eval1 = evaluate_slow(kx, ky)
eval2 = evaluate_fast(kx, ky)

print(np.all(np.isclose(eval1, eval2.reshape((nr, nt)))))

1.69 s ± 73.5 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
262 ms ± 2.86 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
True

在 Python 中评估极坐标网格上的插值函数的最快方法是什么？

What is the quickest way to evaluate an interpolated function on a polar grid in Python?

python

arrays

performance

interpolation

numpy