Repository: sktime/sktime
Pull request: [ENH] forecasting `evaluate`: new `return_model` parameter to return fitted estimator states
#### Reference Issues/PRs

Fixes #7636 

#### What does this implement/fix? Explain your changes.

This draft PR enhances the `evaluate` function in sktime to optionally return fitted models for each fold in the cross-validation. This is implemented through a new `return_model` parameter which, when set to `True`, adds a `fitted_forecaster` column to the results DataFrame containing the fitted forecaster for each fold.

The enhancement is useful for:
- Analyzing model behavior across different folds
- Reusing fitted models from specific folds for further analysis or prediction
- Debugging model performance variations between folds

Changes made:
1. Added `return_model` parameter to `evaluate` function (default=False)
2. Modified `_get_column_order_and_datatype` to include fitted model column
3. Updated `_evaluate_window` to store fitted models in results
4. Added documentation for the new parameter and return value
5. Added example usage in docstring
6. Increased coverage of current tests to account for the new feature

#### Does your contribution introduce a new dependency? If yes, which one?

No new dependencies were introduced.

#### What should a reviewer concentrate their feedback on?

* [] Correctness of model storage and retrieval in parallel execution modes
* [] Memory implications of storing fitted models
* [] Documentation clarity and completeness
* [] Handling of fitted models in error cases
* [] Review updated tests

#### Did you add any tests for the change?

Yes but they could benefit from some review.

#### Any other comments?

The implementation maintains backward compatibility as the `return_model` parameter defaults to False. When using parallel backends (especially dask), users should be aware of potential memory implications when storing fitted models for many folds.

#### PR checklist

##### For all contributions
- [x] I've added myself to the [list of contributors](https://github.com/sktime/sktime/blob/main/CONTRIBUTORS.md) with any new badges I've earned :-)
- [x] The PR title starts with [ENH] as this is a feature enhancement

Problem statement:
[ENH] Improve output of `evaluate` when using grid search forecasters
**Is your feature request related to a problem? Please describe.**
When using `evaluate` with a grid search forecaster (e.g. ) it is currently not possible to retrieve and inspect the best forecaster (or even the different version of the forecaster) fitted in each CV fold. This is relevant if the output of evaluate is to define an overall best forecaster and/or to evaluate the stability of the best forecaster params.

Edit: actually, on reflecting about this, this could also be interesting for any "auto" style forecaster where the fitted model has parameters that could be interesting/informative/worth analysing. So the scope is even bigger.

**Describe the solution you'd like**
Ideally, `evaluate` would return the fitted model for inspection if set through some kind of `return_model: bool = False` param. This would enable the different workflows defined above.

**Describe alternatives you've considered**
Currently, I'm implementing my own version of `evaluate` to be able to access the fitted grid search object and analyse the different best forecasters for each fold.

**Additional context**
A minimal code example to replicate this behaviour would be:
```python
# %% Import necessary libraries
from sktime.datasets import load_airline
from sktime.forecasting.model_evaluation import evaluate
from sktime.forecasting.model_selection import ForecastingGridSearchCV
from sktime.forecasting.naive import NaiveForecaster
from sktime.split import ExpandingWindowSplitter
from sktime.utils.plotting import plot_windows

# %% Load the airline dataset
y = load_airline()
y.plot()
print(f"Length of the time series: {len(y)}")

# %% Define the expanding window splitter
sp = 12  # Seasonal periodicity
expanding_window_splitter = ExpandingWindowSplitter(
    initial_window=sp * 8,
    step_length=sp,
    fh=range(1, 13),  # Forecast horizon of 12 months
)
plot_windows(expanding_window_splitter, y)

# %% Set up the naive forecaster with mean method
naive_forecaster = NaiveForecaster(strategy="mean")

# %% Define a simple parameter grid for grid search
param_grid = {
    "window_length": [1, 12, 24, 36],  # Different windows to test
    "sp": [1],  # Seasonal periodicity
}

# %% Create a grid search forecaster using the naive model and the splitter
grid_cv = ExpandingWindowSplitter(
    initial_window=sp * 5, step_length=sp, fh=range(1, 13)
)
grid_search_forecaster = ForecastingGridSearchCV(
    forecaster=naive_forecaster,
    param_grid=param_grid,
    cv=grid_cv,
)

# %% Evaluate the forecaster using the evaluate method with the defined splitter
results = evaluate(
    grid_search_forecaster, y=y, cv=expanding_window_splitter, return_data=True
)

# %% Display evaluation results
results.head()  # <-- How can I retrieve the best forecaster of the grid search for each fold?
```

Hints from issue discussion:
I see - `get_fitted_params` is not sufficient in the `strategy="update"` case, right? Because intermediate states will be overwritten.
> I see - `get_fitted_params` is not sufficient in the `strategy="update"` case, right? Because intermediate states will be overwritten.

But this is an issue regardless of the strategy, no? There is currently no way to retrieve the fitted models per fold, which I think is the main issue. Thoughts?
oh, I see! You want the folds in `evaluate`. I get it now.

Commit message: Use explicit kwargs when calling functions

Patch to review:
diff --git a/sktime/forecasting/model_evaluation/_functions.py b/sktime/forecasting/model_evaluation/_functions.py
index 8ae6e774190..d8afccf796c 100644
--- a/sktime/forecasting/model_evaluation/_functions.py
+++ b/sktime/forecasting/model_evaluation/_functions.py
@@ -7,6 +7,7 @@
 
 import time
 import warnings
+from copy import deepcopy
 from typing import Optional, Union
 
 import numpy as np
@@ -74,7 +75,11 @@ def _check_scores(metrics) -> dict:
 
 
 def _get_column_order_and_datatype(
-    metric_types: dict, return_data: bool = True, cutoff_dtype=None, old_naming=True
+    metric_types: dict,
+    return_data: bool = True,
+    cutoff_dtype=None,
+    old_naming=True,
+    return_model: bool = False,
 ) -> dict:
     """Get the ordered column name and input datatype of results."""
     others_metadata = {
@@ -105,6 +110,8 @@ def _get_column_order_and_datatype(
     fit_metadata.update(others_metadata)
     if return_data:
         fit_metadata.update(y_metadata)
+    if return_model:
+        fit_metadata["fitted_forecaster"] = "object"
     metrics_metadata.update(fit_metadata)
     return metrics_metadata.copy()
 
@@ -191,6 +198,7 @@ def _evaluate_window(x, meta):
     strategy = meta["strategy"]
     scoring = meta["scoring"]
     return_data = meta["return_data"]
+    return_model = meta["return_model"]
     error_score = meta["error_score"]
     cutoff_dtype = meta["cutoff_dtype"]
 
@@ -302,12 +310,18 @@ def _evaluate_window(x, meta):
         temp_result["y_train"] = [y_train]
         temp_result["y_test"] = [y_test]
         temp_result.update(y_preds_cache)
+    if return_model:
+        temp_result["fitted_forecaster"] = [deepcopy(forecaster)]
     result = pd.DataFrame(temp_result)
     result = result.astype({"len_train_window": int, "cutoff": cutoff_dtype})
     if old_naming:
         result = result.rename(columns=old_name_mapping)
     column_order = _get_column_order_and_datatype(
-        scoring, return_data, cutoff_dtype, old_naming=old_naming
+        metric_types=scoring,
+        return_data=return_data,
+        cutoff_dtype=cutoff_dtype,
+        old_naming=old_naming,
+        return_model=return_model,
     )
     result = result.reindex(columns=column_order.keys())
 
@@ -326,6 +340,7 @@ def evaluate(
     strategy: str = "refit",
     scoring: Optional[Union[callable, list[callable]]] = None,
     return_data: bool = False,
+    return_model: bool = False,
     error_score: Union[str, int, float] = np.nan,
     backend: Optional[str] = None,
     cv_X=None,
@@ -371,6 +386,7 @@ def evaluate(
     * runtimes for fitting and/or predicting, from 2, 3, 7, in the ``i``-th loop
     * cutoff state of ``forecaster``, at 3, in the ``i``-th loop
     * :math:`y_{train, i}`, :math:`y_{test, i}`, ``y_pred`` (optional)
+    * fitted forecaster for each fold (optional)
 
     A distributed and-or parallel back-end can be chosen via the ``backend`` parameter.
 
@@ -400,6 +416,9 @@ def evaluate(
         Returns three additional columns in the DataFrame, by default False.
         The cells of the columns contain each a pd.Series for y_train,
         y_pred, y_test.
+    return_model : bool, default=False
+        If True, returns an additional column 'fitted_forecaster' containing the fitted
+        forecaster for each fold.
     error_score : "raise" or numeric, default=np.nan
         Value to assign to the score if an exception occurs in estimator fitting. If set
         to "raise", the exception is raised. If a numeric value is given,
@@ -468,6 +487,9 @@ def evaluate(
         - y_test: (pd.Series) present if see ``return_data=True``
         testing fold of the i-th split in ``cv``, used to compute the metric.
 
+        - fitted_forecaster: (BaseForecaster) present if see ``return_model=True``
+        fitted forecaster for the i-th split in ``cv``.
+
     Examples
     --------
     The type of evaluation that is done by ``evaluate`` depends on metrics in
@@ -503,15 +525,16 @@ def evaluate(
     ...     scoring=[MeanSquaredError(square_root=True), MeanAbsoluteError()],
     ... )
 
-    An example of an interval metric is the ``PinballLoss``.
-    It can be used with all probabilistic forecasters.
+    To get the fitted models for each fold, set return_model=True:
 
-    >>> from sktime.forecasting.naive import NaiveVariance
-    >>> from sktime.performance_metrics.forecasting.probabilistic import PinballLoss
-    >>> loss = PinballLoss()
-    >>> forecaster = NaiveForecaster(strategy="drift")
-    >>> results = evaluate(forecaster=NaiveVariance(forecaster),
-    ... y=y, cv=cv, scoring=loss)
+    >>> results = evaluate(
+    ...     forecaster=forecaster,
+    ...     y=y,
+    ...     cv=cv,
+    ...     scoring=loss,
+    ...     return_model=True
+    ... )
+    >>> fitted_forecaster = results.iloc[0]["fitted_forecaster"]
     """
     if backend in ["dask", "dask_lazy"]:
         if not _check_soft_dependencies("dask", severity="none"):
@@ -551,6 +574,7 @@ def evaluate(
         "scoring": scoring,
         "strategy": strategy,
         "return_data": return_data,
+        "return_model": return_model,
         "error_score": error_score,
         "cutoff_dtype": cutoff_dtype,
     }
@@ -621,7 +645,9 @@ def gen_y_X_train_test(y, X, cv, cv_X):
     if backend in ["dask", "dask_lazy"] and not not_parallel:
         import dask.dataframe as dd
 
-        metadata = _get_column_order_and_datatype(scoring, return_data, cutoff_dtype)
+        metadata = _get_column_order_and_datatype(
+            scoring, return_data, cutoff_dtype, return_model=return_model
+        )
 
         results = dd.from_delayed(results, meta=metadata)
         if backend == "dask":
diff --git a/sktime/forecasting/model_evaluation/tests/test_evaluate.py b/sktime/forecasting/model_evaluation/tests/test_evaluate.py
index 8579cab72e1..c29c2c7f200 100644
--- a/sktime/forecasting/model_evaluation/tests/test_evaluate.py
+++ b/sktime/forecasting/model_evaluation/tests/test_evaluate.py
@@ -23,6 +23,7 @@
 # from sktime.exceptions import FitFailedWarning
 # commented out until bugs are resolved, see test_evaluate_error_score
 from sktime.forecasting.arima import ARIMA, AutoARIMA
+from sktime.forecasting.base._base import BaseForecaster
 from sktime.forecasting.compose._reduce import DirectReductionForecaster
 from sktime.forecasting.exp_smoothing import ExponentialSmoothing
 from sktime.forecasting.model_evaluation import evaluate
@@ -66,11 +67,13 @@
 BACKENDS = _get_parallel_test_fixtures("estimator")
 
 
-def _check_evaluate_output(out, cv, y, scoring, return_data):
+def _check_evaluate_output(out, cv, y, scoring, return_data, return_model):
     assert isinstance(out, pd.DataFrame)
     # Check column names.
     scoring = _check_scores(scoring)
-    columns = _get_column_order_and_datatype(scoring, return_data)
+    columns = _get_column_order_and_datatype(
+        metric_types=scoring, return_data=return_data, return_model=return_model
+    )
     assert set(out.columns) == columns.keys(), "Columns are not identical"
 
     # Check number of rows against number of splits.
@@ -103,6 +106,14 @@ def _check_evaluate_output(out, cv, y, scoring, return_data):
     else:
         assert np.all(out.loc[:, "len_train_window"] == cv.window_length)
 
+    # Check fitted models
+    if return_model:
+        assert "fitted_forecaster" in out.columns
+        assert all(
+            isinstance(f, (BaseForecaster, type(None)))
+            for f in out["fitted_forecaster"].values
+        )
+
 
 @pytest.mark.skipif(
     not run_test_for_class(evaluate),
@@ -137,7 +148,9 @@ def test_evaluate_common_configs(
         scoring=scoring,
         **backend,
     )
-    _check_evaluate_output(out, cv, y, scoring, False)
+    _check_evaluate_output(
+        out=out, cv=cv, y=y, scoring=scoring, return_data=False, return_model=False
+    )
 
     # check scoring
     actual = out.loc[:, f"test_{scoring.name}"]
@@ -158,8 +171,9 @@ def test_evaluate_common_configs(
     reason="run test only if softdeps are present and incrementally (if requested)",
 )
 @pytest.mark.parametrize("return_data", [True, False])
+@pytest.mark.parametrize("return_model", [True, False])
 @pytest.mark.parametrize("scores", [METRICS, PROBA_METRICS, METRICS + PROBA_METRICS])
-def test_scoring_list(return_data, scores):
+def test_scoring_list(return_data, return_model, 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)
@@ -170,8 +184,16 @@ def test_scoring_list(return_data, scores):
         cv=cv,
         scoring=scores,
         return_data=return_data,
+        return_model=return_model,
+    )
+    _check_evaluate_output(
+        out=out,
+        cv=cv,
+        y=y,
+        scoring=scores,
+        return_data=return_data,
+        return_model=return_model,
     )
-    _check_evaluate_output(out, cv, y, scores, return_data)
 
 
 @pytest.mark.skipif(
@@ -189,7 +211,9 @@ def test_evaluate_initial_window():
     out = evaluate(
         forecaster=forecaster, y=y, cv=cv, strategy="update", scoring=scoring
     )
-    _check_evaluate_output(out, cv, y, scoring, False)
+    _check_evaluate_output(
+        out=out, cv=cv, y=y, scoring=scoring, return_data=False, return_model=False
+    )
     assert out.loc[0, "len_train_window"] == initial_window
 
     # check scoring
@@ -229,10 +253,13 @@ def test_evaluate_no_exog_against_with_exog():
 )
 @pytest.mark.parametrize("error_score", [np.nan, "raise", 1000])
 @pytest.mark.parametrize("return_data", [True, False])
+@pytest.mark.parametrize("return_model", [True, False])
 @pytest.mark.parametrize("strategy", ["refit", "update", "no-update_params"])
 @pytest.mark.parametrize("backend", BACKENDS)
 @pytest.mark.parametrize("scores", [[MeanAbsolutePercentageError()], METRICS])
-def test_evaluate_error_score(error_score, return_data, strategy, backend, scores):
+def test_evaluate_error_score(
+    error_score, return_data, return_model, strategy, backend, scores
+):
     """Test evaluate to raise warnings and exceptions according to error_score value."""
     forecaster = ExponentialSmoothing(sp=12)
     y = load_airline()
@@ -248,6 +275,7 @@ def test_evaluate_error_score(error_score, return_data, strategy, backend, score
         "cv": cv,
         "scoring": scores,
         "return_data": return_data,
+        "return_model": return_model,
         "error_score": error_score,
         "strategy": strategy,
     }


Code review comments:
diff_hunk: @@ -326,6 +340,7 @@ def evaluate(
     strategy: str = "refit",
     scoring: Optional[Union[callable, list[callable]]] = None,
     return_data: bool = False,
+    return_model: bool = False,
path: sktime/forecasting/model_evaluation/_functions.py
text: @user1:
I think this is considered a breaking change, since we are not using keywords only arguments. Just moving to emd should be fine I think. 

@user2:
I fixed it