Repository: aeon-toolkit/aeon
Pull request #2589 (2025-03-05 05:51:21)

Clarify documentation regarding unequal length series limitation
<!--
Thanks for contributing a pull request! Please ensure you have taken a look at our
contribution guide: https://www.aeon-toolkit.org/en/latest/contributing.html.

Feel free to delete sections of this template if they do not apply to your PR,
avoid submitting a blank template or empty sections.
If you are a new contributor, do not delete this template without a suitable
replacement or reason. If in doubt, ask for help. We're here to help!

Please be aware that we are a team of volunteers so patience is
necessary when waiting for a review or reply. There may not be a quick turnaround for
reviews during slow periods. While we value all contributions big or small, pull
requests which do not follow our guidelines may be closed.
-->

#### Reference Issues/PRs

Fixes: #2444
<!--
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 the documentation to clarify the supported input data formats for the collection transformers module
<!--
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.
-->

#### Any other comments?

<!--
Any other information that is important to this PR or helpful for reviewers.
-->

### PR checklist

<!--
Please go through the checklist below. Please feel free to remove points if they are
not applicable. To check a box, replace the space inside the square brackets with an
'x' i.e. [x].
-->

##### For all contributions
- [ ] I've added myself to the [list of contributors](https://github.com/aeon-toolkit/aeon/blob/main/.all-contributorsrc). Alternatively, you can use the [@all-contributors](https://allcontributors.org/docs/en/bot/usage) bot to do this for you after the PR has been merged.
- [x] The PR title starts with either [ENH], [MNT], [DOC], [BUG], [REF], [DEP] or [GOV] indicating whether the PR topic is related to enhancement, maintenance, documentation, bugs, refactoring, deprecation or governance.

##### For new estimators and functions
- [ ] I've added the estimator/function to the online [API documentation](https://www.aeon-toolkit.org/en/latest/api_reference.html).
- [ ] (OPTIONAL) I've added myself as a `__maintainer__` at the top of relevant files and want to be contacted regarding its maintenance. Unmaintained files may be removed. This is for the full file, and you should not add yourself if you are just making minor changes or do not want to help maintain its contents.

##### For developers with write access
- [ ] (OPTIONAL) I've updated aeon's [CODEOWNERS](https://github.com/aeon-toolkit/aeon/blob/main/CODEOWNERS) to receive notifications about future changes to these files.


<!--
Thanks for contributing!
-->


aeon-actions-bot[bot]: 
## Thank you for contributing to `aeon`

I did not find any labels to add based on the title. Please add the [ENH], [MNT], [BUG], [DOC], [REF], [DEP] and/or [GOV] tags to your pull requests titles. For now you can add the labels manually.
I have added the following labels to this PR based on the changes made: **[ $\\color{#41A8F6}{\\textsf{transformations}}$ ]**. Feel free to change these if they do not properly represent the PR.

The [Checks](https://github.com/aeon-toolkit/aeon/pull/2589/checks) tab will show the status of our automated tests. You can click on individual test runs in the tab or "Details" in the panel below to see more information if there is a failure.

If our `pre-commit` code quality check fails, any trivial fixes will automatically be pushed to your PR unless it is a draft.

Don't hesitate to ask questions on the `aeon` [Slack](https://join.slack.com/t/aeon-toolkit/shared_invite/zt-22vwvut29-HDpCu~7VBUozyfL_8j3dLA) channel if you have any.

### PR CI actions

These checkboxes will add labels to enable/disable CI functionality for this PR. This may not take effect immediately, and a new commit may be required to run the new configuration.

- [ ] Run `pre-commit` checks for all files
- [ ] Run `mypy` typecheck tests
- [ ] Run all `pytest` tests and configurations
- [ ] Run all notebook example tests
- [ ] Run numba-disabled `codecov` tests
- [ ] Stop automatic `pre-commit` fixes (always disabled for drafts)
- [ ] Disable numba cache loading
- [ ] Push an empty commit to re-run CI checks
    
MatthewMiddlehurst: My bad, see latest issue comment.
MatthewMiddlehurst: The first PR was more accurate, but we do not need a large warning IMO. Perhaps just tidy the first paragraph and change the part about np-list to reflect that it does not have to be unequal length (but that is its primary usage). The capabilities paragraph below that is sufficient but if you think it can be improved go ahead.
Kaustbh: @MatthewMiddlehurst I have made the changes please review and let me know if any changes are needed.
Kaustbh: kept this `(number of channels, series length)` so that the user can know the shape of the inside array.
MatthewMiddlehurst: Please try to avoid force pushing like this if you can.
MatthewMiddlehurst: or give us a hand in finding out why doing so breaks pre-commit testing i guess 🙂 

Patch:
diff --git a/.github/workflows/pr_precommit.yml b/.github/workflows/pr_precommit.yml
index 2f63ef2b..a74e8b12 100644
--- a/.github/workflows/pr_precommit.yml
+++ b/.github/workflows/pr_precommit.yml
@@ -37,7 +37,7 @@ jobs:
           python-version: "3.10"
 
       - name: Get changed files
-        uses: tj-actions/changed-files@v45
+        uses: tj-actions/changed-files@v46.0.1
         id: changed-files
 
       - name: List changed files
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 9f455ec6..854da8df 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -29,7 +29,7 @@ repos:
         args: [ "--create", "--python-folders", "aeon" ]
 
   - repo: https://github.com/astral-sh/ruff-pre-commit
-    rev: v0.9.10
+    rev: v0.11.0
     hooks:
       - id: ruff
         args: [ "--fix"]
diff --git a/aeon/benchmarking/metrics/anomaly_detection/__init__.py b/aeon/benchmarking/metrics/anomaly_detection/__init__.py
index fdbf13ce..cf6ccac4 100644
--- a/aeon/benchmarking/metrics/anomaly_detection/__init__.py
+++ b/aeon/benchmarking/metrics/anomaly_detection/__init__.py
@@ -14,6 +14,9 @@ __all__ = [
     "range_pr_auc_score",
     "range_pr_vus_score",
     "range_roc_vus_score",
+    "ts_precision",
+    "ts_recall",
+    "ts_fscore",
 ]
 
 from aeon.benchmarking.metrics.anomaly_detection._binary import (
@@ -35,3 +38,8 @@ from aeon.benchmarking.metrics.anomaly_detection._vus_metrics import (
     range_roc_auc_score,
     range_roc_vus_score,
 )
+from aeon.benchmarking.metrics.anomaly_detection.range_metrics import (
+    ts_fscore,
+    ts_precision,
+    ts_recall,
+)
diff --git a/aeon/benchmarking/metrics/anomaly_detection/_binary.py b/aeon/benchmarking/metrics/anomaly_detection/_binary.py
index 85d54a5c..085d7c04 100644
--- a/aeon/benchmarking/metrics/anomaly_detection/_binary.py
+++ b/aeon/benchmarking/metrics/anomaly_detection/_binary.py
@@ -6,11 +6,19 @@ __all__ = ["range_precision", "range_recall", "range_f_score"]
 import warnings
 
 import numpy as np
+from deprecated.sphinx import deprecated
 
 from aeon.benchmarking.metrics.anomaly_detection._util import check_y
 from aeon.utils.validation._dependencies import _check_soft_dependencies
 
 
+# TODO: Remove in v1.2.0
+@deprecated(
+    version="1.1.0",
+    reason="range_precision is deprecated and will be removed in v1.2.0. "
+    "Please use ts_precision from the range_metrics module instead.",
+    category=FutureWarning,
+)
 def range_precision(
     y_true: np.ndarray,
     y_pred: np.ndarray,
@@ -70,6 +78,13 @@ def range_precision(
     return ts_precision(y_true, y_pred, alpha=alpha, cardinality=cardinality, bias=bias)
 
 
+# TODO: Remove in v1.2.0
+@deprecated(
+    version="1.1.0",
+    reason="range_recall is deprecated and will be removed in v1.2.0. "
+    "Please use ts_recall from the range_metrics module instead.",
+    category=FutureWarning,
+)
 def range_recall(
     y_true: np.ndarray,
     y_pred: np.ndarray,
@@ -131,6 +146,13 @@ def range_recall(
     return ts_recall(y_true, y_pred, alpha=alpha, cardinality=cardinality, bias=bias)
 
 
+# TODO: Remove in v1.2.0
+@deprecated(
+    version="1.1.0",
+    reason="range_f_score is deprecated and will be removed in v1.2.0. "
+    "Please use ts_fscore from the range_metrics module instead.",
+    category=FutureWarning,
+)
 def range_f_score(
     y_true: np.ndarray,
     y_pred: np.ndarray,
diff --git a/aeon/benchmarking/metrics/anomaly_detection/range_metrics.py b/aeon/benchmarking/metrics/anomaly_detection/range_metrics.py
new file mode 100644
index 00000000..9084188f
--- /dev/null
+++ b/aeon/benchmarking/metrics/anomaly_detection/range_metrics.py
@@ -0,0 +1,521 @@
+"""Calculate Precision, Recall, and F1-Score for time series anomaly detection."""
+
+__maintainer__ = []
+__all__ = ["ts_precision", "ts_recall", "ts_fscore"]
+
+import numpy as np
+
+
+def _flatten_ranges(ranges):
+    """
+    If the input is a list of lists, it flattens it into a single list.
+
+    Parameters
+    ----------
+    ranges : list of tuples or list of lists of tuples
+        The ranges to flatten. each tuple shoulod be in the format of (start, end).
+
+    Returns
+    -------
+    list of tuples
+        A flattened list of ranges.
+
+    Examples
+    --------
+    >>> _flatten_ranges([[(1, 5), (10, 15)], [(20, 25)]])
+    [(1, 5), (10, 15), (20, 25)]
+    """
+    if not ranges:
+        return []
+    if isinstance(ranges[0], list):
+        flat = []
+        for sublist in ranges:
+            for pred in sublist:
+                flat.append(pred)
+        return flat
+    return ranges
+
+
+def udf_gamma_def(overlap_count):
+    """User-defined gamma function. Should return a gamma value > 1.
+
+    Parameters
+    ----------
+    overlap_count : int
+        The number of overlapping ranges.
+
+    Returns
+    -------
+    float
+        The user-defined gamma value (>1).
+    """
+    return_val = 1 + 0.1 * overlap_count  # modify this function as needed
+
+    return return_val
+
+
+def _calculate_bias(position, length, bias_type="flat"):
+    """Calculate bias value based on position and length.
+
+    Parameters
+    ----------
+    position : int
+        Current position in the range
+    length : int
+        Total length of the range
+    bias_type : str, default="flat"
+        Type of bias to apply, Should be one of ["flat", "front", "middle", "back"].
+    """
+    if bias_type == "flat":
+        return 1.0
+    elif bias_type == "front":
+        return 1.0 - (position - 1) / length
+    elif bias_type == "middle":
+        if length / 2 == 0:
+            return 1.0
+        if position <= length / 2:
+            return position / (length / 2)
+        else:
+            return (length - position + 1) / (length / 2)
+    elif bias_type == "back":
+        return position / length
+    else:
+        raise ValueError(f"Invalid bias type: {bias_type}")
+
+
+def _gamma_select(cardinality, gamma):
+    """Select a gamma value based on the cardinality type.
+
+    Parameters
+    ----------
+    cardinality : int
+        The number of overlapping ranges.
+    gamma : str
+        Gamma to use. Should be one of ["one", "reciprocal", "udf_gamma"].
+
+    Returns
+    -------
+    float
+        The selected gamma value.
+
+    Raises
+    ------
+    ValueError
+        If an invalid `gamma` type is provided or if `udf_gamma` is required
+        but not provided.
+    """
+    if gamma == "one":
+        return 1.0
+    elif gamma == "reciprocal":
+        return 1 / cardinality if cardinality > 1 else 1.0
+    elif gamma == "udf_gamma":
+        if udf_gamma_def(cardinality) is not None:
+            return 1.0 / udf_gamma_def(cardinality)
+        else:
+            raise ValueError("udf_gamma must be provided for 'udf_gamma' gamma type.")
+    else:
+        raise ValueError(
+            "Invalid gamma type. Choose from ['one', 'reciprocal', 'udf_gamma']."
+        )
+
+
+def _calculate_overlap_reward_precision(pred_range, overlap_set, bias_type):
+    """Overlap Reward for y_pred.
+
+    Parameters
+    ----------
+    pred_range : tuple
+        The predicted range.
+    overlap_set : set
+        The set of overlapping positions.
+    bias_type : str
+        Type of bias to apply, Should be one of ["flat", "front", "middle", "back"].
+
+    Returns
+    -------
+    float
+        The weighted value for overlapping positions only.
+    """
+    start, end = pred_range
+    length = end - start + 1
+
+    max_value = 0  # Total possible weighted value for all positions.
+    my_value = 0  # Weighted value for overlapping positions only.
+
+    for i in range(1, length + 1):
+        global_position = start + i - 1
+        bias_value = _calculate_bias(i, length, bias_type)
+        max_value += bias_value
+
+        if global_position in overlap_set:
+            my_value += bias_value
+
+    return my_value / max_value if max_value > 0 else 0.0
+
+
+def _calculate_overlap_reward_recall(real_range, overlap_set, bias_type):
+    """Overlap Reward for y_real.
+
+    Parameters
+    ----------
+    real_range : tuple
+        The real range.
+    overlap_set : set
+        The set of overlapping positions.
+    bias_type : str
+        Type of bias to apply, Should be one of ["flat", "front", "middle", "back"].
+
+    Returns
+    -------
+    float
+        The weighted value for overlapping positions only.
+    """
+    start, end = real_range
+    length = end - start + 1
+
+    max_value = 0.0  # Total possible weighted value for all positions.
+    my_value = 0.0  # Weighted value for overlapping positions only.
+
+    for i in range(1, length + 1):
+        global_position = start + i - 1
+        bias_value = _calculate_bias(i, length, bias_type)
+        max_value += bias_value
+
+        if global_position in overlap_set:
+            my_value += bias_value
+
+    return my_value / max_value if max_value > 0 else 0.0
+
+
+def _binary_to_ranges(binary_sequence):
+    """
+    Convert a binary sequence to a list of anomaly ranges.
+
+    Parameters
+    ----------
+    binary_sequence : list
+        Binary sequence where 1 indicates anomaly and 0 indicates normal.
+
+    Returns
+    -------
+    list of tuples
+        List of anomaly ranges as (start, end) tuples.
+
+    """
+    ranges = []
+    start = None
+
+    for i, val in enumerate(binary_sequence):
+        if val and start is None:
+            start = i
+        elif not val and start is not None:
+            ranges.append((start, i - 1))
+            start = None
+
+    if start is not None:
+        ranges.append((start, len(binary_sequence) - 1))
+
+    return ranges
+
+
+def ts_precision(y_pred, y_real, gamma="one", bias_type="flat"):
+    """
+    Calculate Precision for time series anomaly detection.
+
+    Precision measures the proportion of correctly predicted anomaly positions
+    out of all all the predicted anomaly positions, aggregated across the entire time
+    series.
+
+    Parameters
+    ----------
+    y_pred : list of tuples or binary sequence
+        The predicted anomaly ranges.
+        - For range-based input, each tuple represents a range (start, end) of the
+          anomaly where start is starting index (inclusive) and end is ending index
+          (inclusive).
+        - For binary inputs, the sequence should contain integers (0 or 1), where 1
+          indicates an anomaly. In this case, y_pred and y_real must be of same length.
+    y_real : list of tuples, list of lists of tuples or binary sequence
+        The real/actual (ground truth) ranges.
+        - For range-based input, each tuple represents a range (start, end) of the
+          anomaly where start is starting index (inclusive) and end is ending index
+          (inclusive).
+        - If y_real is in the format of list of lists, they will be flattened into a
+          single list of tuples bringing it to the above format.
+        - For binary inputs, the sequence should contain integers (0 or 1), where 1
+          indicates an anomaly. In this case, y_pred and y_real must be of same length.
+    bias_type : str, default="flat"
+        Type of bias to apply. Should be one of ["flat", "front", "middle", "back"].
+    gamma : str, default="one"
+        Cardinality type. Should be one of ["reciprocal", "one"].
+
+    Returns
+    -------
+    float
+        Precision
+
+    Raises
+    ------
+    ValueError
+        If an invalid `gamma` type is provided.
+    ValueError
+        If input sequence is binary and y_real and y_pred are of different lengths.
+
+    References
+    ----------
+    .. [1] Tatbul, Nesime, Tae Jun Lee, Stan Zdonik, Mejbah Alam,and Justin Gottschlich.
+       "Precision and Recall for Time Series." 32nd Conference on Neural Information
+       Processing Systems (NeurIPS 2018), Montréal, Canada.
+       http://papers.nips.cc/paper/7462-precision-and-recall-for-time-series.pdf
+    """
+    # Check if inputs are binary or range-based
+    is_binary = False
+    if isinstance(y_pred, (list, tuple, np.ndarray)) and isinstance(
+        y_pred[0], (int, np.integer)
+    ):
+        is_binary = True
+    elif isinstance(y_real, (list, tuple, np.ndarray)) and isinstance(
+        y_real[0], (int, np.integer)
+    ):
+        is_binary = True
+
+    if is_binary:
+        if not isinstance(y_pred, (list, tuple, np.ndarray)) or not isinstance(
+            y_real, (list, tuple, np.ndarray)
+        ):
+            raise ValueError(
+                "For binary inputs, y_pred and y_real should be list or tuple, "
+                "or numpy array of integers."
+            )
+        if len(y_pred) != len(y_real):
+            raise ValueError(
+                "For binary inputs, y_pred and y_real must be of the same length."
+            )
+
+        y_pred_ranges = _binary_to_ranges(y_pred)
+        y_real_ranges = _binary_to_ranges(y_real)
+    else:
+        y_pred_ranges = y_pred
+        y_real_ranges = y_real
+
+    if gamma not in ["reciprocal", "one"]:
+        raise ValueError("Invalid gamma type for precision. Use 'reciprocal' or 'one'.")
+
+    # Flattening y_pred and y_real to resolve nested lists
+    flat_y_pred = _flatten_ranges(y_pred_ranges)
+    flat_y_real = _flatten_ranges(y_real_ranges)
+
+    total_overlap_reward = 0.0
+    total_cardinality = 0
+
+    for pred_range in flat_y_pred:
+        overlap_set = set()
+        cardinality = 0
+
+        for real_start, real_end in flat_y_real:
+            overlap_start = max(pred_range[0], real_start)
+            overlap_end = min(pred_range[1], real_end)
+
+            if overlap_start <= overlap_end:
+                overlap_set.update(range(overlap_start, overlap_end + 1))
+                cardinality += 1
+
+        overlap_reward = _calculate_overlap_reward_precision(
+            pred_range, overlap_set, bias_type
+        )
+        gamma_value = _gamma_select(cardinality, gamma)
+        total_overlap_reward += gamma_value * overlap_reward
+        total_cardinality += 1
+
+    precision = (
+        total_overlap_reward / total_cardinality if total_cardinality > 0 else 0.0
+    )
+    return precision
+
+
+def ts_recall(y_pred, y_real, gamma="one", bias_type="flat", alpha=0.0):
+    """
+    Calculate Recall for time series anomaly detection.
+
+    Recall measures the proportion of correctly predicted anomaly positions
+    out of all the real/actual (ground truth) anomaly positions, aggregated across the
+    entire time series.
+
+    Parameters
+    ----------
+    y_pred : list of tuples or binary sequence
+        The predicted anomaly ranges.
+        - For range-based input, each tuple represents a range (start, end) of the
+          anomaly where start is starting index (inclusive) and end is ending index
+          (inclusive).
+        - For binary inputs, the sequence should contain integers (0 or 1), where 1
+          indicates an anomaly. In this case, y_pred and y_real must be of same length.
+    y_real : list of tuples, list of lists of tuples or binary sequence
+        The real/actual (ground truth) ranges.
+        - For range-based input, each tuple represents a range (start, end) of the
+          anomaly where start is starting index (inclusive) and end is ending index
+          (inclusive).
+        - If y_real is in the format of list of lists, they will be flattened into a
+          single list of tuples bringing it to the above format.
+        - For binary inputs, the sequence should contain integers (0 or 1), where 1
+          indicates an anomaly. In this case, y_pred and y_real must be of same length.
+    gamma : str, default="one"
+        Cardinality type. Should be one of ["reciprocal", "one", "udf_gamma"].
+    bias_type : str, default="flat"
+        Type of bias to apply. Should be one of ["flat", "front", "middle", "back"].
+    alpha : float, default: 0.0
+        Weight for existence reward in recall calculation.
+
+    Returns
+    -------
+    float
+        Recall
+
+    Raises
+    ------
+    ValueError
+        If input sequence is binary and y_real and y_pred are of different lengths.
+
+    References
+    ----------
+    .. [1] Tatbul, Nesime, Tae Jun Lee, Stan Zdonik, Mejbah Alam,and Justin Gottschlich.
+       "Precision and Recall for Time Series." 32nd Conference on Neural Information
+       Processing Systems (NeurIPS 2018), Montréal, Canada.
+       http://papers.nips.cc/paper/7462-precision-and-recall-for-time-series.pdf
+    """
+    is_binary = False
+    if isinstance(y_pred, (list, tuple, np.ndarray)) and isinstance(
+        y_pred[0], (int, np.integer)
+    ):
+        is_binary = True
+    elif isinstance(y_real, (list, tuple, np.ndarray)) and isinstance(
+        y_real[0], (int, np.integer)
+    ):
+        is_binary = True
+
+    if is_binary:
+        if not isinstance(y_pred, (list, tuple, np.ndarray)) or not isinstance(
+            y_real, (list, tuple, np.ndarray)
+        ):
+            raise ValueError(
+                "For binary inputs, y_pred and y_real should be list or tuple, "
+                "or numpy array of integers."
+            )
+        if len(y_pred) != len(y_real):
+            raise ValueError(
+                "For binary inputs, y_pred and y_real must be of the same length."
+            )
+
+        y_pred_ranges = _binary_to_ranges(y_pred)
+        y_real_ranges = _binary_to_ranges(y_real)
+    else:
+        y_pred_ranges = y_pred
+        y_real_ranges = y_real
+
+    # Flattening y_pred and y_real to resolve nested lists
+    flat_y_pred = _flatten_ranges(y_pred_ranges)
+    flat_y_real = _flatten_ranges(y_real_ranges)
+
+    total_overlap_reward = 0.0
+
+    for real_range in flat_y_real:
+        overlap_set = set()
+        cardinality = 0
+
+        for pred_range in flat_y_pred:
+            overlap_start = max(real_range[0], pred_range[0])
+            overlap_end = min(real_range[1], pred_range[1])
+
+            if overlap_start <= overlap_end:
+                overlap_set.update(range(overlap_start, overlap_end + 1))
+                cardinality += 1
+
+        existence_reward = 1.0 if overlap_set else 0.0
+
+        if overlap_set:
+            overlap_reward = _calculate_overlap_reward_recall(
+                real_range, overlap_set, bias_type
+            )
+            gamma_value = _gamma_select(cardinality, gamma)
+            overlap_reward *= gamma_value
+        else:
+            overlap_reward = 0.0
+
+        recall_score = alpha * existence_reward + (1 - alpha) * overlap_reward
+        total_overlap_reward += recall_score
+
+    recall = total_overlap_reward / len(flat_y_real) if flat_y_real else 0.0
+    return recall
+
+
+def ts_fscore(
+    y_pred,
+    y_real,
+    gamma="one",
+    p_bias="flat",
+    r_bias="flat",
+    p_alpha=0.0,
+    r_alpha=0.0,
+    beta=1.0,
+):
+    """
+    Calculate F1-Score for time series anomaly detection.
+
+    F-1 Score is the harmonic mean of Precision and Recall, providing
+    a single metric to evaluate the performance of an anomaly detection model.
+
+    Parameters
+    ----------
+    y_pred : list of tuples or binary sequence
+        The predicted anomaly ranges.
+        - For range-based input, each tuple represents a range (start, end) of the
+          anomaly where start is starting index (inclusive) and end is ending index
+          (inclusive).
+        - For binary inputs, the sequence should contain integers (0 or 1), where 1
+          indicates an anomaly. In this case, y_pred and y_real must be of same length.
+    y_real : list of tuples, list of lists of tuples or binary sequence
+        The real/actual (ground truth) ranges.
+        - For range-based input, each tuple represents a range (start, end) of the
+          anomaly where start is starting index (inclusive) and end is ending index
+          (inclusive).
+        - If y_real is in the format of list of lists, they will be flattened into a
+          single list of tuples bringing it to the above format.
+        - For binary inputs, the sequence should contain integers (0 or 1), where 1
+          indicates an anomaly. In this case, y_pred and y_real must be of same length.
+    gamma : str, default="one"
+        Cardinality type. Should be one of ["reciprocal", "one", "udf_gamma"].
+    p_bias : str, default="flat"
+        Type of bias to apply for precision.
+        Should be one of ["flat", "front", "middle", "back"].
+    r_bias : str, default="flat"
+        Type of bias to apply for recall.
+        Should be one of ["flat", "front", "middle", "back"].
+    p_alpha : float, default=0.0
+        Weight for existence reward in Precision calculation.
+    r_alpha : float, default=0.0
+        Weight for existence reward in Recall calculation.
+    beta : float, default=1.0
+        F-score beta determines the weight of recall in the combined score.
+        beta < 1 lends more weight to precision, while beta > 1 favors recall.
+
+    Returns
+    -------
+    float
+        F1-Score
+
+    References
+    ----------
+    .. [1] Tatbul, Nesime, Tae Jun Lee, Stan Zdonik, Mejbah Alam,and Justin Gottschlich.
+       "Precision and Recall for Time Series." 32nd Conference on Neural Information
+       Processing Systems (NeurIPS 2018), Montréal, Canada.
+       http://papers.nips.cc/paper/7462-precision-and-recall-for-time-series.pdf
+    """
+    precision = ts_precision(y_pred, y_real, gamma, p_bias)
+    recall = ts_recall(y_pred, y_real, gamma, r_bias, r_alpha)
+
+    if precision + recall > 0:
+        fscore = ((1 + beta**2) * (precision * recall)) / (beta**2 * precision + recall)
+    else:
+        fscore = 0.0
+
+    return fscore
diff --git a/aeon/classification/dictionary_based/_redcomets.py b/aeon/classification/dictionary_based/_redcomets.py
index f2864674..601f573b 100644
--- a/aeon/classification/dictionary_based/_redcomets.py
+++ b/aeon/classification/dictionary_based/_redcomets.py
@@ -255,8 +255,7 @@ class REDCOMETS(BaseClassifier):
         sfa_clfs = []
         for sfa in sfa_transforms:
             sfa.fit(X_smote, y_smote)
-            sfa_dics = sfa.transform_words(X_smote)
-            X_sfa = sfa_dics[:, 0, :]
+            X_sfa = sfa.transform_words(X_smote)[0]
 
             rf = RandomForestClassifier(
                 n_estimators=self.n_trees,
@@ -417,8 +416,7 @@ class REDCOMETS(BaseClassifier):
         pred_mat = np.zeros((X.shape[0], self.n_classes_))
 
         for sfa, (rf, weight) in zip(self.sfa_transforms, self.sfa_clfs):
-            sfa_dics = sfa.transform_words(X)
-            X_sfa = sfa_dics[:, 0, :]
+            X_sfa = sfa.transform_words(X)[0]
 
             rf_pred_mat = rf.predict_proba(X_sfa)
 
diff --git a/aeon/clustering/_kernel_k_means.py b/aeon/clustering/_kernel_k_means.py
index 6aab712d..98cac931 100644
--- a/aeon/clustering/_kernel_k_means.py
+++ b/aeon/clustering/_kernel_k_means.py
@@ -3,9 +3,122 @@
 from typing import Optional, Union
 
 import numpy as np
+from numba import njit
 from numpy.random import RandomState
 
 from aeon.clustering.base import BaseClusterer
+from aeon.distances.pointwise._squared import squared_pairwise_distance
+
+
+@njit(cache=True, fastmath=True)
+def _kdtw_lk(x, y, local_kernel):
+    channels = np.shape(x)[1]
+    padding_vector = np.zeros((1, channels))
+
+    x = np.concatenate((padding_vector, x), axis=0)
+    y = np.concatenate((padding_vector, y), axis=0)
+
+    x_timepoints, _ = np.shape(x)
+    y_timepoints, _ = np.shape(y)
+
+    cost_matrix = np.zeros((x_timepoints, y_timepoints))
+    cumulative_dp_diag = np.zeros((x_timepoints, y_timepoints))
+    diagonal_weights = np.zeros(max(x_timepoints, y_timepoints))
+
+    min_timepoints = min(x_timepoints, y_timepoints)
+    diagonal_weights[1] = 1.0
+    for i in range(1, min_timepoints):
+        diagonal_weights[i] = local_kernel[i - 1, i - 1]
+
+    cost_matrix[0, 0] = 1
+    cumulative_dp_diag[0, 0] = 1
+
+    for i in range(1, x_timepoints):
+        cost_matrix[i, 1] = cost_matrix[i - 1, 1] * local_kernel[i - 1, 2]
+        cumulative_dp_diag[i, 1] = cumulative_dp_diag[i - 1, 1] * diagonal_weights[i]
+
+    for j in range(1, y_timepoints):
+        cost_matrix[1, j] = cost_matrix[1, j - 1] * local_kernel[2, j - 1]
+        cumulative_dp_diag[1, j] = cumulative_dp_diag[1, j - 1] * diagonal_weights[j]
+
+    for i in range(1, x_timepoints):
+        for j in range(1, y_timepoints):
+            local_cost = local_kernel[i - 1, j - 1]
+            cost_matrix[i, j] = (
+                cost_matrix[i - 1, j]
+                + cost_matrix[i, j - 1]
+                + cost_matrix[i - 1, j - 1]
+            ) * local_cost
+            if i == j:
+                cumulative_dp_diag[i, j] = (
+                    cumulative_dp_diag[i - 1, j - 1] * local_cost
+                    + cumulative_dp_diag[i - 1, j] * diagonal_weights[i]
+                    + cumulative_dp_diag[i, j - 1] * diagonal_weights[j]
+                )
+            else:
+                cumulative_dp_diag[i, j] = (
+                    cumulative_dp_diag[i - 1, j] * diagonal_weights[i]
+                    + cumulative_dp_diag[i, j - 1] * diagonal_weights[j]
+                )
+    cost_matrix = cost_matrix + cumulative_dp_diag
+    return cost_matrix[x_timepoints - 1, y_timepoints - 1]
+
+
+def kdtw(x, y, sigma=1.0, epsilon=1e-3):
+    """
+    Callable kernel function for KernelKMeans.
+
+    Parameters
+    ----------
+    X: np.ndarray, of shape (n_timepoints, n_channels)
+            First time series sample.
+    y: np.ndarray, of shape (n_timepoints, n_channels)
+            Second time series sample.
+    sigma : float, default=1.0
+        Parameter controlling the width of the exponential local kernel. Smaller sigma
+        values lead to a sharper decay of similarity with increasing distance.
+    epsilon : float, default=1e-3
+        A small constant added for numerical stability to avoid zero values in the
+        local kernel matrix.
+
+    Returns
+    -------
+    similarity : float
+        A scalar value representing the computed KDTW similarity between the two time
+        series. Higher values indicate greater similarity.
+    """
+    distance = squared_pairwise_distance(x, y)
+    local_kernel = (np.exp(-distance / sigma) + epsilon) / (3 * (1 + epsilon))
+    return _kdtw_lk(x, y, local_kernel)
+
+
+def factory_kdtw_kernel(channels: int):
+    """
+    Return a kdtw kernel callable function that flattened samples to (T, channels).
+
+    Parameters
+    ----------
+    channels: int
+        Number of channels per timepoint.
+
+    Returns
+    -------
+    kdtw_kernel : callable
+        A callable kernel function that computes the KDTW similarity between two
+        time series samples. The function signature is the same as the kdtw
+        function.
+    """
+
+    def kdtw_kernel(x, y, sigma=1.0, epsilon=1e-3):
+        if x.ndim == 1:
+            T = x.size // channels
+            x = x.reshape(T, channels)
+        if y.ndim == 1:
+            T = y.size // channels
+            y = y.reshape(T, channels)
+        return kdtw(x, y, sigma=sigma, epsilon=epsilon)
+
+    return kdtw_kernel
 
 
 class TimeSeriesKernelKMeans(BaseClusterer):
@@ -141,6 +254,9 @@ class TimeSeriesKernelKMeans(BaseClusterer):
         if self.verbose is True:
             verbose = 1
 
+        if self.kernel == "kdtw":
+            self.kernel = factory_kdtw_kernel(channels=X.shape[1])
+
         self._tslearn_kernel_k_means = TsLearnKernelKMeans(
             n_clusters=self.n_clusters,
             kernel=self.kernel,
diff --git a/aeon/distances/__init__.py b/aeon/distances/__init__.py
index e1d3205e..d6ff3f77 100644
--- a/aeon/distances/__init__.py
+++ b/aeon/distances/__init__.py
@@ -18,6 +18,10 @@ __all__ = [
     "dtw_pairwise_distance",
     "dtw_cost_matrix",
     "dtw_alignment_path",
+    "dtw_gi_distance",
+    "dtw_gi_pairwise_distance",
+    "dtw_gi_cost_matrix",
+    "dtw_gi_alignment_path",
     "ddtw_distance",
     "ddtw_pairwise_distance",
     "ddtw_alignment_path",
@@ -111,6 +115,10 @@ from aeon.distances.elastic import (
     dtw_alignment_path,
     dtw_cost_matrix,
     dtw_distance,
+    dtw_gi_alignment_path,
+    dtw_gi_cost_matrix,
+    dtw_gi_distance,
+    dtw_gi_pairwise_distance,
     dtw_pairwise_distance,
     edr_alignment_path,
     edr_cost_matrix,
diff --git a/aeon/distances/_distance.py b/aeon/distances/_distance.py
index 1263e11c..33f91414 100644
--- a/aeon/distances/_distance.py
+++ b/aeon/distances/_distance.py
@@ -24,6 +24,10 @@ from aeon.distances.elastic import (
     dtw_alignment_path,
     dtw_cost_matrix,
     dtw_distance,
+    dtw_gi_alignment_path,
+    dtw_gi_cost_matrix,
+    dtw_gi_distance,
+    dtw_gi_pairwise_distance,
     dtw_pairwise_distance,
     edr_alignment_path,
     edr_cost_matrix,
@@ -447,6 +451,7 @@ def get_distance_function(method: Union[str, DistanceFunction]) -> DistanceFunct
     method          Distance Function
     =============== ========================================
     'dtw'           distances.dtw_distance
+    'dtw_gi'        distances.dtw_gi_distance
     'shape_dtw'     distances.shape_dtw_distance
     'ddtw'          distances.ddtw_distance
     'wdtw'          distances.wdtw_distance
@@ -728,6 +733,16 @@ DISTANCES = [
         "symmetric": True,
         "unequal_support": True,
     },
+    {
+        "name": "dtw_gi",
+        "distance": dtw_gi_distance,
+        "pairwise_distance": dtw_gi_pairwise_distance,
+        "cost_matrix": dtw_gi_cost_matrix,
+        "alignment_path": dtw_gi_alignment_path,
+        "type": DistanceType.ELASTIC,
+        "symmetric": False,
+        "unequal_support": True,
+    },
     {
         "name": "ddtw",
         "distance": ddtw_distance,
diff --git a/aeon/distances/elastic/__init__.py b/aeon/distances/elastic/__init__.py
index 0c386245..8e5d1aa9 100644
--- a/aeon/distances/elastic/__init__.py
+++ b/aeon/distances/elastic/__init__.py
@@ -10,6 +10,10 @@ __all__ = [
     "dtw_pairwise_distance",
     "dtw_cost_matrix",
     "dtw_alignment_path",
+    "dtw_gi_distance",
+    "dtw_gi_pairwise_distance",
+    "dtw_gi_cost_matrix",
+    "dtw_gi_alignment_path",
     "ddtw_distance",
     "ddtw_pairwise_distance",
     "ddtw_alignment_path",
@@ -71,6 +75,12 @@ from aeon.distances.elastic._dtw import (
     dtw_distance,
     dtw_pairwise_distance,
 )
+from aeon.distances.elastic._dtw_gi import (
+    dtw_gi_alignment_path,
+    dtw_gi_cost_matrix,
+    dtw_gi_distance,
+    dtw_gi_pairwise_distance,
+)
 from aeon.distances.elastic._edr import (
     edr_alignment_path,
     edr_cost_matrix,
diff --git a/aeon/distances/elastic/_bounding_matrix.py b/aeon/distances/elastic/_bounding_matrix.py
index 2b710b97..509dff6e 100644
--- a/aeon/distances/elastic/_bounding_matrix.py
+++ b/aeon/distances/elastic/_bounding_matrix.py
@@ -63,44 +63,61 @@ def create_bounding_matrix(
 def _itakura_parallelogram(x_size: int, y_size: int, max_slope_percent: float):
     """Itakura parallelogram bounding matrix.
 
-    This code was adapted from tslearn. This link to the original code line 974:
-    https://github.com/tslearn-team/tslearn/blob/main/tslearn/metrics/dtw_variants.py
+    This code was adapted from pyts. This link to the original code:
+    https://pyts.readthedocs.io/en/latest/_modules/pyts/metrics/dtw.html#itakura_parallelogram
     """
-    if x_size != y_size:
-        raise ValueError(
-            """Itakura parallelogram does not support unequal length time series.
-Please consider using a full bounding matrix or a sakoe chiba bounding matrix
-instead."""
-        )
     one_percent = min(x_size, y_size) / 100
     max_slope = math.floor((max_slope_percent * one_percent) * 100)
     min_slope = 1 / float(max_slope)
-    max_slope *= float(x_size) / float(y_size)
-    min_slope *= float(x_size) / float(y_size)
-
-    lower_bound = np.empty((2, y_size))
-    lower_bound[0] = min_slope * np.arange(y_size)
-    lower_bound[1] = (
-        (x_size - 1) - max_slope * (y_size - 1) + max_slope * np.arange(y_size)
-    )
-    lower_bound_ = np.empty(y_size)
-    for i in range(y_size):
-        lower_bound_[i] = max(round(lower_bound[0, i], 2), round(lower_bound[1, i], 2))
-    lower_bound_ = np.ceil(lower_bound_)
-
-    upper_bound = np.empty((2, y_size))
-    upper_bound[0] = max_slope * np.arange(y_size)
-    upper_bound[1] = (
-        (x_size - 1) - min_slope * (y_size - 1) + min_slope * np.arange(y_size)
-    )
-    upper_bound_ = np.empty(y_size)
-    for i in range(y_size):
-        upper_bound_[i] = min(round(upper_bound[0, i], 2), round(upper_bound[1, i], 2))
-    upper_bound_ = np.floor(upper_bound_ + 1)
-
-    bounding_matrix = np.full((x_size, y_size), False)
-    for i in range(y_size):
-        bounding_matrix[int(lower_bound_[i]) : int(upper_bound_[i]), i] = True
+    max_slope *= float(y_size - 1) / float(x_size - 2)
+    max_slope = max(max_slope, 1.0)
+
+    min_slope *= float(y_size - 2) / float(x_size - 1)
+    min_slope = min(min_slope, 1.0)
+
+    centered_scale = np.arange(x_size) - x_size + 1
+
+    lower_bound = np.empty(x_size, dtype=np.float64)
+    upper_bound = np.empty(x_size, dtype=np.float64)
+
+    for i in range(x_size):
+        lb0 = min_slope * i
+        lb1 = max_slope * centered_scale[i] + y_size - 1
+        lower_bound[i] = math.ceil(max(round(lb0, 2), round(lb1, 2)))
+
+        ub0 = max_slope * i + 1
+        ub1 = min_slope * centered_scale[i] + y_size
+        upper_bound[i] = math.floor(min(round(ub0, 2), round(ub1, 2)))
+
+    if max_slope == 1.0:
+        if y_size > x_size:
+            for i in range(x_size - 1):
+                upper_bound[i] = lower_bound[i + 1]
+        else:
+            for i in range(x_size):
+                upper_bound[i] = lower_bound[i] + 1
+
+    for i in range(x_size):
+        if lower_bound[i] < 0:
+            lower_bound[i] = 0
+        if lower_bound[i] > y_size:
+            lower_bound[i] = y_size
+        if upper_bound[i] < 0:
+            upper_bound[i] = 0
+        if upper_bound[i] > y_size:
+            upper_bound[i] = y_size
+
+    bounding_matrix = np.empty((x_size, y_size), dtype=np.bool_)
+    for i in range(x_size):
+        for j in range(y_size):
+            bounding_matrix[i, j] = False
+
+    for i in range(x_size):
+        start = int(lower_bound[i])
+        end = int(upper_bound[i])
+        for j in range(start, end):
+            bounding_matrix[i, j] = True
+
     return bounding_matrix
 
 
diff --git a/aeon/distances/elastic/_dtw_gi.py b/aeon/distances/elastic/_dtw_gi.py
new file mode 100644
index 00000000..0026c377
--- /dev/null
+++ b/aeon/distances/elastic/_dtw_gi.py
@@ -0,0 +1,551 @@
+r"""Dynamic time warping with Global Invariances (DTW-GI) between two time series."""
+
+__maintainer__ = []
+
+from typing import Optional, Union
+
+import numpy as np
+from numba import njit
+from numba.typed import List as NumbaList
+
+from aeon.distances.elastic._dtw import dtw_alignment_path, dtw_cost_matrix
+from aeon.utils.conversion._convert_collection import _convert_collection_to_numba_list
+from aeon.utils.validation.collection import _is_numpy_list_multivariate
+
+
+@njit(cache=True, fastmath=True)
+def _path2mat(
+    path: list[tuple[int, int]],
+    x_timepoints: int,
+    y_timepoints: int,
+):
+    r"""Convert a warping alignment path to a binary warping matrix."""
+    w = np.zeros((x_timepoints, y_timepoints))
+    for i, j in path:
+        w[i, j] = 1
+    return w
+
+
+@njit(cache=True, fastmath=True)
+def _dtw_gi(
+    x: np.ndarray,
+    y: np.ndarray,
+    window: Optional[float] = None,
+    itakura_max_slope: Optional[float] = None,
+    init_p: np.ndarray = None,
+    max_iter: int = 20,
+    use_bias: bool = False,
+):
+    r"""
+    Compute Dynamic Time Warping with Global Invariance between the two time series.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        First time series, either univariate, shape ``(n_timepoints,)``, or
+        multivariate, shape ``(n_channels, n_timepoints)``.
+    y : np.ndarray
+        Second time series, either univariate, shape ``(n_timepoints,)``, or
+        multivariate, shape ``(n_channels, n_timepoints)``.
+    window : float or None, default=None
+        The window to use for the bounding matrix. If None, no bounding matrix
+        is used. window is a percentage deviation, so if ``window = 0.1`` then
+        10% of the series length is the max warping allowed.
+        is used.
+    itakura_max_slope : float, default=None
+        Maximum slope as a proportion of the number of time points used to create
+        Itakura parallelogram on the bounding matrix. Must be between 0. and 1.
+    init_p : array-like of shape (x_channels, y_channels), default=None
+        Initial linear transformation. If None, the identity matrix is used.
+    max_iter : int, default=20
+        Maximum number of iterations for the iterative optimization.
+    use_bias : bool, default=False
+        If True, the feature space map is affine (with a bias term).
+
+    Returns
+    -------
+    - w_pi: binary warping matrix of shape (n0, n1)
+    - p: the final linear (Stiefel) matrix of shape (x_channels, y_channels)
+    - cost: final DTW cost considering global invariances
+
+    If use_bias is True, also returns:
+      - bias
+
+    """
+    if x.ndim == 1 and y.ndim == 1:
+        x_ = x.reshape((1, x.shape[0]))
+        y_ = y.reshape((1, y.shape[0]))
+    elif x.ndim == 2 and y.ndim == 2:
+        x_ = x
+        y_ = y
+    else:
+        raise ValueError("x and y must be 1D or 2D")
+
+    x_channels, x_timepoints = x_.shape
+    y_channels, y_timepoints = y_.shape
+
+    x_m = np.sum(x_, axis=1) / x_.shape[1]
+    x_m = x_m.reshape((-1, 1))
+    y_m = np.sum(y_, axis=1) / y_.shape[1]
+    y_m = y_m.reshape((-1, 1))
+
+    w_pi = np.zeros((x_timepoints, y_timepoints))
+    if init_p is None:
+        p = np.eye(x_channels, y_channels, dtype=np.float64)
+    else:
+        p = init_p
+
+    y_ = y_.astype(np.float64)
+    x_ = x_.astype(np.float64)
+
+    bias = np.zeros((x_channels, 1))
+
+    for _ in range(max_iter):
+        w_pi_old = w_pi.copy()
+        y_transformed = p.dot(y_) + bias
+
+        path, cost = dtw_alignment_path(x_, y_transformed, window, itakura_max_slope)
+        w_pi = _path2mat(path, x_timepoints, y_timepoints)
+
+        if np.allclose(w_pi, w_pi_old):
+            break
+
+        if use_bias:
+            m = (x_ - x_m).dot(w_pi).dot((y_ - y_m).T)
+        else:
+            m = x_.dot(w_pi).dot(y_.T)
+
+        u, sigma, vt = np.linalg.svd(m, full_matrices=False)
+        p = u.dot(vt)
+        if use_bias:
+            bias = x_m - p.dot(y_m)
+
+    y_trans = p.dot(y_) + bias
+    path, cost = dtw_alignment_path(x_, y_trans, window, itakura_max_slope)
+
+    if use_bias:
+        return w_pi, p, bias, cost, x_, y_trans
+    else:
+        dummy_bias = np.zeros((x_channels, 1), dtype=np.float64)
+        return w_pi, p, dummy_bias, cost, x_, y_trans
+
+
+@njit(cache=True, fastmath=True)
+def dtw_gi_distance(
+    x: np.ndarray,
+    y: np.ndarray,
+    window: Optional[float] = None,
+    itakura_max_slope: Optional[float] = None,
+    init_p: np.ndarray = None,
+    max_iter: int = 20,
+    use_bias: bool = False,
+) -> float:
+    r"""Compute the DTW_GI distance between two time series x and y.
+
+    The DTW_gi distance between 2 timeseries x and y is the similarity
+    measure that estimates both temporal alignment and does feature space
+    transformation between time series simultaneously. This means that the
+    time series do not have to lie in the same ambient space.
+    A good background into DTW with global invariances can be found in [1]_.
+    This implementation is inspired by [2]_.
+
+    For example, if we have two time series x and y of lengths n and m
+    respectively, and we assume that the time series do not lie in the
+    same ambient space. Lets assume that features of x lie in :math:`\mathbb{R}^p`
+    and features of y lie in :math:`\mathbb{R}^q`. To compare the two time series,
+    we need to find an optimum mapping from the feature space of y to the feature space
+    where features of x lie. So think of it as optimizing on a family of functions F
+    that map features from y onto the feature space in which features of x
+    lie. (This is just one way to do it, the mapping can be
+    in the opposite direction as well. But this function assumes the former way).
+
+    More formally, we define Dynamic Time Warping with Global Invariances (DTW-GI)
+    as the solution of the following joint optimization problem:
+
+    :math:`\text{DTW-GI}(\mathbf{x}, \mathbf{y}) =
+    \min_{f \in \mathcal{F}, \pi \in \mathcal{A}(\mathbf{x}, \mathbf{y})}
+    \sqrt{\sum_{(i,j) \in \pi} d(x_i, f(y_j))^2},`
+
+    This similarity measure estimates temporal alignment
+    with feature space transformation between time series
+    allowing the alignment of time series that do not exist in the
+    same ambient space.
+
+    Note: The optimal warping path :math:`P^*` for a given cost matrix can be found
+    exactly through a dynamic programming formulation. This can be a time consuming
+    operation, and it is common to put a restriction on the amount of warping allowed.
+    This is implemented through the bounding_matrix structure, that supplies a mask for
+    allowable warpings. The most common bounding strategies include the
+    Sakoe-Chiba band [3]_. The width of the allowed warping is controlled through the
+     ``window`` parameter which sets the maximum proportion of warping allowed.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        First time series, either univariate, shape ``(n_timepoints,)``, or
+        multivariate, shape ``(n_channels, n_timepoints)``.
+    y : np.ndarray
+        Second time series, either univariate, shape ``(n_timepoints,)``, or
+        multivariate, shape ``(n_channels, n_timepoints)``.
+    window : float, default=None
+        The window to use for the bounding matrix. If None, no bounding matrix
+        is used. window is a percentage deviation, so if ``window = 0.1``,
+        10% of the series length is the max warping allowed.
+        is used.
+    itakura_max_slope : float, default=None
+        Maximum slope as a proportion of the number of time points used to create
+        Itakura parallelogram on the bounding matrix. Must be between 0. and 1.
+    init_p : array-like of shape (x_channels, y_channels), default=None
+        Initial linear transformation. If None, the identity matrix is used.
+    max_iter : int, default=20
+        Maximum number of iterations for the iterative optimization.
+    use_bias : bool, default=False
+        If True, the feature space map is affine (with a bias term).
+
+    Returns
+    -------
+    float
+        DTW_GI distance between x and y, minimum value 0.
+
+    Raises
+    ------
+    ValueError
+        If x and y are not 1D or 2D arrays.
+
+    References
+    ----------
+    .. [1] T. Vayer, R. Tavenard, L. Chapel, N. Courty, R. Flamary, and Y. Soullard,
+    “Time Series Alignment with Global Invariances,” arXiv.org, 2020.
+    https://arxiv.org/abs/2002.03848
+
+    .. [2] Romain Tavenard, “DTW with Global Invariances,” Github.io, Dec. 17, 2020.
+    https://rtavenar.github.io/hdr/parts/01/dtw/dtw_gi.html
+
+    .. [3] Sakoe H. and Chiba S.: Dynamic programming algorithm optimization for
+    spoken word recognition. IEEE Transactions on Acoustics, Speech, and Signal
+    Processing 26(1):43-49, 1978.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aeon.distances import dtw_gi_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])
+    >>> dtw_gi_distance(x, y) # 1D series
+    768.0
+    >>> x = np.array([[1, 2, 3, 4, 5], [6, 7, 8, 9, 10], [0, 1, 0, 2, 0]])
+    >>> y = np.array([[11, 12, 13, 14],[7, 8, 9, 20],[1, 3, 4, 5]] )
+    >>> round(dtw_gi_distance(x, y), 1) # 2D series with 3 channels, unequal length
+    359.2
+    """
+    if x.ndim == 1 and y.ndim == 1:
+        _x = x.reshape((1, x.shape[0]))
+        _y = y.reshape((1, y.shape[0]))
+        return _dtw_gi(_x, _y, window, itakura_max_slope, init_p, max_iter, use_bias)[3]
+    if x.ndim == 2 and y.ndim == 2:
+        return _dtw_gi(x, y, window, itakura_max_slope, init_p, max_iter, use_bias)[3]
+    raise ValueError("x and y must be 1D or 2D")
+
+
+@njit(cache=True, fastmath=True)
+def dtw_gi_cost_matrix(
+    x: np.ndarray,
+    y: np.ndarray,
+    window: Optional[float] = None,
+    itakura_max_slope: Optional[float] = None,
+    init_p: np.ndarray = None,
+    max_iter: int = 20,
+    use_bias: bool = False,
+) -> np.ndarray:
+    r"""Compute the DTW_GI cost matrix between two time series.
+
+    The cost matrix is the pairwise Euclidean distance between all points
+    :math:`M_{i,j}=(x_i-y_{\text{trans},j})^2`. Where `y_trans` is the time
+    series we get by finding the optimal mapping from feature space of y
+    to feature space where features of x lie. It is used in the DTW GI
+    path calculations.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        First time series, either univariate, shape ``(n_timepoints,)``, or
+        multivariate, shape ``(n_channels, n_timepoints)``.
+    y : np.ndarray
+        Second time series, either univariate, shape ``(n_timepoints,)``, or
+        multivariate, shape ``(n_channels, n_timepoints)``.
+    window : float, default=None
+        The window to use for the bounding matrix. If None, no bounding matrix
+        is used. window is a percentage deviation, so if ``window = 0.1``,
+        10% of the series length is the max warping allowed.
+        is used.
+    itakura_max_slope : float, default=None
+        Maximum slope as a proportion of the number of time points used to create
+        Itakura parallelogram on the bounding matrix. Must be between 0. and 1.
+    init_p : array-like of shape (x_channels, y_channels), default=None
+        Initial linear transformation. If None, the identity matrix is used.
+    max_iter : int, default=20
+        Maximum number of iterations for the iterative optimization.
+    use_bias : bool, default=False
+        If True, the feature space map is affine (with a bias term).
+
+    Returns
+    -------
+    np.ndarray (n_timepoints, m_timepoints)
+        dtw gi cost matrix between x and y.
+
+    Raises
+    ------
+    ValueError
+        If x and y are not 1D or 2D arrays.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aeon.distances import dtw_gi_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]])
+    >>> dtw_gi_cost_matrix(x, y)
+    array([[  0.,   1.,   5.,  14.,  30.,  55.,  91., 140., 204., 285.],
+           [  1.,   0.,   1.,   5.,  14.,  30.,  55.,  91., 140., 204.],
+           [  5.,   1.,   0.,   1.,   5.,  14.,  30.,  55.,  91., 140.],
+           [ 14.,   5.,   1.,   0.,   1.,   5.,  14.,  30.,  55.,  91.],
+           [ 30.,  14.,   5.,   1.,   0.,   1.,   5.,  14.,  30.,  55.],
+           [ 55.,  30.,  14.,   5.,   1.,   0.,   1.,   5.,  14.,  30.],
+           [ 91.,  55.,  30.,  14.,   5.,   1.,   0.,   1.,   5.,  14.],
+           [140.,  91.,  55.,  30.,  14.,   5.,   1.,   0.,   1.,   5.],
+           [204., 140.,  91.,  55.,  30.,  14.,   5.,   1.,   0.,   1.],
+           [285., 204., 140.,  91.,  55.,  30.,  14.,   5.,   1.,   0.]])
+    """
+    _, _, _, _, xnew, y_trans = _dtw_gi(
+        x, y, window, itakura_max_slope, init_p, max_iter, use_bias
+    )
+
+    return dtw_cost_matrix(xnew, y_trans, window, itakura_max_slope)
+
+
+def dtw_gi_pairwise_distance(
+    X: Union[np.ndarray, list[np.ndarray]],
+    y: Optional[Union[np.ndarray, list[np.ndarray]]] = None,
+    window: Optional[float] = None,
+    itakura_max_slope: Optional[float] = None,
+    unequal_length: bool = None,
+    init_p: np.ndarray = None,
+    max_iter: int = 20,
+    use_bias: bool = False,
+) -> np.ndarray:
+    r"""Compute the DTW_GI pairwise distance between a set of time series.
+
+    By default, this takes a collection of :math:`n` time series :math:`X` and returns a
+    matrix
+    :math:`D` where :math:`D_{i,j}` is the DTW_GI distance between the :math:`i^{th}`
+    and the :math:`j^{th}` series in :math:`X`. If :math:`X` is 2 dimensional,
+    it is assumed to be a collection of univariate series with shape ``(n_cases,
+    n_timepoints)``. If it is 3 dimensional, it is assumed to be shape ``(n_cases,
+    n_channels, n_timepoints)``.
+
+    This function has an optional argument, :math:`y`, to allow calculation of the
+    distance matrix between :math:`X` and one or more series stored in :math:`y`. If
+    :math:`y` is 1 dimensional, we assume it is a single univariate series and the
+    distance matrix returned is shape ``(n_cases,1)``. If it is 2D, we assume it
+    is a collection of univariate series with shape ``(m_cases, m_timepoints)``
+    and the distance ``(n_cases,m_cases)``. If it is 3 dimensional,
+    it is assumed to be shape ``(m_cases, m_channels, m_timepoints)``.
+
+    Parameters
+    ----------
+    X : np.ndarray or List of np.ndarray
+        A collection of time series instances  of shape ``(n_cases, n_timepoints)``
+        or ``(n_cases, n_channels, n_timepoints)``.
+    y : np.ndarray or List of np.ndarray or None, default=None
+        A single series or a collection of time series of shape ``(m_timepoints,)`` or
+        ``(m_cases, m_timepoints)`` or ``(m_cases, m_channels, m_timepoints)``.
+        If None, then the dtw gi pairwise distance between the instances of X is
+        calculated.
+    window : float or None, default=None
+        The window to use for the bounding matrix. If None, no bounding matrix
+        is used.
+    itakura_max_slope : float, default=None
+        Maximum slope as a proportion of the number of time points used to create
+        Itakura parallelogram on the bounding matrix. Must be between 0. and 1.
+    init_p : array-like of shape (x_channels, y_channels), default=None
+        Initial linear transformation. If None, the identity matrix is used.
+    max_iter : int, default=20
+        Maximum number of iterations for the iterative optimization.
+    use_bias : bool, default=False
+        If True, the feature space map is affine (with a bias term).
+
+    Returns
+    -------
+    np.ndarray
+        DTW_GI pairwise matrix between the instances of X of shape
+        ``(n_cases, n_cases)`` or between X and y of shape ``(n_cases,
+        n_cases)``.
+
+    Raises
+    ------
+    ValueError
+        If X is not 2D or 3D array and if y is not 1D, 2D or 3D arrays when passing y.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aeon.distances import dtw_gi_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]]])
+    >>> dtw_gi_pairwise_distance(X)
+    array([[  0.,  26., 108.],
+           [ 26.,   0.,  26.],
+           [108.,  26.,   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]]])
+    >>> dtw_gi_pairwise_distance(X, y)
+    array([[300., 507., 768.],
+           [147., 300., 507.],
+           [ 48., 147., 300.]])
+
+    >>> X = np.array([[[1, 2, 3]],[[4, 5, 6]], [[7, 8, 9]]])
+    >>> y_univariate = np.array([11, 12, 13])
+    >>> dtw_gi_pairwise_distance(X, y_univariate)
+    array([[300.],
+           [147.],
+           [ 48.]])
+
+    >>> # Distance between each TS in a collection of unequal-length time series
+    >>> X = [np.array([1, 2, 3]), np.array([4, 5, 6, 7]), np.array([8, 9, 10, 11, 12])]
+    >>> dtw_gi_pairwise_distance(X)
+    array([[  0.,  42., 292.],
+           [ 42.,   0.,  83.],
+           [292.,  83.,   0.]])
+    """
+    multivariate_conversion = _is_numpy_list_multivariate(X, y)
+    _X, unequal_length = _convert_collection_to_numba_list(
+        X, "X", multivariate_conversion
+    )
+
+    if y is None:
+        # To self
+        return _dtw_gi_pairwise_distance(
+            _X, window, itakura_max_slope, unequal_length, init_p, max_iter, use_bias
+        )
+    _y, unequal_length = _convert_collection_to_numba_list(
+        y, "y", multivariate_conversion
+    )
+    return _dtw_gi_from_multiple_to_multiple_distance(
+        _X, _y, window, itakura_max_slope, unequal_length, init_p, max_iter, use_bias
+    )
+
+
+@njit(cache=True, fastmath=True)
+def _dtw_gi_from_multiple_to_multiple_distance(
+    x: NumbaList[np.ndarray],
+    y: NumbaList[np.ndarray],
+    window: Optional[float] = None,
+    itakura_max_slope: Optional[float] = None,
+    unequal_length: bool = None,
+    init_p: np.ndarray = None,
+    max_iter: int = 20,
+    use_bias: bool = False,
+) -> np.ndarray:
+    n_cases = len(x)
+    m_cases = len(y)
+    distances = np.zeros((n_cases, m_cases))
+
+    for i in range(n_cases):
+        for j in range(m_cases):
+            x1, y1 = x[i], y[j]
+            distances[i, j] = dtw_gi_distance(
+                x1, y1, window, itakura_max_slope, init_p, max_iter, use_bias
+            )
+    return distances
+
+
+@njit(cache=True, fastmath=True)
+def _dtw_gi_pairwise_distance(
+    X: NumbaList[np.ndarray],
+    window: Optional[float] = None,
+    itakura_max_slope: Optional[float] = None,
+    unequal_length: bool = None,
+    init_p: np.ndarray = None,
+    max_iter: int = 20,
+    use_bias: bool = False,
+) -> np.ndarray:
+    n_cases = len(X)
+    distances = np.zeros((n_cases, n_cases))
+
+    for i in range(n_cases):
+        for j in range(i + 1, n_cases):
+            x1, x2 = X[i], X[j]
+            distances[i, j] = dtw_gi_distance(
+                x1, x2, window, itakura_max_slope, init_p, max_iter, use_bias
+            )
+            distances[j, i] = distances[i, j]
+
+    return distances
+
+
+@njit(cache=True, fastmath=True)
+def dtw_gi_alignment_path(
+    x: np.ndarray,
+    y: np.ndarray,
+    window: Optional[float] = None,
+    itakura_max_slope: Optional[float] = None,
+    init_p: np.ndarray = None,
+    max_iter: int = 20,
+    use_bias: bool = False,
+) -> tuple[list[tuple[int, int]], float]:
+    """Compute the DTW_GI alignment path between two time series.
+
+    Parameters
+    ----------
+    x : np.ndarray
+        First time series, shape ``(n_channels, n_timepoints)`` or ``(n_timepoints,)``.
+    y : np.ndarray
+        Second time series, shape ``(m_channels, m_timepoints)`` or ``(m_timepoints,)``.
+    window : float, default=None
+        The window to use for the bounding matrix. If None, no bounding matrix
+        is used.
+    itakura_max_slope : float, default=None
+        Maximum slope as a proportion of the number of time points used to create
+        Itakura parallelogram on the bounding matrix. Must be between 0. and 1.
+    init_p : array-like of shape (x_channels, y_channels), default=None
+        Initial linear transformation. If None, the identity matrix is used.
+    max_iter : int, default=20
+        Maximum number of iterations for the iterative optimization.
+    use_bias : bool, default=False
+        If True, the feature space map is affine (with a bias term).
+
+    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 DTW_GI distance betweeen 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 dtw_gi_alignment_path
+    >>> x = np.array([[1, 2, 3, 6]])
+    >>> y = np.array([[1, 2, 3, 4]])
+    >>> dtw_gi_alignment_path(x, y)
+    ([(0, 0), (1, 1), (2, 2), (3, 3)], 4.0)
+    """
+    w_pi, _, _, cost, _, _ = _dtw_gi(
+        x, y, window, itakura_max_slope, init_p, max_iter, use_bias
+    )
+    min_alignment_path = []
+    for i in range(len(w_pi)):
+        for j in range(len(w_pi[0])):
+            if w_pi[i, j] == 1:
+                min_alignment_path.append((i, j))
+
+    return min_alignment_path, cost
diff --git a/aeon/transformations/collection/base.py b/aeon/transformations/collection/base.py
index 013001d8..a3664eb4 100644
--- a/aeon/transformations/collection/base.py
+++ b/aeon/transformations/collection/base.py
@@ -64,8 +64,8 @@ class BaseCollectionTransformer(BaseCollectionEstimator, BaseTransformer):
         X : np.ndarray or list
             Data to fit transform to, of valid collection type. Input data,
             any number of channels, equal length series of shape ``(
-            n_cases, n_channels, n_timepoints)`` or list of numpy arrays (any number
-            of channels, unequal length series) of shape ``[n_cases]``, 2D np.array
+            n_cases, n_channels, n_timepoints)`` or list of numpy arrays (number
+            of channels, series length) of shape ``[n_cases]``, 2D np.array
             ``(n_channels, n_timepoints_i)``, where ``n_timepoints_i`` is length of
             series ``i``. Other types are allowed and converted into one of the above.
 
@@ -118,8 +118,8 @@ class BaseCollectionTransformer(BaseCollectionEstimator, BaseTransformer):
         X : np.ndarray or list
             Data to fit transform to, of valid collection type. Input data,
             any number of channels, equal length series of shape ``(
-            n_cases, n_channels, n_timepoints)`` or list of numpy arrays (any number
-            of channels, unequal length series) of shape ``[n_cases]``, 2D np.array
+            n_cases, n_channels, n_timepoints)`` or list of numpy arrays (number
+            of channels, series length) of shape ``[n_cases]``, 2D np.array
             ``(n_channels, n_timepoints_i)``, where ``n_timepoints_i`` is length of
             series ``i``. Other types are allowed and converted into one of the above.
 
@@ -171,10 +171,10 @@ class BaseCollectionTransformer(BaseCollectionEstimator, BaseTransformer):
         ----------
         X : np.ndarray or list
             Data to fit transform to, of valid collection type. Input data,
-            any number of channels, equal length series of shape ``(n_cases,
-            n_channels, n_timepoints)`` or list of numpy arrays (any number of
-            channels, unequal length series) of shape ``[n_cases]``, 2D np.array ``(
-            n_channels, n_timepoints_i)``, where ``n_timepoints_i`` is length of
+            any number of channels, equal length series of shape ``(
+            n_cases, n_channels, n_timepoints)`` or list of numpy arrays (number
+            of channels, series length) of shape ``[n_cases]``, 2D np.array
+            ``(n_channels, n_timepoints_i)``, where ``n_timepoints_i`` is length of
             series ``i``. Other types are allowed and converted into one of the above.
 
             Different estimators have different capabilities to handle different
@@ -222,8 +222,8 @@ class BaseCollectionTransformer(BaseCollectionEstimator, BaseTransformer):
         X : np.ndarray or list
             Data to fit transform to, of valid collection type. Input data,
             any number of channels, equal length series of shape ``(
-            n_cases, n_channels, n_timepoints)`` or list of numpy arrays (any number
-            of channels, unequal length series)  of shape ``[n_cases]``, 2D np.array
+            n_cases, n_channels, n_timepoints)`` or list of numpy arrays (number
+            of channels, series length) of shape ``[n_cases]``, 2D np.array
             ``(n_channels, n_timepoints_i)``, where ``n_timepoints_i`` is length of
             series ``i``. Other types are allowed and converted into one of the above.
 
diff --git a/pyproject.toml b/pyproject.toml
index 44f35416..1882fdf5 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -47,7 +47,7 @@ requires-python = ">=3.9,<3.13"
 dependencies = [
     "deprecated>=1.2.13",
     "numba>=0.55,<0.61.0",
-    "numpy>=1.21.0,<2.1.0",
+    "numpy>=1.21.0,<2.3.0",
     "packaging>=20.0",
     "pandas>=2.0.0,<2.3.0",
     "scikit-learn>=1.0.0,<1.7.0",


Test patch:
diff --git a/aeon/benchmarking/metrics/anomaly_detection/tests/test_metrics.py b/aeon/benchmarking/metrics/anomaly_detection/tests/test_metrics.py
new file mode 100644
index 00000000..61db5ded
--- /dev/null
+++ b/aeon/benchmarking/metrics/anomaly_detection/tests/test_metrics.py
@@ -0,0 +1,537 @@
+"""Test cases for the range-based anomaly detection metrics."""
+
+import numpy as np
+import pytest
+
+from aeon.benchmarking.metrics.anomaly_detection import (
+    range_f_score,
+    range_precision,
+    range_recall,
+)
+from aeon.benchmarking.metrics.anomaly_detection.range_metrics import (
+    _binary_to_ranges,
+    ts_fscore,
+    ts_precision,
+    ts_recall,
+)
+from aeon.utils.validation._dependencies import _check_soft_dependencies
+
+
+@pytest.mark.skipif(
+    not _check_soft_dependencies("prts", severity="none"),
+    reason="required soft dependency prts not available",
+)
+def test_single_overlapping_range():
+    """Test for single overlapping range."""
+    y_pred_bin = np.array([0, 1, 1, 1, 1, 0, 0])
+    y_real_bin = np.array([0, 0, 1, 1, 1, 1, 1])
+
+    precision = ts_precision(y_pred_bin, y_real_bin, gamma="one", bias_type="flat")
+    recall = ts_recall(y_pred_bin, y_real_bin, gamma="one", bias_type="flat", alpha=0.0)
+    f1_score = ts_fscore(
+        y_pred_bin,
+        y_real_bin,
+        gamma="one",
+        p_bias="flat",
+        r_bias="flat",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    rb_prec = range_precision(y_real_bin, y_pred_bin, cardinality="one", bias="flat")
+    rb_rec = range_recall(y_real_bin, y_pred_bin, cardinality="one", bias="flat")
+    rb_fsc = range_f_score(
+        y_real_bin,
+        y_pred_bin,
+        beta=1,
+        cardinality="one",
+        p_bias="flat",
+        r_bias="flat",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    np.testing.assert_almost_equal(
+        precision,
+        rb_prec,
+        decimal=6,
+        err_msg=f"Precision mismatch: ts_precision={precision} vs prts={rb_prec}",
+    )
+    np.testing.assert_almost_equal(
+        recall,
+        rb_rec,
+        decimal=6,
+        err_msg=f"Recall mismatch: ts_recall={recall} vs prts={rb_rec}",
+    )
+    np.testing.assert_almost_equal(
+        f1_score,
+        rb_fsc,
+        decimal=6,
+        err_msg=f"F1-Score mismatch: ts_fscore={f1_score} vs prts={rb_fsc}",
+    )
+
+
+@pytest.mark.skipif(
+    not _check_soft_dependencies("prts", severity="none"),
+    reason="required soft dependency prts not available",
+)
+def test_multiple_non_overlapping_ranges():
+    """Test for multiple non-overlapping ranges."""
+    y_pred_bin = np.array([0, 1, 1, 0, 0, 0, 0, 1, 1, 0, 0])
+    y_real_bin = np.array([0, 0, 0, 1, 1, 0, 0, 0, 0, 1, 1])
+
+    precision = ts_precision(y_pred_bin, y_real_bin, gamma="one", bias_type="flat")
+    recall = ts_recall(y_pred_bin, y_real_bin, gamma="one", bias_type="flat", alpha=0.0)
+    f1_score = ts_fscore(
+        y_pred_bin,
+        y_real_bin,
+        gamma="one",
+        beta=1,
+        p_bias="flat",
+        r_bias="flat",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    rb_prec = range_precision(y_real_bin, y_pred_bin, cardinality="one", bias="flat")
+    rb_rec = range_recall(y_real_bin, y_pred_bin, cardinality="one", bias="flat")
+    rb_fsc = range_f_score(
+        y_real_bin,
+        y_pred_bin,
+        beta=1,
+        cardinality="one",
+        p_bias="flat",
+        r_bias="flat",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    np.testing.assert_almost_equal(
+        precision,
+        rb_prec,
+        decimal=6,
+        err_msg=f"Precision mismatch: ts_precision={precision} vs prts={rb_prec}",
+    )
+    np.testing.assert_almost_equal(
+        recall,
+        rb_rec,
+        decimal=6,
+        err_msg=f"Recall mismatch: ts_recall={recall} vs prts={rb_rec}",
+    )
+    np.testing.assert_almost_equal(
+        f1_score,
+        rb_fsc,
+        decimal=6,
+        err_msg=f"F1-Score mismatch: ts_fscore={f1_score} vs prts={rb_fsc}",
+    )
+
+
+@pytest.mark.skipif(
+    not _check_soft_dependencies("prts", severity="none"),
+    reason="required soft dependency prts not available",
+)
+def test_multiple_overlapping_ranges():
+    """Test for multiple overlapping ranges."""
+    y_pred_bin = np.array([0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0])
+    y_real_bin = np.array([0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1])
+
+    precision = ts_precision(y_pred_bin, y_real_bin, gamma="one", bias_type="flat")
+    recall = ts_recall(y_pred_bin, y_real_bin, gamma="one", bias_type="flat", alpha=0.0)
+    f1_score = ts_fscore(
+        y_pred_bin,
+        y_real_bin,
+        gamma="one",
+        beta=1,
+        p_bias="flat",
+        r_bias="flat",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    rb_prec = range_precision(y_real_bin, y_pred_bin, cardinality="one", bias="flat")
+    rb_rec = range_recall(y_real_bin, y_pred_bin, cardinality="one", bias="flat")
+    rb_fsc = range_f_score(
+        y_real_bin,
+        y_pred_bin,
+        beta=1,
+        cardinality="one",
+        p_bias="flat",
+        r_bias="flat",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    np.testing.assert_almost_equal(
+        precision,
+        rb_prec,
+        decimal=6,
+        err_msg=f"Precision mismatch: ts_precision={precision} vs prts={rb_prec}",
+    )
+    np.testing.assert_almost_equal(
+        recall,
+        rb_rec,
+        decimal=6,
+        err_msg=f"Recall mismatch: ts_recall={recall} vs prts={rb_rec}",
+    )
+    np.testing.assert_almost_equal(
+        f1_score,
+        rb_fsc,
+        decimal=6,
+        err_msg=f"F1-Score mismatch: ts_fscore={f1_score} vs prts={rb_fsc}",
+    )
+
+
+@pytest.mark.skipif(
+    not _check_soft_dependencies("prts", severity="none"),
+    reason="required soft dependency prts not available",
+)
+def test_nested_lists_of_predictions():
+    """Test for nested lists of predictions."""
+    y_pred_bin = np.array([0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 1, 1, 1])
+    y_real_bin = np.array([0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1, 0, 0])
+
+    precision = ts_precision(y_pred_bin, y_real_bin, gamma="one", bias_type="flat")
+    recall = ts_recall(y_pred_bin, y_real_bin, gamma="one", bias_type="flat", alpha=0.0)
+    f1_score = ts_fscore(
+        y_pred_bin,
+        y_real_bin,
+        gamma="one",
+        beta=1,
+        p_bias="flat",
+        r_bias="flat",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    rb_prec = range_precision(y_real_bin, y_pred_bin, cardinality="one", bias="flat")
+    rb_rec = range_recall(y_real_bin, y_pred_bin, cardinality="one", bias="flat")
+    rb_fsc = range_f_score(
+        y_real_bin,
+        y_pred_bin,
+        beta=1,
+        cardinality="one",
+        p_bias="flat",
+        r_bias="flat",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    np.testing.assert_almost_equal(
+        precision,
+        rb_prec,
+        decimal=6,
+        err_msg=f"Precision mismatch: ts_precision={precision} vs prts={rb_prec}",
+    )
+    np.testing.assert_almost_equal(
+        recall,
+        rb_rec,
+        decimal=6,
+        err_msg=f"Recall mismatch: ts_recall={recall} vs prts={rb_rec}",
+    )
+    np.testing.assert_almost_equal(
+        f1_score,
+        rb_fsc,
+        decimal=6,
+        err_msg=f"F1-Score mismatch: ts_fscore={f1_score} vs prts={rb_fsc}",
+    )
+
+
+@pytest.mark.skipif(
+    not _check_soft_dependencies("prts", severity="none"),
+    reason="required soft dependency prts not available",
+)
+def test_all_encompassing_range():
+    """Test for all encompassing range."""
+    y_pred_bin = np.array([0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1])
+    y_real_bin = np.array([0, 0, 1, 1, 0, 1, 1, 0, 1, 1, 0])
+
+    precision = ts_precision(y_pred_bin, y_real_bin, gamma="one", bias_type="flat")
+    recall = ts_recall(y_pred_bin, y_real_bin, gamma="one", bias_type="flat", alpha=0.0)
+    f1_score = ts_fscore(
+        y_pred_bin,
+        y_real_bin,
+        gamma="one",
+        beta=1,
+        p_bias="flat",
+        r_bias="flat",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    rb_prec = range_precision(y_real_bin, y_pred_bin, cardinality="one", bias="flat")
+    rb_rec = range_recall(y_real_bin, y_pred_bin, cardinality="one", bias="flat")
+    rb_fsc = range_f_score(
+        y_real_bin,
+        y_pred_bin,
+        beta=1,
+        cardinality="one",
+        p_bias="flat",
+        r_bias="flat",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    np.testing.assert_almost_equal(
+        precision,
+        rb_prec,
+        decimal=6,
+        err_msg=f"Precision mismatch: ts_precision={precision} vs prts={rb_prec}",
+    )
+    np.testing.assert_almost_equal(
+        recall,
+        rb_rec,
+        decimal=6,
+        err_msg=f"Recall mismatch: ts_recall={recall} vs prts={rb_rec}",
+    )
+    np.testing.assert_almost_equal(
+        f1_score,
+        rb_fsc,
+        decimal=6,
+        err_msg=f"F1-Score mismatch: ts_fscore={f1_score} vs prts={rb_fsc}",
+    )
+
+
+@pytest.mark.skipif(
+    not _check_soft_dependencies("prts", severity="none"),
+    reason="required soft dependency prts not available",
+)
+def test_binary_input_example():
+    """Comparing prts package(_binary.py) and range-based metrics(range_metrics.py)."""
+    y_pred_binary = np.array([0, 1, 1, 0])
+    y_true_binary = np.array([0, 1, 0, 0])
+
+    precision = ts_precision(
+        y_pred_binary, y_true_binary, gamma="reciprocal", bias_type="flat"
+    )
+    recall = ts_recall(
+        y_pred_binary, y_true_binary, gamma="reciprocal", bias_type="flat", alpha=0.0
+    )
+    f1_score = ts_fscore(
+        y_pred_binary,
+        y_true_binary,
+        beta=1,
+        gamma="reciprocal",
+        p_bias="flat",
+        r_bias="flat",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    rb_prec = range_precision(
+        y_true_binary, y_pred_binary, cardinality="reciprocal", bias="flat"
+    )
+    rb_rec = range_recall(
+        y_true_binary, y_pred_binary, cardinality="reciprocal", bias="flat"
+    )
+    rb_fsc = range_f_score(
+        y_true_binary,
+        y_pred_binary,
+        beta=1,
+        cardinality="reciprocal",
+        p_bias="flat",
+        r_bias="flat",
+    )
+
+    np.testing.assert_almost_equal(
+        precision,
+        rb_prec,
+        decimal=6,
+        err_msg=f"Precision mismatch: ts_precision={precision} vs prts={rb_prec}",
+    )
+    np.testing.assert_almost_equal(
+        recall,
+        rb_rec,
+        decimal=6,
+        err_msg=f"Recall mismatch: ts_recall={recall} vs prts={rb_rec}",
+    )
+    np.testing.assert_almost_equal(
+        f1_score,
+        rb_fsc,
+        decimal=6,
+        err_msg=f"F1-Score mismatch: ts_fscore={f1_score} vs prts={rb_fsc}",
+    )
+
+
+@pytest.mark.skipif(
+    not _check_soft_dependencies("prts", severity="none"),
+    reason="required soft dependency prts not available",
+)
+def test_multiple_overlapping_ranges_with_gamma_reciprocal():
+    """Test for multiple overlapping ranges with gamma=reciprocal."""
+    y_pred_bin = np.array([0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0])
+    y_real_bin = np.array([0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1])
+
+    precision = ts_precision(
+        y_pred_bin, y_real_bin, gamma="reciprocal", bias_type="flat"
+    )
+    recall = ts_recall(
+        y_pred_bin, y_real_bin, gamma="reciprocal", bias_type="flat", alpha=0.0
+    )
+    f1_score = ts_fscore(
+        y_pred_bin,
+        y_real_bin,
+        gamma="reciprocal",
+        beta=1,
+        p_bias="flat",
+        r_bias="flat",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    rb_prec = range_precision(
+        y_real_bin, y_pred_bin, cardinality="reciprocal", bias="flat"
+    )
+    rb_rec = range_recall(y_real_bin, y_pred_bin, cardinality="reciprocal", bias="flat")
+    rb_fsc = range_f_score(
+        y_real_bin,
+        y_pred_bin,
+        beta=1,
+        cardinality="reciprocal",
+        p_bias="flat",
+        r_bias="flat",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    np.testing.assert_almost_equal(
+        precision,
+        rb_prec,
+        decimal=6,
+        err_msg=f"Precision mismatch: ts_precision={precision} vs prts={rb_prec}",
+    )
+    np.testing.assert_almost_equal(
+        recall,
+        rb_rec,
+        decimal=6,
+        err_msg=f"Recall mismatch: ts_recall={recall} vs prts={rb_rec}",
+    )
+    np.testing.assert_almost_equal(
+        f1_score,
+        rb_fsc,
+        decimal=6,
+        err_msg=f"F1-Score mismatch: ts_fscore={f1_score} vs prts={rb_fsc}",
+    )
+
+
+@pytest.mark.skipif(
+    not _check_soft_dependencies("prts", severity="none"),
+    reason="required soft dependency prts not available",
+)
+def test_multiple_overlapping_ranges_with_bias_middle():
+    """Test for multiple overlapping ranges with bias_type=middle using range-binary conversion."""  # noqa E501
+    y_pred_bin = np.array([0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0])
+    y_real_bin = np.array([0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1])
+
+    # metric functions can also handle range-based input (internal binary conversion)
+    y_pred = _binary_to_ranges(
+        y_pred_bin
+    )  # Convert binary to range to show compatibility
+    y_real = _binary_to_ranges(y_real_bin)
+
+    precision = ts_precision(y_pred, y_real, gamma="one", bias_type="middle")
+    recall = ts_recall(y_pred, y_real, gamma="one", bias_type="middle", alpha=0.0)
+    f1_score = ts_fscore(
+        y_pred,
+        y_real,
+        gamma="one",
+        beta=1,
+        p_bias="middle",
+        r_bias="middle",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    rb_prec = range_precision(y_real_bin, y_pred_bin, cardinality="one", bias="middle")
+    rb_rec = range_recall(y_real_bin, y_pred_bin, cardinality="one", bias="middle")
+    rb_fsc = range_f_score(
+        y_real_bin,
+        y_pred_bin,
+        beta=1,
+        cardinality="one",
+        p_bias="middle",
+        r_bias="middle",
+        p_alpha=0.0,
+        r_alpha=0.0,
+    )
+
+    np.testing.assert_almost_equal(
+        precision,
+        rb_prec,
+        decimal=6,
+        err_msg=f"Precision mismatch: ts_precision={precision} vs prts={rb_prec}",
+    )
+    np.testing.assert_almost_equal(
+        recall,
+        rb_rec,
+        decimal=6,
+        err_msg=f"Recall mismatch: ts_recall={recall} vs prts={rb_rec}",
+    )
+    np.testing.assert_almost_equal(
+        f1_score,
+        rb_fsc,
+        decimal=6,
+        err_msg=f"F1-Score mismatch: ts_fscore={f1_score} vs prts={rb_fsc}",
+    )
+
+
+@pytest.mark.skipif(
+    not _check_soft_dependencies("prts", severity="none"),
+    reason="required soft dependency prts not available",
+)
+def test_multiple_overlapping_ranges_with_bias_middle_gamma_reciprocal():
+    """Test for multiple overlapping ranges with bias_type=middle, gamma=reciprocal."""
+    y_pred_bin = np.array([0, 1, 1, 1, 0, 1, 1, 1, 0, 0, 0])
+    y_real_bin = np.array([0, 0, 1, 1, 1, 1, 1, 0, 1, 1, 1])
+
+    precision = ts_precision(
+        y_pred_bin, y_real_bin, gamma="reciprocal", bias_type="middle"
+    )
+    recall = ts_recall(
+        y_pred_bin, y_real_bin, gamma="reciprocal", bias_type="middle", alpha=1.0
+    )
+    f1_score = ts_fscore(
+        y_pred_bin,
+        y_real_bin,
+        gamma="reciprocal",
+        beta=1,
+        p_bias="middle",
+        r_bias="middle",
+        p_alpha=0.0,
+        r_alpha=1.0,
+    )
+
+    rb_prec = range_precision(
+        y_real_bin, y_pred_bin, cardinality="reciprocal", bias="middle"
+    )
+    rb_rec = range_recall(
+        y_real_bin, y_pred_bin, cardinality="reciprocal", bias="middle", alpha=1.0
+    )
+    rb_fsc = range_f_score(
+        y_real_bin,
+        y_pred_bin,
+        beta=1,
+        cardinality="reciprocal",
+        p_bias="middle",
+        r_bias="middle",
+        p_alpha=0.0,
+        r_alpha=1.0,
+    )
+
+    np.testing.assert_almost_equal(
+        precision,
+        rb_prec,
+        decimal=6,
+        err_msg=f"Precision mismatch: ts_precision={precision} vs prts={rb_prec}",
+    )
+    np.testing.assert_almost_equal(
+        recall,
+        rb_rec,
+        decimal=6,
+        err_msg=f"Recall mismatch: ts_recall={recall} vs prts={rb_rec}",
+    )
+    np.testing.assert_almost_equal(
+        f1_score,
+        rb_fsc,
+        decimal=6,
+        err_msg=f"F1-Score mismatch: ts_fscore={f1_score} vs prts={rb_fsc}",
+    )
diff --git a/aeon/clustering/tests/test_kernel_k_means.py b/aeon/clustering/tests/test_kernel_k_means.py
index f4af21f4..36a761a4 100644
--- a/aeon/clustering/tests/test_kernel_k_means.py
+++ b/aeon/clustering/tests/test_kernel_k_means.py
@@ -13,6 +13,12 @@ expected_iters = 2
 
 expected_results = [0, 0, 0, 0, 0]
 
+expected_labels_kdtw = [0, 0, 0, 1, 2]
+
+expected_iters_kdtw = 2
+
+expected_results_kdtw = [0, 2, 0, 0, 0]
+
 
 @pytest.mark.skipif(
     not _check_estimator_deps(TimeSeriesKernelKMeans, severity="none"),
@@ -37,3 +43,21 @@ def test_kernel_k_means():
 
     for val in proba:
         assert np.count_nonzero(val == 1.0) == 1
+
+    kernel_kmeans_kdtw = TimeSeriesKernelKMeans(
+        kernel="kdtw",
+        random_state=1,
+        n_clusters=3,
+        kernel_params={"sigma": 2.0, "epsilon": 1e-4},
+    )
+    kernel_kmeans_kdtw.fit(X_train[0:max_train])
+    kdtw_results = kernel_kmeans_kdtw.predict(X_test[0:max_train])
+    kdtw_proba = kernel_kmeans_kdtw.predict_proba(X_test[0:max_train])
+
+    assert np.array_equal(kdtw_results, expected_results_kdtw)
+    assert kernel_kmeans_kdtw.n_iter_ == expected_iters_kdtw
+    assert np.array_equal(kernel_kmeans_kdtw.labels_, expected_labels_kdtw)
+    assert kdtw_proba.shape == (max_train, 3)
+
+    for val in kdtw_proba:
+        assert np.count_nonzero(val == 1.0) == 1
diff --git a/aeon/distances/elastic/tests/test_bounding.py b/aeon/distances/elastic/tests/test_bounding.py
index 32ac1987..f0c2d837 100644
--- a/aeon/distances/elastic/tests/test_bounding.py
+++ b/aeon/distances/elastic/tests/test_bounding.py
@@ -1,7 +1,6 @@
 """Test for bounding matrix."""
 
 import numpy as np
-import pytest
 
 from aeon.distances import create_bounding_matrix
 
@@ -37,10 +36,34 @@ def test_itakura_parallelogram():
     matrix = create_bounding_matrix(10, 10, itakura_max_slope=0.2)
     assert isinstance(matrix, np.ndarray)
 
-    with pytest.raises(
-        ValueError,
-        match="""Itakura parallelogram does not support unequal length time series.
-Please consider using a full bounding matrix or a sakoe chiba bounding matrix
-instead.""",
-    ):
-        create_bounding_matrix(5, 10, itakura_max_slope=0.2)
+    expected_result_5_7 = np.array(
+        [
+            [True, False, False, False, False, False, False],
+            [False, True, True, True, True, False, False],
+            [False, False, True, True, True, False, False],
+            [False, False, True, True, True, True, False],
+            [False, False, False, False, False, False, True],
+        ]
+    )
+
+    expected_result_7_5 = np.array(
+        [
+            [True, False, False, False, False],
+            [False, True, False, False, False],
+            [False, True, True, True, False],
+            [False, True, True, True, False],
+            [False, True, True, True, False],
+            [False, False, False, True, False],
+            [False, False, False, False, True],
+        ]
+    )
+
+    matrix = create_bounding_matrix(5, 7, itakura_max_slope=0.5)
+    assert isinstance(matrix, np.ndarray)
+    assert matrix.shape == (5, 7)
+    assert np.array_equal(matrix, expected_result_5_7)
+
+    matrix = create_bounding_matrix(7, 5, itakura_max_slope=0.5)
+    assert isinstance(matrix, np.ndarray)
+    assert matrix.shape == (7, 5)
+    assert np.array_equal(matrix, expected_result_7_5)
diff --git a/aeon/distances/elastic/tests/test_distance_correctness.py b/aeon/distances/elastic/tests/test_distance_correctness.py
index f33ca088..2a0c203b 100644
--- a/aeon/distances/elastic/tests/test_distance_correctness.py
+++ b/aeon/distances/elastic/tests/test_distance_correctness.py
@@ -10,6 +10,7 @@ from aeon.datasets import load_basic_motions, load_unit_test
 from aeon.distances import (
     ddtw_distance,
     dtw_distance,
+    dtw_gi_distance,
     edr_distance,
     erp_distance,
     euclidean_distance,
@@ -23,6 +24,7 @@ from aeon.distances import (
 
 distances = [
     "dtw",
+    "dtw_gi",
     "wdtw",
     "lcss",
     "msm",
@@ -36,6 +38,7 @@ distances = [
 
 distance_parameters = {
     "dtw": [0.0, 0.1, 1.0],  # window
+    "dtw_gi": [0.0, 0.1, 1.0],  # window
     "wdtw": [0.0, 0.1, 1.0],  # parameter g
     "wddtw": [0.0, 0.1, 1.0],  # parameter g
     "erp": [0.0, 0.1, 1.0],  # window
@@ -63,6 +66,7 @@ basic_motions_distances = {
     "euclidean": 27.51835240,
     "squared": 757.25971908652,
     "dtw": [757.259719, 330.834497, 330.834497],
+    "dtw_gi": [259.5333502342899, 310.10738471013804, 310.10738471013804],
     "wdtw": [165.41724, 3.308425, 0],
     "msm": [70.014828, 89.814828, 268.014828],
     "erp": [169.3715, 102.0979, 102.097904],
@@ -90,6 +94,8 @@ def test_multivariate_correctness():
     for j in range(0, 3):
         d = dtw_distance(case1, case2, window=distance_parameters["dtw"][j])
         assert_almost_equal(d, basic_motions_distances["dtw"][j], 4)
+        d = dtw_gi_distance(case1, case2, window=distance_parameters["dtw_gi"][j])
+        assert_almost_equal(d, basic_motions_distances["dtw_gi"][j], 4)
         d = wdtw_distance(case1, case2, g=distance_parameters["wdtw"][j])
         assert_almost_equal(d, basic_motions_distances["wdtw"][j], 4)
         d = lcss_distance(case1, case2, epsilon=distance_parameters["lcss"][j] / 50.0)
diff --git a/aeon/distances/tests/test_distances.py b/aeon/distances/tests/test_distances.py
index 3633c9ca..4efe396a 100644
--- a/aeon/distances/tests/test_distances.py
+++ b/aeon/distances/tests/test_distances.py
@@ -57,6 +57,7 @@ def _validate_distance_result(
     assert isinstance(dist_result_to_self, float)
 
     # If unequal length swap where x and y are to ensure it works both ways around
+
     if symmetric and original_x.shape[-1] != original_y.shape[-1] and check_xy_permuted:
         _validate_distance_result(
             original_y,
diff --git a/aeon/testing/expected_results/expected_distance_results.py b/aeon/testing/expected_results/expected_distance_results.py
index 8e7c6873..7126c5c6 100644
--- a/aeon/testing/expected_results/expected_distance_results.py
+++ b/aeon/testing/expected_results/expected_distance_results.py
@@ -41,6 +41,13 @@ _expected_distance_results = {
         4.0997661869195205,
         25.0,
     ],
+    "dtw_gi": [
+        0.344520787316184,
+        0.344520787316184,
+        0.2998607605839068,
+        5.893210968537887,
+        25.0,
+    ],
     "ddtw": [
         0.2963709096971962,
         0.2963709096971962,
diff --git a/aeon/testing/testing_config.py b/aeon/testing/testing_config.py
index 3d05d667..be7d04f0 100644
--- a/aeon/testing/testing_config.py
+++ b/aeon/testing/testing_config.py
@@ -24,7 +24,6 @@ NUMBA_DISABLED = os.environ.get("NUMBA_DISABLE_JIT") == "1"
 
 # exclude estimators here for short term fixes
 EXCLUDE_ESTIMATORS = [
-    "REDCOMETS",
     "HydraTransformer",  # returns a pytorch Tensor
 ]
 
@@ -52,7 +51,7 @@ EXCLUDED_TESTS = {
     "RSASTClassifier": ["check_fit_deterministic"],
     "SAST": ["check_fit_deterministic"],
     "RSAST": ["check_fit_deterministic"],
-    "MatrixProfile": ["check_persistence_via_pickle"],
+    "MatrixProfile": ["check_fit_deterministic", "check_persistence_via_pickle"],
     # missed in legacy testing, changes state in predict/transform
     "FLUSSSegmenter": ["check_non_state_changing_method"],
     "InformationGainSegmenter": ["check_non_state_changing_method"],
@@ -73,6 +72,10 @@ EXCLUDED_TESTS = {
 EXCLUDED_TESTS_NO_NUMBA = {
     # See issue #622
     "HIVECOTEV2": ["check_classifier_against_expected_results"],
+    # Other failures
+    "TemporalDictionaryEnsemble": ["check_classifier_against_expected_results"],
+    "OrdinalTDE": ["check_classifier_against_expected_results"],
+    "CanonicalIntervalForestRegressor": ["check_regressor_against_expected_results"],
 }