Repository: sktime/sktime
Pull request #7615 (2025-01-06 20:03:50)

[ENH] Implement efficient `_evaluate_by_index` for MedianSquaredError
<!--
Welcome to sktime, and thanks for contributing!
Please have a look at our contribution guide:
https://www.sktime.net/en/latest/get_involved/contributing.html
-->

#### Reference Issues/PRs
Towards #4304 

#### What does this implement/fix? Explain your changes.
Implemented efficient `_evaluate_by_index` for MedianSquaredError

#### PR checklist
<!--
Please go through the checklist below. Please feel free to remove points if they are not applicable.
-->

##### For all contributions
- [x] The PR title starts with either [ENH], [MNT], [DOC], or [BUG]. [BUG] - bugfix, [MNT] - CI, test framework, [ENH] - adding or improving code, [DOC] - writing or improving documentation or docstrings.


satvshr: > Thanks.
> 
> Non-blocking, but it would be good to say in the docstring what `evaluate_by_index` does. See `MeanAbsoluteError` for an example.

Done!
satvshr: I  am so sorry I mistakenly pushed to the wrong branch, I sent the PR for this issue here -> #7616 
satvshr: Errors dont seem to stem from my commit, kindly take a look @fkiraly 
satvshr: > Thanks - I fixed the math in the docstring.
Thank you, looks all good!

satvshr: @fkiraly could you please merge and close this PR? Thank you!

Patch:
diff --git a/sktime/classification/sklearn/_rotation_forest.py b/sktime/classification/sklearn/_rotation_forest.py
index 645e590d7..338100a19 100644
--- a/sktime/classification/sklearn/_rotation_forest.py
+++ b/sktime/classification/sklearn/_rotation_forest.py
@@ -10,18 +10,17 @@ __all__ = ["RotationForest"]
 import time
 
 import numpy as np
-import pandas as pd
-from sklearn.base import BaseEstimator
+from sklearn.base import BaseEstimator, ClassifierMixin
 from sklearn.decomposition import PCA
 from sklearn.tree import DecisionTreeClassifier
 from sklearn.utils import check_random_state
+from sklearn.utils.validation import check_is_fitted
 
 from sktime.base._base import _clone_estimator
-from sktime.exceptions import NotFittedError
 from sktime.utils.validation import check_n_jobs
 
 
-class RotationForest(BaseEstimator):
+class RotationForest(BaseEstimator, ClassifierMixin):
     """A rotation forest (RotF) vector classifier.
 
     Implementation of the Rotation Forest classifier described in Rodriguez et al
@@ -160,18 +159,14 @@ class RotationForest(BaseEstimator):
         """
         from joblib import Parallel, delayed
 
-        if isinstance(X, np.ndarray) and len(X.shape) == 3 and X.shape[1] == 1:
-            X = np.reshape(X, (X.shape[0], -1))
-        elif isinstance(X, pd.DataFrame) and len(X.shape) == 2:
-            X = X.to_numpy()
-        elif not isinstance(X, np.ndarray) or len(X.shape) > 2:
-            raise ValueError(
-                "RotationForest is not a time series classifier. "
-                "A valid sklearn input such as a 2d numpy array is required."
-                "Sparse input formats are currently not supported."
-            )
-        X, y = self._validate_data(X=X, y=y, ensure_min_samples=2)
-
+        X, y = self._validate_data(
+            X,
+            y,
+            dtype=[np.float32, np.float64],
+            ensure_2d=True,
+            allow_nd=True,
+            force_all_finite=True,
+        )
         self._n_jobs = check_n_jobs(self.n_jobs)
 
         self.n_instances_, self.n_atts_ = X.shape
@@ -269,6 +264,7 @@ class RotationForest(BaseEstimator):
         y : array-like, shape = [n_instances]
             Predicted class labels.
         """
+        check_is_fitted(self)
         rng = check_random_state(self.random_state)
         return np.array(
             [
@@ -292,27 +288,20 @@ class RotationForest(BaseEstimator):
         """
         from joblib import Parallel, delayed
 
-        if not self._is_fitted:
-            raise NotFittedError(
-                f"This instance of {self.__class__.__name__} has not "
-                f"been fitted yet; please call `fit` first."
-            )
+        check_is_fitted(self)
 
         # treat case of single class seen in fit
         if self.n_classes_ == 1:
             return np.repeat([[1]], X.shape[0], axis=0)
 
-        if isinstance(X, np.ndarray) and len(X.shape) == 3 and X.shape[1] == 1:
-            X = np.reshape(X, (X.shape[0], -1))
-        elif isinstance(X, pd.DataFrame) and len(X.shape) == 2:
-            X = X.to_numpy()
-        elif not isinstance(X, np.ndarray) or len(X.shape) > 2:
-            raise ValueError(
-                "RotationForest is not a time series classifier. "
-                "A valid sklearn input such as a 2d numpy array is required."
-                "Sparse input formats are currently not supported."
-            )
-        X = self._validate_data(X=X, reset=False)
+        X = self._validate_data(
+            X,
+            dtype=[np.float32, np.float64],
+            reset=False,
+            ensure_2d=True,
+            allow_nd=True,
+            force_all_finite=True,
+        )
 
         # replace missing values with 0 and remove useless attributes
         X = X[:, self._useful_atts]
@@ -338,21 +327,15 @@ class RotationForest(BaseEstimator):
     def _get_train_probs(self, X, y):
         from joblib import Parallel, delayed
 
-        if not self._is_fitted:
-            raise NotFittedError(
-                f"This instance of {self.__class__.__name__} has not "
-                f"been fitted yet; please call `fit` first."
-            )
-        if isinstance(X, np.ndarray) and len(X.shape) == 3 and X.shape[1] == 1:
-            X = np.reshape(X, (X.shape[0], -1))
-        elif isinstance(X, pd.DataFrame) and len(X.shape) == 2:
-            X = X.to_numpy()
-        elif not isinstance(X, np.ndarray) or len(X.shape) > 2:
-            raise ValueError(
-                "RotationForest is not a time series classifier. "
-                "A valid sklearn input such as a 2d numpy array is required."
-                "Sparse input formats are currently not supported."
-            )
+        check_is_fitted(self)
+        X = self._validate_data(
+            X,
+            dtype=[np.float32, np.float64],
+            reset=False,
+            ensure_2d=True,
+            allow_nd=True,
+            force_all_finite=True,
+        )
         X = self._validate_data(X=X, reset=False)
 
         # handle the single-class-label case
diff --git a/sktime/performance_metrics/forecasting/_classes.py b/sktime/performance_metrics/forecasting/_classes.py
index c5986cb0f..b23d1d8d2 100644
--- a/sktime/performance_metrics/forecasting/_classes.py
+++ b/sktime/performance_metrics/forecasting/_classes.py
@@ -2139,22 +2139,48 @@ class MeanSquaredErrorPercentage(BaseForecastingErrorMetricFunc):
 
 
 class MedianSquaredError(BaseForecastingErrorMetricFunc):
-    """Median squared error (MdSE) or root median squared error (RMdSE).
+    r"""Median squared error (MdSE) or root median squared error (RMdSE).
 
     If ``square_root`` is False then calculates MdSE and if ``square_root`` is True
     then RMdSE is calculated. Both MdSE and RMdSE return non-negative floating
     point. The best value is 0.0.
 
+    The Median Squared Error (MdSE) is calculated as:
+
+    .. math::
+        \\text{MdSE} = \text{median}\left((y_{i} - \\widehat{y}_{i})^2\right)_{i=1}^{n}
+
+    where:
+    - \\( y_i \\) are the true values,
+    - \\( \\hat{y}_i \\) are the predicted values,
+    - \\( n \\) is the number of observations.
+
+    If ``square_root`` is set to True,
+    the Root Median Squared Error Percentage (RMdSE) is computed:
+
+    .. math::
+        \\text{RMdSE} = \\sqrt{ \\text{MdSE} }
+
     Like MSE, MdSE is measured in squared units of the input data. RMdSE is
-    on the same scale as the input data like RMSE. Because MdSE and RMdSE
-    square the forecast error rather than taking the absolute value, they
-    penalize large errors more than MAE or MdAE.
+    on the same scale as the input data, like RMSE.
 
     Taking the median instead of the mean of the squared errors makes
     this metric more robust to error outliers relative to a meean based metric
     since the median tends to be a more robust measure of central tendency in
     the presence of outliers.
 
+    Because the median commutes with monotonic transformations, MdSE and RMdSE
+    do not penalize large errors more than the MdAE.
+
+    ``evaluate_by_index`` returns, at a time index :math:`t_i`:
+
+    * if ``square_root`` is False, the squared error at that time index,
+        :math:`(y_i - \widehat{y}_i)^2`,
+    * if ``square_root`` is True, the absolute error at that time index,
+      the absolute error at that time index, :math:`|y_i - \widehat{y}_i|`,
+
+    for all time indices :math:`t_1, \dots, t_n` in the input.
+
     Parameters
     ----------
     square_root : bool, default = False
@@ -2234,6 +2260,53 @@ class MedianSquaredError(BaseForecastingErrorMetricFunc):
         self.square_root = square_root
         super().__init__(multioutput=multioutput, multilevel=multilevel)
 
+    def _evaluate_by_index(self, y_true, y_pred, **kwargs):
+        """Return the metric evaluated at each time point.
+
+        private _evaluate_by_index containing core logic, called from evaluate_by_index
+
+        Parameters
+        ----------
+        y_true : time series in sktime compatible pandas based data container format
+            Ground truth (correct) target values
+            y can be in one of the following formats:
+            Series scitype: pd.DataFrame
+            Panel scitype: pd.DataFrame with 2-level row MultiIndex
+            Hierarchical scitype: pd.DataFrame with 3 or more level row MultiIndex
+        y_pred : time series in sktime compatible data container format
+            Forecasted values to evaluate
+            must be of same format as y_true, same indices and columns if indexed
+
+        Returns
+        -------
+        loss : pd.Series or pd.DataFrame
+            Calculated metric, by time point.
+            pd.Series if self.multioutput="uniform_average" or array-like
+                index is equal to index of y_true
+                entry at index i is metric at time i, averaged over variables
+            pd.DataFrame if self.multioutput="raw_values"
+                index and columns equal to those of y_true
+                i,j-th entry is metric at time i, at variable j
+        """
+        multioutput = self.multioutput
+
+        raw_values = (y_true - y_pred) ** 2
+
+        if self.square_root:
+            raw_values = np.sqrt(raw_values)
+
+        raw_values = self._get_weighted_df(raw_values, **kwargs)
+
+        if isinstance(multioutput, str):
+            if multioutput == "raw_values":
+                return raw_values
+
+            if multioutput == "uniform_average":
+                return raw_values.median(axis=1)
+
+        # else, we expect multioutput to be array-like
+        return raw_values.dot(multioutput)
+
     @classmethod
     def get_test_params(cls, parameter_set="default"):
         """Return testing parameter settings for the estimator.
diff --git a/sktime/performance_metrics/forecasting/_functions.py b/sktime/performance_metrics/forecasting/_functions.py
index 28fa97357..b109c833d 100644
--- a/sktime/performance_metrics/forecasting/_functions.py
+++ b/sktime/performance_metrics/forecasting/_functions.py
@@ -1122,7 +1122,7 @@ def median_squared_error(
     penalize large errors more than MAE or MdAE.
 
     Taking the median instead of the mean of the squared errors makes
-    this metric more robust to error outliers relative to a meean based metric
+    this metric more robust to error outliers relative to a mean based metric
     since the median tends to be a more robust measure of central tendency in
     the presence of outliers.
 


Test patch:
diff --git a/sktime/classification/sklearn/tests/test_rotation_forest.py b/sktime/classification/sklearn/tests/test_rotation_forest.py
index 9186e71e6..1ba351de6 100644
--- a/sktime/classification/sklearn/tests/test_rotation_forest.py
+++ b/sktime/classification/sklearn/tests/test_rotation_forest.py
@@ -2,12 +2,23 @@
 
 import numpy as np
 import pytest
+from sklearn.utils.estimator_checks import parametrize_with_checks
 
 from sktime.classification.sklearn import RotationForest
 from sktime.datasets import load_unit_test
 from sktime.tests.test_switch import run_test_for_class
 
 
+@pytest.mark.skipif(
+    not run_test_for_class(RotationForest),
+    reason="run test if softdeps are present and incrementally (if requested)",
+)
+@parametrize_with_checks([RotationForest()])
+def test_sklearn_compatible_estimator(estimator, check):
+    """Run sklearn estimator compatibility checks."""
+    check(estimator)
+
+
 @pytest.mark.skipif(
     not run_test_for_class(RotationForest),
     reason="run test only if softdeps are present and incrementally (if requested)",