Repository: sktime/sktime

Issue:
[BUG] fix logic for multiple scorers in `evaluate`
I looked at the logic for multiple scorers in `evaluate` and it is in essence broken - this seems to date back to the PR that originally claimed to add functionality to compute multiple scores, see https://github.com/sktime/sktime/pull/3883.

The nature of the breakage is as follows:

* if a list of scorers are passed, then only the first one is run through the entire machinery (!) - this includes crucial logic such as selecting the right predict method based on the scorer type, e.g., probabilistic, point prediction etc. The current logic inserts `scoring[0]` everywhere in the list case and does not touch the code otherwise.
* for the remaining scorers, evaluation happens separately, in the part that says "final formatting of results DataFrame". This part, in fact, does not just reformat, but runs the scoring loop for all other scorers - but omits crucial logic as mentioned above, so it will break for probabilistic scores or those requiring special arguments.

As a consequence, this also breaks logic of the tuning estimators that internally use `evaluate`.

What should happen is that `_evaluate_window` computes multiple scores, and takes multiple `scoring` arguments.

The final part of the code that does results formatting should do that and only that.

FYI @hazrulakmal due to the benchmarking topic and your recent work on this.

Discussion:
hazrulakmal: Also, since `evaluate` accepts any combination of metrics types, the column names that store `pred_time` and `y_pred` for respective predict types are slightly different. e.g probabilistic metrics will have "pred_time_proba", "y_pred_proba".  interval prediction will have "pred_time_interval" and "y_pred_interval" etc in order to help differentiate them.
hazrulakmal: For clarity, if it is not too difficult, would you kindly be able to produce the simplest example in which there is a substantial difference, with pre/post ordering?

1. pre-ordering looks like
metrics_1, fit_time, pred_time, len_train_window, cutoff, y_train, y_pred, y_test, metrics_2, metric_3 ...

2. post-ordering is sorted alphabetically, 
cutoff, fit_time, len_train_window, pred_time, metrics_1, metric_2, metric_3, ... , y_pred, y_test, y_train

    when there are two or more metric scitypes, there will be additional columns for pred_time, and y_pred. 

    the reason for changing the order is so that we can handle different cases consistently and also easier to handle metadata for dask.

hazrulakmal: When I think about it again, I think the nicest/proper solution (maybe I'm being nickpicking here) will be to maintain the order of the columns based on the status quo but bundle up the same thing next to each other. For example, all metrics coulmns are in the first few left columns, not like the current status quo where the first metrics at the very left while others are at the right of the dataframe.

Nevertheless this requires depereciation, informing the users that all metrics will be shown in the first few columns instead being left at the right.
fkiraly: > When I think about it again, I think the nicest/proper solution (maybe I'm being nitpicking here) will be to maintain the order of the columns based on the status quo but bundle up the same thing next to each other. For example, all metrics columns are in the first few left columns, not like the current status quo where the first metrics at the very left while others are at the right of the dataframe. I think is the closest to the original ordering with a lil improvements. what do you think?

I think this is a matter of taste and would be ok with anything that "clusters" same kind of columns or column groups next to each other.

From a user interface design perspective, the typical advice is to order columns in a way that the most important information (i.e., what the user will most likely want to look at) is furthest to the left. Similar, for rows or documents, the more important the further to the top it should be. Just a rule of thumb which has additional caveats if you are looking at a nested or structured presentation (like multiple tables, sections, etc).

My personal estimate of what is important would be, in decreasing order: metric, cutoff, timings. Or, it could be cutoff, metric, timings. Perhaps better to put cutoff to the left, as that's only one column and there are many metrics.
hazrulakmal: > From a user interface design perspective, the typical advice is to order columns in a way that the most important information 

I believe that the status quo is already arranged according to importance. However, it may not be properly ordered for multiple metrics. so I decided to follow the following arrangement, ordered by importance that will introduce minimal change to the existing user interface. 

1. In the case of multiple metrics of the same `.predict` method,
metrics_1, metrics_2, metric_3 ..., fit_time, pred_time, len_train_window, cutoff, y_train, y_test, y_pred, 

2. In the case of multiple metrics of different `predict` methods, say a combo of point, proba and interval forecast,
metrics_1, metrics_2, metric_3, fit_time, pred_time, pred_proba_time, pred_interval_time, len_train_window, cutoff, y_train, y_test, y_pred, y_pred_proba, y_pred_interval

I added a depreciation warning stating that this new arrangement will take place in v0.23.0. In the meantime, it will return the original ordering when only a single metric is given. The warning will only be raised when multiple metrics are given. Kindly review to ensure proper implementation. Thank you.

fkiraly: > I added a depreciation warning stating that this new arrangement will take place in v0.23.0.

We need to give users more warning, as 0.23.0 is due next week - 0.24.0 at least.

For the change, I think the message should say "please ensure you use `loc` references when addressing the columns", that should take care of the changes without needing to make changes to the `evaluate` function itself.

Also, please add todo comments `# todo 0.24.0` with explicit instructions what needs to be changed or removed at 0.24.0, so the release manager just needs to naively execute your instructions (rather than trying to understand what's going on).
fkiraly: (excellent work overall, it´s just deprecation is fiddly. All for our users, so they don´t suddenly get hit with a breaking change...)
hazrulakmal: > to do deprecation, I would suggest: the _get_column_order should return the current order, until 0.24.0. There should be an easy switch with a 0.24.0 release manager instruction on deleting a code block and replacing it with sth else.

I added two todo comments of deleting a code block. 

> the warning should state explicitly version of change, what will change where and how, before giving the "ensure" message

mentioned this too in the revised warning msg.

Sorry for the fiddly work. it wasnt quite clear abt how to do this depreciation before. 😓
fkiraly: (let us know when you are done and we should merge)
hazrulakmal: > (let us know when you are done and we should merge)

All good now to be merged!

Patch:
diff --git a/sktime/forecasting/model_evaluation/_functions.py b/sktime/forecasting/model_evaluation/_functions.py
index 8f2d2e236..cf71c62e9 100644
--- a/sktime/forecasting/model_evaluation/_functions.py
+++ b/sktime/forecasting/model_evaluation/_functions.py
@@ -7,7 +7,7 @@ __all__ = ["evaluate"]
 
 import time
 import warnings
-from typing import List, Optional, Union
+from typing import Dict, List, Optional, Union
 
 import numpy as np
 import pandas as pd
@@ -40,6 +40,65 @@ def _check_strategy(strategy):
         raise ValueError(f"`strategy` must be one of {valid_strategies}")
 
 
+def _check_scores(metrics) -> Dict:
+    """Validate and coerce to BaseMetric and segregate them based on predict type.
+
+    Parameters
+    ----------
+    metrics : sktime accepted metrics object or a list of them or None
+
+    Return
+    ------
+    metrics_type : Dict
+        The key is metric types and its value is a list of its corresponding metrics.
+    """
+    if not isinstance(metrics, List):
+        metrics = [metrics]
+
+    metrics_type = {}
+    for metric in metrics:
+        metric = check_scoring(metric)
+        # collect predict type
+        if hasattr(metric, "get_tag"):
+            scitype = metric.get_tag(
+                "scitype:y_pred", raise_error=False, tag_value_default="pred"
+            )
+        else:  # If no scitype exists then metric is a point forecast type
+            scitype = "pred"
+        if scitype not in metrics_type.keys():
+            metrics_type[scitype] = [metric]
+        else:
+            metrics_type[scitype].append(metric)
+    return metrics_type
+
+
+def _get_column_order_and_datatype(
+    metric_types: Dict, return_data: bool = True, cutoff_dtype=None
+) -> Dict:
+    """Get the ordered column name and input datatype of results."""
+    others_metadata = {
+        "len_train_window": "int",
+        "cutoff": cutoff_dtype,
+    }
+    y_metadata = {
+        "y_train": "object",
+        "y_test": "object",
+    }
+    fit_metadata, metrics_metadata = {"fit_time": "float"}, {}
+    for scitype in metric_types:
+        fit_metadata[f"{scitype}_time"] = "float"
+        if return_data:
+            y_metadata[f"y_{scitype}"] = "object"
+        for metric in metric_types.get(scitype):
+            metrics_metadata[f"test_{metric.name}"] = "float"
+    fit_metadata.update(others_metadata)
+    if return_data:
+        fit_metadata.update(y_metadata)
+    metrics_metadata.update(fit_metadata)
+    return metrics_metadata.copy()
+
+
+# should we remove _split since this is no longer being used?
 def _split(
     y,
     X,
@@ -112,7 +171,6 @@ def _evaluate_window(
     strategy,
     scoring,
     return_data,
-    score_name,
     error_score,
     cutoff_dtype,
 ):
@@ -122,6 +180,7 @@ def _evaluate_window(
     pred_time = np.nan
     cutoff = pd.Period(pd.NaT) if cutoff_dtype.startswith("period") else pd.NA
     y_pred = pd.NA
+    temp_result = dict()
 
     if fh is None:
         fh = _select_fh_from_y(y_test)
@@ -137,40 +196,38 @@ def _evaluate_window(
             forecaster.update(y_train, X_train, update_params=update_params)
         fit_time = time.perf_counter() - start_fit
 
+        # predict based on metrics
         pred_type = {
             "pred_quantiles": "predict_quantiles",
             "pred_interval": "predict_interval",
             "pred_proba": "predict_proba",
-            None: "predict",
+            "pred": "predict",
         }
-        # predict
-        start_pred = time.perf_counter()
-
-        if hasattr(scoring, "metric_args"):
-            metric_args = scoring.metric_args
-        else:
-            metric_args = {}
-
-        if hasattr(scoring, "get_tag"):
-            scitype = scoring.get_tag("scitype:y_pred", raise_error=False)
-        else:
-            # If no scitype exists then metric is not proba and no args needed
-            scitype = None
-
-        methodname = pred_type[scitype]
-        method = getattr(forecaster, methodname)
-
-        y_pred = method(fh=fh, X=X_test, **metric_args)
-        pred_time = time.perf_counter() - start_pred
-        # score
-        score = scoring(y_test, y_pred, y_train=y_train)
+        # cache prediction from the first scitype and reuse it to compute other metrics
+        for scitype in scoring:
+            method = getattr(forecaster, pred_type[scitype])
+            start_pred = time.perf_counter()
+            y_pred = method(fh, X_test)
+            pred_time = time.perf_counter() - start_pred
+            temp_result[f"{scitype}_time"] = [pred_time]
+            for metric in scoring.get(scitype):
+                score = metric(y_test, y_pred, y_train=y_train)
+                temp_result[f"test_{metric.name}"] = [score]
+            if return_data:
+                temp_result[f"y_{scitype}"] = [y_pred]
         # get cutoff
         cutoff = forecaster.cutoff
 
     except Exception as e:
         if error_score == "raise":
             raise e
-        else:
+        else:  # assign default value when fitting failed
+            for scitype in scoring:
+                temp_result[f"{scitype}_time"] = [pred_time]
+                if return_data:
+                    temp_result[f"y_{scitype}"] = [y_pred]
+                for metric in scoring.get(scitype):
+                    temp_result[f"test_{metric.name}"] = [score]
             warnings.warn(
                 f"""
                 In evaluate, fitting of forecaster {type(forecaster).__name__} failed,
@@ -190,18 +247,17 @@ def _evaluate_window(
     else:
         cutoff_ind = cutoff[0]
 
-    result = pd.DataFrame(
-        {
-            score_name: [score],
-            "fit_time": [fit_time],
-            "pred_time": [pred_time],
-            "len_train_window": [len(y_train)],
-            "cutoff": [cutoff_ind],
-            "y_train": [y_train if return_data else pd.NA],
-            "y_test": [y_test if return_data else pd.NA],
-            "y_pred": [y_pred if return_data else pd.NA],
-        }
-    ).astype({"cutoff": cutoff_dtype})
+    # Storing the remaining evaluate detail
+    temp_result["fit_time"] = [fit_time]
+    temp_result["len_train_window"] = [len(y_train)]
+    temp_result["cutoff"] = [cutoff_ind]
+    if return_data:
+        temp_result["y_train"] = [y_train]
+        temp_result["y_test"] = [y_test]
+    result = pd.DataFrame(temp_result)
+    result = result.astype({"len_train_window": int, "cutoff": cutoff_dtype})
+    column_order = _get_column_order_and_datatype(scoring, return_data, cutoff_dtype)
+    result = result.reindex(columns=column_order.keys())
 
     # Return forecaster if "update"
     if strategy == "update" or (strategy == "no-update_params" and i == 0):
@@ -394,10 +450,13 @@ def evaluate(
 
     _check_strategy(strategy)
     cv = check_cv(cv, enforce_start_with_window=True)
-    if isinstance(scoring, List):
-        scoring = [check_scoring(s) for s in scoring]
+    # TODO: remove lines(four lines below) and 599-612 in v0.25.0
+    if isinstance(scoring, list):
+        raise_warn, num = True, len(scoring)
     else:
-        scoring = check_scoring(scoring)
+        raise_warn, num = False, 1
+    # removal until here
+    scoring = _check_scores(scoring)
 
     ALLOWED_SCITYPES = ["Series", "Panel", "Hierarchical"]
 
@@ -419,20 +478,14 @@ def evaluate(
             )
         X = convert_to(X, to_type=PANDAS_MTYPES)
 
-    score_name = (
-        f"test_{scoring.name}"
-        if not isinstance(scoring, List)
-        else f"test_{scoring[0].name}"
-    )
     cutoff_dtype = str(y.index.dtype)
     _evaluate_window_kwargs = {
         "fh": cv.fh,
         "forecaster": forecaster,
-        "scoring": scoring if not isinstance(scoring, List) else scoring[0],
+        "scoring": scoring,
         "strategy": strategy,
-        "return_data": True,
+        "return_data": return_data,
         "error_score": error_score,
-        "score_name": score_name,
         "cutoff_dtype": cutoff_dtype,
     }
 
@@ -504,6 +557,7 @@ def evaluate(
         from dask import delayed as dask_delayed
 
         results = []
+        metadata = _get_column_order_and_datatype(scoring, return_data, cutoff_dtype)
         for i, (y_train, y_test, X_train, X_test) in enumerate(yx_splits):
             results.append(
                 dask_delayed(_evaluate_window)(
@@ -517,16 +571,7 @@ def evaluate(
             )
         results = dd.from_delayed(
             results,
-            meta={
-                score_name: "float",
-                "fit_time": "float",
-                "pred_time": "float",
-                "len_train_window": "int",
-                "cutoff": cutoff_dtype,
-                "y_train": "object",
-                "y_test": "object",
-                "y_pred": "object",
-            },
+            meta=metadata,
         )
         if compute:
             results = results.compute()
@@ -550,19 +595,22 @@ def evaluate(
 
     # final formatting of results DataFrame
     results = results.reset_index(drop=True)
-    if isinstance(scoring, List):
-        for s in scoring[1:]:
-            results[f"test_{s.name}"] = np.nan
-            for row in results.index:
-                results.loc[row, f"test_{s.name}"] = s(
-                    results["y_test"].loc[row],
-                    results["y_pred"].loc[row],
-                    y_train=results["y_train"].loc[row],
-                )
-
-    # drop pointer to data if not requested
-    if not return_data:
-        results = results.drop(columns=["y_train", "y_test", "y_pred"])
-    results = results.astype({"len_train_window": int})
 
+    # TODO: remove 16 lines below and 451-455 in v0.25.0
+    if raise_warn:
+        warnings.warn(
+            "Starting v0.25.0 model_evaluation.evaluate module will rearrange "
+            "all metric columns to the left of its output result DataFrame. "
+            "Please use loc references when addressing the columns. You can "
+            "safely ignore this warning if you don't use evaluate function directly.",
+            DeprecationWarning,
+            stacklevel=2,
+        )
+        columns = results.columns.to_list()
+        non_first_metrics = []
+        for _ in range(1, num):
+            metric = columns.pop(1)
+            non_first_metrics.append(metric)
+        results = results.reindex(columns=columns + non_first_metrics)
+    #  removal until here
     return results


Test patch:
diff --git a/sktime/forecasting/model_evaluation/tests/test_evaluate.py b/sktime/forecasting/model_evaluation/tests/test_evaluate.py
index 523e2e418..07278adaf 100644
--- a/sktime/forecasting/model_evaluation/tests/test_evaluate.py
+++ b/sktime/forecasting/model_evaluation/tests/test_evaluate.py
@@ -24,6 +24,10 @@ from sktime.forecasting.arima import ARIMA, AutoARIMA
 from sktime.forecasting.compose._reduce import DirectReductionForecaster
 from sktime.forecasting.exp_smoothing import ExponentialSmoothing
 from sktime.forecasting.model_evaluation import evaluate
+from sktime.forecasting.model_evaluation._functions import (
+    _check_scores,
+    _get_column_order_and_datatype,
+)
 from sktime.forecasting.naive import NaiveForecaster
 from sktime.forecasting.tests._config import TEST_FHS, TEST_STEP_LENGTHS_INT
 from sktime.performance_metrics.forecasting import (
@@ -45,18 +49,16 @@ from sktime.utils._testing.hierarchical import _make_hierarchical
 from sktime.utils._testing.series import _make_series
 from sktime.utils.validation._dependencies import _check_soft_dependencies
 
+METRICS = [MeanAbsolutePercentageError(symmetric=True), MeanAbsoluteScaledError()]
+PROBA_METRICS = [CRPS(), EmpiricalCoverage(), LogLoss(), PinballLoss()]
 
-def _check_evaluate_output(out, cv, y, scoring):
-    assert isinstance(out, pd.DataFrame)
 
+def _check_evaluate_output(out, cv, y, scoring, return_data):
+    assert isinstance(out, pd.DataFrame)
     # Check column names.
-    assert set(out.columns) == {
-        "cutoff",
-        "fit_time",
-        "len_train_window",
-        "pred_time",
-        f"test_{scoring.name}",
-    }
+    scoring = _check_scores(scoring)
+    columns = _get_column_order_and_datatype(scoring, return_data)
+    assert set(out.columns) == columns.keys(), "Columns are not identical"
 
     # Check number of rows against number of splits.
     n_splits = cv.get_n_splits(y)
@@ -100,13 +102,7 @@ def _check_evaluate_output(out, cv, y, scoring):
 @pytest.mark.parametrize("window_length", [7, 10])
 @pytest.mark.parametrize("step_length", TEST_STEP_LENGTHS_INT)
 @pytest.mark.parametrize("strategy", ["refit", "update", "no-update_params"])
-@pytest.mark.parametrize(
-    "scoring",
-    [
-        MeanAbsolutePercentageError(symmetric=True),
-        MeanAbsoluteScaledError(),
-    ],
-)
+@pytest.mark.parametrize("scoring", METRICS)
 @pytest.mark.parametrize("backend", [None, "dask", "loky", "threading"])
 def test_evaluate_common_configs(
     CV, fh, window_length, step_length, strategy, scoring, backend
@@ -128,7 +124,7 @@ def test_evaluate_common_configs(
         scoring=scoring,
         backend=backend,
     )
-    _check_evaluate_output(out, cv, y, scoring)
+    _check_evaluate_output(out, cv, y, scoring, False)
 
     # check scoring
     actual = out.loc[:, f"test_{scoring.name}"]
@@ -144,11 +140,12 @@ def test_evaluate_common_configs(
 
 
 @pytest.mark.skipif(
-    not run_test_for_class(evaluate),
+    not run_test_for_class([evaluate] + PROBA_METRICS),
     reason="run test only if softdeps are present and incrementally (if requested)",
 )
 @pytest.mark.parametrize("return_data", [True, False])
-def test_scoring_list(return_data):
+@pytest.mark.parametrize("scores", [METRICS, PROBA_METRICS, METRICS + PROBA_METRICS])
+def test_scoring_list(return_data, scores):
     y = make_forecasting_problem(n_timepoints=30, index_type="int")
     forecaster = NaiveForecaster()
     cv = SlidingWindowSplitter(fh=[1, 2, 3], initial_window=15, step_length=5)
@@ -157,22 +154,10 @@ def test_scoring_list(return_data):
         forecaster=forecaster,
         y=y,
         cv=cv,
-        scoring=[
-            MeanAbsolutePercentageError(symmetric=True),
-            MeanAbsoluteScaledError(),
-        ],
+        scoring=scores,
         return_data=return_data,
     )
-    assert "test_MeanAbsolutePercentageError" in out.columns
-    assert "test_MeanAbsoluteScaledError" in out.columns
-    if return_data:
-        assert "y_pred" in out.columns
-        assert "y_train" in out.columns
-        assert "y_test" in out.columns
-    else:
-        assert "y_pred" not in out.columns
-        assert "y_train" not in out.columns
-        assert "y_test" not in out.columns
+    _check_evaluate_output(out, cv, y, scores, return_data)
 
 
 @pytest.mark.skipif(
@@ -190,7 +175,7 @@ def test_evaluate_initial_window():
     out = evaluate(
         forecaster=forecaster, y=y, cv=cv, strategy="update", scoring=scoring
     )
-    _check_evaluate_output(out, cv, y, scoring)
+    _check_evaluate_output(out, cv, y, scoring, False)
     assert out.loc[0, "len_train_window"] == initial_window
 
     # check scoring
@@ -232,7 +217,8 @@ def test_evaluate_no_exog_against_with_exog():
 @pytest.mark.parametrize("return_data", [True, False])
 @pytest.mark.parametrize("strategy", ["refit", "update", "no-update_params"])
 @pytest.mark.parametrize("backend", [None, "dask", "loky", "threading"])
-def test_evaluate_error_score(error_score, return_data, strategy, backend):
+@pytest.mark.parametrize("scores", [[MeanAbsolutePercentageError()], METRICS])
+def test_evaluate_error_score(error_score, return_data, strategy, backend, scores):
     """Test evaluate to raise warnings and exceptions according to error_score value."""
     # skip test for dask backend if dask is not installed
     if backend == "dask" and not _check_soft_dependencies("dask", severity="none"):
@@ -243,22 +229,24 @@ def test_evaluate_error_score(error_score, return_data, strategy, backend):
     # add NaN to make ExponentialSmoothing fail
     y.iloc[1] = np.nan
     fh = [1, 2, 3]
-    cv = ExpandingWindowSplitter(step_length=48, initial_window=12, fh=fh)
+    cv = SlidingWindowSplitter(step_length=33, initial_window=36, fh=fh)
+    scoring_name = [f"test_{score.name}" for score in scores]
     if error_score in [np.nan, 1000]:
         with pytest.warns(FitFailedWarning):
             results = evaluate(
                 forecaster=forecaster,
                 y=y,
                 cv=cv,
+                scoring=scores,
                 return_data=return_data,
                 error_score=error_score,
                 strategy=strategy,
                 backend=backend,
             )
         if isinstance(error_score, type(np.nan)):
-            assert results["test_MeanAbsolutePercentageError"].isna().sum() > 0
+            assert all(results[scoring_name].isna().sum() > 0)
         if error_score == 1000:
-            assert results["test_MeanAbsolutePercentageError"].max() == 1000
+            assert all(results[scoring_name].max() == 1000)
     if error_score == "raise":
         with pytest.raises(Exception):  # noqa: B017
             evaluate(
@@ -331,22 +319,18 @@ def test_evaluate_bigger_X(cls):
     evaluate(forecaster=f, y=y, X=X, cv=cv, error_score="raise", scoring=loss)
 
 
-PROBA_METRICS = [CRPS, EmpiricalCoverage, LogLoss, PinballLoss]
-
-
 @pytest.mark.skipif(
     not run_test_for_class([evaluate] + PROBA_METRICS),
     reason="run test only if softdeps are present and incrementally (if requested)",
 )
 @pytest.mark.parametrize("n_columns", [1, 2])
-@pytest.mark.parametrize("metric", PROBA_METRICS)
-def test_evaluate_probabilistic(n_columns, metric):
+@pytest.mark.parametrize("scoring", PROBA_METRICS)
+def test_evaluate_probabilistic(n_columns, scoring):
     """Check that evaluate works with interval, quantile, and distribution forecasts."""
     y = _make_series(n_columns=n_columns)
 
     forecaster = NaiveForecaster()
     cv = SlidingWindowSplitter()
-    scoring = metric()
     try:
         out = evaluate(
             forecaster,
@@ -397,4 +381,5 @@ def test_evaluate_hierarchical_unequal_X_y():
     assert res.shape == (5, 5)
 
     expected_cols = np.array([1 / 2, 1 / 3, 1 / 4, 1 / 5, 1 / 6])
-    _assert_array_almost_equal(res.iloc[:, 0].values, expected_cols)
+    output_metrics = res.loc[:, "test_MeanAbsolutePercentageError"].values
+    _assert_array_almost_equal(output_metrics, expected_cols)