perf(ordinal-patterns): serial path for small data and bypass decorator overhead

delaynet/connectivities/continuous_ordinal_patterns.py

@@ -4,7 +4,10 @@
 from numba import njit, prange
 
 from ..decorators import connectivity
-from .granger import gt_multi_lag
+from ..utils.lag_steps import find_optimal_lag
+from .granger import gt_multi_lag, gt_single_lag
+
+_SERIAL_THRESHOLD = 200
 
 
 @connectivity
@@ -53,7 +56,9 @@ def random_patterns(
     t_ts2 = pattern_transform(np.copy(ts2), rnd_patterns)
 
     for i in range(num_rnd_patterns):
-        p_v, pv_idx = gt_multi_lag(t_ts1[i, :], t_ts2[i, :], lag_steps=lag_steps)
+        p_v, pv_idx = find_optimal_lag(
+            gt_single_lag, t_ts1[i, :], t_ts2[i, :], lag_steps=lag_steps
+        )
         if best_pv > p_v:
             best_pv = p_v
             best_lag = pv_idx
@@ -99,7 +104,6 @@ def pattern_transform(ts: np.ndarray, patterns: np.ndarray) -> np.ndarray:
     return pattern_transform_2d(ts, patterns).squeeze()
 
 
-@njit(nogil=True, parallel=True)
 def pattern_transform_2d(
     ts: np.ndarray, patterns: np.ndarray
 ) -> np.ndarray:  # pragma: no cover
@@ -128,34 +132,74 @@ def pattern_transform_2d(
     return transformed
 
 
-@njit(nogil=True, parallel=True)
-def norm_windows(ts: np.ndarray, window_size: int) -> np.ndarray:  # pragma: no cover
-    """Normalise sliding windows of a time series to values between -1 and 1.
+@njit(nogil=True, cache=True)
+def _norm_windows_serial(
+    ts: np.ndarray, window_size: int
+) -> np.ndarray:  # pragma: no cover
+    windows = np.lib.stride_tricks.as_strided(
+        x=ts,
+        strides=(ts.strides[0], ts.strides[0]),
+        shape=(ts.shape[0] - window_size + 1, window_size),
+    )
+    normed_windows = np.zeros_like(windows)
+    for i in range(windows.shape[0]):
+        normed_windows[i] = norm_window(windows[i])
+    return normed_windows
 
-    :param ts: Time series.
-    :type ts: numpy.ndarray
-    :param window_size: Size of the window.
-    :type window_size: int
-    :return: Normalised windows.
-    :rtype: numpy.ndarray
-    """
-    # Create a sliding window view of the input array
-    # windows = np.lib.stride_tricks.sliding_window_view(
-    #     x=ts, window_shape=window_size, writeable=False
-    # )
+
+@njit(nogil=True, parallel=True, cache=True)
+def _norm_windows_parallel(
+    ts: np.ndarray, window_size: int
+) -> np.ndarray:  # pragma: no cover
     windows = np.lib.stride_tricks.as_strided(
         x=ts,
         strides=(ts.strides[0], ts.strides[0]),
         shape=(ts.shape[0] - window_size + 1, window_size),
     )
     normed_windows = np.zeros_like(windows)
-    # Normalise each window to [-1, 1]
     for i in prange(windows.shape[0]):
         normed_windows[i] = norm_window(windows[i])
     return normed_windows
 
 
-@njit(nogil=True, parallel=True)
+def norm_windows(ts: np.ndarray, window_size: int) -> np.ndarray:
+    """Normalise sliding windows of a time series to values between -1 and 1.
+
+    :param ts: Time series.
+    :type ts: numpy.ndarray
+    :param window_size: Size of the window.
+    :type window_size: int
+    :return: Normalised windows.
+    :rtype: numpy.ndarray
+    """
+    n_windows = ts.shape[0] - window_size + 1
+    if n_windows < _SERIAL_THRESHOLD:
+        return _norm_windows_serial(ts, window_size)
+    return _norm_windows_parallel(ts, window_size)
+
+
+@njit(nogil=True, cache=True)
+def _pattern_distance_serial(
+    windows: np.ndarray, pattern: np.ndarray
+) -> np.ndarray:  # pragma: no cover
+    distances = np.zeros(windows.shape[0])
+    for i in range(windows.shape[0]):
+        for j in range(pattern.shape[0]):
+            distances[i] += np.abs(windows[i, j] - pattern[j])
+    return distances / pattern.shape[0] / 2.0
+
+
+@njit(nogil=True, parallel=True, cache=True)
+def _pattern_distance_parallel(
+    windows: np.ndarray, pattern: np.ndarray
+) -> np.ndarray:  # pragma: no cover
+    distances = np.zeros(windows.shape[0])
+    for i in prange(windows.shape[0]):
+        for j in prange(pattern.shape[0]):
+            distances[i] += np.abs(windows[i, j] - pattern[j])
+    return distances / pattern.shape[0] / 2.0
+
+
 def pattern_distance(
     windows: np.ndarray, pattern: np.ndarray
 ) -> np.ndarray:  # pragma: no cover
@@ -168,9 +212,6 @@ def pattern_distance(
     :return: Distance between the windows and the pattern.
     :rtype: numpy.ndarray
     """
-    distances = np.zeros(windows.shape[0])
-    for i in prange(windows.shape[0]):
-        for j in prange(pattern.shape[0]):
-            distances[i] += np.abs(windows[i, j] - pattern[j])
-    return distances / pattern.shape[0] / 2.0
-    # equiv. to np.sum(np.abs(windows - pattern), axis=1) / pattern.shape[0] / 2.0
+    if windows.shape[0] < _SERIAL_THRESHOLD:
+        return _pattern_distance_serial(windows, pattern)
+    return _pattern_distance_parallel(windows, pattern)

tests/connectivities/test_continuous_ordinal_patterns.py

@@ -1,17 +1,53 @@
 """Tests for the continuous ordinal patterns connectivity measure."""
 
+import random as _random
+import numpy as np
 import pytest
 from numpy import array, allclose, linspace, random, roll
+from numpy.lib.stride_tricks import as_strided
 
 from delaynet import connectivity
 from delaynet.connectivities.continuous_ordinal_patterns import (
+    _norm_windows_parallel,
+    _norm_windows_serial,
+    _pattern_distance_parallel,
+    _pattern_distance_serial,
     norm_window,
     norm_windows,
     pattern_distance,
     pattern_transform,
     random_patterns,
 )
 
+_SEED = _random.randint(0, 10000)  # deterministic per-session
+
+
+def _norm_windows_expected(ts, window_size):
+    """Pure-numpy reference for sliding window normalisation (no numba)."""
+    windows = as_strided(
+        ts,
+        strides=(ts.strides[0], ts.strides[0]),
+        shape=(ts.shape[0] - window_size + 1, window_size),
+    )
+    result = np.zeros_like(windows)
+    for i in range(windows.shape[0]):
+        w = windows[i]
+        nw = w - w.min()
+        nw = nw / nw.max()
+        nw = (nw - 0.5) * 2.0
+        nw[np.isnan(nw)] = 0.0
+        result[i] = nw
+    return result
+
+
+def _pattern_distance_expected(windows, pattern):
+    """Pure-numpy reference for pattern distance (no numba)."""
+    distances = np.zeros(windows.shape[0])
+    for i in range(windows.shape[0]):
+        for j in range(pattern.shape[0]):
+            distances[i] += np.abs(windows[i, j] - pattern[j])
+    return distances / pattern.shape[0] / 2.0
+
 
 @pytest.mark.parametrize(
     "ts, expected",
@@ -58,6 +94,26 @@ def test_norm_windows(ts, window_size, expected):
     assert allclose(norm_windows(array(ts), window_size), array(expected))
 
 
+@pytest.mark.parametrize("window_size", [2, 5, 10])
+def test_norm_windows_random_noise(window_size):
+    """norm_windows on non-monotonic random data matches independent reference."""
+    rng = np.random.default_rng(_SEED)
+    ts = rng.normal(0, 1, size=200)
+    assert allclose(
+        norm_windows(ts, window_size), _norm_windows_expected(ts, window_size)
+    )
+
+
+def test_norm_windows_serial_parallel_equivalence():
+    """Serial and parallel paths produce identical output on large data."""
+    rng = np.random.default_rng(_SEED + 1)
+    ts = rng.normal(0, 1, size=500)
+    window_size = 2  # 499 windows > _SERIAL_THRESHOLD = 200
+    serial = _norm_windows_serial(ts, window_size)
+    parallel = _norm_windows_parallel(ts, window_size)
+    assert allclose(serial, parallel)
+
+
 @pytest.mark.parametrize(
     "windows, pattern, expected",
     [
@@ -73,6 +129,25 @@ def test_pattern_distance(windows, pattern, expected):
     assert allclose(pattern_distance(array(windows), array(pattern)), array(expected))
 
 
+def test_pattern_distance_random_noise():
+    """pattern_distance on random data matches independent reference."""
+    rng = np.random.default_rng(_SEED + 2)
+    windows = rng.uniform(-1, 1, (50, 5))
+    pattern = rng.uniform(-1, 1, (5,))
+    expected = _pattern_distance_expected(windows, pattern)
+    assert allclose(pattern_distance(windows, pattern), expected)
+
+
+def test_pattern_distance_serial_parallel_equivalence():
+    """Serial and parallel paths produce identical output on large data."""
+    rng = np.random.default_rng(_SEED + 3)
+    windows = rng.uniform(-1, 1, (500, 3))
+    pattern = rng.uniform(-1, 1, (3,))
+    serial = _pattern_distance_serial(windows, pattern)
+    parallel = _pattern_distance_parallel(windows, pattern)
+    assert allclose(serial, parallel)
+
+
 @pytest.mark.parametrize(
     "ts, patterns, expected",
     [