diff --git a/aeon/transformations/collection/channel_selection/__init__.py b/aeon/transformations/collection/channel_selection/__init__.py index d5e18c4e26..e9c0ca387a 100644 --- a/aeon/transformations/collection/channel_selection/__init__.py +++ b/aeon/transformations/collection/channel_selection/__init__.py @@ -9,6 +9,7 @@ "ChannelScorer", "ElbowClassPairwise", "ElbowClassSum", + "MSTS", "RandomChannelSelector", "TSelect", ] @@ -21,6 +22,7 @@ ElbowClassPairwise, ElbowClassSum, ) +from aeon.transformations.collection.channel_selection._msts import MSTS from aeon.transformations.collection.channel_selection._random import ( RandomChannelSelector, ) diff --git a/aeon/transformations/collection/channel_selection/_msts.py b/aeon/transformations/collection/channel_selection/_msts.py new file mode 100644 index 0000000000..bb187b1ea4 --- /dev/null +++ b/aeon/transformations/collection/channel_selection/_msts.py @@ -0,0 +1,230 @@ +"""Merit Score for Time Series channel selection.""" + +__maintainer__ = ["aeon developers"] +__all__ = ["MSTS"] + +from itertools import combinations +from typing import Any + +import numpy as np +from sklearn.metrics import adjusted_mutual_info_score +from sklearn.model_selection import StratifiedKFold + +from aeon.classification.distance_based import KNeighborsTimeSeriesClassifier +from aeon.transformations.collection.channel_selection.base import BaseChannelSelector +from aeon.utils.validation import check_n_jobs + + +def _merit_score( + class_scores: np.ndarray, + pair_scores: np.ndarray, + channels: tuple[int, ...], +) -> float: + """Calculate the MSTS merit score for a channel subset.""" + k = len(channels) + dc = float(np.mean(class_scores[list(channels)])) + if k == 1: + return dc + + pair_values = pair_scores[np.ix_(channels, channels)] + upper_values = pair_values[np.triu_indices(k, 1)] + dd = float(np.mean(upper_values)) + denominator_squared = k + k * (k - 1) * dd + denominator = np.sqrt(max(denominator_squared, np.finfo(float).eps)) + return float(k * dc / denominator) + + +class MSTS(BaseChannelSelector): + """Select channels using the Merit Score for Time Series (MSTS). + + MSTS evaluates each channel with out-of-fold predictions from a 1-nearest + neighbour dynamic time warping classifier. It measures channel-to-class and + channel-to-channel agreement with adjusted mutual information, then uses a + greedy forward search to select the subset with the highest merit score. + + Parameters + ---------- + n_splits : int, default=3 + Number of stratified folds used to generate single-channel predictions. + n_jobs : int, default=1 + Number of parallel jobs used by the 1-NN DTW classifiers. + random_state : int or None, default=None + Controls shuffling of the stratified folds. If None, folds are not shuffled. + + Attributes + ---------- + channels_selected_ : list[int] + Indices of the selected channels, in their original channel order. + channel_predictions_ : np.ndarray + Out-of-fold predictions for each training case and channel. + channel_class_scores_ : np.ndarray + Adjusted mutual information between each channel's predictions and ``y``. + channel_pair_scores_ : np.ndarray + Pairwise adjusted mutual information between channel predictions. + selection_history_ : list[tuple[tuple[int, ...], float]] + Selected subsets and their merit scores at each forward-search step. + selected_merit_score_ : float + Merit score of the selected channel subset. + + Notes + ----- + MSTS is implemented as a channel selector, so it does not fit a final + classifier. The selected channels are applied by ``transform`` and can be + passed to a downstream multivariate time-series estimator. + + References + ---------- + .. [1] Kathirgamanathan, B. and Cunningham, P. "Correlation Based Feature + Subset Selection for Multivariate Time-Series Data." 2021. + https://arxiv.org/abs/2112.03705 + .. [2] aeon issue #1481, "Implement Merit Score Function channel selection + algorithm." + https://github.com/aeon-toolkit/aeon/issues/1481 + + Examples + -------- + >>> from aeon.testing.data_generation import make_example_3d_numpy + >>> from aeon.transformations.collection.channel_selection import MSTS + >>> X, y = make_example_3d_numpy( + ... n_cases=12, n_channels=3, n_timepoints=10, random_state=0 + ... ) + >>> selector = MSTS(n_splits=2, random_state=0).fit(X, y) + >>> X_selected = selector.transform(X) + >>> X_selected.shape[1] <= X.shape[1] + True + """ + + _tags = { + "capability:multivariate": True, + "requires_y": True, + "X_inner_type": "numpy3D", + } + + def __init__( + self, + n_splits: int = 3, + n_jobs: int = 1, + random_state: int | None = None, + ): + self.n_splits = n_splits + self.n_jobs = n_jobs + self.random_state = random_state + super().__init__() + + def _fit(self, X: np.ndarray, y: np.ndarray): + """Fit channel predictions, agreement scores, and the MSTS subset.""" + self._validate_parameters() + X = np.asarray(X, dtype=float) + y = np.asarray(y) + n_cases, n_channels, _ = X.shape + self.n_cases_ = n_cases + self.n_channels_ = n_channels + + class_counts = np.unique(y, return_counts=True)[1] + if class_counts.size == 0 or np.min(class_counts) < self.n_splits: + raise ValueError( + "Each class must contain at least n_splits cases for MSTS." + ) + + splitter = StratifiedKFold( + n_splits=self.n_splits, + shuffle=self.random_state is not None, + random_state=self.random_state, + ) + fold_indices = tuple(splitter.split(np.zeros(n_cases), y)) + + predictions = np.empty((n_cases, n_channels), dtype=y.dtype) + for channel in range(n_channels): + channel_data = X[:, channel : channel + 1, :] + for train_index, validation_index in fold_indices: + classifier = KNeighborsTimeSeriesClassifier( + n_neighbors=1, + distance="dtw", + n_jobs=self._n_jobs, + ) + classifier.fit(channel_data[train_index], y[train_index]) + predictions[validation_index, channel] = classifier.predict( + channel_data[validation_index] + ) + + self.channel_predictions_ = predictions + self.channel_class_scores_ = np.asarray( + [ + adjusted_mutual_info_score(y, predictions[:, channel]) + for channel in range(n_channels) + ], + dtype=float, + ) + self.channel_pair_scores_ = np.eye(n_channels, dtype=float) + for first, second in combinations(range(n_channels), 2): + score = adjusted_mutual_info_score( + predictions[:, first], predictions[:, second] + ) + self.channel_pair_scores_[first, second] = score + self.channel_pair_scores_[second, first] = score + + selected, history = self._forward_search() + self.channels_selected_ = list(selected) + self.selection_history_ = history + self.selected_merit_score_ = history[-1][1] + return self + + def _forward_search( + self, + ) -> tuple[tuple[int, ...], list[tuple[tuple[int, ...], float]]]: + """Select a channel subset by greedy forward merit maximisation.""" + if self.n_channels_ == 1: + subset = (0,) + return subset, [ + ( + subset, + _merit_score( + self.channel_class_scores_, self.channel_pair_scores_, subset + ), + ) + ] + + best_subset = (0, 1) + best_score = _merit_score( + self.channel_class_scores_, self.channel_pair_scores_, best_subset + ) + for subset in combinations(range(self.n_channels_), 2): + score = _merit_score( + self.channel_class_scores_, self.channel_pair_scores_, subset + ) + if score > best_score: + best_subset = subset + best_score = score + + history = [(best_subset, best_score)] + while len(best_subset) < self.n_channels_: + candidate_subset = None + candidate_score = best_score + for channel in range(self.n_channels_): + if channel in best_subset: + continue + subset = tuple(sorted((*best_subset, channel))) + score = _merit_score( + self.channel_class_scores_, self.channel_pair_scores_, subset + ) + if score > candidate_score: + candidate_subset = subset + candidate_score = score + if candidate_subset is None: + break + best_subset = candidate_subset + best_score = candidate_score + history.append((best_subset, best_score)) + + return best_subset, history + + def _validate_parameters(self): + """Validate constructor parameters.""" + if not isinstance(self.n_splits, (int, np.integer)) or self.n_splits < 2: + raise ValueError("n_splits must be an integer >= 2.") + self._n_jobs = check_n_jobs(self.n_jobs) + + @classmethod + def _get_test_params(cls, parameter_set: str = "default") -> dict[str, Any]: + """Return a small parameter set for estimator checks.""" + return {"n_splits": 2, "n_jobs": 1, "random_state": 0} diff --git a/aeon/transformations/collection/channel_selection/tests/test_msts.py b/aeon/transformations/collection/channel_selection/tests/test_msts.py new file mode 100644 index 0000000000..1fef540946 --- /dev/null +++ b/aeon/transformations/collection/channel_selection/tests/test_msts.py @@ -0,0 +1,74 @@ +"""Tests for the MSTS channel selector.""" + +import numpy as np +import pytest + +from aeon.transformations.collection.channel_selection import MSTS +from aeon.transformations.collection.channel_selection._msts import _merit_score + + +def test_msts_merit_score_matches_definition(): + """The merit score uses mean class agreement and pair agreement.""" + class_scores = np.array([0.8, 0.6, 0.2]) + pair_scores = np.eye(3) + pair_scores[0, 1] = pair_scores[1, 0] = 0.5 + + expected = 2 * np.mean([0.8, 0.6]) / np.sqrt(2 + 2 * 0.5) + assert _merit_score(class_scores, pair_scores, (0, 1)) == pytest.approx(expected) + + +def test_msts_selects_informative_channel_and_subsets_transform(): + """MSTS retains an informative channel and preserves collection shape.""" + rng = np.random.RandomState(0) + n_cases, n_channels, n_timepoints = 18, 3, 12 + y = np.array([0, 1] * (n_cases // 2)) + X = rng.normal(scale=0.05, size=(n_cases, n_channels, n_timepoints)) + X[y == 1, 0, :] += 2.0 + + selector = MSTS(n_splits=3, random_state=0) + Xt = selector.fit_transform(X, y) + + assert 0 in selector.channels_selected_ + assert Xt.shape == ( + n_cases, + len(selector.channels_selected_), + n_timepoints, + ) + np.testing.assert_array_equal(Xt, X[:, selector.channels_selected_, :]) + + +def test_msts_stores_prediction_and_agreement_metadata(): + """MSTS exposes training predictions and square agreement matrices.""" + rng = np.random.RandomState(1) + n_cases, n_channels, n_timepoints = 12, 4, 8 + X = rng.normal(size=(n_cases, n_channels, n_timepoints)) + y = np.array(["a", "b"] * (n_cases // 2)) + + selector = MSTS(n_splits=2, random_state=0).fit(X, y) + + assert selector.channel_predictions_.shape == (n_cases, n_channels) + assert selector.channel_class_scores_.shape == (n_channels,) + assert selector.channel_pair_scores_.shape == (n_channels, n_channels) + np.testing.assert_allclose(np.diag(selector.channel_pair_scores_), 1.0) + + +@pytest.mark.parametrize( + "kwargs, message", + [ + ({"n_splits": 1}, "n_splits must be an integer"), + ({"n_splits": 3, "n_jobs": 1.5}, "n_jobs"), + ], +) +def test_msts_validates_parameters(kwargs, message): + """MSTS rejects invalid cross-validation and threading parameters.""" + with pytest.raises(ValueError, match=message): + MSTS(**kwargs)._validate_parameters() + + +def test_msts_requires_enough_cases_per_class(): + """MSTS requires every class to appear in each validation-fold scheme.""" + X = np.zeros((5, 2, 8)) + y = np.array([0, 0, 0, 1, 1]) + + with pytest.raises(ValueError, match="Each class must contain"): + MSTS(n_splits=3).fit(X, y) diff --git a/docs/api_reference/transformations.md b/docs/api_reference/transformations.md index 4a2e40cd7e..71bf3e8d8e 100644 --- a/docs/api_reference/transformations.md +++ b/docs/api_reference/transformations.md @@ -47,6 +47,7 @@ all_tags_for_estimator`` function with the argument ``"transformer"``. ChannelScorer ElbowClassPairwise ElbowClassSum + MSTS RandomChannelSelector TSelect ```