Repository: aeon-toolkit/aeon
Pull request #333 (2023-04-25 01:23:27)

[ENH] Update erp 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 #332

#### What does this implement/fix? Explain your changes.
This pr updates erp to the new format for distances. In addition in doing so it makes optimisations for computing the distances between multiple time series. 

In addition, this pr updates the 'g' parameter allowing a array to be passed that allows for a different value of g per channel (the suggested way in the original paper to set g is the standard deviation of each channel for each channel).

<!--
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 753170a1..9ae07ae4 100644
--- a/aeon/distances/__init__.py
+++ b/aeon/distances/__init__.py
@@ -27,6 +27,10 @@ __all__ = [
     "lcss_pairwise_distance",
     "lcss_alignment_path",
     "lcss_cost_matrix",
+    "erp_distance",
+    "erp_pairwise_distance",
+    "erp_alignment_path",
+    "erp_cost_matrix",
     "distance",
     "distance_factory",
     "pairwise_distance",
@@ -58,8 +62,6 @@ from aeon.distances._distance import (
     distance_factory,
     edr_alignment_path,
     edr_distance,
-    erp_alignment_path,
-    erp_distance,
     msm_alignment_path,
     msm_distance,
     pairwise_distance,
@@ -72,6 +74,12 @@ from aeon.distances._dtw import (
     dtw_distance,
     dtw_pairwise_distance,
 )
+from aeon.distances._erp import (
+    erp_alignment_path,
+    erp_cost_matrix,
+    erp_distance,
+    erp_pairwise_distance,
+)
 from aeon.distances._euclidean import euclidean_distance, euclidean_pairwise_distance
 from aeon.distances._lcss import (
     lcss_alignment_path,
diff --git a/aeon/distances/_distance.py b/aeon/distances/_distance.py
index 9591ba3e..5d3f7d6a 100644
--- a/aeon/distances/_distance.py
+++ b/aeon/distances/_distance.py
@@ -13,7 +13,7 @@ from aeon.distances._ddtw import (
 )
 from aeon.distances._dtw import dtw_alignment_path, dtw_distance, dtw_pairwise_distance
 from aeon.distances._edr import _EdrDistance
-from aeon.distances._erp import _ErpDistance
+from aeon.distances._erp import erp_alignment_path, erp_distance, erp_pairwise_distance
 from aeon.distances._euclidean import euclidean_distance, euclidean_pairwise_distance
 from aeon.distances._lcss import (
     lcss_alignment_path,
@@ -52,79 +52,6 @@ from aeon.distances.base import (
 )
 
 
-def erp_distance(
-    x: np.ndarray,
-    y: np.ndarray,
-    window: Union[float, None] = None,
-    g: float = 0.0,
-    **kwargs: Any,
-) -> float:
-    """Compute the Edit distance for real penalty (ERP) distance between two series.
-
-    ERP, first proposed in [1]_, attempts align time series
-    by better considering how indexes are carried forward through the cost matrix.
-    Usually in the dtw cost matrix, if an alignment can't be found the previous value
-    is carried forward. Erp instead proposes the idea of gaps or sequences of points
-    that have no matches. These gaps are then punished based on their distance from 'g'.
-
-    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.
-        The reference value to penalise gaps.
-    kwargs: Any
-        Extra kwargs.
-
-    Returns
-    -------
-    float
-        ERP distance between 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 3 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
-        If g is not a float.
-
-    Examples
-    --------
-    >>> x_1d = np.array([1, 2, 3, 4])  # 1d array
-    >>> y_1d = np.array([5, 6, 7, 8])  # 1d array
-    >>> erp_distance(x_1d, y_1d)
-    16.0
-
-    >>> 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
-    >>> erp_distance(x_2d, y_2d)
-    45.254833995939045
-
-    References
-    ----------
-    .. [1] Lei Chen and Raymond Ng. 2004. On the marriage of Lp-norms and edit distance.
-    In Proceedings of the Thirtieth international conference on Very large data bases
-     - Volume 30 (VLDB '04). VLDB Endowment, 792–803.
-    """
-    format_kwargs = {
-        "window": window,
-        "g": g,
-    }
-    format_kwargs = {**format_kwargs, **kwargs}
-
-    return distance(x, y, metric="erp", **format_kwargs)
-
-
 def edr_distance(
     x: np.ndarray,
     y: np.ndarray,
@@ -427,77 +354,6 @@ def edr_alignment_path(
     )
 
 
-def erp_alignment_path(
-    x: np.ndarray,
-    y: np.ndarray,
-    return_cost_matrix: bool = False,
-    window: Union[float, None] = None,
-    g: float = 0.0,
-    **kwargs: Any,
-) -> AlignmentPathReturn:
-    """Compute the Edit distance for real penalty (ERP) alignment path.
-
-    ERP, first proposed in [1]_, attempts align time series
-    by better considering how indexes are carried forward through the cost matrix.
-    Usually in the dtw cost matrix, if an alignment can't be found the previous value
-    is carried forward. Erp instead proposes the idea of gaps or sequences of points
-    that have no matches. These gaps are then punished based on their distance from 'g'.
-
-    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.
-        The reference value to penalise gaps.
-    kwargs: Any
-        Extra kwargs.
-
-    Returns
-    -------
-    list[tuple]
-        List of tuples containing the erp alignment path.
-    float
-        Erp 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 3 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
-        If g is not a float.
-
-    References
-    ----------
-    .. [1] Lei Chen and Raymond Ng. 2004. On the marriage of Lp-norms and edit distance.
-    In Proceedings of the Thirtieth international conference on Very large data bases
-     - Volume 30 (VLDB '04). VLDB Endowment, 792–803.
-    """
-    format_kwargs = {
-        "window": window,
-        "g": g,
-    }
-    format_kwargs = {**format_kwargs, **kwargs}
-
-    return distance_alignment_path(
-        x, y, metric="erp", return_cost_matrix=return_cost_matrix, **format_kwargs
-    )
-
-
 def msm_alignment_path(
     x: np.ndarray,
     y: np.ndarray,
@@ -650,7 +506,7 @@ def twe_alignment_path(
     )
 
 
-NEW_DISTANCES = ["squared", "euclidean", "dtw", "ddtw", "wdtw", "wddtw", "lcss"]
+NEW_DISTANCES = ["squared", "euclidean", "dtw", "ddtw", "wdtw", "wddtw", "lcss", "erp"]
 
 
 def distance(
@@ -747,6 +603,8 @@ def distance(
             return wddtw_distance(x, y, **kwargs)
         elif metric == "lcss":
             return lcss_distance(x, y, **kwargs)
+        elif metric == "erp":
+            return erp_distance(x, y, **kwargs)
     _x = to_numba_timeseries(x)
     _y = to_numba_timeseries(y)
 
@@ -828,6 +686,8 @@ def distance_factory(
             return wddtw_distance
         elif metric == "lcss":
             return lcss_distance
+        elif metric == "erp":
+            return erp_distance
     global dist_callable
 
     if x is None:
@@ -958,6 +818,8 @@ def pairwise_distance(
             return wddtw_pairwise_distance(_x, _y, **kwargs)
         elif metric == "lcss":
             return lcss_pairwise_distance(_x, _y, **kwargs)
+        elif metric == "erp":
+            return erp_pairwise_distance(_x, _y, **kwargs)
 
     symmetric = np.array_equal(_x, _y)
     _metric_callable = _resolve_metric_to_factory(
@@ -1046,6 +908,8 @@ def distance_alignment_path(
             return wddtw_alignment_path(x, y, **kwargs)
         elif metric == "lcss":
             return lcss_alignment_path(x, y, **kwargs)
+        elif metric == "erp":
+            return erp_alignment_path(x, y, **kwargs)
     _x = to_numba_timeseries(x)
     _y = to_numba_timeseries(y)
 
@@ -1123,7 +987,7 @@ def distance_alignment_path_factory(
     if metric in NEW_DISTANCES:
         if metric == "dtw":
             return dtw_alignment_path
-        if metric == "ddtw":
+        elif metric == "ddtw":
             return ddtw_alignment_path
         elif metric == "wdtw":
             return wdtw_alignment_path
@@ -1131,6 +995,8 @@ def distance_alignment_path_factory(
             return wddtw_alignment_path
         elif metric == "lcss":
             return lcss_alignment_path
+        elif metric == "erp":
+            return erp_alignment_path
     if x is None:
         x = np.zeros((1, 10))
     if y is None:
@@ -1153,13 +1019,6 @@ def distance_alignment_path_factory(
 
 
 _METRIC_INFOS = [
-    MetricInfo(
-        canonical_name="erp",
-        aka={"erp", "edit distance with real penalty"},
-        dist_func=erp_distance,
-        dist_instance=_ErpDistance(),
-        dist_alignment_path_func=erp_alignment_path,
-    ),
     MetricInfo(
         canonical_name="edr",
         aka={"edr", "edit distance for real sequences"},
@@ -1192,7 +1051,6 @@ _METRICS_NAMES = list(_METRICS.keys())
 
 ALL_DISTANCES = (
     edr_distance,
-    erp_distance,
     msm_distance,
     twe_distance,
 )
diff --git a/aeon/distances/_erp.py b/aeon/distances/_erp.py
index 967cc3e4..c5a24bef 100644
--- a/aeon/distances/_erp.py
+++ b/aeon/distances/_erp.py
@@ -1,213 +1,446 @@
 # -*- coding: utf-8 -*-
+r"""Edit real penalty (erp) distance between two time series.
+
+ERP, first proposed in [1]_, attempts to align time series
+by better considering how indexes are carried forward through the cost matrix.
+Usually in the dtw cost matrix, if an alignment can't be found the previous value
+is carried forward. Erp instead proposes the idea of gaps or sequences of points
+that have no matches. These gaps are then punished based on their distance from 'g'.
+
+References
+----------
+.. [1] Lei Chen and Raymond Ng. 2004. On the marriage of Lp-norms and edit distance.
+In Proceedings of the Thirtieth international conference on Very large data bases
+ - Volume 30 (VLDB '04). VLDB Endowment, 792–803.
+"""
 __author__ = ["chrisholder", "TonyBagnall"]
 
-import warnings
-from typing import Any, List, Tuple
+from typing import List, Tuple, Union
 
 import numpy as np
 from numba import njit
-from numba.core.errors import NumbaWarning
 
 from aeon.distances._alignment_paths import (
     _add_inf_to_out_of_bounds_cost_matrix,
     compute_min_return_path,
 )
 from aeon.distances._bounding_matrix import create_bounding_matrix
-from aeon.distances.base import (
-    DistanceAlignmentPathCallable,
-    DistanceCallable,
-    NumbaDistance,
-)
+from aeon.distances._euclidean import _univariate_euclidean_distance
 
-# Warning occurs when using large time series (i.e. 1000x1000)
-warnings.simplefilter("ignore", category=NumbaWarning)
-
-
-class _ErpDistance(NumbaDistance):
-    """Edit distance with real penalty (erp) between two time series."""
-
-    def _distance_alignment_path_factory(
-        self,
-        x: np.ndarray,
-        y: np.ndarray,
-        return_cost_matrix: bool = False,
-        window: float = None,
-        g: float = 0.0,
-        **kwargs: Any,
-    ) -> DistanceAlignmentPathCallable:
-        """Create a no_python compiled erp distance alignment path callable.
-
-        Similar to LCSS with a different penalty.
-        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.
-            The reference value to penalise gaps.
-        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 times eries is not a numpy array.
-            If the input time series doesn't have exactly 2 dimensions.
-            If g is not a float.
-        """
-        _bounding_matrix = create_bounding_matrix(x.shape[1], y.shape[1], window)
-        if not isinstance(g, float):
-            raise ValueError("The value of g must be a float.")
-
-        if return_cost_matrix is True:
-
-            @njit(cache=True)
-            def numba_erp_distance_alignment_path(
-                _x: np.ndarray, _y: np.ndarray
-            ) -> Tuple[List, float, np.ndarray]:
-                cost_matrix = _erp_cost_matrix(_x, _y, _bounding_matrix, g)
-                temp_cm = _add_inf_to_out_of_bounds_cost_matrix(
-                    cost_matrix, _bounding_matrix
-                )
-                path = compute_min_return_path(temp_cm)
-                return path, cost_matrix[-1, -1], cost_matrix
 
-        else:
+@njit(cache=True, fastmath=True)
+def erp_distance(
+    x: np.ndarray,
+    y: np.ndarray,
+    window: float = None,
+    g: Union[float, np.ndarray] = 0.0,
+) -> float:
+    """Compute the ERP distance between two time series.
 
-            @njit(cache=True)
-            def numba_erp_distance_alignment_path(
-                _x: np.ndarray, _y: np.ndarray
-            ) -> Tuple[List, float]:
-                cost_matrix = _erp_cost_matrix(_x, _y, _bounding_matrix, g)
-                temp_cm = _add_inf_to_out_of_bounds_cost_matrix(
-                    cost_matrix, _bounding_matrix
-                )
-                path = compute_min_return_path(temp_cm)
-                return path, cost_matrix[-1, -1]
-
-        return numba_erp_distance_alignment_path
-
-    def _distance_factory(
-        self,
-        x: np.ndarray,
-        y: np.ndarray,
-        window: float = None,
-        g: float = 0.0,
-        **kwargs: Any,
-    ) -> DistanceCallable:
-        """Create a no_python compiled erp distance callable.
-
-        Similar to LCSS with a different penalty.
-        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.
-        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.
-            The reference value to penalise gaps.
-        kwargs: Any
-            Extra kwargs.
-
-        Returns
-        -------
-        Callable[[np.ndarray, np.ndarray], float]
-            No_python compiled erp distance callable.
-
-        Raises
-        ------
-        ValueError
-            If the input time series is not a numpy array.
-            If the input time series doesn't have exactly 2 dimensions.
-            If g is not a float.
-        """
-        _bounding_matrix = create_bounding_matrix(x.shape[1], y.shape[1], window)
-
-        if not isinstance(g, float):
-            raise ValueError("The value of g must be a float.")
-
-        @njit(cache=True)
-        def numba_erp_distance(_x: np.ndarray, _y: np.ndarray) -> float:
-            cost_matrix = _erp_cost_matrix(_x, _y, _bounding_matrix, g)
-
-            return cost_matrix[-1, -1]
-
-        return numba_erp_distance
-
-
-@njit(cache=True)
-def _erp_cost_matrix(
-    x: np.ndarray, y: np.ndarray, bounding_matrix: np.ndarray, g: float
-):
-    """Compute the erp cost matrix between two time series.
+    ERP, first proposed in [1]_, attempts to align time series
+    by better considering how indexes are carried forward through the cost matrix.
+    Usually in the dtw cost matrix, if an alignment can't be found the previous value
+    is carried forward. Erp instead proposes the idea of gaps or sequences of points
+    that have no matches. These gaps are then punished based on their distance from 'g'.
+
+    The optimal value of g is selected from the range [σ/5, σ], where σ is the
+    standard deviation of the training data. When there is > 1 channel, g should
+    be a np.ndarray where the nth value is the standard deviation of the nth
+    channel.
+
+    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 or np.ndarray of shape (n_channels), defaults=0.
+        The reference value to penalise gaps. The default is 0. If it is an array
+        then it must be the length of the number of channels in x and y. If a single
+        value is provided then that value is used across each channel
+
+    Returns
+    -------
+    float
+        ERP 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 erp_distance
+    >>> 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]])
+    >>> dist = erp_distance(x, y)
+
+    References
+    ----------
+    .. [1] Lei Chen and Raymond Ng. 2004. On the marriage of Lp-norms and edit distance.
+    In Proceedings of the Thirtieth international conference on Very large data bases
+     - Volume 30 (VLDB '04). VLDB Endowment, 792–803.
+    """
+    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 _erp_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 _erp_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 += _erp_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 erp_cost_matrix(
+    x: np.ndarray,
+    y: np.ndarray,
+    window: float = None,
+    g: Union[float, np.ndarray] = 0.0,
+) -> np.ndarray:
+    """Compute the ERP cost matrix between two time series.
+
+    The optimal value of g is selected from the range [σ/5, σ], where σ is the
+    standard deviation of the training data. When there is > 1 channel, g should
+    be a np.ndarray where the nth value is the standard deviation of the nth
+    channel.
 
     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
-        The reference value to penalise gaps ('gap' defined when an alignment to
-        the next value (in x) in value can't be found).
+    window: float, defaults=None
+        The window to use for the bounding matrix. If None, no bounding matrix
+        is used.
+    g: float or np.ndarray of shape (n_channels), defaults=0.
+        The reference value to penalise gaps. The default is 0. If it is an array
+        then it must be the length of the number of channels in x and y. If a single
+        value is provided then that value is used across each channel.
 
     Returns
     -------
-    np.ndarray (2d of size mxn where m is len(x) and n is len(y))
-        Erp cost matrix between x and y.
+    np.ndarray (n_timepoints_x, n_timepoints_y)
+        ERP 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 erp_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]])
+    >>> erp_cost_matrix(x, y)
+    array([[ 0.,  2.,  5.,  9., 14., 20., 27., 35., 44., 54.],
+           [ 2.,  0.,  3.,  7., 12., 18., 25., 33., 42., 52.],
+           [ 5.,  3.,  0.,  4.,  9., 15., 22., 30., 39., 49.],
+           [ 9.,  7.,  4.,  0.,  5., 11., 18., 26., 35., 45.],
+           [14., 12.,  9.,  5.,  0.,  6., 13., 21., 30., 40.],
+           [20., 18., 15., 11.,  6.,  0.,  7., 15., 24., 34.],
+           [27., 25., 22., 18., 13.,  7.,  0.,  8., 17., 27.],
+           [35., 33., 30., 26., 21., 15.,  8.,  0.,  9., 19.],
+           [44., 42., 39., 35., 30., 24., 17.,  9.,  0., 10.],
+           [54., 52., 49., 45., 40., 34., 27., 19., 10.,  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 _erp_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 _erp_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, _erp_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 _erp_distance(
+    x: np.ndarray,
+    y: np.ndarray,
+    bounding_matrix: np.ndarray,
+    g: Union[float, np.ndarray],
+) -> float:
+    return _erp_cost_matrix(x, y, bounding_matrix, g)[x.shape[1] - 1, y.shape[1] - 1]
+
+
+@njit(cache=True, fastmath=True)
+def _erp_cost_matrix(
+    x: np.ndarray,
+    y: np.ndarray,
+    bounding_matrix: np.ndarray,
+    g: Union[float, np.ndarray],
+) -> np.ndarray:
     x_size = x.shape[1]
     y_size = y.shape[1]
+
     cost_matrix = np.zeros((x_size + 1, y_size + 1))
-    gx_distance = np.zeros(x_size)
-    gy_distance = np.zeros(y_size)
-    for j in range(x_size):
-        for i in range(dimensions):
-            gx_distance[j] += (x[i][j] - g) * (x[i][j] - g)
-        gx_distance[j] = np.sqrt(gx_distance[j])
-    for j in range(y_size):
-        for i in range(dimensions):
-            gy_distance[j] += (y[i][j] - g) * (y[i][j] - g)
-        gy_distance[j] = np.sqrt(gy_distance[j])
-    cost_matrix[1:, 0] = np.sum(gx_distance)
-    cost_matrix[0, 1:] = np.sum(gy_distance)
+
+    gx_distance, x_sum = _precompute_g(x, g)
+    gy_distance, y_sum = _precompute_g(y, g)
+
+    cost_matrix[1:, 0] = x_sum
+    cost_matrix[0, 1:] = y_sum
 
     for i in range(1, x_size + 1):
         for j in range(1, y_size + 1):
             if bounding_matrix[i - 1, j - 1]:
-                curr_dist = 0
-                for k in range(dimensions):
-                    curr_dist += (x[k][i - 1] - y[k][j - 1]) * (
-                        x[k][i - 1] - y[k][j - 1]
-                    )
-                curr_dist = np.sqrt(curr_dist)
                 cost_matrix[i, j] = min(
-                    cost_matrix[i - 1, j - 1] + curr_dist,
+                    cost_matrix[i - 1, j - 1]
+                    + _univariate_euclidean_distance(x[:, i - 1], y[:, j - 1]),
                     cost_matrix[i - 1, j] + gx_distance[i - 1],
                     cost_matrix[i, j - 1] + gy_distance[j - 1],
                 )
+
     return cost_matrix[1:, 1:]
+
+
+@njit(cache=True, fastmath=True)
+def _precompute_g(
+    x: np.ndarray, g: Union[float, np.ndarray]
+) -> Tuple[np.ndarray, float]:
+    gx_distance = np.zeros(x.shape[1])
+    if isinstance(g, float):
+        g_arr = np.full(x.shape[0], g)
+    else:
+        if g.shape[0] != x.shape[0]:
+            raise ValueError("g must be a float or an array with shape (x.shape[0],)")
+        g_arr = g
+    x_sum = 0
+
+    for i in range(x.shape[1]):
+        temp = _univariate_euclidean_distance(x[:, i], g_arr)
+        gx_distance[i] = temp
+        x_sum += temp
+    return gx_distance, x_sum
+
+
+@njit(cache=True, fastmath=True)
+def erp_pairwise_distance(
+    X: np.ndarray,
+    y: np.ndarray = None,
+    window: float = None,
+    g: Union[float, np.ndarray] = 0.0,
+) -> np.ndarray:
+    """Compute the erp pairwise distance between a set of time series.
+
+    The optimal value of g is selected from the range [σ/5, σ], where σ is the
+    standard deviation of the training data. When there is > 1 channel, g should
+    be a np.ndarray where the nth value is the standard deviation of the nth
+    channel.
+
+    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 or np.ndarray of shape (n_channels), defaults=0
+        The reference value to penalise gaps. The default is 0. If it is an array
+        then it must be the length of the number of channels in x and y. If a single
+        value is provided then that value is used across each channel.
+
+    Returns
+    -------
+    np.ndarray (n_instances, n_instances)
+        erp 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 erp_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]]])
+    >>> erp_pairwise_distance(X)
+    array([[ 0.,  9., 18.],
+           [ 9.,  0.,  9.],
+           [18.,  9.,  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]]])
+    >>> erp_pairwise_distance(X, y)
+    array([[30., 39., 48.],
+           [21., 30., 39.],
+           [12., 21., 30.]])
+
+    >>> 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]])
+    >>> erp_pairwise_distance(X, y_univariate)
+    array([[30.],
+           [21.],
+           [12.]])
+    """
+    if y is None:
+        # To self
+        if X.ndim == 3:
+            return _erp_pairwise_distance(X, window, g)
+        if X.ndim == 2:
+            _X = X.reshape((X.shape[0], 1, X.shape[1]))
+            return _erp_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 _erp_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 _erp_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 _erp_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 _erp_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 _erp_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 _erp_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 _erp_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 _erp_pairwise_distance(
+    X: np.ndarray, window: float, g: Union[float, np.ndarray]
+) -> 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] = _erp_distance(X[i], X[j], bounding_matrix, g)
+            distances[j, i] = distances[i, j]
+
+    return distances
+
+
+@njit(cache=True, fastmath=True)
+def _erp_from_multiple_to_multiple_distance(
+    x: np.ndarray, y: np.ndarray, window: float, g: Union[float, np.ndarray]
+) -> 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] = _erp_distance(x[i], y[j], bounding_matrix, g)
+    return distances
+
+
+@njit(cache=True, fastmath=True)
+def erp_alignment_path(
+    x: np.ndarray,
+    y: np.ndarray,
+    window: float = None,
+    g: Union[float, np.ndarray] = 0.0,
+) -> Tuple[List[Tuple[int, int]], float]:
+    """Compute the erp alignment path between two time series.
+
+    The optimal value of g is selected from the range [σ/5, σ], where σ is the
+    standard deviation of the training data. When there is > 1 channel, g should
+    be a np.ndarray where the nth value is the standard deviation of the nth
+    channel.
+
+    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 or np.ndarray of shape (n_channels), defaults=0.
+        The reference value to penalise gaps. The default is 0. If it is an array
+        then it must be the length of the number of channels in x and y. If a single
+        value is provided then that value is used across each channel.
+
+    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 erp 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 erp_alignment_path
+    >>> x = np.array([[1, 2, 3, 6]])
+    >>> y = np.array([[1, 2, 3, 4]])
+    >>> erp_alignment_path(x, y)
+    ([(0, 0), (1, 1), (2, 2), (3, 3)], 2.0)
+    """
+    bounding_matrix = create_bounding_matrix(x.shape[-1], y.shape[-1], window)
+    cost_matrix = _add_inf_to_out_of_bounds_cost_matrix(
+        erp_cost_matrix(x, y, window, g), bounding_matrix
+    )
+    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 65122a90..2f40d3b2 100644
--- a/aeon/distances/tests/_expected_results.py
+++ b/aeon/distances/tests/_expected_results.py
@@ -14,7 +14,15 @@ _expected_distance_results = {
         4.938969178714248,
         49.694052115721675,
     ],
-    "erp": [5.0, 5.037672883414786, 20.724482073800456],
+    "erp": [
+        5.0,
+        5.037672883414786,
+        20.724482073800456,
+        208.797697936059,
+        5.16666987535642,
+        18.353616712062177,
+        171.1112853298397,
+    ],
     "edr": [1.0, 0.6, 1.0],
     "lcss": [1.0, 0.1, 1.0, 9.8, 0.0, 1.0, 9.4],
     "squared": [
@@ -75,6 +83,7 @@ _expected_distance_results_params = {
     "erp": [
         [0.6648081862148058, 4.365472428062562],
         [5.279748833082764, 20.71830043391765],
+        [7.399565602170839, 21.90092421200079],
     ],
     "edr": [
         [0.3, 0.3],
diff --git a/aeon/distances/tests/test_new_distances.py b/aeon/distances/tests/test_new_distances.py
index 83a9f789..3d17e838 100644
--- a/aeon/distances/tests/test_new_distances.py
+++ b/aeon/distances/tests/test_new_distances.py
@@ -15,6 +15,10 @@ from aeon.distances import (
     dtw_cost_matrix,
     dtw_distance,
     dtw_pairwise_distance,
+    erp_alignment_path,
+    erp_cost_matrix,
+    erp_distance,
+    erp_pairwise_distance,
     euclidean_distance,
     euclidean_pairwise_distance,
     lcss_alignment_path,
@@ -81,6 +85,13 @@ DISTANCES = [
         "cost_matrix": lcss_cost_matrix,
         "alignment_path": lcss_alignment_path,
     },
+    {
+        "name": "erp",
+        "distance": erp_distance,
+        "pairwise_distance": erp_pairwise_distance,
+        "cost_matrix": erp_cost_matrix,
+        "alignment_path": erp_alignment_path,
+    },
 ]
 
 
@@ -102,7 +113,6 @@ def _validate_distance_result(
 @pytest.mark.parametrize("dist", DISTANCES)
 def test_new_distances(dist):
     # Test univariate
-
     if dist["name"] != "ddtw" and dist["name"] != "wddtw":
         _validate_distance_result(
             np.array([10.0]),
diff --git a/aeon/distances/tests/test_numba_distance_parameters.py b/aeon/distances/tests/test_numba_distance_parameters.py
index 29d4a7f3..466df36b 100644
--- a/aeon/distances/tests/test_numba_distance_parameters.py
+++ b/aeon/distances/tests/test_numba_distance_parameters.py
@@ -29,9 +29,15 @@ def _test_distance_params(
 
     i = 0
     for param_dict in param_list:
+        g_none = False
+        if distance_str == "erp" and "g" in param_dict and param_dict["g"] is None:
+            g_none = True
+
         j = 0
         curr_results = []
         for x, y in test_ts:
+            if g_none:
+                param_dict["g"] = np.std([x, y], axis=0).sum(axis=1)
             results = []
             curr_dist_fact = distance_factory(x, y, metric=distance_str, **param_dict)
             results.append(distance_func(x, y, **param_dict))
@@ -51,6 +57,16 @@ def _test_distance_params(
         i += 1
         results_to_fill.append(curr_results)
 
+    # if distance_str == "erp":
+    #     param_dict = {"g": np.array([0.5])}
+    #     first_uni = distance_func(x_univ, y_univ, **param_dict)
+    #     second_uni = distance(x_univ, y_univ, metric=distance_str, **param_dict)
+    #     param_dict = {"g": np.array(list(range(10)))}
+    #     first_multi = distance_func(x_multi, y_multi, **param_dict)
+    #     second_multi = distance(x_multi, y_multi, metric=distance_str, **param_dict)
+    #     assert first_uni == pytest.approx(second_uni)
+    #     assert first_multi == pytest.approx(second_multi)
+
 
 BASIC_BOUNDING_PARAMS = [
     {"window": 0.2},
@@ -64,7 +80,7 @@ def _test_derivative(q: np.ndarray):
 
 DIST_PARAMS = {
     "dtw": BASIC_BOUNDING_PARAMS,
-    "erp": BASIC_BOUNDING_PARAMS + [{"g": 0.5}],
+    "erp": BASIC_BOUNDING_PARAMS + [{"g": 0.5}, {"g": None}],
     "edr": BASIC_BOUNDING_PARAMS + [{"epsilon": 0.5}],
     "lcss": BASIC_BOUNDING_PARAMS + [{"epsilon": 0.5}],
     "ddtw": BASIC_BOUNDING_PARAMS,