在（x，y）相空间中高效计算轨道的符号变化

Question

我有使用时间积分计算出的轨道，每个位置的坐标对应着 (x, y)。下图展示了这样一条轨道的示例： 现在需要使用 x 和 y 的数组来构造以下序列: 这实际上表示一个序列，其中 0 表示 x 在时间 t_n 和 t_n+1 之间改变了符号，而 1 则表示 y 在时间 t_n 和 t_n+1 之间改变了符号。 对于上述图片，该序列将如下所示：Pxy[n] = 1 0 1 0 1 0 ... 目前我已经找到多种解决方案，但它们在我使用的大型数据集上都非常慢。上述示例轨道的 x 和 y 坐标数组共有 10,000,000 个元素。由于为了足够精确，积分算法中的时间步长需要非常小，所以需要这么多的元素。我能想到的最快算法需要大约 5.3 秒才能运行。然而，这个相同的算法需要针对不同的轨道运行 100 多次。这意味着它将需要很长时间。

我实际上有两种计算方式： 第一种算法基于这个问题：“Efficiently detect sign-changes in python”，并使用了numpy：

def calc_Pxy():
    Pxy = np.full(x.shape, -1)
    Pxy[np.where(np.diff(np.signbit(x)))] = 0
    Pxy[np.where(np.diff(np.signbit(y)))] = 1
    return Pxy[Pxy >= 0]

第二个解决方案只是使用了一个简单的for循环：

def calc_Pxy():
    Pxy = []
    for n, [xn, yn] in enumerate(zip(x[:-1], y[:-1])):
        if xn * x[n + 1] < 0:
            Pxy.append(0)
        if yn * y[n + 1] < 0:
            Pxy.append(1)
    return np.array(Pxy)

!请看编辑! 经过对两种方法性能的测量，使用numpy的第一种方法似乎比普通的for循环慢2.5倍。

我的第一个问题是，为什么会这样？是因为需要重新调整数组大小，因此需要创建新的数组吗？

第二个问题是，我该如何才能比for循环更快呢？

最后，第一种numpy方法的另一个问题是，当x和y在同一时间步骤中都改变符号时，在第一个示例中只记录1次交叉，而在第二个示例中它将始终先记录x，然后是y的交叉。我更喜欢这样做，因为这样不会丢失任何数据。

如果您有任何改进性能的想法，将不胜感激。

编辑：在大型数据集上进行测试后，如果增加数据集的大小，numpy方法变得更快。我的初始测试基于20个值的小数据集，但对于100万个值，numpy方法更快。但是，我仍在寻找使该方法的行为与for循环相同的方法，并且任何进一步的速度提升也是受欢迎的。

- Viktor VN

Answer 1

结合方法1和方法2的思路：

def calc_Pxy3():
    Pxy = []
    for [xn, yn] in zip(np.diff(np.signbit(x)), np.diff(np.signbit(y))):
        if xn:
            Pxy.append(0)
        if yn:
            Pxy.append(1)
    return np.array(Pxy)

与方法2相比，我们的速度提高了约5倍（但仍然比方法1慢5倍）。

Results for n = 10_000_000 random x/y values:

%timeit calc_Pxy1()  # numpy
253 ms ± 3.08 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit calc_Pxy2()  # simple loop
10.8 s ± 124 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit calc_Pxy3()
2.18 s ± 17.2 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

---

Update:

def calc_Pxy4():
    z = np.ravel((np.where(np.diff(np.signbit(x)), 0, -1), np.where(np.diff(np.signbit(y)), 1, -1)), order='F')
    return np.compress(z >= 0, z)

使用

%timeit calc_Pxy4()
361 ms ± 5.27 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

这个方法比Dominik Stańczak的直接使用numba解决方案快（422毫秒），但比I'mahdi的优化后的numba解决方案慢（175毫秒）。

- Stef

Answer 2

您可以简单地使用 numba：

import numba
@numba.njit
def calc_Pxy_numba(x, y):
    Pxy = []
    for n, [xn, yn] in enumerate(zip(x[:-1], y[:-1])):
        if xn * x[n + 1] < 0:
            Pxy.append(0)
        if yn * y[n + 1] < 0:
            Pxy.append(1)
    return np.array(Pxy)

请注意，我在这种方法中避免使用全局变量。

- Dominik Stańczak

Answer 3

为了获得更好的性能，最好不要使用append，你可以尝试这样做：

import numpy as np
import numba as nb

x = np.random.uniform(-2, 2, 10_000_000)
y = np.random.uniform(-2, 2, 10_000_000)

@nb.njit('float32[:](float64[:], float64[:])')
def calc_Pxy5(x, y):
    Pxy = np.empty(len(x)+len(y), dtype=np.float32)
    idx = 0
    for xn, yn in zip(x[:-1] * x[1:] , y[:-1] * y[1:]):
        if xn < 0:
            Pxy[idx] = 0
            idx += 1
        if yn < 0:
            Pxy[idx] = 1
            idx += 1
    return Pxy[:idx]


@nb.njit('float32[:](float64[:], float64[:])')
def calc_Pxy6(x, y):
    Pxy = np.empty(len(x)+len(y), dtype=np.float32)
    idx = 0
    for xn, yn in zip(np.diff(np.signbit(x)) , np.diff(np.signbit(y))):
        if xn :
            Pxy[idx] = 0
            idx += 1
        if yn :
            Pxy[idx] = 1
            idx += 1
    return Pxy[:idx]

在 Colab 上进行基准测试：

>>> %timeit calc_Pxy5(x,y)
10 loops, best of 5: 176.1 ms per loop

>>> %timeit calc_Pxy6(x,y)
10 loops, best of 5: 178 ms per loop

>>> np.all(calc_Pxy5(x,y) == calc_Pxy6(x,y))
True

- I'mahdi