Repository: aeon-toolkit/aeon
Pull request #1231 (2024-02-23 16:09:43)

[ENH] Add implementation for sbd distance (from k-Shape)
#### What does this implement/fix? Explain your changes.

Introduces the shape-based distance (SBD) from the k-Shape paper (https://doi.org/10.1145/2949741.2949758) to aeon.

- new `sbd_distance` and `sbd_pairwise_distance` functions with documentation
- updated distance dispatcher functions
- Numba-optimized code

> [!WARNING]
> The current implemenetation computes independent distances for multivariate TS (mean over all channels). 

I compared the implementation of SBD using FFT in Numba's object mode and `np.convolve`:

- time series length < 100: `np.convolve` is faster than FFT
- time series length ~= 100: both are equally fast
- time series length > 100: FFT is significantly faster than `np.convolve` (the speed-up scales with the TS length)

➡️ **always use FFT**

```python
for i in range(1, 6):
    ts1 = np.random.random(10**i)
    ts2 = np.random.random(10**i)
    print(f"With time series size: {10**i}: FFT vs. convolve:")
    %timeit sbd_distance_correlate(ts1, ts2, True)
    %timeit sbd_distance_conv(ts1, ts2, True)
```

```
With time series size: 10: FFT vs. convolve:
90.4 µs ± 5.99 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
5.6 µs ± 519 ns per loop (mean ± std. dev. of 7 runs, 100,000 loops each)

With time series size: 100: FFT vs. convolve:
94.5 µs ± 11.9 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
64.4 µs ± 6.17 µs per loop (mean ± std. dev. of 7 runs, 10,000 loops each)

With time series size: 1000: FFT vs. convolve:
130 µs ± 706 ns per loop (mean ± std. dev. of 7 runs, 10,000 loops each)
1.3 ms ± 9.04 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)

With time series size: 10000: FFT vs. convolve:
659 µs ± 4.93 µs per loop (mean ± std. dev. of 7 runs, 1,000 loops each)
89 ms ± 260 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

With time series size: 100000: FFT vs. convolve:
20.8 ms ± 383 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
31.2 s ± 6.54 s per loop (mean ± std. dev. of 7 runs, 1 loop each)
```

#### Does your contribution introduce a new dependency? If yes, which one?

No.

### PR checklist

##### For all contributions

- [x] I've added myself to the [list of contributors](https://github.com/aeon-toolkit/aeon/blob/main/.all-contributorsrc). Alternatively, you can use the [@all-contributors](https://allcontributors.org/docs/en/bot/usage) bot to do this for you.
- [x] The PR title starts with either [ENH], [MNT], [DOC], [BUG], [REF], [DEP] or [GOV] indicating whether the PR topic is related to enhancement, maintenance, documentation, bugs, refactoring, deprecation or governance.

##### For new estimators and functions

- [x] I've added the estimator to the online [API documentation](https://www.aeon-toolkit.org/en/latest/api_reference.html).



Patch:
diff --git a/.all-contributorsrc b/.all-contributorsrc
index a6898f15..34be435c 100644
--- a/.all-contributorsrc
+++ b/.all-contributorsrc
@@ -2338,7 +2338,9 @@
       "profile": "https://github.com/codelionx",
       "contributions": [
         "bug",
-        "code"
+        "code",
+        "doc",
+        "research"
       ]
     }
   ],
diff --git a/aeon/distances/__init__.py b/aeon/distances/__init__.py
index f3604163..e8915766 100644
--- a/aeon/distances/__init__.py
+++ b/aeon/distances/__init__.py
@@ -63,6 +63,8 @@ __all__ = [
     "shape_dtw_cost_matrix",
     "shape_dtw_alignment_path",
     "shape_dtw_pairwise_distance",
+    "sbd_distance",
+    "sbd_pairwise_distance",
 ]
 
 
@@ -123,6 +125,7 @@ from aeon.distances._msm import (
     msm_distance,
     msm_pairwise_distance,
 )
+from aeon.distances._sbd import sbd_distance, sbd_pairwise_distance
 from aeon.distances._shape_dtw import (
     shape_dtw_alignment_path,
     shape_dtw_cost_matrix,
diff --git a/aeon/distances/_distance.py b/aeon/distances/_distance.py
index f4113e25..4935d69c 100644
--- a/aeon/distances/_distance.py
+++ b/aeon/distances/_distance.py
@@ -49,6 +49,7 @@ from aeon.distances._msm import (
     msm_distance,
     msm_pairwise_distance,
 )
+from aeon.distances._sbd import sbd_distance, sbd_pairwise_distance
 from aeon.distances._shape_dtw import (
     shape_dtw_alignment_path,
     shape_dtw_cost_matrix,
@@ -222,6 +223,8 @@ def distance(
             kwargs.get("itakura_max_slope"),
             kwargs.get("warp_penalty", 1.0),
         )
+    elif metric == "sbd":
+        return sbd_distance(x, y, kwargs.get("standardize", True))
     else:
         if isinstance(metric, Callable):
             return metric(x, y, **kwargs)
@@ -385,6 +388,8 @@ def pairwise_distance(
             kwargs.get("itakura_max_slope"),
             kwargs.get("warp_penalty", 1.0),
         )
+    elif metric == "sbd":
+        return sbd_pairwise_distance(x, y, kwargs.get("standardize", True))
     else:
         if isinstance(metric, Callable):
             return _custom_func_pairwise(x, y, metric, **kwargs)
@@ -758,6 +763,7 @@ def get_distance_function(metric: Union[str, DistanceFunction]) -> DistanceFunct
     'squared'       distances.squared_distance
     'manhattan'     distances.manhattan_distance
     'minkowski'     distances.minkowski_distance
+    'sbd'           distances.sbd_distance
     =============== ========================================
 
     Parameters
@@ -814,6 +820,7 @@ def get_pairwise_distance_function(
     'squared'       distances.squared_pairwise_distance
     'manhattan'     distances.manhattan_pairwise_distance
     'minkowski'     distances.minkowski_pairwise_distance
+    'sbd'           distances.sbd_pairwise_distance
     =============== ========================================
 
     Parameters
@@ -1062,6 +1069,11 @@ DISTANCES = [
         "cost_matrix": shape_dtw_cost_matrix,
         "alignment_path": shape_dtw_alignment_path,
     },
+    {
+        "name": "sbd",
+        "distance": sbd_distance,
+        "pairwise_distance": sbd_pairwise_distance,
+    },
 ]
 
 DISTANCES_DICT = {d["name"]: d for d in DISTANCES}
diff --git a/aeon/distances/_sbd.py b/aeon/distances/_sbd.py
new file mode 100644
index 00000000..c3db181d
--- /dev/null
+++ b/aeon/distances/_sbd.py
@@ -0,0 +1,232 @@
+"""Shape-based distance (SBD) between two time series."""
+
+__maintainer__ = ["codelionx"]
+
+from typing import Optional
+
+import numpy as np
+from numba import njit, objmode
+from scipy.signal import correlate
+
+from aeon.distances._utils import reshape_pairwise_to_multiple
+
+
+@njit(cache=True, fastmath=True)
+def sbd_distance(x: np.ndarray, y: np.ndarray, standardize: bool = True) -> float:
+    r"""
+    Compute the shape-based distance (SBD) between two time series.
+
+    Shape-based distance (SBD) [1]_ is a normalized version of
+    cross-correlation (CC) that is shifting and scaling (if standardization
+    is used) invariant.
+
+    For two series, possibly of unequal length, :math:`\mathbf{x}=\{x_1,x_2,\ldots,
+    x_n\}` and :math:`\mathbf{y}=\{y_1,y_2, \ldots,y_m\}`, SBD works by (optionally)
+    first standardizing both time series using the z-score
+    (:math:`x' = \frac{x - \mu}{\sigma}`), then computing the cross-correlation
+    between x and y (:math:`CC(\mathbf{x}, \mathbf{y})`), then deviding it by the
+    geometric mean of both autocorrelations of the individual sequences to normalize
+    it to :math:`[-1, 1]` (coefficient normalization), and finally detecting the
+    position with the maximum normalized cross-correlation:
+
+    .. math::
+        SBD(\mathbf{x}, \mathbf{y}) = 1 - max_w\left( \frac{
+            CC_w(\mathbf{x}, \mathbf{y})
+        }{
+            \sqrt{ (\mathbf{x} \cdot \mathbf{x}) * (\mathbf{y} \cdot \mathbf{y}) }
+        }\right)
+
+    This distance measure has values between 0 and 2; 0 is perfect similarity.
+
+    The computation of the cross-correlation :math:`CC(\mathbf{x}, \mathbf{y})` for
+    all values of w requires :math:`O(m^2)` time, where m is the maximum time-series
+    length. We can however use the convolution theorem to our advantage, and use the
+    fast (inverse) fourier transform (FFT) to perform the computation of
+    :math:`CC(\mathbf{x}, \mathbf{y})` in :math:`O(m \cdot log(m))`:
+
+    .. math::
+        CC(x, y) = \mathcal{F}^{-1}\{
+            \mathcal{F}(\mathbf{x}) * \mathcal{F}(\mathbf{y})
+        \}
+
+    For multivariate time series, SBD is computed independently for each channel and
+    then averaged.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        First time series, either univariate, shape ``(n_timepoints,)``, or
+        multivariate, shape ``(n_channels, n_timepoints)``.
+    y : np.ndarray
+        Second time series, either univariate, shape ``(n_timepoints,)``, or
+        multivariate, shape ``(n_channels, n_timepoints)``.
+    standardize : bool, default=True
+        Apply z-score to both input time series for standardization before
+        computing the distance. This makes SBD scaling invariant. Default is True.
+
+    Returns
+    -------
+    float
+        SBD distance between x and y.
+
+    Raises
+    ------
+    ValueError
+        If x and y are not 1D or 2D arrays.
+
+    See Also
+    --------
+    :func:`~aeon.distances.sbd_pairwise_distance` : Compute the shape-based
+        distance (SBD) between all pairs of time series.
+
+    References
+    ----------
+    .. [1] Paparrizos, John, and Luis Gravano: Fast and Accurate Time-Series
+           Clustering. ACM Transactions on Database Systems 42, no. 2 (2017):
+           8:1-8:49. https://doi.org/10.1145/3044711.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aeon.distances import sbd_distance
+    >>> x = np.array([[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]])
+    >>> y = np.array([[11, 12, 13, 14, 15, 16, 17, 18, 19, 20]])
+    >>> dist = sbd_distance(x, y)
+    """
+    if x.ndim == 1 and y.ndim == 1:
+        return _univariate_sbd_distance(x, y, standardize)
+    if x.ndim == 2 and y.ndim == 2:
+        if x.shape[0] == y.shape[0] == 1:
+            _x = x.ravel()
+            _y = y.ravel()
+            return _univariate_sbd_distance(_x, _y, standardize)
+        else:
+            # independent
+            nchannels = min(x.shape[0], y.shape[0])
+            distance = 0.0
+            for i in range(nchannels):
+                distance += _univariate_sbd_distance(x[i], y[i], standardize)
+            return distance / nchannels
+
+    raise ValueError("x and y must be 1D or 2D")
+
+
+@njit(cache=True, fastmath=True)
+def sbd_pairwise_distance(
+    x: np.ndarray, y: Optional[np.ndarray] = None, standardize: bool = True
+) -> np.ndarray:
+    """
+    Compute the shape-based distance (SBD) between all pairs of time series.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        A collection of time series instances  of shape ``(n_instances, n_timepoints)``
+        or ``(n_instances, n_channels, n_timepoints)``.
+    y : np.ndarray or None, default=None
+        A single series or a collection of time series of shape ``(m_timepoints,)`` or
+        ``(m_instances, m_timepoints)`` or ``(m_instances, m_channels, m_timepoints)``.
+        If None, then the SBD is calculated between pairwise instances of x.
+    standardize : bool, default=True
+        Apply z-score to both input time series for standardization before
+        computing the distance. This makes SBD scaling invariant. Default is True.
+
+    Returns
+    -------
+    np.ndarray (n_instances, n_instances)
+        SBD matrix between the instances of x (and y).
+
+    Raises
+    ------
+    ValueError
+        If x is not 2D or 3D array when only passing x.
+        If x and y are not 1D, 2D or 3D arrays when passing both x and y.
+
+    See Also
+    --------
+    :func:`~aeon.distances.sbd_distance` : Compute the shape-based distance between
+        two time series.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aeon.distances import sbd_pairwise_distance
+    >>> # Distance between each time series in a collection of time series
+    >>> X = np.array([[[1, 2, 3]],[[4, 5, 6]], [[7, 8, 9]]])
+    >>> sbd_pairwise_distance(X)
+    array([[0., 0., 0.],
+           [0., 0., 0.],
+           [0., 0., 0.]])
+
+    >>> # Distance between two collections of time series
+    >>> X = np.array([[[1, 2, 3]],[[4, 5, 6]], [[7, 8, 9]]])
+    >>> y = np.array([[[11, 12, 13]],[[14, 15, 16]], [[17, 18, 19]]])
+    >>> sbd_pairwise_distance(X, y)
+    array([[0., 0., 0.],
+           [0., 0., 0.],
+           [0., 0., 0.]])
+
+    >>> X = np.array([[[1, 2, 3]],[[4, 5, 6]], [[7, 8, 9]]])
+    >>> y_univariate = np.array([[11, 12, 13],[14, 15, 16], [17, 18, 19]])
+    >>> sbd_pairwise_distance(X, y_univariate)
+    array([[0.],
+           [0.],
+           [0.]])
+    """
+    if y is None:
+        # To self
+        if x.ndim == 3:
+            return _sbd_pairwise_distance_single(x, standardize)
+        elif x.ndim == 2:
+            _X = x.reshape((x.shape[0], 1, x.shape[1]))
+            return _sbd_pairwise_distance_single(_X, standardize)
+        raise ValueError("X must be 2D or 3D")
+
+    _x, _y = reshape_pairwise_to_multiple(x, y)
+    return _sbd_pairwise_distance(_x, _y, standardize)
+
+
+@njit(cache=True, fastmath=True)
+def _sbd_pairwise_distance_single(x: np.ndarray, standardize: bool) -> np.ndarray:
+    n_instances = x.shape[0]
+    distances = np.zeros((n_instances, n_instances))
+
+    for i in range(n_instances):
+        for j in range(i + 1, n_instances):
+            distances[i, j] = sbd_distance(x[i], x[j], standardize)
+            distances[j, i] = distances[i, j]
+
+    return distances
+
+
+@njit(cache=True, fastmath=True)
+def _sbd_pairwise_distance(
+    x: np.ndarray, y: np.ndarray, standardize: bool
+) -> np.ndarray:
+    n_instances = x.shape[0]
+    m_instances = y.shape[0]
+    distances = np.zeros((n_instances, m_instances))
+
+    for i in range(n_instances):
+        for j in range(m_instances):
+            distances[i, j] = sbd_distance(x[i], y[j], standardize)
+    return distances
+
+
+@njit(cache=True, fastmath=True)
+def _univariate_sbd_distance(x: np.ndarray, y: np.ndarray, standardize: bool) -> float:
+    x = x.astype(np.float64)
+    y = y.astype(np.float64)
+
+    if standardize:
+        if x.size == 1 or y.size == 1:
+            return 0.0
+
+        x = (x - np.mean(x)) / np.std(x)
+        y = (y - np.mean(y)) / np.std(y)
+
+    with objmode(a="float64[:]"):
+        a = correlate(x, y, method="fft")
+
+    b = np.sqrt(np.dot(x, x) * np.dot(y, y))
+    return np.abs(1.0 - np.max(a / b))
diff --git a/docs/api_reference/distances.rst b/docs/api_reference/distances.rst
index 9780e827..ef799123 100644
--- a/docs/api_reference/distances.rst
+++ b/docs/api_reference/distances.rst
@@ -128,6 +128,18 @@ Move-Split-Merge (MSM)
     msm_cost_matrix
     msm_alignment_path
 
+Shape-based Distance (SBD)
+--------------------------
+
+.. currentmodule:: aeon.distances
+
+.. autosummary::
+    :toctree: auto_generated/
+    :template: function.rst
+
+    sbd_distance
+    sbd_pairwise_distance
+
 Shape Dynamic Time Warping (Shape DTW)
 --------------------------------------
 


Test patch:
diff --git a/aeon/distances/tests/test_expected_results.py b/aeon/distances/tests/test_expected_results.py
index 2c20ec80..24de5ab2 100644
--- a/aeon/distances/tests/test_expected_results.py
+++ b/aeon/distances/tests/test_expected_results.py
@@ -97,6 +97,13 @@ _expected_distance_results = {
         4.788349892989078,
         28.544576071968976,
     ],
+    "sbd": [
+        0.0,
+        0.33449276,
+        0.53488311,
+        0.59674837,
+        0.51688041,
+    ],
 }
 
 _expected_distance_results_params = {
diff --git a/aeon/distances/tests/test_numba_distance_parameters.py b/aeon/distances/tests/test_numba_distance_parameters.py
index 2eabbe40..e7cd0238 100644
--- a/aeon/distances/tests/test_numba_distance_parameters.py
+++ b/aeon/distances/tests/test_numba_distance_parameters.py
@@ -69,6 +69,7 @@ DIST_PARAMS = {
     "msm": BASIC_BOUNDING_PARAMS + [{"independent": False}, {"c": 0.1}],
     "adtw": BASIC_BOUNDING_PARAMS + [{"warp_penalty": 5.0}],
     "minkowski": [{"p": 1.0}, {"p": 2.0}],
+    "sbd": [{"standardize": False}],
 }