Repository: aeon-toolkit/aeon
Pull request #322 (2023-04-23 01:55:01)

[ENH] Update wdtw distance
<!--
Thanks for contributing a pull request! Please ensure you have taken a look at our contribution guide: https://github.com/aeon-toolkit/aeon/blob/main/CONTRIBUTING.md

Feel free to delete sections of this template if they do not apply to your PR, avoid submitting a blank template or empty sections.
-->

#### Reference Issues/PRs
This pr depends on pr #319 

<!--
Example: Fixes #1234. See also #3456.

Please use keywords (e.g., Fixes) to create link to the issues or pull requests you resolved, so that they will automatically be closed when your pull request is merged. See https://github.com/blog/1506-closing-issues-via-pull-requests
-->

#### What does this implement/fix? Explain your changes.
This pr updates wdtw to the new format for distances. In addition in doing so it makes optimisations for computing the distances between multiple time series. 

<!--
A clear and concise description of what you have implemented.
-->

#### Does your contribution introduce a new dependency? If yes, which one?

<!--
If your contribution does add a dependency, we may suggest adding it as an optional/soft dependency to keep external dependencies of the core aeon package to a minimum.
-->

#### What should a reviewer concentrate their feedback on?

<!--
This section is particularly useful if you have a pull request that is still in development. You can guide the reviews to focus on the parts that are ready for their comments. We suggest using bullets (indicated by * or -) and filled checkboxes [x] here
-->

#### Did you add any tests for the change?

<!--
This section is useful if you have added a test in addition to the existing ones. This will ensure that further changes to these files won't introduce the same kind of bug. It is considered good practice to add tests with newly added code to enforce the fact that the code actually works. This will reduce the chance of introducing logical bugs.
-->

#### Any other comments?

<!--
Please be aware that we are a loose team of volunteers so patience is necessary; assistance handling other issues is very welcome. We value all user contributions, no matter how minor they are. If we are slow to review, either the pull request needs some benchmarking, tinkering, convincing, etc. or more likely the reviewers are simply busy. In either case, we ask for your understanding during the review process.
-->

### PR checklist

<!--
Please go through the checklist below. Please feel free to remove points if they are not applicable.
-->

##### For all contributions
- [x] I've added myself to the [list of contributors](https://github.com/aeon-toolkit/aeon/blob/main/.all-contributorsrc).
- [x] Optionally, I've updated aeon's [CODEOWNERS](https://github.com/aeon-toolkit/aeon/blob/main/CODEOWNERS) to receive notifications about future changes to these files.
- [x] The PR title starts with either [ENH], [MNT], [DOC], or [BUG] indicating whether the PR topic is related to enhancement, maintenance, documentation, or bug.
- [x] For new estimators, I've added the estimator to the online documentation.


<!--
Thanks for contributing!
-->



Patch:
diff --git a/aeon/distances/__init__.py b/aeon/distances/__init__.py
index 03149473..62e99cf5 100644
--- a/aeon/distances/__init__.py
+++ b/aeon/distances/__init__.py
@@ -15,19 +15,21 @@ __all__ = [
     "ddtw_pairwise_distance",
     "ddtw_alignment_path",
     "ddtw_cost_matrix",
+    "wdtw_distance",
+    "wdtw_pairwise_distance",
+    "wdtw_cost_matrix",
+    "wdtw_alignment_path",
     "distance",
     "distance_factory",
     "pairwise_distance",
     "euclidean_distance",
     "squared_distance",
-    "wdtw_distance",
     "wddtw_distance",
     "edr_distance",
     "erp_distance",
     "msm_distance",
     "lcss_distance",
     "twe_distance",
-    "wdtw_alignment_path",
     "lcss_alignment_path",
     "msm_alignment_path",
     "erp_alignment_path",
@@ -63,8 +65,6 @@ from aeon.distances._distance import (
     twe_distance,
     wddtw_alignment_path,
     wddtw_distance,
-    wdtw_alignment_path,
-    wdtw_distance,
 )
 from aeon.distances._dtw import (
     dtw_alignment_path,
@@ -74,3 +74,9 @@ from aeon.distances._dtw import (
 )
 from aeon.distances._euclidean import euclidean_distance, euclidean_pairwise_distance
 from aeon.distances._squared import squared_distance, squared_pairwise_distance
+from aeon.distances._wdtw import (
+    wdtw_alignment_path,
+    wdtw_cost_matrix,
+    wdtw_distance,
+    wdtw_pairwise_distance,
+)
diff --git a/aeon/distances/_distance.py b/aeon/distances/_distance.py
index 685ee0ed..3b8e7248 100644
--- a/aeon/distances/_distance.py
+++ b/aeon/distances/_distance.py
@@ -31,7 +31,11 @@ from aeon.distances._resolve_metric import (
 from aeon.distances._squared import squared_distance, squared_pairwise_distance
 from aeon.distances._twe import _TweDistance
 from aeon.distances._wddtw import _WddtwDistance
-from aeon.distances._wdtw import _WdtwDistance
+from aeon.distances._wdtw import (
+    wdtw_alignment_path,
+    wdtw_distance,
+    wdtw_pairwise_distance,
+)
 from aeon.distances.base import (
     AlignmentPathReturn,
     DistanceAlignmentPathCallable,
@@ -350,88 +354,6 @@ def wddtw_distance(
     return distance(x, y, metric="wddtw", **format_kwargs)
 
 
-def wdtw_distance(
-    x: np.ndarray,
-    y: np.ndarray,
-    window: Union[float, None] = None,
-    g: float = 0.05,
-    **kwargs: Any,
-) -> float:
-    """Compute the weighted dynamic time warping (WDTW) distance between time series.
-
-    First proposed in [1]_, WDTW adds a  adds a multiplicative weight penalty based on
-    the warping distance. This means that time series with lower phase difference have
-    a smaller weight imposed (i.e less penalty imposed) and time series with larger
-    phase difference have a larger weight imposed (i.e. larger penalty imposed).
-
-    Formally this can be described as:
-
-    .. math::
-        d_{w}(x_{i}, y_{j}) = ||w_{|i-j|}(x_{i} - y_{j})||
-
-    Where d_w is the distance with a the weight applied to it for points i, j, where
-    w(|i-j|) is a positive weight between the two points x_i and y_j.
-
-
-    Parameters
-    ----------
-    x: np.ndarray (1d or 2d array)
-        First time series.
-    y: np.ndarray (1d or 2d array)
-        Second time series.
-    window: float, defaults = None
-        Float that is the radius of the sakoe chiba window (if using Sakoe-Chiba
-        lower bounding). Value must be between 0. and 1.
-    g: float, defaults = 0.
-        Constant that controls the curvature (slope) of the function; that is, g
-        controls the level of penalisation for the points with larger phase
-        difference.
-    kwargs: Any
-        Extra kwargs.
-
-    Returns
-    -------
-    float
-        Wdtw distance between the x and y.
-
-    Raises
-    ------
-    ValueError
-        If the sakoe_chiba_window_radius is not a float.
-        If the value of x or y provided is not a numpy array.
-        If the value of x or y has more than 2 dimensions.
-        If a metric string provided, and is not a defined valid string.
-        If a metric object (instance of class) is provided and doesn't inherit from
-        NumbaDistance.
-        If the metric type cannot be determined
-
-    Examples
-    --------
-    >>> x_1d = np.array([1, 2, 3, 4])  # 1d array
-    >>> y_1d = np.array([5, 6, 7, 8])  # 1d array
-    >>> wdtw_distance(x_1d, y_1d)
-    27.975712863958133
-
-    >>> x_2d = np.array([[1, 2, 3, 4], [5, 6, 7, 8]])  # 2d array
-    >>> y_2d = np.array([[9, 10, 11, 12], [13, 14, 15, 16]])  # 2d array
-    >>> wdtw_distance(x_2d, y_2d)
-    243.2106560107827
-
-    References
-    ----------
-    .. [1] Young-Seon Jeong, Myong K. Jeong, Olufemi A. Omitaomu, Weighted dynamic time
-    warping for time series classification, Pattern Recognition, Volume 44, Issue 9,
-    2011, Pages 2231-2240, ISSN 0031-3203, https://doi.org/10.1016/j.patcog.2010.09.022.
-    """
-    format_kwargs = {
-        "window": window,
-        "g": g,
-    }
-    format_kwargs = {**format_kwargs, **kwargs}
-
-    return distance(x, y, metric="wdtw", **format_kwargs)
-
-
 def msm_distance(
     x: np.ndarray,
     y: np.ndarray,
@@ -577,86 +499,6 @@ def twe_distance(
     return distance(x, y, metric="twe", **format_kwargs)
 
 
-def wdtw_alignment_path(
-    x: np.ndarray,
-    y: np.ndarray,
-    return_cost_matrix: bool = False,
-    window: Union[float, None] = None,
-    g: float = 0.05,
-    **kwargs: Any,
-) -> AlignmentPathReturn:
-    """Compute the weighted dynamic time warping (wdtw) alignment path.
-
-    First proposed in [1]_, WDTW adds a  adds a multiplicative weight penalty based on
-    the warping distance. This means that time series with lower phase difference have
-    a smaller weight imposed (i.e less penalty imposed) and time series with larger
-    phase difference have a larger weight imposed (i.e. larger penalty imposed).
-
-    Formally this can be described as:
-
-    .. math::
-        d_{w}(x_{i}, y_{j}) = ||w_{|i-j|}(x_{i} - y_{j})||
-
-    Where d_w is the distance with a the weight applied to it for points i, j, where
-    w(|i-j|) is a positive weight between the two points x_i and y_j.
-
-
-    Parameters
-    ----------
-    x: np.ndarray (1d or 2d array)
-        First time series.
-    y: np.ndarray (1d or 2d array)
-        Second time series.
-    return_cost_matrix: bool, defaults = False
-        Boolean that when true will also return the cost matrix.
-    window: float, defaults = None
-        Float that is the radius of the sakoe chiba window (if using Sakoe-Chiba
-        lower bounding). Value must be between 0. and 1.
-    g: float, defaults = 0.
-        Constant that controls the curvature (slope) of the function; that is, g
-        controls the level of penalisation for the points with larger phase
-        difference.
-    kwargs: Any
-        Extra kwargs.
-
-    Returns
-    -------
-    list[tuple]
-        List of tuples containing the wdtw alignment path.
-    float
-        Wdtw distance between x and y.
-    np.ndarray (of shape (len(x), len(y)).
-        Optional return only given if return_cost_matrix = True.
-        Cost matrix used to compute the distance.
-
-    Raises
-    ------
-    ValueError
-        If the sakoe_chiba_window_radius is not a float.
-        If the value of x or y provided is not a numpy array.
-        If the value of x or y has more than 2 dimensions.
-        If a metric string provided, and is not a defined valid string.
-        If a metric object (instance of class) is provided and doesn't inherit from
-        NumbaDistance.
-        If the metric type cannot be determined
-
-    References
-    ----------
-    .. [1] Young-Seon Jeong, Myong K. Jeong, Olufemi A. Omitaomu, Weighted dynamic time
-    warping for time series classification, Pattern Recognition, Volume 44, Issue 9,
-    2011, Pages 2231-2240, ISSN 0031-3203, https://doi.org/10.1016/j.patcog.2010.09.022.
-    """
-    format_kwargs = {
-        "window": window,
-        "g": g,
-    }
-    format_kwargs = {**format_kwargs, **kwargs}
-
-    return distance_alignment_path(
-        x, y, metric="wdtw", return_cost_matrix=return_cost_matrix, **format_kwargs
-    )
-
-
 def wddtw_alignment_path(
     x: np.ndarray,
     y: np.ndarray,
@@ -1121,7 +963,7 @@ def twe_alignment_path(
     )
 
 
-NEW_DISTANCES = ["squared", "euclidean", "dtw", "ddtw"]
+NEW_DISTANCES = ["squared", "euclidean", "dtw", "ddtw", "wdtw"]
 
 
 def distance(
@@ -1212,6 +1054,8 @@ def distance(
             return dtw_distance(x, y, **kwargs)
         elif metric == "ddtw":
             return ddtw_distance(x, y, **kwargs)
+        elif metric == "wdtw":
+            return wdtw_distance(x, y, **kwargs)
     _x = to_numba_timeseries(x)
     _y = to_numba_timeseries(y)
 
@@ -1287,6 +1131,8 @@ def distance_factory(
             return dtw_distance
         elif metric == "ddtw":
             return ddtw_distance
+        elif metric == "wdtw":
+            return wdtw_distance
     global dist_callable
 
     if x is None:
@@ -1411,6 +1257,8 @@ def pairwise_distance(
             return dtw_pairwise_distance(_x, _y, **kwargs)
         elif metric == "ddtw":
             return ddtw_pairwise_distance(_x, _y, **kwargs)
+        elif metric == "wdtw":
+            return wdtw_pairwise_distance(_x, _y, **kwargs)
 
     symmetric = np.array_equal(_x, _y)
     _metric_callable = _resolve_metric_to_factory(
@@ -1493,6 +1341,8 @@ def distance_alignment_path(
             return dtw_alignment_path(x, y, **kwargs)
         elif metric == "ddtw":
             return ddtw_alignment_path(x, y, **kwargs)
+        elif metric == "wdtw":
+            return wdtw_alignment_path(x, y, **kwargs)
     _x = to_numba_timeseries(x)
     _y = to_numba_timeseries(y)
 
@@ -1572,6 +1422,8 @@ def distance_alignment_path_factory(
             return dtw_alignment_path
         if metric == "ddtw":
             return ddtw_alignment_path
+        elif metric == "wdtw":
+            return wdtw_alignment_path
     if x is None:
         x = np.zeros((1, 10))
     if y is None:
@@ -1615,13 +1467,6 @@ _METRIC_INFOS = [
         dist_instance=_LcssDistance(),
         dist_alignment_path_func=lcss_alignment_path,
     ),
-    MetricInfo(
-        canonical_name="wdtw",
-        aka={"wdtw", "weighted dynamic time warping"},
-        dist_func=wdtw_distance,
-        dist_instance=_WdtwDistance(),
-        dist_alignment_path_func=wdtw_alignment_path,
-    ),
     MetricInfo(
         canonical_name="wddtw",
         aka={"wddtw", "weighted derivative dynamic time warping"},
@@ -1658,6 +1503,5 @@ ALL_DISTANCES = (
     lcss_distance,
     msm_distance,
     wddtw_distance,
-    wdtw_distance,
     twe_distance,
 )
diff --git a/aeon/distances/_dtw.py b/aeon/distances/_dtw.py
index 2f7116a1..2545f466 100644
--- a/aeon/distances/_dtw.py
+++ b/aeon/distances/_dtw.py
@@ -229,6 +229,8 @@ def dtw_pairwise_distance(
     np.ndarray (n_instances, n_instances)
         dtw pairwise matrix between the instances of X.
 
+    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.
diff --git a/aeon/distances/_wddtw.py b/aeon/distances/_wddtw.py
index 2c9589d8..6a79c12b 100644
--- a/aeon/distances/_wddtw.py
+++ b/aeon/distances/_wddtw.py
@@ -11,7 +11,7 @@ from numba.core.errors import NumbaWarning
 from aeon.distances._alignment_paths import compute_min_return_path
 from aeon.distances._bounding_matrix import create_bounding_matrix
 from aeon.distances._ddtw import average_of_slope
-from aeon.distances._wdtw import _weighted_cost_matrix
+from aeon.distances._wdtw import _wdtw_cost_matrix
 from aeon.distances.base import (
     DistanceAlignmentPathCallable,
     DistanceCallable,
@@ -95,7 +95,7 @@ class _WddtwDistance(NumbaDistance):
             ) -> Tuple[List, float, np.ndarray]:
                 _x = compute_derivative(_x)
                 _y = compute_derivative(_y)
-                cost_matrix = _weighted_cost_matrix(_x, _y, _bounding_matrix, g)
+                cost_matrix = _wdtw_cost_matrix(_x, _y, _bounding_matrix, g)
                 path = compute_min_return_path(cost_matrix)
                 return path, cost_matrix[-1, -1], cost_matrix
 
@@ -108,7 +108,7 @@ class _WddtwDistance(NumbaDistance):
             ) -> Tuple[List, float]:
                 _x = compute_derivative(_x)
                 _y = compute_derivative(_y)
-                cost_matrix = _weighted_cost_matrix(_x, _y, _bounding_matrix, g)
+                cost_matrix = _wdtw_cost_matrix(_x, _y, _bounding_matrix, g)
                 path = compute_min_return_path(cost_matrix)
                 return path, cost_matrix[-1, -1]
 
@@ -175,7 +175,7 @@ class _WddtwDistance(NumbaDistance):
         ) -> float:
             _x = compute_derivative(_x)
             _y = compute_derivative(_y)
-            cost_matrix = _weighted_cost_matrix(_x, _y, _bounding_matrix, g)
+            cost_matrix = _wdtw_cost_matrix(_x, _y, _bounding_matrix, g)
             return cost_matrix[-1, -1]
 
         return numba_wddtw_distance
diff --git a/aeon/distances/_wdtw.py b/aeon/distances/_wdtw.py
index 9540fd0e..c42b5e6b 100644
--- a/aeon/distances/_wdtw.py
+++ b/aeon/distances/_wdtw.py
@@ -1,225 +1,412 @@
 # -*- coding: utf-8 -*-
+r"""Weighted dynamic time warping (WDTW) distance between two time series.
+
+WDTW uses DTW with a weighted pairwise distance matrix rather than a window. When
+creating the distance matrix :math:'M', a weight penalty  :math:'w_{|i-j|}' for a
+warping distance of :math:'|i-j|' is applied, so that for series
+:math:'a = <a_1, ..., a_m>' and :math:'b=<b_1,...,b_m>',
+.. math:: M_{i,j}=  w(|i-j|) (a_i-b_j)^2.
+A logistic weight function, proposed in [1] is used, so that a warping of :math:'x'
+places imposes a weighting of
+.. math:: w(x)=\frac{w_{max}}{1+e^{-g(x-m/2)}},
+where :math:'w_{max}' is an upper bound on the weight (set to 1), :math:'m' is
+the series length and :math:'g' is a parameter that controls the penalty level
+for larger warpings. The greater :math:'g' is, the greater the penalty for warping.
+Once :math:'M' is found, standard dynamic time warping is applied.
+
+WDTW is set up so you can use it with a bounding box in addition to the weight
+function is so desired. This is for consistency with the other distance functions.
+
+References
+----------
+.. [1] Jeong, Y., Jeong, M., Omitaomu, O.: Weighted dynamic time warping for time
+series classification. Pattern Recognition 44, 2231–2240 (2011)
+"""
 __author__ = ["chrisholder", "TonyBagnall"]
 
-import warnings
-from typing import Any, List, Tuple
+from typing import List, Tuple
 
 import numpy as np
 from numba import njit
-from numba.core.errors import NumbaWarning
 
 from aeon.distances._alignment_paths import compute_min_return_path
 from aeon.distances._bounding_matrix import create_bounding_matrix
-from aeon.distances.base import (
-    DistanceAlignmentPathCallable,
-    DistanceCallable,
-    NumbaDistance,
-)
-
-# Warning occurs when using large time series (i.e. 1000x1000)
-warnings.simplefilter("ignore", category=NumbaWarning)
-
-
-class _WdtwDistance(NumbaDistance):
-    r"""Weighted dynamic time warping (WDTW) distance between two time series.
-
-    WDTW uses DTW with a weighted pairwise distance matrix rather than a window. When
-    creating the distance matrix :math:'M', a weight penalty  :math:'w_{|i-j|}' for a
-    warping distance of :math:'|i-j|' is applied, so that for series
-    :math:'a = <a_1, ..., a_m>' and :math:'b=<b_1,...,b_m>',
-    .. math:: M_{i,j}=  w(|i-j|) (a_i-b_j)^2.
-    A logistic weight function, proposed in [1] is used, so that a warping of :math:'x'
-    places imposes a weighting of
-    .. math:: w(x)=\frac{w_{max}}{1+e^{-g(x-m/2)}},
-    where :math:'w_{max}' is an upper bound on the weight (set to 1), :math:'m' is
-    the series length and :math:'g' is a parameter that controls the penalty level
-    for larger warpings. The greater :math:'g' is, the greater the penalty for warping.
-    Once :math:'M' is found, standard dynamic time warping is applied.
-
-    WDTW is set up so you can use it with a bounding box in addition to the weight
-    function is so desired. This is for consistency with the other distance functions.
+from aeon.distances._squared import _univariate_squared_distance
 
-    References
-    ----------
-    .. [1] Jeong, Y., Jeong, M., Omitaomu, O.: Weighted dynamic time warping for time
-    series classification. Pattern Recognition 44, 2231–2240 (2011)
-    """
 
-    def _distance_alignment_path_factory(
-        self,
-        x: np.ndarray,
-        y: np.ndarray,
-        return_cost_matrix: bool = False,
-        window: int = None,
-        g: float = 0.05,
-        **kwargs: Any,
-    ) -> DistanceAlignmentPathCallable:
-        """Create a no_python compiled wdtw distance alignment path callable.
-
-        Series should be shape (d, m), where d is the number of dimensions, m the series
-        length. Series can be different lengths.
-
-        Parameters
-        ----------
-        x: np.ndarray (2d array of shape (d,m1)).
-            First time series.
-        y: np.ndarray (2d array of shape (d,m2)).
-            Second time series.
-        return_cost_matrix: bool, defaults = False
-            Boolean that when true will also return the cost matrix.
-        window: Float, defaults = None
-            Float that is the radius of the sakoe chiba window (if using Sakoe-Chiba
-            lower bounding). Must be between 0 and 1.
-        g: float, defaults = 0.
-            Constant that controls the curvature (slope) of the function; that is, g
-            controls the level of penalisation for the points with larger phase
-            difference.
-        kwargs: any
-            extra kwargs.
-
-        Returns
-        -------
-        Callable[[np.ndarray, np.ndarray], tuple[np.ndarray, float]]
-            No_python compiled wdtw distance path callable.
-
-        Raises
-        ------
-        ValueError
-            If the input time series are not numpy array.
-            If the input time series do not have exactly 2 dimensions.
-            If the value of g is not a float
-        """
-        if not isinstance(g, float):
-            raise ValueError(
-                f"The value of g must be a float. The current value is {g}"
-            )
+@njit(cache=True, fastmath=True)
+def wdtw_distance(
+    x: np.ndarray, y: np.ndarray, window: float = None, g: float = 0.05
+) -> float:
+    """Compute the wdtw distance between two time series.
 
-        _bounding_matrix = create_bounding_matrix(x.shape[1], y.shape[1], window)
+    First proposed in [1]_, WDTW adds a  adds a multiplicative weight penalty based on
+    the warping distance. This means that time series with lower phase difference have
+    a smaller weight imposed (i.e less penalty imposed) and time series with larger
+    phase difference have a larger weight imposed (i.e. larger penalty imposed).
 
-        if return_cost_matrix is True:
-
-            @njit(cache=True)
-            def numba_wdtw_distance_alignment_path(
-                _x: np.ndarray,
-                _y: np.ndarray,
-            ) -> Tuple[List, float, np.ndarray]:
-                cost_matrix = _weighted_cost_matrix(_x, _y, _bounding_matrix, g)
-                path = compute_min_return_path(cost_matrix)
-                return path, cost_matrix[-1, -1], cost_matrix
-
-        else:
+    Formally this can be described as:
 
-            @njit(cache=True)
-            def numba_wdtw_distance_alignment_path(
-                _x: np.ndarray,
-                _y: np.ndarray,
-            ) -> Tuple[List, float]:
-                cost_matrix = _weighted_cost_matrix(_x, _y, _bounding_matrix, g)
-                path = compute_min_return_path(cost_matrix)
-                return path, cost_matrix[-1, -1]
-
-        return numba_wdtw_distance_alignment_path
-
-    def _distance_factory(
-        self,
-        x: np.ndarray,
-        y: np.ndarray,
-        window: int = None,
-        g: float = 0.05,
-        **kwargs: Any,
-    ) -> DistanceCallable:
-        """Create a no_python compiled wdtw distance callable.
-
-        Parameters
-        ----------
-        x: np.ndarray (2d array of shape (d,m1)).
-            First time series.
-        y: np.ndarray (2d array of shape (d,m2)).
-            Second time series.
-        window: float, defaults = None
-            Integer that is the radius of the sakoe chiba window (if using Sakoe-Chiba
-            lower bounding). Must be between 0 and 1.
-        g: float, defaults = 0.
-            Constant that controls the curvature (slope) of the function; that is, g
-            controls the level of penalisation for the points with larger phase
-            difference.
-        kwargs: Any
-            Extra kwargs.
-
-
-        Returns
-        -------
-        Callable[[np.ndarray, np.ndarray], float]
-            No_python compiled wdtw distance callable.
-
-        Raises
-        ------
-        ValueError
-            If the input time series are not numpy array.
-            If the input time series do not have exactly 2 dimensions.
-            If the value of g is not a float
-        """
-        _bounding_matrix = create_bounding_matrix(x.shape[1], y.shape[1], window)
-
-        if not isinstance(g, float):
-            raise ValueError(
-                f"The value of g must be a float. The current value is {g}"
-            )
+    .. math::
+        d_{w}(x_{i}, y_{j}) = ||w_{|i-j|}(x_{i} - y_{j})||
 
-        @njit(cache=True)
-        def numba_wdtw_distance(
-            _x: np.ndarray,
-            _y: np.ndarray,
-        ) -> float:
-            cost_matrix = _weighted_cost_matrix(_x, _y, _bounding_matrix, g)
-            return cost_matrix[-1, -1]
+    Where d_w is the distance with a the weight applied to it for points i, j, where
+    w(|i-j|) is a positive weight between the two points x_i and y_j.
 
-        return numba_wdtw_distance
+    Parameters
+    ----------
+    x: np.ndarray, of shape (n_channels, n_timepoints) or (n_timepoints,) or
+            (n_instances, n_channels, n_timepoints)
+        First time series.
+    y: np.ndarray, of shape (m_channels, m_timepoints) or (m_timepoints,) or
+            (m_instances, m_channels, m_timepoints)
+        Second time series.
+    window: float, defaults=None
+        The window to use for the bounding matrix. If None, no bounding matrix
+        is used.
+    g: float, defaults=0.05
+        Constant that controls the level of penalisation for the points with larger
+        phase difference. Default is 0.05.
 
+    Returns
+    -------
+    float
+        wdtw distance between x and y.
+
+    Raises
+    ------
+    ValueError
+        If x and y are not 1D, 2D, or 3D arrays.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aeon.distances import wdtw_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 = wdtw_distance(x, y)
 
-@njit(cache=True)
-def _weighted_cost_matrix(
-    x: np.ndarray, y: np.ndarray, bounding_matrix: np.ndarray, g: float
-):
+    References
+    ----------
+    .. [1] Young-Seon Jeong, Myong K. Jeong, Olufemi A. Omitaomu, Weighted dynamic time
+    warping for time series classification, Pattern Recognition, Volume 44, Issue 9,
+    2011, Pages 2231-2240, ISSN 0031-3203, https://doi.org/10.1016/j.patcog.2010.09.022.
+    """
+    if x.ndim == 1 and y.ndim == 1:
+        _x = x.reshape((1, x.shape[0]))
+        _y = y.reshape((1, y.shape[0]))
+        bounding_matrix = create_bounding_matrix(_x.shape[1], _y.shape[1], window)
+        return _wdtw_distance(_x, _y, bounding_matrix, g)
+    if x.ndim == 2 and y.ndim == 2:
+        bounding_matrix = create_bounding_matrix(x.shape[1], y.shape[1], window)
+        return _wdtw_distance(x, y, bounding_matrix, g)
+    if x.ndim == 3 and y.ndim == 3:
+        distance = 0
+        bounding_matrix = create_bounding_matrix(x.shape[2], y.shape[2], window)
+        for curr_x, curr_y in zip(x, y):
+            distance += _wdtw_distance(curr_x, curr_y, bounding_matrix, g)
+        return distance
+    raise ValueError("x and y must be 1D, 2D, or 3D arrays")
+
+
+@njit(cache=True, fastmath=True)
+def wdtw_cost_matrix(
+    x: np.ndarray, y: np.ndarray, window: float = None, g: float = 0.05
+) -> np.ndarray:
     """Compute the wdtw cost matrix between two time series.
 
     Parameters
     ----------
-    x: np.ndarray (2d array)
+    x: np.ndarray, of shape (n_channels, n_timepoints) or (n_timepoints,) or
+            (n_instances, n_channels, n_timepoints)
         First time series.
-    y: np.ndarray (2d array)
+    y: np.ndarray, of shape (m_channels, m_timepoints) or (m_timepoints,) or
+            (m_instances, m_channels, m_timepoints)
         Second time series.
-    bounding_matrix: np.ndarray (2d of size mxn where m is len(x) and n is len(y))
-        Bounding matrix where the values in bound are marked by finite values and
-        outside bound points are infinite values.
-    g: float
-        Constant that controls the curvature (slope) of the function; that is, g
-        controls the level of penalisation for the points with larger phase difference.
+    window: float, defaults=None
+        The window to use for the bounding matrix. If None, no bounding matrix
+        is used.
+    g: float, defaults=0.05
+        Constant that controls the level of penalisation for the points with larger
+        phase difference. Default is 0.05.
 
     Returns
     -------
-    np.ndarray
-        Weighted cost matrix between x and y time series.
+    np.ndarray (n_timepoints_x, n_timepoints_y)
+        wdtw cost matrix between x and y.
+
+    Raises
+    ------
+    ValueError
+        If x and y are not 1D, 2D, or 3D arrays.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aeon.distances import wdtw_cost_matrix
+    >>> x = np.array([[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]])
+    >>> y = np.array([[1, 2, 3, 4, 5, 6, 7, 8, 9, 10]])
+    >>> wdtw_cost_matrix(x, y)
+    array([[  0.        ,   0.450166  ,   2.30044662,   6.57563393,
+             14.37567559,  26.87567559,  45.32558186,  71.04956205,
+            105.44507215, 149.98162593],
+           [  0.450166  ,   0.        ,   0.450166  ,   2.30044662,
+              6.57563393,  14.37567559,  26.87567559,  45.32558186,
+             71.04956205, 105.44507215],
+           [  2.30044662,   0.450166  ,   0.        ,   0.450166  ,
+              2.30044662,   6.57563393,  14.37567559,  26.87567559,
+             45.32558186,  71.04956205],
+           [  6.57563393,   2.30044662,   0.450166  ,   0.        ,
+              0.450166  ,   2.30044662,   6.57563393,  14.37567559,
+             26.87567559,  45.32558186],
+           [ 14.37567559,   6.57563393,   2.30044662,   0.450166  ,
+              0.        ,   0.450166  ,   2.30044662,   6.57563393,
+             14.37567559,  26.87567559],
+           [ 26.87567559,  14.37567559,   6.57563393,   2.30044662,
+              0.450166  ,   0.        ,   0.450166  ,   2.30044662,
+              6.57563393,  14.37567559],
+           [ 45.32558186,  26.87567559,  14.37567559,   6.57563393,
+              2.30044662,   0.450166  ,   0.        ,   0.450166  ,
+              2.30044662,   6.57563393],
+           [ 71.04956205,  45.32558186,  26.87567559,  14.37567559,
+              6.57563393,   2.30044662,   0.450166  ,   0.        ,
+              0.450166  ,   2.30044662],
+           [105.44507215,  71.04956205,  45.32558186,  26.87567559,
+             14.37567559,   6.57563393,   2.30044662,   0.450166  ,
+              0.        ,   0.450166  ],
+           [149.98162593, 105.44507215,  71.04956205,  45.32558186,
+             26.87567559,  14.37567559,   6.57563393,   2.30044662,
+              0.450166  ,   0.        ]])
     """
-    dimensions = x.shape[0]
+    if x.ndim == 1 and y.ndim == 1:
+        _x = x.reshape((1, x.shape[0]))
+        _y = y.reshape((1, y.shape[0]))
+        bounding_matrix = create_bounding_matrix(_x.shape[1], _y.shape[1], window)
+        return _wdtw_cost_matrix(_x, _y, bounding_matrix, g)
+    if x.ndim == 2 and y.ndim == 2:
+        bounding_matrix = create_bounding_matrix(x.shape[1], y.shape[1], window)
+        return _wdtw_cost_matrix(x, y, bounding_matrix, g)
+    if x.ndim == 3 and y.ndim == 3:
+        bounding_matrix = create_bounding_matrix(x.shape[2], y.shape[2], window)
+        cost_matrix = np.zeros((x.shape[2], y.shape[2]))
+        for curr_x, curr_y in zip(x, y):
+            cost_matrix = np.add(
+                cost_matrix, _wdtw_cost_matrix(curr_x, curr_y, bounding_matrix, g)
+            )
+        return cost_matrix
+    raise ValueError("x and y must be 1D, 2D, or 3D arrays")
+
+
+@njit(cache=True, fastmath=True)
+def _wdtw_distance(
+    x: np.ndarray, y: np.ndarray, bounding_matrix: np.ndarray, g: float
+) -> float:
+    return _wdtw_cost_matrix(x, y, bounding_matrix, g)[x.shape[1] - 1, y.shape[1] - 1]
+
+
+@njit(cache=True, fastmath=True)
+def _wdtw_cost_matrix(
+    x: np.ndarray, y: np.ndarray, bounding_matrix: np.ndarray, g: float
+) -> np.ndarray:
     x_size = x.shape[1]
     y_size = y.shape[1]
     cost_matrix = np.full((x_size + 1, y_size + 1), np.inf)
     cost_matrix[0, 0] = 0.0
 
+    max_size = max(x_size, y_size)
     weight_vector = np.array(
-        [1 / (1 + np.exp(-g * (i - x_size / 2))) for i in range(0, x_size)]
+        [1 / (1 + np.exp(-g * (i - max_size / 2))) for i in range(0, max_size)]
     )
 
     for i in range(x_size):
         for j in range(y_size):
             if bounding_matrix[i, j]:
-                sum = 0
-                for k in range(dimensions):
-                    sum += (x[k][i] - y[k][j]) * (x[k][i] - y[k][j])
-                cost_matrix[i + 1, j + 1] = (
-                    min(cost_matrix[i, j + 1], cost_matrix[i + 1, j], cost_matrix[i, j])
-                    + weight_vector[np.abs(i - j)] * sum
+                cost_matrix[i + 1, j + 1] = _univariate_squared_distance(
+                    x[:, i], y[:, j]
+                ) * weight_vector[abs(i - j)] + min(
+                    cost_matrix[i, j + 1],
+                    cost_matrix[i + 1, j],
+                    cost_matrix[i, j],
                 )
 
     return cost_matrix[1:, 1:]
+
+
+@njit(cache=True, fastmath=True)
+def wdtw_pairwise_distance(
+    X: np.ndarray, y: np.ndarray = None, window: float = None, g: float = 0.05
+) -> np.ndarray:
+    """Compute the wdtw pairwise distance between a set of time series.
+
+    Parameters
+    ----------
+    X: np.ndarray, of shape (n_instances, n_channels, n_timepoints) or
+            (n_instances, n_timepoints)
+        A collection of time series instances.
+    y: np.ndarray, of shape (m_instances, m_channels, m_timepoints) or
+            (m_instances, m_timepoints) or (m_timepoints,), default=None
+        A collection of time series instances.
+    window: float, default=None
+        The window to use for the bounding matrix. If None, no bounding matrix
+        is used.
+    g: float, defaults=0.05
+        Constant that controls the level of penalisation for the points with larger
+        phase difference. Default is 0.05.
+
+    Returns
+    -------
+    np.ndarray (n_instances, n_instances)
+        wdtw pairwise matrix between the instances of X.
+
+    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.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aeon.distances import wdtw_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]]])
+    >>> wdtw_pairwise_distance(X)
+    array([[ 0.        , 12.61266072, 51.97594869],
+           [12.61266072,  0.        , 12.61266072],
+           [51.97594869, 12.61266072,  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]]])
+    >>> wdtw_pairwise_distance(X, y)
+    array([[144.37763524, 243.99820355, 369.60674621],
+           [ 70.74504127, 144.37763524, 243.99820355],
+           [ 23.10042164,  70.74504127, 144.37763524]])
+
+    >>> 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]])
+    >>> wdtw_pairwise_distance(X, y_univariate)
+    array([[144.37763524],
+           [ 70.74504127],
+           [ 23.10042164]])
+    """
+    if y is None:
+        # To self
+        if X.ndim == 3:
+            return _wdtw_pairwise_distance(X, window, g)
+        if X.ndim == 2:
+            _X = X.reshape((X.shape[0], 1, X.shape[1]))
+            return _wdtw_pairwise_distance(_X, window, g)
+        raise ValueError("x and y must be 2D or 3D arrays")
+    elif y.ndim == X.ndim:
+        # Multiple to multiple
+        if y.ndim == 3 and X.ndim == 3:
+            return _wdtw_from_multiple_to_multiple_distance(X, y, window, g)
+        if y.ndim == 2 and X.ndim == 2:
+            _x = X.reshape((X.shape[0], 1, X.shape[1]))
+            _y = y.reshape((y.shape[0], 1, y.shape[1]))
+            return _wdtw_from_multiple_to_multiple_distance(_x, _y, window, g)
+        if y.ndim == 1 and X.ndim == 1:
+            _x = X.reshape((1, 1, X.shape[0]))
+            _y = y.reshape((1, 1, y.shape[0]))
+            return _wdtw_from_multiple_to_multiple_distance(_x, _y, window, g)
+        raise ValueError("x and y must be 1D, 2D, or 3D arrays")
+    else:
+        # Single to multiple
+        if X.ndim == 3 and y.ndim == 2:
+            _y = y.reshape((1, y.shape[0], y.shape[1]))
+            return _wdtw_from_multiple_to_multiple_distance(X, _y, window, g)
+        if y.ndim == 3 and X.ndim == 2:
+            _x = X.reshape((1, X.shape[0], X.shape[1]))
+            return _wdtw_from_multiple_to_multiple_distance(_x, y, window, g)
+        if X.ndim == 2 and y.ndim == 1:
+            _x = X.reshape((X.shape[0], 1, X.shape[1]))
+            _y = y.reshape((1, 1, y.shape[0]))
+            return _wdtw_from_multiple_to_multiple_distance(_x, _y, window, g)
+        if y.ndim == 2 and X.ndim == 1:
+            _x = X.reshape((1, 1, X.shape[0]))
+            _y = y.reshape((y.shape[0], 1, y.shape[1]))
+            return _wdtw_from_multiple_to_multiple_distance(_x, _y, window, g)
+        else:
+            raise ValueError("x and y must be 2D or 3D arrays")
+
+
+@njit(cache=True, fastmath=True)
+def _wdtw_pairwise_distance(X: np.ndarray, window: float, g: float) -> np.ndarray:
+    n_instances = X.shape[0]
+    distances = np.zeros((n_instances, n_instances))
+    bounding_matrix = create_bounding_matrix(X.shape[2], X.shape[2], window)
+
+    for i in range(n_instances):
+        for j in range(i + 1, n_instances):
+            distances[i, j] = _wdtw_distance(X[i], X[j], bounding_matrix, g)
+            distances[j, i] = distances[i, j]
+
+    return distances
+
+
+@njit(cache=True, fastmath=True)
+def _wdtw_from_multiple_to_multiple_distance(
+    x: np.ndarray, y: np.ndarray, window: float, g: float
+) -> np.ndarray:
+    n_instances = x.shape[0]
+    m_instances = y.shape[0]
+    distances = np.zeros((n_instances, m_instances))
+    bounding_matrix = create_bounding_matrix(x.shape[2], y.shape[2], window)
+
+    for i in range(n_instances):
+        for j in range(m_instances):
+            distances[i, j] = _wdtw_distance(x[i], y[j], bounding_matrix, g)
+    return distances
+
+
+@njit(cache=True, fastmath=True)
+def wdtw_alignment_path(
+    x: np.ndarray, y: np.ndarray, window: float = None, g: float = 0.05
+) -> Tuple[List[Tuple[int, int]], float]:
+    """Compute the wdtw alignment path between two time series.
+
+    Parameters
+    ----------
+    x: np.ndarray, of shape (n_channels, n_timepoints) or (n_timepoints,) or
+            (n_instances, n_channels, n_timepoints)
+        First time series.
+    y: np.ndarray, of shape (m_channels, m_timepoints) or (m_timepoints,) or
+            (m_instances, m_channels, m_timepoints)
+        Second time series.
+    window: float, default=None
+        The window to use for the bounding matrix. If None, no bounding matrix
+        is used.
+    g: float, defaults=0.05
+        Constant that controls the level of penalisation for the points with larger
+        phase difference. Default is 0.05.
+
+    Returns
+    -------
+    List[Tuple[int, int]]
+        The alignment path between the two time series where each element is a tuple
+        of the index in x and the index in y that have the best alignment according
+        to the cost matrix.
+    float
+        The wdtw distance betweeen the two time series.
+
+    Raises
+    ------
+    ValueError
+        If x and y are not 1D, 2D, or 3D arrays.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aeon.distances import wdtw_alignment_path
+    >>> x = np.array([[1, 2, 3, 6]])
+    >>> y = np.array([[1, 2, 3, 4]])
+    >>> path, dist = wdtw_alignment_path(x, y)
+    >>> path
+    [(0, 0), (1, 1), (2, 2), (3, 3)]
+    """
+    cost_matrix = wdtw_cost_matrix(x, y, window, g)
+    return (
+        compute_min_return_path(cost_matrix),
+        cost_matrix[x.shape[-1] - 1, y.shape[-1] - 1],
+    )


Test patch:
diff --git a/aeon/distances/tests/_expected_results.py b/aeon/distances/tests/_expected_results.py
index cd886914..ef0a3b5c 100644
--- a/aeon/distances/tests/_expected_results.py
+++ b/aeon/distances/tests/_expected_results.py
@@ -44,7 +44,15 @@ _expected_distance_results = {
         35.916981095128804,
         281.92058999075414,
     ],
-    "wdtw": [12.343758137512241, 0.985380547171357, 21.265839226825413],
+    "wdtw": [
+        12.343758137512241,
+        0.985380547171357,
+        21.265839226825413,
+        205.30721939046768,
+        2.0040890166976926,
+        20.795690703034445,
+        189.3408257880001,
+    ],
     "wddtw": [0.0, 1.0442409418611003, 17.418900020282003],
     "msm": [5.0, 6.828557434224288, None],
     "twe": [5.0, 11.548698748091073, 39.87793560457224],
@@ -73,7 +81,7 @@ _expected_distance_results_params = {
     ],
     "wdtw": [
         [1.364177898415523, 21.265839226825413],
-        [0.2643712293510365, 3.8171163209131693],
+        [0.02752598656586074, 0.33677832093219406],
     ],
     "wddtw": [
         [1.3419794992173555, 17.418900020282003],
diff --git a/aeon/distances/tests/test_new_distances.py b/aeon/distances/tests/test_new_distances.py
index fe7e87a0..b15a2b92 100644
--- a/aeon/distances/tests/test_new_distances.py
+++ b/aeon/distances/tests/test_new_distances.py
@@ -19,6 +19,10 @@ from aeon.distances import (
     euclidean_pairwise_distance,
     squared_distance,
     squared_pairwise_distance,
+    wdtw_alignment_path,
+    wdtw_cost_matrix,
+    wdtw_distance,
+    wdtw_pairwise_distance,
 )
 from aeon.distances.tests._expected_results import _expected_distance_results
 from aeon.distances.tests._utils import create_test_distance_numpy
@@ -48,6 +52,13 @@ DISTANCES = [
         "cost_matrix": ddtw_cost_matrix,
         "alignment_path": ddtw_alignment_path,
     },
+    {
+        "name": "wdtw",
+        "distance": wdtw_distance,
+        "pairwise_distance": wdtw_pairwise_distance,
+        "cost_matrix": wdtw_cost_matrix,
+        "alignment_path": wdtw_alignment_path,
+    },
 ]
 
 
diff --git a/aeon/distances/tests/test_numba_distance_parameters.py b/aeon/distances/tests/test_numba_distance_parameters.py
index d3f04220..d495fca0 100644
--- a/aeon/distances/tests/test_numba_distance_parameters.py
+++ b/aeon/distances/tests/test_numba_distance_parameters.py
@@ -68,7 +68,7 @@ DIST_PARAMS = {
     "edr": BASIC_BOUNDING_PARAMS + [{"epsilon": 0.5}],
     "lcss": BASIC_BOUNDING_PARAMS + [{"epsilon": 0.5}],
     "ddtw": BASIC_BOUNDING_PARAMS,
-    "wdtw": BASIC_BOUNDING_PARAMS + [{"g": 0.5}],
+    "wdtw": BASIC_BOUNDING_PARAMS + [{"g": 1.0}],
     "wddtw": BASIC_BOUNDING_PARAMS
     + [{"compute_derivative": _test_derivative}]
     + [{"g": 0.5}],