Repository: aeon-toolkit/aeon
Pull request #1580 (2024-05-29 12:46:33)

[ENH] DWT-MLEAD anomaly detection method
#### Reference Issues/PRs

Depends on #1579 being merged first!

#### What does this implement/fix? Explain your changes.

Adds the anomaly detection method DWT-MLEAD to the `aeon.anomaly_detection` package (https://blogs.gm.fh-koeln.de/ciop/files/2019/01/thillwavelet.pdf).

#### Does your contribution introduce a new dependency? If yes, which one?

no

### PR checklist

##### For all contributions
- [x] I've added myself to the [list of contributors](https://github.com/aeon-toolkit/aeon/blob/main/.all-contributorsrc). Alternatively, you can use the [@all-contributors](https://allcontributors.org/docs/en/bot/usage) bot to do this for you.
- [x] The PR title starts with either [ENH], [MNT], [DOC], [BUG], [REF], [DEP] or [GOV] indicating whether the PR topic is related to enhancement, maintenance, documentation, bugs, refactoring, deprecation or governance.

##### For new estimators and functions
- [x] I've added the estimator to the online [API documentation](https://www.aeon-toolkit.org/en/latest/api_reference.html).
- [x] (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
- [x] (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 would have added the following labels to this PR based on the changes made: **[ $\\color{#6F6E8D}{\\textsf{anomaly detection}}$ ]**, however some package labels are already present.

The [Checks](https://github.com/aeon-toolkit/aeon/pull/1580/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.
    
MatthewMiddlehurst: Ah, didnt realise this was not being merged into main. CI does not run properly currently if it is not, we can probably make some changes so that it does though.

Patch:
diff --git a/aeon/anomaly_detection/__init__.py b/aeon/anomaly_detection/__init__.py
index edaf9ec2..ddbc2c12 100644
--- a/aeon/anomaly_detection/__init__.py
+++ b/aeon/anomaly_detection/__init__.py
@@ -1,9 +1,11 @@
 """Time Series Anomaly Detection."""
 
 __all__ = [
+    "DWT_MLEAD",
     "MERLIN",
     "STRAY",
 ]
 
+from aeon.anomaly_detection._dwt_mlead import DWT_MLEAD
 from aeon.anomaly_detection._merlin import MERLIN
 from aeon.anomaly_detection._stray import STRAY
diff --git a/aeon/anomaly_detection/_dwt_mlead.py b/aeon/anomaly_detection/_dwt_mlead.py
new file mode 100644
index 00000000..22c5d57c
--- /dev/null
+++ b/aeon/anomaly_detection/_dwt_mlead.py
@@ -0,0 +1,250 @@
+"""DWT-MLEAD anomaly detector."""
+
+__maintainer__ = ["CodeLionX"]
+__all__ = ["DWT_MLEAD"]
+
+import warnings
+from typing import Any, Iterable, List, Tuple
+
+import numpy as np
+from numpy.lib.stride_tricks import sliding_window_view
+from sklearn.covariance import EmpiricalCovariance
+
+from aeon.anomaly_detection.base import BaseAnomalyDetector
+from aeon.utils.numba.wavelets import multilevel_haar_transform
+
+
+def _pad_series(x: np.ndarray) -> Tuple[np.ndarray, int, int]:
+    """Pad input signal to the next power of 2 using periodic padding mode."""
+    n = x.shape[0]
+    exp = np.ceil(np.log2(n))
+    m = int(np.power(2, exp))
+    return np.pad(x, (0, m - n), mode="wrap"), n, m
+
+
+def _combine_alternating(xs: List[Any], ys: List[Any]) -> Iterable[Any]:
+    """Combine two lists by alternating their elements."""
+    for x, y in zip(xs, ys):
+        yield x
+        yield y
+
+
+class DWT_MLEAD(BaseAnomalyDetector):
+    """DWT-MLEAD anomaly detector.
+
+    DWT-MLEAD is an anomaly detection algorithm that uses the Discrete Wavelet Transform
+    (DWT) and Maximum Likelihood Estimation (MLE) to detect anomalies in univariate
+    time series. The algorithm performs mutli-level DWT using the Haar wavelet, slides
+    windows over the DWT coefficients, and estimates the likelihood of each window
+    using a Gaussian distribution. Anomalies are detected by comparing the likelihoods
+    to a quantile boundary in each level and passing down the anomaly counts to the
+    individual time points, which we use as anomaly scores. The original paper [1]_
+    subsequently clusters the anomalies to determine the anomaly centers. This step is
+    not implemented in this version.
+
+    Notes
+    -----
+    This implementation does not exactly match the original paper [1]_. We make the
+    following changes:
+
+    - We use window sizes for the DWT coefficients that decrease with the level number
+      because otherwise we would have too few items to slide the window over.
+    - We exclude the highest level coefficients because they contain only a single entry
+      and are, thus, not suitable for sliding a window of length 2 over it.
+    - We have not implemented the Monte Carlo quantile boundary type yet.
+    - We do not perform the anomaly clustering step to determine the anomaly centers.
+      Instead, we return the anomaly scores for each timestep in the original time
+      series.
+
+    Parameters
+    ----------
+    start_level : int, default=3
+        The level at which to start the anomaly detection. Must be >= 0 and less than
+        log_2(n_timepoints).
+    quantile_boundary_type : str, default='percentile'
+        The type of boundary to use for the quantile. Must be 'percentile',
+        'monte-carlo' is not implemented yet.
+    quantile_epsilon : float, default=0.01
+        The epsilon value for the quantile boundary. Must be in [0, 1].
+
+    References
+    ----------
+    .. [1] Thill, Markus, Wolfgang Konen, and Thomas Bäck. "Time Series Anomaly
+           Detection with Discrete Wavelet Transforms and Maximum Likelihood
+           Estimation." In Proceedings of the International Conference on Time Series
+           (ITISE). Granada, Spain, 2017.
+
+    Examples
+    --------
+    >>> import numpy as np
+    >>> from aeon.anomaly_detection import DWT_MLEAD
+    >>> X = np.array([1, 2, 3, 4, 1, 2, 3, 3, 2, 8, 9, 8, 1, 2, 3, 4], dtype=np.float_)
+    >>> detector = DWT_MLEAD(
+    ...    start_level=1, quantile_boundary_type='percentile', quantile_epsilon=0.01
+    ... )
+    >>> detector.fit_predict(X)
+    array([0., 0., 0., 0., 0., 0., 0., 0., 2., 2., 2., 2., 0., 0., 0., 0.])
+    """
+
+    _tags = {
+        "capability:univariate": True,
+        "capability:multivariate": False,
+        "capability:missing_values": False,
+    }
+
+    def __init__(
+        self,
+        start_level: int = 3,
+        quantile_boundary_type: str = "percentile",
+        quantile_epsilon: float = 0.01,
+    ):
+        self.start_level = start_level
+        self.quantile_boundary_type = quantile_boundary_type
+        self.quantile_epsilon = quantile_epsilon
+
+        super().__init__(axis=0)
+
+    def _predict(self, X) -> np.ndarray:
+        X = X.squeeze()
+
+        if self.start_level < 0:
+            raise ValueError("start_level must be >= 0")
+        if self.quantile_boundary_type != "percentile":
+            if self.quantile_boundary_type not in ["percentile", "monte-carlo"]:
+                raise ValueError(
+                    "quantile_boundary_type must be 'percentile' or 'monte-carlo', "
+                    f"but is {self.quantile_boundary_type}"
+                )
+            else:
+                raise NotImplementedError(
+                    f"The quantile boundary type '{self.quantile_boundary_type}' "
+                    "is not implemented yet!"
+                )
+        if self.quantile_epsilon < 0 or self.quantile_epsilon > 1:
+            raise ValueError("quantile_epsilon must be in [0, 1]")
+
+        X, self.N_, self.M_ = _pad_series(X)
+        max_level = int(np.log2(self.M_))
+
+        if self.start_level >= max_level:
+            raise ValueError(
+                f"start_level ({self.start_level}) must be less than "
+                f"log_2(n_timepoints) ({max_level})"
+            )
+
+        # perform multilevel DWT and capture coefficients
+        levels, approx_coeffs, detail_coeffs = self._multilevel_dwt(X, max_level)
+
+        # extract anomalies in each level
+        window_sizes = np.array(
+            [
+                max(2, max_level - level - self.start_level + 1)
+                for level in range(max_level)
+            ],
+            dtype=np.int_,
+        )
+        coef_anomaly_counts = []
+        for x, level in zip(
+            _combine_alternating(detail_coeffs, approx_coeffs), levels.repeat(2, axis=0)
+        ):
+            w = window_sizes[level]
+            windows = sliding_window_view(x, w)
+
+            p = self._estimate_gaussian_likelihoods(windows)
+            a = self._mark_anomalous_windows(p)
+            xa = self._reverse_windowing(a, window_length=w, full_length=x.shape[0])
+            coef_anomaly_counts.append(xa)
+
+        # aggregate anomaly counts (leaf counters)
+        point_anomaly_scores = self._push_anomaly_counts_down_to_points(
+            coef_anomaly_counts
+        )
+        return point_anomaly_scores
+
+    def _multilevel_dwt(
+        self, X: np.ndarray, max_level: int
+    ) -> Tuple[np.ndarray, List[np.ndarray], List[np.ndarray]]:
+        ls_ = np.arange(self.start_level - 1, max_level - 1, dtype=np.int_) + 1
+        as_, ds_ = multilevel_haar_transform(X, max_level - 1)
+        as_ = as_[self.start_level :]
+        ds_ = ds_[self.start_level - 1 :]
+        return ls_, as_, ds_
+
+    @staticmethod
+    def _estimate_gaussian_likelihoods(x_windows: np.ndarray) -> np.ndarray:
+        with warnings.catch_warnings():
+            warnings.filterwarnings(action="ignore", category=UserWarning)
+
+            # fit gaussion distribution with mean and covariance
+            estimator = EmpiricalCovariance(assume_centered=False)
+            estimator.fit(x_windows)
+
+            # compute log likelihood for each window x in x_view
+            n_windows = x_windows.shape[0]
+            p = np.empty(shape=n_windows)
+            for i in range(n_windows):
+                p[i] = estimator.score(x_windows[i].reshape(1, -1))
+        return p
+
+    def _mark_anomalous_windows(self, p: np.ndarray) -> np.ndarray:
+        if self.quantile_boundary_type == "percentile":
+            z_eps = np.percentile(p, self.quantile_epsilon * 100)
+        else:  # self.quantile_boundary_type == "monte-carlo"
+            raise ValueError(
+                f"The quantile boundary type '{self.quantile_boundary_type}' "
+                "is not implemented yet!"
+            )
+
+        return p < z_eps
+
+    @staticmethod
+    def _reverse_windowing(
+        x: np.ndarray, window_length: int, full_length: int
+    ) -> np.ndarray:
+        mapped = np.full(shape=(full_length, window_length), fill_value=0)
+        mapped[: x.shape[0], 0] = x
+
+        for w in range(1, window_length):
+            mapped[:, w] = np.roll(mapped[:, 0], w)
+
+        return np.sum(mapped, axis=1)
+
+    def _push_anomaly_counts_down_to_points(
+        self, coef_anomaly_counts: List[np.ndarray]
+    ) -> np.ndarray:
+        # sum up counters of detail coeffs (orig. D^l) and approx coeffs (orig. C^l)
+        anomaly_counts = coef_anomaly_counts[0::2] + coef_anomaly_counts[1::2]
+
+        # extrapolate anomaly counts to the original series' points
+        counter = np.zeros(self.M_)
+        for ac in anomaly_counts:
+            counter += ac.repeat(self.M_ // ac.shape[0], axis=0)
+        # set event counters with count < 2 to 0
+        counter[counter < 2] = 0
+        return counter[: self.N_]
+
+    @classmethod
+    def get_test_params(cls, parameter_set="default"):
+        """Return testing parameter settings for the estimator.
+
+        Only supports 'default'-parameter set.
+
+        Parameters
+        ----------
+        parameter_set : str, default="default"
+            Name of the set of test parameters to return, for use in tests. If no
+            special parameters are defined for a value, will return `"default"` set.
+
+        Returns
+        -------
+        params : dict or list of dict, default={}
+            Parameters to create testing instances of the class.
+            Each dict are parameters to construct an "interesting" test instance, i.e.,
+            `MyClass(**params)` or `MyClass(**params[i])` creates a valid test instance.
+            `create_test_instance` uses the first (or only) dictionary in `params`.
+        """
+        return {
+            "start_level": 2,
+            "quantile_boundary_type": "percentile",
+            "quantile_epsilon": 0.01,
+        }
diff --git a/aeon/base/_base.py b/aeon/base/_base.py
index 5132e564..7b046e58 100644
--- a/aeon/base/_base.py
+++ b/aeon/base/_base.py
@@ -332,7 +332,7 @@ class BaseObject(_BaseEstimator):
         See Also
         --------
         get_tag : Get a single tag from an object.
-        get_clas_tags : Get all tags from a class.
+        get_class_tags : Get all tags from a class.
         get_class_tag : Get a single tag from a class.
 
         Examples
diff --git a/aeon/forecasting/base/_fh.py b/aeon/forecasting/base/_fh.py
index ea14d0d4..940e4369 100644
--- a/aeon/forecasting/base/_fh.py
+++ b/aeon/forecasting/base/_fh.py
@@ -770,7 +770,7 @@ def _check_cutoff(cutoff, index):
         assert isinstance(cutoff, pd.Timestamp)
 
 
-def _coerce_to_period(x, freq=None):
+def _coerce_to_period(x, freq):
     """Coerce pandas time index to a alternative pandas time index.
 
     This coerces pd.Timestamp to pd.Period or pd.DatetimeIndex to
@@ -788,10 +788,6 @@ def _coerce_to_period(x, freq=None):
     index : pd.Period or pd.PeriodIndex
         Index or index element coerced to period based format.
     """
-    if isinstance(x, pd.Timestamp) and freq is None:
-        raise ValueError(
-            "_coerce_to_period requires freq argument to be passed if x is pd.Timestamp"
-        )
     try:
         return x.to_period(freq)
     except (ValueError, AttributeError) as e:
diff --git a/docs/api_reference/anomaly_detection.rst b/docs/api_reference/anomaly_detection.rst
index 10d92e51..b5286be9 100644
--- a/docs/api_reference/anomaly_detection.rst
+++ b/docs/api_reference/anomaly_detection.rst
@@ -6,7 +6,7 @@ Anomaly Detection
 Time Series Anomaly Detection aims at discovering regions of a time series that in
 some way not representative of the underlying generative process.
 The :mod:`aeon.anomaly_detection` module contains algorithms and tools
-for time series segmentation.
+for time series anomaly detection.
 
 .. currentmodule:: aeon.anomaly_detection
 
@@ -15,3 +15,5 @@ for time series segmentation.
     :template: class.rst
 
     STRAY
+    MERLIN
+    DWT_MLEAD
diff --git a/docs/api_reference/transformations.rst b/docs/api_reference/transformations.rst
index 042ee5ea..99546d95 100644
--- a/docs/api_reference/transformations.rst
+++ b/docs/api_reference/transformations.rst
@@ -355,6 +355,7 @@ Binning and segmentation
 
     IntervalSegmenter
     RandomIntervalSegmenter
+    SlidingWindowSegmenter
 
 Missing value imputation
 ~~~~~~~~~~~~~~~~~~~~~~~~


Test patch:
diff --git a/aeon/anomaly_detection/tests/test_dwt_mlead.py b/aeon/anomaly_detection/tests/test_dwt_mlead.py
new file mode 100644
index 00000000..ba055bed
--- /dev/null
+++ b/aeon/anomaly_detection/tests/test_dwt_mlead.py
@@ -0,0 +1,32 @@
+"""Tests for the DWT_MLEAD class."""
+
+__maintainer__ = ["CodeLionX"]
+
+import numpy as np
+import pytest
+
+from aeon.anomaly_detection import DWT_MLEAD
+from aeon.testing.utils.data_gen import make_series
+
+
+def test_dwt_mlead_output():
+    """Test DWT_MLEAD output."""
+    series = make_series(n_timepoints=100, return_numpy=True, random_state=42)
+    series[50:58] -= 5
+
+    ad = DWT_MLEAD(start_level=2)
+    pred = ad.predict(series)
+
+    assert pred.shape == (100,)
+    assert pred.dtype == np.float_
+    assert 50 <= np.argmax(pred) <= 58
+
+
+def test_dwt_mlead_monte_carlo_unimplemented():
+    """Test that the monte-carlo method is not implemented."""
+    ad = DWT_MLEAD(quantile_boundary_type="monte-carlo")
+    with pytest.raises(
+        NotImplementedError,
+        match=".*quantile boundary type 'monte-carlo' is not implemented.*",
+    ):
+        ad.predict(np.array([0.5, 1.0, 0.8]))