Repository: aeon-toolkit/aeon
Pull request #332 (2023-04-25 00:12:08)

[ENH] Update lcss 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 #323 

#### What does this implement/fix? Explain your changes.
This pr updates lcss 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 19deb1b7..753170a1 100644
--- a/aeon/distances/__init__.py
+++ b/aeon/distances/__init__.py
@@ -23,6 +23,10 @@ __all__ = [
     "wddtw_pairwise_distance",
     "wddtw_alignment_path",
     "wddtw_cost_matrix",
+    "lcss_distance",
+    "lcss_pairwise_distance",
+    "lcss_alignment_path",
+    "lcss_cost_matrix",
     "distance",
     "distance_factory",
     "pairwise_distance",
@@ -31,9 +35,7 @@ __all__ = [
     "edr_distance",
     "erp_distance",
     "msm_distance",
-    "lcss_distance",
     "twe_distance",
-    "lcss_alignment_path",
     "msm_alignment_path",
     "erp_alignment_path",
     "edr_alignment_path",
@@ -58,8 +60,6 @@ from aeon.distances._distance import (
     edr_distance,
     erp_alignment_path,
     erp_distance,
-    lcss_alignment_path,
-    lcss_distance,
     msm_alignment_path,
     msm_distance,
     pairwise_distance,
@@ -73,6 +73,12 @@ from aeon.distances._dtw import (
     dtw_pairwise_distance,
 )
 from aeon.distances._euclidean import euclidean_distance, euclidean_pairwise_distance
+from aeon.distances._lcss import (
+    lcss_alignment_path,
+    lcss_cost_matrix,
+    lcss_distance,
+    lcss_pairwise_distance,
+)
 from aeon.distances._squared import squared_distance, squared_pairwise_distance
 from aeon.distances._wddtw import (
     wddtw_alignment_path,
diff --git a/aeon/distances/_alignment_paths.py b/aeon/distances/_alignment_paths.py
index a18ba2fe..0e8a380a 100644
--- a/aeon/distances/_alignment_paths.py
+++ b/aeon/distances/_alignment_paths.py
@@ -6,10 +6,10 @@ from typing import List, Tuple
 import numpy as np
 from numba import njit
 
-from aeon.distances._squared import _univariate_squared_distance
+from aeon.distances._euclidean import _univariate_euclidean_distance
 
 
-@njit(cache=True)
+@njit(cache=True, fastmath=True)
 def compute_min_return_path(cost_matrix: np.ndarray) -> List[Tuple]:
     """Compute the minimum return path through a cost matrix.
 
@@ -56,7 +56,7 @@ def compute_min_return_path(cost_matrix: np.ndarray) -> List[Tuple]:
     return alignment[::-1]
 
 
-@njit(cache=True)
+@njit(cache=True, fastmath=True)
 def compute_lcss_return_path(
     x: np.ndarray,
     y: np.ndarray,
@@ -84,24 +84,25 @@ def compute_lcss_return_path(
     List[Tuple]
         List of indices that make up the return path.
     """
-    x_size, y_size = cost_matrix.shape
+    x_size = x.shape[1]
+    y_size = y.shape[1]
 
     i, j = (x_size, y_size)
     path = []
 
     while i > 0 and j > 0:
         if bounding_matrix[i - 1, j - 1]:
-            if _univariate_squared_distance(x[:, i - 1], y[:, j - 1]) <= epsilon:
+            if _univariate_euclidean_distance(x[:, i - 1], y[:, j - 1]) <= epsilon:
                 path.append((i - 1, j - 1))
                 i, j = (i - 1, j - 1)
-            elif cost_matrix[i - 1][j] > cost_matrix[i][j - 1]:
+            elif cost_matrix[i - 1, j] > cost_matrix[i, j - 1]:
                 i = i - 1
             else:
                 j = j - 1
     return path[::-1]
 
 
-@njit(cache=True)
+@njit(cache=True, fastmath=True)
 def _add_inf_to_out_of_bounds_cost_matrix(
     cost_matrix: np.ndarray, bounding_matrix: np.ndarray
 ) -> np.ndarray:
diff --git a/aeon/distances/_distance.py b/aeon/distances/_distance.py
index 32cb5036..9591ba3e 100644
--- a/aeon/distances/_distance.py
+++ b/aeon/distances/_distance.py
@@ -15,7 +15,11 @@ from aeon.distances._dtw import dtw_alignment_path, dtw_distance, dtw_pairwise_d
 from aeon.distances._edr import _EdrDistance
 from aeon.distances._erp import _ErpDistance
 from aeon.distances._euclidean import euclidean_distance, euclidean_pairwise_distance
-from aeon.distances._lcss import _LcssDistance
+from aeon.distances._lcss import (
+    lcss_alignment_path,
+    lcss_distance,
+    lcss_pairwise_distance,
+)
 from aeon.distances._msm import _MsmDistance
 from aeon.distances._numba_utils import (
     _compute_pairwise_distance,
@@ -201,72 +205,6 @@ def edr_distance(
     return distance(x, y, metric="edr", **format_kwargs)
 
 
-def lcss_distance(
-    x: np.ndarray,
-    y: np.ndarray,
-    window: Union[float, None] = None,
-    epsilon: float = 1.0,
-    **kwargs: Any,
-) -> float:
-    """Compute the longest common subsequence (LCSS) score between two time series.
-
-    LCSS attempts to find the longest common sequence between two time series and
-    returns a value that is the percentage that longest common sequence assumes.
-    Originally present in [1]_, LCSS is computed by matching indexes that are
-    similar up until a defined threshold (epsilon).
-
-    The value returned will be between 0.0 and 1.0, where 0.0 means the two time series
-    are exactly the same and 1.0 means they are complete opposites.
-
-    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.
-    epsilon : float, defaults = 1.
-        Matching threshold to determine if two subsequences are considered close
-        enough to be considered 'common'.
-    kwargs: Any
-        Extra kwargs.
-
-    Returns
-    -------
-    float
-        Lcss distance between x and y. The value returned will be between 0.0 and 1.0,
-        where 0.0 means the two time series are exactly the same and 1.0 means they
-        are complete opposites.
-
-    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] M. Vlachos, D. Gunopoulos, and G. Kollios. 2002. "Discovering
-        Similar Multidimensional Trajectories", In Proceedings of the
-        18th International Conference on Data Engineering (ICDE '02).
-        IEEE Computer Society, USA, 673.
-    """
-    format_kwargs = {
-        "window": window,
-        "epsilon": epsilon,
-    }
-    format_kwargs = {**format_kwargs, **kwargs}
-
-    return distance(x, y, metric="lcss", **format_kwargs)
-
-
 def msm_distance(
     x: np.ndarray,
     y: np.ndarray,
@@ -560,82 +498,6 @@ def erp_alignment_path(
     )
 
 
-def lcss_alignment_path(
-    x: np.ndarray,
-    y: np.ndarray,
-    return_cost_matrix: bool = False,
-    window: Union[float, None] = None,
-    epsilon: float = 1.0,
-    **kwargs: Any,
-) -> AlignmentPathReturn:
-    """Compute the longest common subsequence (LCSS) alignment path.
-
-    LCSS attempts to find the longest common sequence between two time series and
-    returns a value that is the percentage that longest common sequence assumes.
-    Originally present in [1]_, LCSS is computed by matching indexes that are
-    similar up until a defined threshold (epsilon).
-
-    The value returned will be between 0.0 and 1.0, where 0.0 means the two time series
-    are exactly the same and 1.0 means they are complete opposites.
-
-    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.
-    epsilon : float, defaults = 1.
-        Matching threshold to determine if two subsequences are considered close
-        enough to be considered 'common'.
-    kwargs: Any
-        Extra kwargs.
-
-    Returns
-    -------
-    list[tuple]
-        List of tuples containing the lcss alignment path.
-    float
-        Lcss distance between x and y. The value returned will be between 0.0 and 1.0,
-        where 0.0 means the two time series are exactly the same and 1.0 means they
-        are complete opposites.
-    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] M. Vlachos, D. Gunopoulos, and G. Kollios. 2002. "Discovering
-        Similar Multidimensional Trajectories", In Proceedings of the
-        18th International Conference on Data Engineering (ICDE '02).
-        IEEE Computer Society, USA, 673.
-    """
-    format_kwargs = {
-        "window": window,
-        "epsilon": epsilon,
-    }
-    format_kwargs = {**format_kwargs, **kwargs}
-
-    return distance_alignment_path(
-        x, y, metric="lcss", return_cost_matrix=return_cost_matrix, **format_kwargs
-    )
-
-
 def msm_alignment_path(
     x: np.ndarray,
     y: np.ndarray,
@@ -788,7 +650,7 @@ def twe_alignment_path(
     )
 
 
-NEW_DISTANCES = ["squared", "euclidean", "dtw", "ddtw", "wdtw", "wddtw"]
+NEW_DISTANCES = ["squared", "euclidean", "dtw", "ddtw", "wdtw", "wddtw", "lcss"]
 
 
 def distance(
@@ -883,6 +745,8 @@ def distance(
             return wdtw_distance(x, y, **kwargs)
         elif metric == "wddtw":
             return wddtw_distance(x, y, **kwargs)
+        elif metric == "lcss":
+            return lcss_distance(x, y, **kwargs)
     _x = to_numba_timeseries(x)
     _y = to_numba_timeseries(y)
 
@@ -962,6 +826,8 @@ def distance_factory(
             return wdtw_distance
         elif metric == "wddtw":
             return wddtw_distance
+        elif metric == "lcss":
+            return lcss_distance
     global dist_callable
 
     if x is None:
@@ -1090,6 +956,8 @@ def pairwise_distance(
             return wdtw_pairwise_distance(_x, _y, **kwargs)
         elif metric == "wddtw":
             return wddtw_pairwise_distance(_x, _y, **kwargs)
+        elif metric == "lcss":
+            return lcss_pairwise_distance(_x, _y, **kwargs)
 
     symmetric = np.array_equal(_x, _y)
     _metric_callable = _resolve_metric_to_factory(
@@ -1176,6 +1044,8 @@ def distance_alignment_path(
             return wdtw_alignment_path(x, y, **kwargs)
         elif metric == "wddtw":
             return wddtw_alignment_path(x, y, **kwargs)
+        elif metric == "lcss":
+            return lcss_alignment_path(x, y, **kwargs)
     _x = to_numba_timeseries(x)
     _y = to_numba_timeseries(y)
 
@@ -1259,6 +1129,8 @@ def distance_alignment_path_factory(
             return wdtw_alignment_path
         elif metric == "wddtw":
             return wddtw_alignment_path
+        elif metric == "lcss":
+            return lcss_alignment_path
     if x is None:
         x = np.zeros((1, 10))
     if y is None:
@@ -1295,13 +1167,6 @@ _METRIC_INFOS = [
         dist_instance=_EdrDistance(),
         dist_alignment_path_func=edr_alignment_path,
     ),
-    MetricInfo(
-        canonical_name="lcss",
-        aka={"lcss", "longest common subsequence"},
-        dist_func=lcss_distance,
-        dist_instance=_LcssDistance(),
-        dist_alignment_path_func=lcss_alignment_path,
-    ),
     MetricInfo(
         canonical_name="msm",
         aka={"msm", "move-split-merge"},
@@ -1328,7 +1193,6 @@ _METRICS_NAMES = list(_METRICS.keys())
 ALL_DISTANCES = (
     edr_distance,
     erp_distance,
-    lcss_distance,
     msm_distance,
     twe_distance,
 )
diff --git a/aeon/distances/_lcss.py b/aeon/distances/_lcss.py
index 4945efcd..3479c0f8 100644
--- a/aeon/distances/_lcss.py
+++ b/aeon/distances/_lcss.py
@@ -1,258 +1,408 @@
 # -*- coding: utf-8 -*-
+r"""Longest common subsequence (LCSS) between two time series.
+
+The LCSS distance for time series is based on the solution to the
+longest common subsequence problem in pattern matching [1]. The typical problem
+is to find the longest subsequence that is common to two discrete series based on
+the edit distance. This approach can be extended to consider real-valued time series
+by using a distance threshold epsilon, which defines the maximum difference
+between a pair of values that is allowed for them to be considered a match.
+LCSS finds the optimal alignment between two series by find the greatest number
+of matching pairs. The LCSS distance uses a matrix :math:'L' that records the
+sequence of matches over valid warpings. For two series :math:'a = a_1,... a_m
+and b = b_1,... b_m, L' is found by iterating over all valid windows (i.e.
+where bounding_matrix is not infinity, which by default is the constant band
+:math:'|i-j|<w*m', where :math:'w' is the window parameter value and m is series
+length), then calculating
+
+::math
+if(|a_i - b_j| < espilon) \\
+        &L_{i,j} \leftarrow L_{i-1,j-1}+1 \\
+else\\
+        &L_{i,j} \leftarrow \max(L_{i,j-1}, L_{i-1,j})\\
+
+The distance is an inverse function of the final LCSS.
+::math
+d_{LCSS}({\bf a,b}) = 1- \frac{LCSS({\bf a,b})}{m}.\]
+
+Note that series a and b need not be equal length.
+
+References
+----------
+.. [1] D. Hirschberg, Algorithms for the longest common subsequence problem, Journal
+of the ACM 24(4), 664--675, 1977
+"""
 __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 (
-    _add_inf_to_out_of_bounds_cost_matrix,
-    compute_lcss_return_path,
-)
+from aeon.distances._alignment_paths import compute_lcss_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 _LcssDistance(NumbaDistance):
-    r"""Longest common subsequence (LCSS) between two time series.
-
-    The LCSS distance for time series is based on the solution to the
-    longest common subsequence problem in pattern matching [1]. The typical problem
-    is to find the longest subsequence that is common to two discrete series based on
-    the edit distance. This approach can be extended to consider real-valued time series
-    by using a distance threshold epsilon, which defines the maximum difference
-    between a pair of values that is allowed for them to be considered a match.
-    LCSS finds the optimal alignment between two series by find the greatest number
-    of matching pairs. The LCSS distance uses a matrix :math:'L' that records the
-    sequence of matches over valid warpings. For two series :math:'a = a_1,... a_m
-    and b = b_1,... b_m, L' is found by iterating over all valid windows (i.e.
-    where bounding_matrix is not infinity, which by default is the constant band
-    :math:'|i-j|<w*m', where :math:'w' is the window parameter value and m is series
-    length), then calculating
-
-    ::math
-    if(|a_i - b_j| < espilon) \\
-            &L_{i,j} \leftarrow L_{i-1,j-1}+1 \\
-    else\\
-            &L_{i,j} \leftarrow \max(L_{i,j-1}, L_{i-1,j})\\
-
-    The distance is an inverse function of the final LCSS.
-    ::math
-    d_{LCSS}({\bf a,b}) = 1- \frac{LCSS({\bf a,b})}{m}.\]
-
-    Note that series a and b need not be equal length.
+from aeon.distances._euclidean import _univariate_euclidean_distance
 
-    References
-    ----------
-    .. [1] D. Hirschberg, Algorithms for the longest common subsequence problem, Journal
-    of the ACM 24(4), 664--675, 1977
-    """
 
-    def _distance_alignment_path_factory(
-        self,
-        x: np.ndarray,
-        y: np.ndarray,
-        return_cost_matrix: bool = False,
-        epsilon: float = 1.0,
-        window: float = None,
-        **kwargs: Any,
-    ) -> DistanceAlignmentPathCallable:
-        """Create a no_python compiled lcss 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.
-        epsilon : float, default = 1.
-            Matching threshold to determine if two subsequences are considered close
-            enough to be considered 'common'.
-        window: float, default = None, radius of the bounding window (if using
-        Sakoe-Chiba lower bounding). Must be between 0 and 1.
-        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 is not a numpy array.
-            If the input time series doesn't have exactly 2 dimensions.
-          If
-        """
-        _bounding_matrix = create_bounding_matrix(x.shape[1], y.shape[1], window)
-
-        if not isinstance(epsilon, float):
-            raise ValueError("The value of epsilon must be a float.")
-
-        if return_cost_matrix is True:
-
-            @njit(cache=True)
-            def numba_lcss_distance_alignment_path(
-                _x: np.ndarray,
-                _y: np.ndarray,
-            ) -> Tuple[List, float, np.ndarray]:
-                x_size = _x.shape[1]
-                y_size = _y.shape[1]
-                cost_matrix = _sequence_cost_matrix(_x, _y, _bounding_matrix, epsilon)
-                distance = 1 - float(
-                    cost_matrix[x_size - 1, y_size - 1] / min(x_size, y_size)
-                )
-                temp_cm = _add_inf_to_out_of_bounds_cost_matrix(
-                    cost_matrix, _bounding_matrix
-                )
-                path = compute_lcss_return_path(
-                    _x,
-                    _y,
-                    epsilon=epsilon,
-                    bounding_matrix=_bounding_matrix,
-                    cost_matrix=temp_cm,
-                )
-                return path, distance, cost_matrix
+@njit(cache=True, fastmath=True)
+def lcss_distance(
+    x: np.ndarray, y: np.ndarray, window: float = None, epsilon: float = 1.0
+) -> float:
+    r"""Return the lcss distance between x and y.
 
-        else:
+    LCSS attempts to find the longest common sequence between two time series and
+    returns a value that is the percentage that longest common sequence assumes.
+    Originally present in [1]_, LCSS is computed by matching indexes that are
+    similar up until a defined threshold (epsilon).
+
+    The value returned will be between 0.0 and 1.0, where 0.0 means the two time series
+    are exactly the same and 1.0 means they are complete opposites.
+
+    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.
+    epsilon: float, defaults=1.
+        Matching threshold to determine if two subsequences are considered close
+        enough to be considered 'common'. The default is 1.
+
+    Returns
+    -------
+    float
+        The lcss distance between x and y.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aeon.distances import lcss_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 = lcss_distance(x, y)
+
+    Raises
+    ------
+    ValueError
+        If x and y are not 1D, 2D, or 3D arrays.
 
-            @njit(cache=True)
-            def numba_lcss_distance_alignment_path(
-                _x: np.ndarray,
-                _y: np.ndarray,
-            ) -> Tuple[List, float]:
-                x_size = _x.shape[1]
-                y_size = _y.shape[1]
-                cost_matrix = _sequence_cost_matrix(_x, _y, _bounding_matrix, epsilon)
-                distance = 1 - float(
-                    cost_matrix[x_size - 1, y_size - 1] / min(x_size, y_size)
-                )
-                temp_cm = _add_inf_to_out_of_bounds_cost_matrix(
-                    cost_matrix, _bounding_matrix
-                )
-                path = compute_lcss_return_path(
-                    _x,
-                    _y,
-                    epsilon=epsilon,
-                    bounding_matrix=_bounding_matrix,
-                    cost_matrix=temp_cm,
-                )
-                return path, distance
-
-        return numba_lcss_distance_alignment_path
-
-    def _distance_factory(
-        self,
-        x: np.ndarray,
-        y: np.ndarray,
-        epsilon: float = 1.0,
-        window: float = None,
-        **kwargs: Any,
-    ) -> DistanceCallable:
-        """Create a no_python compiled lcss distance 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.
-        epsilon : float, default = 1.
-            Matching threshold to determine if two subsequences are considered close
-            enough to be considered 'common'.
-        window: float, default = None, radius of the bounding window (if using
-        Sakoe-Chiba lower bounding). Must be between 0 and 1.
-        kwargs: Any Extra kwargs.
-
-        Returns
-        -------
-        Callable[[np.ndarray, np.ndarray], float]
-            No_python compiled lcss 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
-        """
-        _bounding_matrix = create_bounding_matrix(x.shape[1], y.shape[1], window)
-
-        if not isinstance(epsilon, float):
-            raise ValueError("The value of epsilon must be a float.")
-
-        @njit(cache=True)
-        def numba_lcss_distance(
-            _x: np.ndarray,
-            _y: np.ndarray,
-        ) -> float:
-            x_size = _x.shape[1]
-            y_size = _y.shape[1]
-            cost_matrix = _sequence_cost_matrix(_x, _y, _bounding_matrix, epsilon)
-            return 1 - float(cost_matrix[x_size - 1, y_size - 1] / min(x_size, y_size))
-
-        return numba_lcss_distance
-
-
-@njit(cache=True)
-def _sequence_cost_matrix(
-    x: np.ndarray,
-    y: np.ndarray,
-    bounding_matrix: np.ndarray,
-    epsilon: float,
-):
-    """Compute the lcss cost matrix between two time series.
+    References
+    ----------
+    .. [1] M. Vlachos, D. Gunopoulos, and G. Kollios. 2002. "Discovering
+        Similar Multidimensional Trajectories", In Proceedings of the
+        18th International Conference on Data Engineering (ICDE '02).
+        IEEE Computer Society, USA, 673.
+    """
+    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 _lcss_distance(_x, _y, bounding_matrix, epsilon)
+    if x.ndim == 2 and y.ndim == 2:
+        bounding_matrix = create_bounding_matrix(x.shape[1], y.shape[1], window)
+        return _lcss_distance(x, y, bounding_matrix, epsilon)
+    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 += _lcss_distance(curr_x, curr_y, bounding_matrix, epsilon)
+        return distance
+    raise ValueError("x and y must be 1D, 2D, or 3D arrays")
+
+
+@njit(cache=True, fastmath=True)
+def lcss_cost_matrix(
+    x: np.ndarray, y: np.ndarray, window: float = None, epsilon: float = 1.0
+) -> np.ndarray:
+    r"""Return the lcss cost matrix between x and y.
 
     Parameters
     ----------
-    x: np.ndarray (2d array), first time series.
-    y: np.ndarray (2d array), 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.
-    epsilon : float
-        Matching threshold to determine if distance between two subsequences are
-        considered similar (similar if distance less than the threshold).
+    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.
+    epsilon: float, defaults=1.
+        Matching threshold to determine if two subsequences are considered close
+        enough to be considered 'common'. The default is 1.
 
     Returns
     -------
-    np.ndarray (2d of size mxn where m is len(x) and n is len(y))
-        Lcss cost matrix between x and y.
+    np.ndarray
+        The lcss 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 lcss_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]])
+    >>> lcss_cost_matrix(x, y)
+    array([[ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
+           [ 0.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.],
+           [ 0.,  1.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.,  2.],
+           [ 0.,  1.,  2.,  3.,  3.,  3.,  3.,  3.,  3.,  3.,  3.],
+           [ 0.,  1.,  2.,  3.,  4.,  4.,  4.,  4.,  4.,  4.,  4.],
+           [ 0.,  1.,  2.,  3.,  4.,  5.,  5.,  5.,  5.,  5.,  5.],
+           [ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  6.,  6.,  6.,  6.],
+           [ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  7.,  7.,  7.],
+           [ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  8.,  8.],
+           [ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9.,  9.],
+           [ 0.,  1.,  2.,  3.,  4.,  5.,  6.,  7.,  8.,  9., 10.]])
     """
-    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 _lcss_cost_matrix(_x, _y, bounding_matrix, epsilon)
+    if x.ndim == 2 and y.ndim == 2:
+        bounding_matrix = create_bounding_matrix(x.shape[1], y.shape[1], window)
+        return _lcss_cost_matrix(x, y, bounding_matrix, epsilon)
+    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] + 1, y.shape[2] + 1))
+        for curr_x, curr_y in zip(x, y):
+            cost_matrix = np.add(
+                cost_matrix, _lcss_cost_matrix(curr_x, curr_y, bounding_matrix, epsilon)
+            )
+        return cost_matrix
+    raise ValueError("x and y must be 1D, 2D, or 3D arrays")
+
+
+@njit(cache=True, fastmath=True)
+def _lcss_distance(
+    x: np.ndarray, y: np.ndarray, bounding_matrix: np.ndarray, epsilon: float
+) -> float:
+    distance = _lcss_cost_matrix(x, y, bounding_matrix, epsilon)[x.shape[1], y.shape[1]]
+    return 1 - (float(distance / min(x.shape[1], y.shape[1])))
+
+
+@njit(cache=True, fastmath=True)
+def _lcss_cost_matrix(
+    x: np.ndarray, y: np.ndarray, bounding_matrix: np.ndarray, epsilon
+) -> np.ndarray:
     x_size = x.shape[1]
     y_size = y.shape[1]
+
     cost_matrix = np.zeros((x_size + 1, y_size + 1))
+
     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]) ** 2
-                curr_dist = np.sqrt(curr_dist)
-                if curr_dist <= epsilon:
+                if _univariate_euclidean_distance(x[:, i - 1], y[:, j - 1]) <= epsilon:
                     cost_matrix[i, j] = 1 + cost_matrix[i - 1, j - 1]
                 else:
                     cost_matrix[i, j] = max(
                         cost_matrix[i, j - 1], cost_matrix[i - 1, j]
                     )
-    return cost_matrix[1:, 1:]
+    return cost_matrix
+
+
+@njit(cache=True, fastmath=True)
+def lcss_pairwise_distance(
+    X: np.ndarray, y: np.ndarray = None, window: float = None, epsilon: float = 1.0
+) -> np.ndarray:
+    """Compute the lcss 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.
+    epsilon: float, defaults=1.
+        Matching threshold to determine if two subsequences are considered close
+        enough to be considered 'common'. The default is 1.
+
+    Returns
+    -------
+    np.ndarray (n_instances, n_instances)
+        lcss 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 lcss_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]]])
+    >>> lcss_pairwise_distance(X)
+    array([[0.        , 0.66666667, 1.        ],
+           [0.66666667, 0.        , 0.66666667],
+           [1.        , 0.66666667, 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]]])
+    >>> lcss_pairwise_distance(X, y)
+    array([[1., 1., 1.],
+           [1., 1., 1.],
+           [1., 1., 1.]])
+
+    >>> 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]])
+    >>> lcss_pairwise_distance(X, y_univariate)
+    array([[1.],
+           [1.],
+           [1.]])
+    """
+    if y is None:
+        # To self
+        if X.ndim == 3:
+            return _lcss_pairwise_distance(X, window, epsilon)
+        if X.ndim == 2:
+            _X = X.reshape((X.shape[0], 1, X.shape[1]))
+            return _lcss_pairwise_distance(_X, window, epsilon)
+        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 _lcss_from_multiple_to_multiple_distance(X, y, window, epsilon)
+        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 _lcss_from_multiple_to_multiple_distance(_x, _y, window, epsilon)
+        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 _lcss_from_multiple_to_multiple_distance(_x, _y, window, epsilon)
+        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 _lcss_from_multiple_to_multiple_distance(X, _y, window, epsilon)
+        if y.ndim == 3 and X.ndim == 2:
+            _x = X.reshape((1, X.shape[0], X.shape[1]))
+            return _lcss_from_multiple_to_multiple_distance(_x, y, window, epsilon)
+        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 _lcss_from_multiple_to_multiple_distance(_x, _y, window, epsilon)
+        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 _lcss_from_multiple_to_multiple_distance(_x, _y, window, epsilon)
+        else:
+            raise ValueError("x and y must be 2D or 3D arrays")
+
+
+@njit(cache=True, fastmath=True)
+def _lcss_pairwise_distance(X: np.ndarray, window: float, epsilon: 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] = _lcss_distance(X[i], X[j], bounding_matrix, epsilon)
+            distances[j, i] = distances[i, j]
+
+    return distances
+
+
+@njit(cache=True, fastmath=True)
+def _lcss_from_multiple_to_multiple_distance(
+    x: np.ndarray, y: np.ndarray, window: float, epsilon: 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] = _lcss_distance(x[i], y[j], bounding_matrix, epsilon)
+    return distances
+
+
+@njit(cache=True, fastmath=True)
+def lcss_alignment_path(
+    x: np.ndarray, y: np.ndarray, window: float = None, epsilon: float = 1.0
+) -> Tuple[List[Tuple[int, int]], float]:
+    """Compute the lcss alignment path between two time series.
+
+    Parameters
+    ----------
+    x: np.ndarray, of shape (n_channels, n_timepoints) or (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.
+    epsilon: float, defaults=1.
+        Matching threshold to determine if two subsequences are considered close
+        enough to be considered 'common'. The default is 1.
+
+
+    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 lcss distance between the two time series.
+
+    Raises
+    ------
+    ValueError
+        If x and y are not 1D, or 2D arrays.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aeon.distances import lcss_alignment_path
+    >>> x = np.array([[1, 2, 3, 6]])
+    >>> y = np.array([[1, 2, 3, 4]])
+    >>> path, dist = lcss_alignment_path(x, y)
+    >>> path
+    [(0, 0), (1, 1), (2, 2)]
+    """
+    x_size = x.shape[-1]
+    y_size = y.shape[-1]
+    bounding_matrix = create_bounding_matrix(x_size, y_size, window)
+    if x.ndim == 1 and y.ndim == 1:
+        _x = x.reshape((1, x.shape[0]))
+        _y = y.reshape((1, y.shape[0]))
+        cost_matrix = _lcss_cost_matrix(_x, _y, bounding_matrix, epsilon)
+        distance = 1 - (float(cost_matrix[x_size, y_size] / min(x_size, y_size)))
+        return (
+            compute_lcss_return_path(_x, _y, epsilon, bounding_matrix, cost_matrix),
+            distance,
+        )
+    if x.ndim == 2 and y.ndim == 2:
+        cost_matrix = _lcss_cost_matrix(x, y, bounding_matrix, epsilon)
+        distance = 1 - (float(cost_matrix[x_size, y_size] / min(x_size, y_size)))
+        return (
+            compute_lcss_return_path(x, y, epsilon, bounding_matrix, cost_matrix),
+            distance,
+        )
+    raise ValueError("x and y must be 1D, or 2D arrays")


Test patch:
diff --git a/aeon/distances/tests/_expected_results.py b/aeon/distances/tests/_expected_results.py
index 90fc9981..65122a90 100644
--- a/aeon/distances/tests/_expected_results.py
+++ b/aeon/distances/tests/_expected_results.py
@@ -16,7 +16,7 @@ _expected_distance_results = {
     ],
     "erp": [5.0, 5.037672883414786, 20.724482073800456],
     "edr": [1.0, 0.6, 1.0],
-    "lcss": [1.0, 0.09999999999999998, 1.0],
+    "lcss": [1.0, 0.1, 1.0, 9.8, 0.0, 1.0, 9.4],
     "squared": [
         25.0,
         6.932617853961479,
diff --git a/aeon/distances/tests/test_alignment_path.py b/aeon/distances/tests/test_alignment_path.py
index 8479ff01..8c992f8a 100644
--- a/aeon/distances/tests/test_alignment_path.py
+++ b/aeon/distances/tests/test_alignment_path.py
@@ -30,7 +30,7 @@ def test_cost_matrix(dist):
         return
 
     # Test univariate
-    if dist["name"] != "ddtw" and dist["name"] != "wddtw":
+    if dist["name"] != "ddtw" and dist["name"] != "wddtw" and dist["name"] != "lcss":
         _validate_cost_matrix_result(
             np.array([10.0]),
             np.array([15.0]),
@@ -64,13 +64,14 @@ def test_cost_matrix(dist):
         dist["alignment_path"],
     )
 
-    _validate_cost_matrix_result(
-        create_test_distance_numpy(10, 10, 10),
-        create_test_distance_numpy(10, 10, 10, random_state=2),
-        dist["name"],
-        dist["distance"],
-        dist["alignment_path"],
-    )
+    if dist["name"] != "lcss":
+        _validate_cost_matrix_result(
+            create_test_distance_numpy(10, 10, 10),
+            create_test_distance_numpy(10, 10, 10, random_state=2),
+            dist["name"],
+            dist["distance"],
+            dist["alignment_path"],
+        )
 
     # Test unequal length
     _validate_cost_matrix_result(
@@ -89,10 +90,11 @@ def test_cost_matrix(dist):
         dist["alignment_path"],
     )
 
-    _validate_cost_matrix_result(
-        create_test_distance_numpy(10, 10, 5),
-        create_test_distance_numpy(10, 10, 10, random_state=2),
-        dist["name"],
-        dist["distance"],
-        dist["alignment_path"],
-    )
+    if dist["name"] != "lcss":
+        _validate_cost_matrix_result(
+            create_test_distance_numpy(10, 10, 5),
+            create_test_distance_numpy(10, 10, 10, random_state=2),
+            dist["name"],
+            dist["distance"],
+            dist["alignment_path"],
+        )
diff --git a/aeon/distances/tests/test_cost_matrix.py b/aeon/distances/tests/test_cost_matrix.py
index 1fd9473f..290bb238 100644
--- a/aeon/distances/tests/test_cost_matrix.py
+++ b/aeon/distances/tests/test_cost_matrix.py
@@ -19,12 +19,22 @@ def _validate_cost_matrix_result(
     assert isinstance(cost_matrix_result, np.ndarray)
     if name == "ddtw" or name == "wddtw":
         assert cost_matrix_result.shape == (x.shape[-1] - 2, y.shape[-1] - 2)
+    elif name == "lcss":
+        # lcss cm is one larger than the input
+        assert cost_matrix_result.shape == (x.shape[-1] + 1, y.shape[-1] + 1)
     else:
         assert cost_matrix_result.shape == (x.shape[-1], y.shape[-1])
 
     distance_result = distance(x, y)
 
-    assert_almost_equal(cost_matrix_result[-1, -1], distance_result)
+    if name == "lcss":
+        if x.ndim != 3:
+            distance = 1 - (
+                float(cost_matrix_result[-1, -1] / min(x.shape[-1], y.shape[-1]))
+            )
+            assert_almost_equal(distance, distance_result)
+    else:
+        assert_almost_equal(cost_matrix_result[-1, -1], distance_result)
 
 
 @pytest.mark.parametrize("dist", DISTANCES)
diff --git a/aeon/distances/tests/test_new_distances.py b/aeon/distances/tests/test_new_distances.py
index 63c86360..83a9f789 100644
--- a/aeon/distances/tests/test_new_distances.py
+++ b/aeon/distances/tests/test_new_distances.py
@@ -17,6 +17,10 @@ from aeon.distances import (
     dtw_pairwise_distance,
     euclidean_distance,
     euclidean_pairwise_distance,
+    lcss_alignment_path,
+    lcss_cost_matrix,
+    lcss_distance,
+    lcss_pairwise_distance,
     squared_distance,
     squared_pairwise_distance,
     wddtw_alignment_path,
@@ -70,6 +74,13 @@ DISTANCES = [
         "cost_matrix": wddtw_cost_matrix,
         "alignment_path": wddtw_alignment_path,
     },
+    {
+        "name": "lcss",
+        "distance": lcss_distance,
+        "pairwise_distance": lcss_pairwise_distance,
+        "cost_matrix": lcss_cost_matrix,
+        "alignment_path": lcss_alignment_path,
+    },
 ]
 
 
diff --git a/aeon/distances/tests/test_pairwise.py b/aeon/distances/tests/test_pairwise.py
index 60686808..32d8b0c5 100644
--- a/aeon/distances/tests/test_pairwise.py
+++ b/aeon/distances/tests/test_pairwise.py
@@ -1,6 +1,7 @@
 # -*- coding: utf-8 -*-
 import numpy as np
 import pytest
+from numpy.testing import assert_almost_equal
 
 from aeon.distances.tests._utils import _make_3d_series, create_test_distance_numpy
 from aeon.distances.tests.test_new_distances import DISTANCES
@@ -83,7 +84,7 @@ def _validate_single_to_multiple_result(
         curr_y = y[i]
         curr = single_to_multiple_result[0, i]
         dist = distance(x, curr_y)
-        assert dist == curr
+        assert_almost_equal(dist, curr)
 
 
 @pytest.mark.parametrize("dist", DISTANCES)