Add variant reduce_accumulate2

This adds a new variant reduce_accumulate2 using the scan with an associative operator that is slower in a single-threaded context but can take full advantage of parallel scan implementations.

Author: avitase

README.md

@@ -18,7 +18,8 @@ More precisely,
 1. [`branchless`](fast_frechet/branchless.py): A variant w/o branches.
 1. [`linear_memory`](fast_frechet/linear_memory.py): This formulation reduces the quadratic memory footprint to a linear one.
 1. [`accumulate`](fast_frechet/accumulate.py): Formulation using a scan operation.
-1. [`reduce_accumulate`](fast_frechet/reduce_accumulate.py): Formulation using a fold operation.
+1. [`reduce_accumulate`](fast_frechet/reduce_accumulate.py): Formulation using scan and fold operations.
+1. [`reduce_accumulate2`](fast_frechet/reduce_accumulate2.py): Alternative formulation using the scan with an associative operator that is slower in a single-threaded context but can take full advantage of parallel scan implementations. 
 1. [`compiled`](fast_frechet/compiled.py): Variant of [`reduce_accumulate`](fast_frechet/reduce_accumulate.py) using the [Numba library](https://numba.pydata.org/) for JIT compilation of the innermost loop.
 
 Implementations of all these variants can be found under [`fast_frechet/`](fast_frechet/) or by simply clicking on the listed names above.
@@ -40,6 +41,7 @@ $ pre-commit install
 The snippet below estimates the Fréchet distance between the polygonal curves `p` and `q` using the Euclidean distance as a metric to measure distances between points:
 
 ```python
+>>> import numpy as np
 >>> from fast_frechet.linear_memory import frechet_distance
 
 >>> p = np.array([[1, 2], [3, 4]])
@@ -55,13 +57,14 @@ For invoking the [benchmark script](fast_frechet/__main__.py), run:
 $ python fast_frechet
 Length of trajectory = 1024
 
-        no_recursion: 2303 ms
-          vectorized:  553 ms
-          branchless:  508 ms
-       linear_memory:  348 ms
-          accumulate:  286 ms
-   reduce_accumulate:  282 ms
-            compiled:   11 ms
+        no_recursion: 1915 ms
+          vectorized:  495 ms
+          branchless:  466 ms
+       linear_memory:  294 ms
+          accumulate:  258 ms
+   reduce_accumulate:  249 ms
+  reduce_accumulate2:  360 ms
+            compiled:    9 ms
 ```
 (Note that we don't even try to benchmark the [`vanilla`](fast_frechet/vanilla.py) version here, as it already crashes for polygonal curves with a few hundred points due to its recursive nature.)
 

fast_frechet/__main__.py

@@ -10,6 +10,7 @@
     linear_memory,
     no_recursion,
     reduce_accumulate,
+    reduce_accumulate2,
     vectorized,
 )
 
@@ -43,6 +44,7 @@ def main(*, n=1024, seed=42):
         linear_memory,
         accumulate,
         reduce_accumulate,
+        reduce_accumulate2,
         compiled,
     ]:
         f = partial(v.frechet_distance, metric=metric)

fast_frechet/reduce_accumulate2.py

@@ -0,0 +1,25 @@
+from functools import partial, reduce
+from itertools import accumulate
+
+import numpy as np
+
+
+def frechet_combine(a, b):
+    v1, d1 = a
+    v2, d2 = b
+    v, d = min(v1, v2), max(d1, d2)
+    return d if d2 > v1 else v, v if d1 > v2 else d
+
+
+def frechet_next(v, d):
+    v[1:] = np.minimum(v[:-1], v[1:])
+    v = np.maximum(v, d)
+
+    d[0] = v[0]
+    return [d for _, d in accumulate(zip(v, d), frechet_combine)]
+
+
+def frechet_distance(p, q, metric):
+    d = map(partial(metric, q), p)
+    init = np.maximum.accumulate(next(d))
+    return reduce(frechet_next, d, init)[-1]

tests/test_frechet.py

@@ -8,6 +8,7 @@
     linear_memory,
     no_recursion,
     reduce_accumulate,
+    reduce_accumulate2,
     vanilla,
     vectorized,
 )
@@ -41,6 +42,7 @@ def generate_trajectory(n, *, dim, rng):
         linear_memory,
         accumulate,
         reduce_accumulate,
+        reduce_accumulate2,
         compiled,
     ],
 )
@@ -63,6 +65,7 @@ def test_simple_example(variant):
         linear_memory,
         accumulate,
         reduce_accumulate,
+        reduce_accumulate2,
         compiled,
     ],
 )
@@ -83,3 +86,17 @@ def test_frechet(variant, P, Q, dim, seed):
     dqp = variant.frechet_distance(q, p, metric=f)
     assert dpq == d_exp
     assert dqp == d_exp
+
+
+@pytest.mark.parametrize("seed", range(1_000))
+def test_frechet_combine_associativity(seed):
+    rng = np.random.default_rng(seed)
+
+    vd1, vd2, vd3 = rng.integers(0, 6, size=(3, 2))
+
+    vd1 = np.sort(vd1)[::-1]
+    vd2 = np.sort(vd2)[::-1]
+    vd3 = np.sort(vd3)[::-1]
+
+    f = reduce_accumulate2.frechet_combine
+    assert f(f(vd1, vd2), vd3) == f(vd1, f(vd2, vd3))