Repository: alan-turing-institute/sktime
Pull request #2601 (2022-05-02 20:37:47)

[ENH] metrics rework part IV - hierarchical metrics
This PR continues rework of the current metrics interface.

This PR uses the tests in https://github.com/alan-turing-institute/sktime/pull/2496 and builds upon #2500 and #2502

The following changes are made relative to #2502:

* extending `BaseForecastingErrorMetricFunc` to accept hierarchical input
* enable averaging options via the `multilevel` parameter


Patch:
diff --git a/docs/source/developer_guide/coding_standards.rst b/docs/source/developer_guide/coding_standards.rst
index 6a4e1f775..fe1b8a4e3 100644
--- a/docs/source/developer_guide/coding_standards.rst
+++ b/docs/source/developer_guide/coding_standards.rst
@@ -12,19 +12,22 @@ Coding style
 
 In coding, we follow:
 
-* the `PEP8 <https://www.python.org/dev/peps/pep-0008/>`__
-coding guidelines. A good example can be found
-`here <https://gist.github.com/nateGeorge/5455d2c57fb33c1ae04706f2dc4fee01>`__.
+*  the `PEP8 <https://www.python.org/dev/peps/pep-0008/>`__ coding guidelines. A good example can be found `here <https://gist.github.com/nateGeorge/5455d2c57fb33c1ae04706f2dc4fee01>`__.
+
 * code formatting according to ``black``, ``flake8``, ``isort``, ``numpydoc``
 
 Code formatting and linting
 ---------------------------
 
 We adhere to the following code formatting standards:
+
 * `black <https://black.readthedocs.io/en/stable/>`__ with default settings
+
 * `flake8 <https://flake8.pycqa.org/en/latest/>`__ with a ``max_line_length=88`` and some exceptions as per ``setup.cfg``
+
 * ``isort`` with default settings
-* ``numpydoc``to enforce numpy docstring standard <https://numpydoc.readthedocs.io/en/latest/format.html#docstring-standard>`_, along with sktime specific conventions described in our :ref:`developer_guide`'s :ref:`documentation section <developer_guide_documentation>`.
+
+* ``numpydoc`` to enforce numpy `docstring standard <https://numpydoc.readthedocs.io/en/latest/format.html#docstring-standard>`_ , along with sktime specific conventions described in our :ref:`developer_guide`'s :ref:`documentation section <developer_guide_documentation>`.
 
 This is enforced through our CI/CD workflows via `pre-commit <https://pre-commit.com/>`_.
 
diff --git a/extension_templates/classification.py b/extension_templates/classification.py
index f10b070d8..f24951996 100644
--- a/extension_templates/classification.py
+++ b/extension_templates/classification.py
@@ -118,20 +118,26 @@ class MyTimeSeriesClassifier(BaseClassifier):
     def _fit(self, X, y):
         """Fit time series classifier to training data.
 
-        core logic
+        private _fit containing the core logic, called from fit
+
+        Writes to self:
+            Sets fitted model attributes ending in "_".
 
         Parameters
         ----------
-        X : Training data of type self.get_tag("X_inner_mtype")
-        y : array-like, shape = [n_instances] - the class labels
+        X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
+            if self.get_tag("X_inner_mtype") = "numpy3D":
+                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "nested_univ":
+                pd.DataFrame with each column a dimension, each cell a pd.Series
+            for list of other mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
+        y : 1D np.array of int, of shape [n_instances] - class labels for fitting
+            indices correspond to instance indices in X
 
         Returns
         -------
-        self : reference to self.
-
-        State change
-        ------------
-        creates fitted model (attributes ending in "_")
+        self : Reference to self.
         """
 
         # implement here
@@ -147,15 +153,28 @@ class MyTimeSeriesClassifier(BaseClassifier):
     def _predict(self, X) -> np.ndarray:
         """Predict labels for sequences in X.
 
-        core logic
+        private _predict containing the core logic, called from predict
+
+        State required:
+            Requires state to be "fitted".
+
+        Accesses in self:
+            Fitted model attributes ending in "_"
 
         Parameters
         ----------
-        X : data not used in training, of type self.get_tag("X_inner_mtype")
+        X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
+            if self.get_tag("X_inner_mtype") = "numpy3D":
+                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "nested_univ":
+                pd.DataFrame with each column a dimension, each cell a pd.Series
+            for list of other mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
 
         Returns
         -------
-        y : predictions of labels for X, np.ndarray
+        y : 1D np.array of int, of shape [n_instances] - predicted class labels
+            indices correspond to instance indices in X
         """
 
         # implement here
@@ -167,17 +186,30 @@ class MyTimeSeriesClassifier(BaseClassifier):
     def _predict_proba(self, X) -> np.ndarray:
         """Predicts labels probabilities for sequences in X.
 
-        Default behaviour is to call _predict and set the predicted class probability
-        to 1, other class probabilities to 0. Override if better estimates are
-        obtainable.
+        private _predict_proba containing the core logic, called from predict_proba
+
+        State required:
+            Requires state to be "fitted".
+
+        Accesses in self:
+            Fitted model attributes ending in "_"
 
         Parameters
         ----------
-        X : data to predict y with, of type self.get_tag("X_inner_mtype")
+        X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
+            if self.get_tag("X_inner_mtype") = "numpy3D":
+                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "nested_univ":
+                pd.DataFrame with each column a dimension, each cell a pd.Series
+            for list of other mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
 
         Returns
         -------
-        y : predictions of probabilities for class values of X, np.ndarray
+        y : 2D array of shape [n_instances, n_classes] - predicted class probabilities
+            1st dimension indices correspond to instance indices in X
+            2nd dimension indices correspond to possible labels (integers)
+            (i, j)-th entry is predictive probability that i-th instance is of class j
         """
 
         # implement here
diff --git a/sktime/classification/_delegate.py b/sktime/classification/_delegate.py
new file mode 100644
index 000000000..cef272bd0
--- /dev/null
+++ b/sktime/classification/_delegate.py
@@ -0,0 +1,124 @@
+# -*- coding: utf-8 -*-
+"""Delegator mixin that delegates all methods to wrapped classifiers.
+
+Useful for building estimators where all but one or a few methods are delegated.
+For that purpose, inherit from this estimator and then override only the methods
+    that are not delegated.
+"""
+# copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
+
+__author__ = ["fkiraly"]
+__all__ = ["_DelegatedClassifier"]
+
+from sktime.classification.base import BaseClassifier
+
+
+class _DelegatedClassifier(BaseClassifier):
+    """Delegator mixin that delegates all methods to wrapped classifier.
+
+    Delegates inner classifier methods to a wrapped estimator.
+        Wrapped estimator is value of attribute with name self._delegate_name.
+        By default, this is "estimator_", i.e., delegates to self.estimator_
+        To override delegation, override _delegate_name attribute in child class.
+
+    Delegates the following inner underscore methods:
+        _fit, _predict, _predict_proba
+
+    Does NOT delegate get_params, set_params.
+        get_params, set_params will hence use one additional nesting level by default.
+
+    Does NOT delegate or copy tags, this should be done in a child class if required.
+    """
+
+    _delegate_name = "estimator_"
+
+    def _get_delegate(self):
+        return getattr(self, self._delegate_name)
+
+    def _fit(self, X, y):
+        """Fit time series classifier to training data.
+
+        private _fit containing the core logic, called from fit
+
+        Writes to self:
+            Sets fitted model attributes ending in "_".
+
+        Parameters
+        ----------
+        X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
+            if self.get_tag("X_inner_mtype") = "numpy3D":
+                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "nested_univ":
+                pd.DataFrame with each column a dimension, each cell a pd.Series
+            for list of other mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
+        y : 1D np.array of int, of shape [n_instances] - class labels for fitting
+            indices correspond to instance indices in X
+
+        Returns
+        -------
+        self : Reference to self.
+        """
+        estimator = self._get_delegate()
+        estimator.fit(X=X, y=y)
+        return self
+
+    def _predict(self, X):
+        """Predict labels for sequences in X.
+
+        private _predict containing the core logic, called from predict
+
+        State required:
+            Requires state to be "fitted".
+
+        Accesses in self:
+            Fitted model attributes ending in "_"
+
+        Parameters
+        ----------
+        X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
+            if self.get_tag("X_inner_mtype") = "numpy3D":
+                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "nested_univ":
+                pd.DataFrame with each column a dimension, each cell a pd.Series
+            for list of other mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
+
+        Returns
+        -------
+        y : 1D np.array of int, of shape [n_instances] - predicted class labels
+            indices correspond to instance indices in X
+        """
+        estimator = self._get_delegate()
+        return estimator.predict(X=X)
+
+    def _predict_proba(self, X):
+        """Predicts labels probabilities for sequences in X.
+
+        private _predict_proba containing the core logic, called from predict_proba
+
+        State required:
+            Requires state to be "fitted".
+
+        Accesses in self:
+            Fitted model attributes ending in "_"
+
+        Parameters
+        ----------
+        X : guaranteed to be of a type in self.get_tag("X_inner_mtype")
+            if self.get_tag("X_inner_mtype") = "numpy3D":
+                3D np.ndarray of shape = [n_instances, n_dimensions, series_length]
+            if self.get_tag("X_inner_mtype") = "nested_univ":
+                pd.DataFrame with each column a dimension, each cell a pd.Series
+            for list of other mtypes, see datatypes.SCITYPE_REGISTER
+            for specifications, see examples/AA_datatypes_and_datasets.ipynb
+
+        Returns
+        -------
+        y : 2D array of shape [n_instances, n_classes] - predicted class probabilities
+            1st dimension indices correspond to instance indices in X
+            2nd dimension indices correspond to possible labels (integers)
+            (i, j)-th entry is predictive probability that i-th instance is of class j
+        """
+        estimator = self._get_delegate()
+        return estimator.predict_proba(X=X)
diff --git a/sktime/classification/base.py b/sktime/classification/base.py
index 208f41d4b..25fd8f2c1 100644
--- a/sktime/classification/base.py
+++ b/sktime/classification/base.py
@@ -67,11 +67,13 @@ class BaseClassifier(BaseEstimator, ABC):
     }
 
     def __init__(self):
-        self.classes_ = []
-        self.n_classes_ = 0
-        self.fit_time_ = 0
+        # reserved attributes written to in fit
+        self.classes_ = []  # classes seen in y, unique labels
+        self.n_classes_ = 0  # number of unique classes in y
+        self.fit_time_ = 0  # time elapsed in last fit call
         self._class_dictionary = {}
         self._threads_to_use = 1
+        self._X_metadata = []  # metadata/properties of X seen in fit
 
         # required for compatability with some sklearn interfaces
         # i.e. CalibratedClassifierCV
@@ -170,6 +172,7 @@ class BaseClassifier(BaseEstimator, ABC):
 
         self.classes_ = np.unique(y)
         self.n_classes_ = self.classes_.shape[0]
+        self._X_metadata = X_metadata
         self._class_dictionary = {}
         for index, class_val in enumerate(self.classes_):
             self._class_dictionary[class_val] = index
@@ -524,9 +527,10 @@ class BaseClassifier(BaseEstimator, ABC):
             2nd dimension indices correspond to possible labels (integers)
             (i, j)-th entry is predictive probability that i-th instance is of class j
         """
-        dists = np.zeros((X.shape[0], self.n_classes_))
         preds = self._predict(X)
-        for i in range(0, X.shape[0]):
+        n_pred = len(preds)
+        dists = np.zeros((n_pred, self.n_classes_))
+        for i in range(n_pred):
             dists[i, self._class_dictionary[preds[i]]] = 1
 
         return dists
diff --git a/sktime/classification/compose/_pipeline.py b/sktime/classification/compose/_pipeline.py
index ad3bd754c..b0d6f50d2 100644
--- a/sktime/classification/compose/_pipeline.py
+++ b/sktime/classification/compose/_pipeline.py
@@ -423,7 +423,7 @@ class SklearnClassifierPipeline(ClassifierPipeline):
 
         # can handle multivariate iff all transformers can
         # sklearn transformers always support multivariate
-        multivariate = self.transformers_.get_tag("univariate-only", True)
+        multivariate = not self.transformers_.get_tag("univariate-only", True)
         # can handle missing values iff transformer chain removes missing data
         # sklearn classifiers might be able to handle missing data (but no tag there)
         # so better set the tag liberally
@@ -546,8 +546,12 @@ class SklearnClassifierPipeline(ClassifierPipeline):
         y : predictions of probabilities for class values of X, np.ndarray
         """
         Xt = self.transformers_.transform(X)
-        Xt_sklearn = self._convert_X_to_sklearn(Xt)
-        return self.classifier_.predict_proba(Xt_sklearn)
+        if hasattr(self.classifier_, "predict_proba"):
+            Xt_sklearn = self._convert_X_to_sklearn(Xt)
+            return self.classifier_.predict_proba(Xt_sklearn)
+        else:
+            # if sklearn classifier does not have predict_proba
+            return BaseClassifier._predict_proba(self, X)
 
     def set_params(self, **kwargs):
         """Set the parameters of estimator in `transformers`.
diff --git a/sktime/classification/kernel_based/_rocket_classifier.py b/sktime/classification/kernel_based/_rocket_classifier.py
index d2ee43c4b..bf990b0f9 100644
--- a/sktime/classification/kernel_based/_rocket_classifier.py
+++ b/sktime/classification/kernel_based/_rocket_classifier.py
@@ -4,15 +4,15 @@
 Pipeline classifier using the ROCKET transformer and RidgeClassifierCV estimator.
 """
 
-__author__ = ["MatthewMiddlehurst", "victordremov"]
+__author__ = ["MatthewMiddlehurst", "victordremov", "fkiraly"]
 __all__ = ["RocketClassifier"]
 
 import numpy as np
 from sklearn.linear_model import RidgeClassifierCV
-from sklearn.pipeline import make_pipeline
 from sklearn.preprocessing import StandardScaler
 
-from sktime.classification.base import BaseClassifier
+from sktime.classification._delegate import _DelegatedClassifier
+from sktime.pipeline import make_pipeline
 from sktime.transformations.panel.rocket import (
     MiniRocket,
     MiniRocketMultivariate,
@@ -22,23 +22,46 @@ from sktime.transformations.panel.rocket import (
 )
 
 
-class RocketClassifier(BaseClassifier):
+class RocketClassifier(_DelegatedClassifier):
     """Classifier wrapped for the Rocket transformer using RidgeClassifierCV.
 
     This classifier simply transforms the input data using the Rocket [1]_
     transformer and builds a RidgeClassifierCV estimator using the transformed data.
 
+    Shorthand for the pipeline
+    `rocket * StandardScaler(with_mean=False) * RidgeClassifierCV(alphas)`
+    where `alphas = np.logspace(-3, 3, 10)`, and
+    where `rocket` depends on params `rocket_transform`, `use_multivariate` as follows:
+
+        | rocket_transform | `use_multivariate` | rocket (class)          |
+        |------------------|--------------------|-------------------------|
+        | "rocket"         | any                | Rocket                  |
+        | "minirocket"     | "yes               | MiniRocketMultivariate  |
+        | "minirocket"     | "no"               | MiniRocket              |
+        | "multirocket"    | "yes"              | MultiRocketMultivariate |
+        | "multirocket"    | "no"               | MultiRocket             |
+
+    classes are sktime classes, other parameters are passed on to the rocket class.
+
+    To build other classifiers with rocket transformers, use `make_pipeline` or the
+    pipeline dunder `*`, and different transformers/classifiers in combination.
+
     Parameters
     ----------
-    num_kernels : int, default=10,000
+    num_kernels : int, optional, default=10,000
         The number of kernels the for Rocket transform.
-    rocket_transform : str, default="rocket"
+    rocket_transform : str, optional, default="rocket"
         The type of Rocket transformer to use.
-        Valid inputs = ["rocket","minirocket","multirocket"]
-    max_dilations_per_kernel : int, default=32
+        Valid inputs = ["rocket", "minirocket", "multirocket"]
+    max_dilations_per_kernel : int, optional, default=32
         MiniRocket and MultiRocket only. The maximum number of dilations per kernel.
-    n_features_per_kernel : int, default=4
+    n_features_per_kernel : int, optional, default=4
         MultiRocket only. The number of features per kernel.
+    use_multivariate : str, ["auto", "yes", "no"], optional, default="auto"
+        whether to use multivariate rocket transforms or univariate ones
+        "auto" = multivariate iff data seen in fit is multivariate, otherwise univariate
+        "yes" = always uses multivariate transformers, native multi/univariate
+        "no" = always univariate transformers, multivariate by framework vectorization
     n_jobs : int, default=1
         The number of jobs to run in parallel for both `fit` and `predict`.
         ``-1`` means using all processors.
@@ -51,6 +74,8 @@ class RocketClassifier(BaseClassifier):
         The number of classes.
     classes_ : list
         The classes labels.
+    estimator_ : ClassifierPipeline
+        RocketClassifier as a ClassifierPipeline, fitted to data internally
 
     See Also
     --------
@@ -86,12 +111,17 @@ class RocketClassifier(BaseClassifier):
         "classifier_type": "kernel",
     }
 
+    # valid rocket strings for input validity checking
+    VALID_ROCKET_STRINGS = ["rocket", "minirocket", "multirocket"]
+    VALID_MULTIVAR_VALUES = ["auto", "yes", "no"]
+
     def __init__(
         self,
         num_kernels=10000,
         rocket_transform="rocket",
         max_dilations_per_kernel=32,
         n_features_per_kernel=4,
+        use_multivariate="auto",
         n_jobs=1,
         random_state=None,
     ):
@@ -99,120 +129,77 @@ class RocketClassifier(BaseClassifier):
         self.rocket_transform = rocket_transform
         self.max_dilations_per_kernel = max_dilations_per_kernel
         self.n_features_per_kernel = n_features_per_kernel
+        self.use_multivariate = use_multivariate
 
         self.n_jobs = n_jobs
         self.random_state = random_state
 
-        self._pipeline = None
-
         super(RocketClassifier, self).__init__()
 
-    def _fit(self, X, y):
-        """Build a pipeline containing the Rocket transformer and RidgeClassifierCV.
+        if use_multivariate not in self.VALID_MULTIVAR_VALUES:
+            raise ValueError(
+                f"Invalid use_multivariate value, must be one of "
+                f"{self.VALID_MULTIVAR_VALUES}, but found {use_multivariate}"
+            )
 
-        Parameters
-        ----------
-        X : 3D np.array of shape = [n_instances, n_dimensions, series_length]
-            The training data.
-        y : array-like, shape = [n_instances]
-            The class labels.
+        common_params = {
+            "num_kernels": self.num_kernels,
+            "random_state": self.random_state,
+            "max_dilations_per_kernel": self.max_dilations_per_kernel,
+            "n_jobs": self._threads_to_use,
+        }
 
-        Returns
-        -------
-        self :
-            Reference to self.
+        if rocket_transform == "rocket":
+            del common_params["max_dilations_per_kernel"]
+            univar_rocket = Rocket(**common_params)
+            multivar_rocket = univar_rocket
 
-        Notes
-        -----
-        Changes state by creating a fitted model that updates attributes
-        ending in "_" and sets is_fitted flag to True.
-        """
-        _, n_dims, _ = X.shape
+        elif rocket_transform == "minirocket":
+            multivar_rocket = MiniRocketMultivariate(**common_params)
+            univar_rocket = MiniRocket(**common_params)
 
-        if self.rocket_transform == "rocket":
-            rocket = Rocket(
-                num_kernels=self.num_kernels,
-                random_state=self.random_state,
-                n_jobs=self._threads_to_use,
-            )
-        elif self.rocket_transform == "minirocket":
-            if n_dims > 1:
-                rocket = MiniRocketMultivariate(
-                    num_kernels=self.num_kernels,
-                    max_dilations_per_kernel=self.max_dilations_per_kernel,
-                    random_state=self.random_state,
-                    n_jobs=self._threads_to_use,
-                )
-            else:
-                rocket = MiniRocket(
-                    num_kernels=self.num_kernels,
-                    max_dilations_per_kernel=self.max_dilations_per_kernel,
-                    random_state=self.random_state,
-                    n_jobs=self._threads_to_use,
-                )
         elif self.rocket_transform == "multirocket":
-            if n_dims > 1:
-                rocket = MultiRocketMultivariate(
-                    num_kernels=self.num_kernels,
-                    max_dilations_per_kernel=self.max_dilations_per_kernel,
-                    n_features_per_kernel=self.n_features_per_kernel,
-                    random_state=self.random_state,
-                    n_jobs=self._threads_to_use,
-                )
-            else:
-                rocket = MultiRocket(
-                    num_kernels=self.num_kernels,
-                    max_dilations_per_kernel=self.max_dilations_per_kernel,
-                    n_features_per_kernel=self.n_features_per_kernel,
-                    random_state=self.random_state,
-                    n_jobs=self._threads_to_use,
-                )
+            common_params["n_features_per_kernel"] = self.n_features_per_kernel
+            multivar_rocket = MultiRocketMultivariate(**common_params)
+            univar_rocket = MultiRocket(**common_params)
+
         else:
-            raise ValueError(f"Invalid Rocket transformer: {self.rocket_transform}")
+            raise ValueError(
+                f"Invalid rocket_transform string, must be one of "
+                f"{self.VALID_ROCKET_STRINGS}, but found {rocket_transform}"
+            )
 
-        self._pipeline = rocket_pipeline = make_pipeline(
-            rocket,
+        self.multivar_rocket_ = make_pipeline(
+            multivar_rocket,
+            StandardScaler(with_mean=False),
+            RidgeClassifierCV(alphas=np.logspace(-3, 3, 10)),
+        )
+        self.univar_rocket_ = make_pipeline(
+            univar_rocket,
             StandardScaler(with_mean=False),
             RidgeClassifierCV(alphas=np.logspace(-3, 3, 10)),
         )
-        rocket_pipeline.fit(X, y)
-
-        return self
-
-    def _predict(self, X) -> np.ndarray:
-        """Predicts labels for sequences in X.
-
-        Parameters
-        ----------
-        X : 3D np.array of shape = [n_instances, n_dimensions, series_length]
-            The data to make predictions for.
 
-        Returns
-        -------
-        y : array-like, shape = [n_instances]
-            Predicted class labels.
-        """
-        return self._pipeline.predict(X)
+        if not use_multivariate:
+            self.set_tags(**{"capability:multivariate": False})
 
-    def _predict_proba(self, X) -> np.ndarray:
-        """Predicts labels probabilities for sequences in X.
+    @property
+    def estimator_(self):
+        """Shorthand for the internal estimator that is fitted."""
+        return self._get_delegate()
 
-        Parameters
-        ----------
-        X : 3D np.array of shape = [n_instances, n_dimensions, series_length]
-            The data to make predict probabilities for.
+    def _get_delegate(self):
+        use_multivariate = self.use_multivariate
+        if use_multivariate == "auto":
+            code_dict = {True: "yes", False: "no"}
+            use_multivariate = code_dict[not self._X_metadata["is_univariate"]]
 
-        Returns
-        -------
-        y : array-like, shape = [n_instances, n_classes_]
-            Predicted probabilities using the ordering in classes_.
-        """
-        dists = np.zeros((X.shape[0], self.n_classes_))
-        preds = self._pipeline.predict(X)
-        for i in range(0, X.shape[0]):
-            dists[i, np.where(self.classes_ == preds[i])] = 1
+        if use_multivariate == "yes":
+            delegate = self.multivar_rocket_
+        else:
+            delegate = self.univar_rocket_
 
-        return dists
+        return delegate
 
     @classmethod
     def get_test_params(cls, parameter_set="default"):
diff --git a/sktime/performance_metrics/forecasting/_classes.py b/sktime/performance_metrics/forecasting/_classes.py
index 3eea4e338..b50effdbd 100644
--- a/sktime/performance_metrics/forecasting/_classes.py
+++ b/sktime/performance_metrics/forecasting/_classes.py
@@ -7,12 +7,14 @@ Classes named as ``*Score`` return a value to maximize: the higher the better.
 Classes named as ``*Error`` or ``*Loss`` return a value to minimize:
 the lower the better.
 """
+from copy import deepcopy
 from warnings import warn
 
+import numpy as np
 import pandas as pd
 from sklearn.utils import check_array
 
-from sktime.datatypes import check_is_scitype, convert_to
+from sktime.datatypes import VectorizedDF, check_is_scitype, convert_to
 from sktime.performance_metrics.base import BaseMetric
 from sktime.performance_metrics.forecasting._functions import (
     geometric_mean_absolute_error,
@@ -65,6 +67,23 @@ __all__ = [
 ]
 
 
+def _coerce_to_scalar(obj):
+    """Coerce obj to scalar, from polymorphic input scalar or pandas."""
+    if isinstance(obj, pd.DataFrame):
+        assert len(obj) == 1
+        assert len(obj.columns) == 1
+        return obj.iloc[0, 0]
+    if isinstance(obj, pd.Series):
+        assert len(obj) == 1
+        return obj.iloc[0]
+    return obj
+
+
+def _coerce_to_df(obj):
+    """Coerce to pd.DataFrame, from polymorphic input scalar or pandas."""
+    return pd.DataFrame(obj)
+
+
 class BaseForecastingErrorMetric(BaseMetric):
     """Base class for defining forecasting error metrics in sktime.
 
@@ -83,9 +102,10 @@ class BaseForecastingErrorMetric(BaseMetric):
         If array-like, errors are weighted averaged across variables, values as weights.
         If 'raw_values', does not average errors across variables, columns are retained.
 
-    multilevel : {'raw_values', 'uniform_average'}
+    multilevel : {'raw_values', 'uniform_average', 'uniform_average_time'}
         Defines how to aggregate metric for hierarchical data (with levels).
         If 'uniform_average' (default), errors are mean-averaged across levels.
+        If 'uniform_average_time', errors are mean-averaged across rows.
         If 'raw_values', does not average errors across levels, hierarchy is retained.
     """
 
@@ -121,12 +141,18 @@ class BaseForecastingErrorMetric(BaseMetric):
 
         Returns
         -------
-        loss : float or np.ndarray
+        loss : float, np.ndarray, or pd.DataFrame
             Calculated metric, averaged or by variable.
             float if self.multioutput="uniform_average" or array-like
-                value is metric averaged over variables (see class docstring)
+                and self.multilevel="uniform_average" or "uniform_average_time"
+                value is metric averaged over variables and levels (see class docstring)
             np.ndarray of shape (y_true.columns,) if self.multioutput="raw_values"
+                and self.multilevel="uniform_average" or "uniform_average_time"
                 i-th entry is metric calculated for i-th variable
+            pd.DataFrame if self.multilevel=raw.values
+                of shape (n_levels, ) if self.multioutput = "uniform_average" or array
+                of shape (n_levels, y_true.columns) if self.multioutput="raw_values"
+                metric is applied per level, row averaging (yes/no) as in multioutput
         """
         return self.evaluate(y_true, y_pred, **kwargs)
 
@@ -148,12 +174,18 @@ class BaseForecastingErrorMetric(BaseMetric):
 
         Returns
         -------
-        loss : float or np.ndarray
+        loss : float, np.ndarray, or pd.DataFrame
             Calculated metric, averaged or by variable.
             float if self.multioutput="uniform_average" or array-like
-                value is metric averaged over variables (see class docstring)
+                and self.multilevel="uniform_average" or "uniform_average_time"
+                value is metric averaged over variables and levels (see class docstring)
             np.ndarray of shape (y_true.columns,) if self.multioutput="raw_values"
+                and self.multilevel="uniform_average" or "uniform_average_time"
                 i-th entry is metric calculated for i-th variable
+            pd.DataFrame if self.multilevel=raw.values
+                of shape (n_levels, ) if self.multioutput = "uniform_average" or array
+                of shape (n_levels, y_true.columns) if self.multioutput="raw_values"
+                metric is applied per level, row averaging (yes/no) as in multioutput
         """
         multioutput = self.multioutput
         multilevel = self.multilevel
@@ -161,8 +193,25 @@ class BaseForecastingErrorMetric(BaseMetric):
         y_true_inner, y_pred_inner, multioutput, multilevel, kwargs = self._check_ys(
             y_true, y_pred, multioutput, multilevel, **kwargs
         )
-        # pass to inner function
-        out_df = self._evaluate(y_true=y_true_inner, y_pred=y_pred_inner, **kwargs)
+
+        requires_vectorization = isinstance(y_true_inner, VectorizedDF)
+        if not requires_vectorization:
+            # pass to inner function
+            out_df = self._evaluate(y_true=y_true_inner, y_pred=y_pred_inner, **kwargs)
+        else:
+            out_df = self._evaluate_vectorized(
+                y_true=y_true_inner, y_pred=y_pred_inner, **kwargs
+            )
+            if multilevel == "uniform_average":
+                out_df = out_df.mean(axis=0)
+                # if level is averaged, but not variables, return numpy
+                if multioutput == "raw_values":
+                    out_df = out_df.values
+
+        if multilevel == "uniform_average" and multioutput == "uniform_average":
+            out_df = _coerce_to_scalar(out_df)
+        if multilevel == "raw_values":
+            out_df = _coerce_to_df(out_df)
 
         return out_df
 
@@ -203,6 +252,41 @@ class BaseForecastingErrorMetric(BaseMetric):
         except RecursionError:
             RecursionError("Must implement one of _evaluate or _evaluate_by_index")
 
+    def _evaluate_vectorized(self, y_true, y_pred, **kwargs):
+        """Vectorized version of _evaluate.
+
+        Runs _evaluate for all instances in y_true, y_pred,
+        and returns results in a hierarchical pandas.DataFrame.
+
+        Parameters
+        ----------
+        y_true : pandas.DataFrame with MultiIndex, last level time-like
+        y_pred : pandas.DataFrame with MultiIndex, last level time-like
+        non-time-like instanceso of y_true, y_pred must be identical
+        """
+        kwargsi = deepcopy(kwargs)
+        n_batches = len(y_true)
+        res = []
+        for i in range(n_batches):
+            if "y_train" in kwargs:
+                kwargsi["y_train"] = kwargs["y_train"][i]
+            if "y_pred_benchmark" in kwargs:
+                kwargsi["y_pred_benchmark"] = kwargs["y_pred_benchmark"][i]
+            resi = self._evaluate(y_true=y_true[i], y_pred=y_pred[i], **kwargsi)
+            if isinstance(resi, float):
+                resi = pd.Series(resi)
+            if self.multioutput == "raw_values":
+                assert isinstance(resi, np.ndarray)
+                df = pd.DataFrame(columns=y_true.X.columns)
+                df.loc[0] = resi
+                resi = df
+            res += [resi]
+        out_df = y_true.reconstruct(res)
+        if out_df.index.nlevels == y_true.X.index.nlevels:
+            out_df.index = out_df.index.droplevel(-1)
+
+        return out_df
+
     def evaluate_by_index(self, y_true, y_pred, **kwargs):
         """Return the metric evaluated at each time point.
 
@@ -296,6 +380,15 @@ class BaseForecastingErrorMetric(BaseMetric):
 
     def _check_consistent_input(self, y_true, y_pred, multioutput, multilevel):
 
+        y_true_orig = y_true
+        y_pred_orig = y_pred
+
+        # unwrap y_true, y_pred, if wrapped in VectorizedDF
+        if isinstance(y_true, VectorizedDF):
+            y_true = y_true.X
+        if isinstance(y_pred, VectorizedDF):
+            y_pred = y_pred.X
+
         # check row and column indices if y_true vs y_pred
         same_rows = y_true.index.equals(y_pred.index)
         same_row_num = len(y_true.index) == len(y_pred.index)
@@ -315,17 +408,22 @@ class BaseForecastingErrorMetric(BaseMetric):
                 "This may indicate incorrect objects passed to the metric. "
                 "Indices of y_true will be used for y_pred."
             )
-            y_pred = y_pred.copy()
-            y_pred.index = y_true.index
+            y_pred_orig = y_pred_orig.copy()
+            if isinstance(y_pred_orig, VectorizedDF):
+                y_pred_orig.X.index = y_true.index
+            else:
+                y_pred_orig.index = y_true.index
         if not same_cols:
             warn(
                 "y_pred and y_true do not have the same column index. "
                 "This may indicate incorrect objects passed to the metric. "
                 "Indices of y_true will be used for y_pred."
             )
-            y_pred = y_pred.copy()
-            y_pred.columns = y_true.columns
-
+            y_pred_orig = y_pred_orig.copy()
+            if isinstance(y_pred_orig, VectorizedDF):
+                y_pred_orig.X.columns = y_true.columns
+            else:
+                y_pred_orig.columns = y_true.columns
         # check multioutput arg
         # todo: add this back when variance_weighted is supported
         # ("raw_values", "uniform_average", "variance_weighted")
@@ -346,7 +444,11 @@ class BaseForecastingErrorMetric(BaseMetric):
                 )
 
         # check multilevel arg
-        allowed_multilevel_str = ("raw_values", "uniform_average")
+        allowed_multilevel_str = (
+            "raw_values",
+            "uniform_average",
+            "uniform_average_time",
+        )
 
         if not isinstance(multilevel, str):
             raise ValueError(f"multilevel must be a str, but found {type(multilevel)}")
@@ -356,7 +458,7 @@ class BaseForecastingErrorMetric(BaseMetric):
                 f"but found multilevel={multilevel}"
             )
 
-        return y_true, y_pred, multioutput, multilevel
+        return y_true_orig, y_pred_orig, multioutput, multilevel
 
     def _check_ys(self, y_true, y_pred, multioutput, multilevel, **kwargs):
 
@@ -364,18 +466,27 @@ class BaseForecastingErrorMetric(BaseMetric):
         INNER_MTYPES = ["pd.DataFrame", "pd-multiindex", "pd_multiindex_hier"]
 
         def _coerce_to_df(y, var_name="y"):
-            valid, msg, _ = check_is_scitype(
+            valid, msg, metadata = check_is_scitype(
                 y, scitype=SCITYPES, return_metadata=True, var_name=var_name
             )
             if not valid:
                 raise TypeError(msg)
             y_inner = convert_to(y, to_type=INNER_MTYPES)
+
+            scitype = metadata["scitype"]
+            ignore_index = multilevel == "uniform_average_time"
+            if scitype in ["Panel", "Hierarchical"] and not ignore_index:
+                y_inner = VectorizedDF(y_inner, is_scitype=scitype)
             return y_inner
 
         y_true = _coerce_to_df(y_true, var_name="y_true")
         y_pred = _coerce_to_df(y_pred, var_name="y_pred")
         if "y_train" in kwargs.keys():
             kwargs["y_train"] = _coerce_to_df(kwargs["y_train"], var_name="y_train")
+        if "y_pred_benchmark" in kwargs.keys():
+            kwargs["y_pred_benchmark"] = _coerce_to_df(
+                kwargs["y_pred_benchmark"], var_name="y_pred_benchmark"
+            )
 
         y_true, y_pred, multioutput, multilevel = self._check_consistent_input(
             y_true, y_pred, multioutput, multilevel


Test patch:
diff --git a/sktime/performance_metrics/tests/test_metrics_classes.py b/sktime/performance_metrics/tests/test_metrics_classes.py
index d4f988f81..82c2b7a12 100644
--- a/sktime/performance_metrics/tests/test_metrics_classes.py
+++ b/sktime/performance_metrics/tests/test_metrics_classes.py
@@ -6,12 +6,13 @@ import numpy as np
 import pandas as pd
 import pytest
 
-from sktime.performance_metrics.forecasting import _classes
+from sktime.performance_metrics.forecasting import MeanSquaredError, _classes
+from sktime.utils._testing.hierarchical import _make_hierarchical
 from sktime.utils._testing.series import _make_series
 
 metric_classes = getmembers(_classes, isclass)
 
-exclude_starts_with = ("_", "Base")
+exclude_starts_with = ("_", "Base", "Vectorized")
 metric_classes = [x for x in metric_classes if not x[0].startswith(exclude_starts_with)]
 
 names, metrics = zip(*metric_classes)
@@ -59,3 +60,39 @@ def test_metric_output_direct(metric, multioutput, n_columns):
 
     # assert results from all options are equal
     assert np.allclose(res[1], res[2])
+
+
+@pytest.mark.parametrize("n_columns", [1, 2])
+@pytest.mark.parametrize(
+    "multilevel", ["uniform_average", "uniform_average_time", "raw_values"]
+)
+@pytest.mark.parametrize("multioutput", ["uniform_average", "raw_values"])
+def test_metric_hierarchical(multioutput, multilevel, n_columns):
+    """Test hierarchical input for metrics."""
+    y_pred = _make_hierarchical(random_state=21, n_columns=n_columns)
+    y_true = _make_hierarchical(random_state=42, n_columns=n_columns)
+
+    metric = MeanSquaredError(multioutput=multioutput, multilevel=multilevel)
+
+    res = metric(
+        y_true=y_true,
+        y_pred=y_pred,
+    )
+
+    if multilevel == "raw_values":
+        assert isinstance(res, (pd.DataFrame, pd.Series))
+        assert isinstance(res.index, pd.MultiIndex)
+
+        expected_index = y_true.index.droplevel(-1).unique()
+        found_index = res.index.unique()
+        assert set(expected_index) == set(found_index)
+        if multioutput == "raw_values" and isinstance(res, pd.DataFrame):
+            assert all(y_true.columns == res.columns)
+    # if multilevel == "uniform_average" or "uniform_average_time"
+    else:
+        if multioutput == "uniform_average":
+            assert isinstance(res, float)
+        elif multioutput == "raw_values":
+            assert isinstance(res, np.ndarray)
+            assert res.ndim == 1
+            assert len(res) == len(y_true.columns)