有效计算 (x, y) 相空间中轨道的符号变化

Question

我有使用时间积分计算出的轨道。(x, y)这些轨道的每个位置都有坐标。这种轨道的一个例子如下图所示：

现在使用这些数组作为 x 和 y，需要构建以下系列：

这实际上意味着一个序列，其中 0 表示 x 在时间t_n和之间改变了符号，1 表示 y 在和t_n+1之间改变了符号。对于上图，该系列将如下所示：t_nt_n+1Pxy[n] = 1 0 1 0 1 0 ...

目前我已经找到了多种解决方案。但它们在我使用的大型数据集上都非常慢。上面显示的示例轨道的 x 和 y 坐标数组有 10,000,000 个元素。它需要这么多元素，因为为了足够准确，积分算法中的时间步需要非常小。我能想到的最快算法运行时间大约为 5.3 秒。然而，相同的算法需要针对不同的轨道运行 100 多次。这意味着这将需要很长时间。

我本质上有两种计算方法：第一种算法基于这个问题：“ EfficientlyDetectsign-changesinpython ”并使用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 = 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]

！请参阅编辑！在测量了两种方法的性能后，似乎第一种使用 numpy 的方法比普通的 for 循环慢 2.5 倍。我的第一个问题是，这是为什么？是否因为需要创建新数组而重塑数组需要时间？第二个问题是，如何才能比 for 循环更快？

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

任何有关如何提高性能的想法将不胜感激。

编辑：在大型数据集上进行测试后，如果数据集的大小增加，numpy 方法似乎会变得更快。我最初的测试基于 20 个值的小数据集，但对于 1,000,000 个值，numpy 方法更快。不过，我仍在寻找一种方法，使该方法的行为与 for 循环相同，当然，任何进一步的速度改进也是受欢迎的。

Answer 1

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

\n

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

\n

与方法 2 相比，我们获得了约 5 倍的加速（但仍比方法 1 慢 5 倍）。

\n

Results for n = 10_000_000 random x/y values:\n\n%timeit calc_Pxy1()  # numpy\n253 ms \xc2\xb1 3.08 ms per loop (mean \xc2\xb1 std. dev. of 7 runs, 1 loop each)\n\n%timeit calc_Pxy2()  # simple loop\n10.8 s \xc2\xb1 124 ms per loop (mean \xc2\xb1 std. dev. of 7 runs, 1 loop each)\n\n%timeit calc_Pxy3()\n2.18 s \xc2\xb1 17.2 ms per loop (mean \xc2\xb1 std. dev. of 7 runs, 1 loop each)\n

\n \n

更新：

\n

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

\n

和

\n

%timeit calc_Pxy4()\n361 ms \xc2\xb1 5.27 ms per loop (mean \xc2\xb1 std. dev. of 7 runs, 1 loop each)\n

\n

比 Dominik Sta\xc5\x84czak\ 的直接 numba 解决方案（422 毫秒）更快，但比 I\'mahdi\ 的优化 numba 解决方案（175 毫秒）慢。

\n