Repository: sktime/sktime
Pull request: [ENH] Enable `EnbPIForecaster` to work with `sktime` native bootstrap algorithms
#### Reference Issues/PRs
fixes #7561

Depends on https://github.com/sktime/sktime/pull/8023 for test fixes, which should be merged first.

#### What does this implement/fix? Explain your changes.
Add the option to bootstrap algorithms to return the indices. Change the Enbpi forecaster so that it works with native bootstrapper by changing the internally expected interface.

#### Does your contribution introduce a new dependency? If yes, which one?

No 
#### What should a reviewer concentrate their feedback on?

* Are the deprecation cycles okay?
* Should a tag be introduced showing the possibility of the transformations to return indices? 

#### 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 :-)
  How to: add yourself to the [all-contributors file](https://github.com/sktime/sktime/blob/main/.all-contributorsrc) in the `sktime` root directory (not the `CONTRIBUTORS.md`). Common badges: `code` - fixing a bug, or adding code logic. `doc` - writing or improving documentation or docstrings. `bug` - reporting or diagnosing a bug (get this plus `code` if you also fixed the bug in the PR).`maintenance` - CI, test framework, release.
  See here for [full badge reference](https://allcontributors.org/docs/en/emoji-key)
- [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.

<!--
Thanks for contributing!
-->


Problem statement:
[ENH] Enable Enbpi with native bootstrap algorithms.
**Is your feature request related to a problem? Please describe.**
Currently enbpi is only usable with bootstrap algorithms from tsbootstrap. Change it so that it uses only bootstrap algorithms implemented in sktime. tsbootstrap would still be possible by using the tsbootstrapadapter.

Requires that the bootstrap algorithms return the original indices that led to the bootstrap times series.

Hints from issue discussion:
Hey @benHeid I have tried this issue ..I hope I am doing the right thing

Commit message: Fix comments

Patch to review:
diff --git a/sktime/forecasting/enbpi.py b/sktime/forecasting/enbpi.py
index 65c140d3b1a..536de213740 100644
--- a/sktime/forecasting/enbpi.py
+++ b/sktime/forecasting/enbpi.py
@@ -1,6 +1,8 @@
 # copyright: sktime developers, BSD-3-Clause License (see LICENSE file)
 """Implements EnbPIForecaster."""
 
+from warnings import warn
+
 import numpy as np
 import pandas as pd
 from sklearn.base import clone
@@ -9,6 +11,10 @@
 from sktime.forecasting.base import BaseForecaster
 from sktime.forecasting.naive import NaiveForecaster
 from sktime.libs._aws_fortuna_enbpi.enbpi import EnbPI
+from sktime.transformations.bootstrap import (
+    MovingBlockBootstrapTransformer,
+    TSBootstrapAdapter,
+)
 
 __all__ = ["EnbPIForecaster"]
 __author__ = ["benheid"]
@@ -62,7 +68,7 @@ class EnbPIForecaster(BaseForecaster):
     Examples
     --------
     >>> from tsbootstrap import MovingBlockBootstrap
-    >>> from sktime.forecasting.compose import EnbPIForecaster
+    >>> from sktime.forecasting.enbpi import EnbPIForecaster
     >>> from sktime.forecasting.naive import NaiveForecaster
     >>> from sktime.datasets import load_airline
     >>> from sktime.transformations.series.difference import Differencer
@@ -71,7 +77,7 @@ class EnbPIForecaster(BaseForecaster):
     >>> forecaster = Differencer(lags=[1]) * Deseasonalizer(sp=12) * EnbPIForecaster(
     ...    forecaster=NaiveForecaster(sp=12),
     ...    bootstrap_transformer=MovingBlockBootstrap(n_bootstraps=10))
-    >>> forecaster.fit(y, fh=range(12))
+    >>> forecaster = forecaster.fit(y, fh=range(1, 12)) # Range is broken
     >>> res = forecaster.predict()
     >>> res_int = forecaster.predict_interval(coverage=[0.5])
 
@@ -128,13 +134,41 @@ def __init__(
 
         super().__init__()
 
-        from tsbootstrap import MovingBlockBootstrap
+        from tsbootstrap import BaseTimeSeriesBootstrap, MovingBlockBootstrap
 
-        self.bootstrap_transformer_ = (
-            clone(bootstrap_transformer)
-            if bootstrap_transformer is not None
-            else MovingBlockBootstrap()
-        )
+        if isinstance(bootstrap_transformer, BaseTimeSeriesBootstrap):
+            warn(
+                "Passing directly bootstrapper from tsbootstrap is deprecated."
+                + "And will change in version 0.37.0"
+                + "Please wrap bootstrapper from tsbootstrap using tsbootstrapAdapter.",
+                DeprecationWarning,
+            )
+            self.bootstrap_transformer_ = TSBootstrapAdapter(
+                bootstrap_transformer, return_indices=True
+            )
+        else:
+            self.bootstrap_transformer_ = bootstrap_transformer
+
+        if self.bootstrap_transformer is None:
+            warn(
+                "The default value for the bootstrap_transformer will change to the"
+                + "sktime MovingBlockBootstrap in version 0.37.0."
+                + "If you want to use the current behavior"
+                + "Please bass  TSBootstrapAdapter(MovingBlockBootstrap())"
+            )
+
+            self.bootstrap_transformer_ = (
+                clone(bootstrap_transformer)
+                if bootstrap_transformer is not None
+                else TSBootstrapAdapter(MovingBlockBootstrap())
+            )
+        if (
+            not self.bootstrap_transformer_.get_tags()["can_return_indices"]
+            or not self.bootstrap_transformer_.return_indices
+        ):
+            raise ValueError(
+                "The bootstrap_transformer needs to be able to return indices"
+            )
 
     def _fit(self, X, y, fh=None):
         self._fh = fh
@@ -144,16 +178,16 @@ def _fit(self, X, y, fh=None):
         self.random_state_ = check_random_state(self.random_state)
 
         # fit/transform the transformer to obtain bootstrap samples
-        bs_ts_index = list(
-            self.bootstrap_transformer_.bootstrap(y, test_ratio=0, return_indices=True)
-        )
-        self.indexes = np.stack(list(map(lambda x: x[1], bs_ts_index)))
-        bootstrapped_ts = list(map(lambda x: x[0], bs_ts_index))
+        bs_ts_index = self.bootstrap_transformer_.fit_transform(y)
+
+        self.indexes = bs_ts_index["resampled_index"].values.reshape((-1, len(y)))
+        bootstrapped_ts = bs_ts_index[y.columns]
 
         self.forecasters = []
         self._preds = []
         # Fit Models per Bootstrap Sample
-        for bs_ts in bootstrapped_ts:
+        for bs_index in bootstrapped_ts.index.get_level_values(0).unique():
+            bs_ts = bootstrapped_ts.loc[bs_index]
             bs_df = pd.DataFrame(bs_ts, index=y.index)
             forecaster = clone(self.forecaster_)
             forecaster.fit(y=bs_df, fh=fh, X=X)
@@ -177,7 +211,7 @@ def _predict(self, X, fh=None):
     def _predict_interval(self, fh, X, coverage):
         preds = []
         for forecaster in self.forecasters:
-            preds.append(forecaster.predict(fh=fh, X=X))
+            preds.append(forecaster.predict(fh=fh, X=X).values)
 
         train_targets = self._y.copy()
         train_targets.index = pd.RangeIndex(len(train_targets))
@@ -233,21 +267,18 @@ def get_test_params(cls):
             instance.
             ``create_test_instance`` uses the first (or only) dictionary in ``params``
         """
-        from sktime.utils.dependencies import _check_soft_dependencies
-
-        deps = cls.get_class_tag("python_dependencies")
-
-        if _check_soft_dependencies(deps, severity="none"):
-            from tsbootstrap import BlockBootstrap
-
-            params = [
-                {},
-                {
-                    "forecaster": NaiveForecaster(),
-                    "bootstrap_transformer": BlockBootstrap(),
-                },
-            ]
-        else:
-            params = {}
+        params = [
+            {
+                "bootstrap_transformer": MovingBlockBootstrapTransformer(
+                    return_indices=True
+                ),
+            },
+            {
+                "forecaster": NaiveForecaster(),
+                "bootstrap_transformer": MovingBlockBootstrapTransformer(
+                    return_indices=True
+                ),
+            },
+        ]
 
         return params
diff --git a/sktime/registry/_tags.py b/sktime/registry/_tags.py
index 0c46bb5b725..c3af8d483c2 100644
--- a/sktime/registry/_tags.py
+++ b/sktime/registry/_tags.py
@@ -2155,6 +2155,12 @@ class visual_block_kind(_BaseTag):
         "dict",
         "deprecated tag for dependency import aliases",
     ),
+    (
+        "can_return_indices",
+        ["transformer"],
+        "bool",
+        "The bootstrapper can return indices of the bootstrapped samples",
+    ),
 ]
 
 # construct the tag register from all classes in this module
diff --git a/sktime/transformations/bootstrap/_mbb.py b/sktime/transformations/bootstrap/_mbb.py
index 1cde691de4a..b901a711d3d 100644
--- a/sktime/transformations/bootstrap/_mbb.py
+++ b/sktime/transformations/bootstrap/_mbb.py
@@ -3,7 +3,6 @@
 
 __author__ = ["ltsaprounis"]
 
-from copy import copy
 from typing import Union
 
 import numpy as np
@@ -119,6 +118,9 @@ class STLBootstrapTransformer(BaseTransformer):
         This param goes into STL.fit() from statsmodels.
     random_state : int, np.random.RandomState or None, by default None
         Controls the randomness of the estimator
+    return_indices : bool, optional
+        If True, the output will contain the resampled indices as extra column,
+        by default False.
 
     See Also
     --------
@@ -189,6 +191,7 @@ class STLBootstrapTransformer(BaseTransformer):
         "enforce_index_type": None,  # index type that needs to be enforced in X/y
         "fit_is_empty": False,  # is fit empty and can be skipped? Yes = True
         "transform-returns-same-time-index": False,
+        "can_return_indices": True,
     }
 
     def __init__(
@@ -213,6 +216,7 @@ def __init__(
         inner_iter: int = None,
         outer_iter: int = None,
         random_state: Union[int, np.random.RandomState] = None,
+        return_indices=False,
     ):
         self.n_series = n_series
         self.sp = sp
@@ -234,6 +238,7 @@ def __init__(
         self.inner_iter = inner_iter
         self.outer_iter = outer_iter
         self.random_state = random_state
+        self.return_indices = return_indices
 
         super().__init__()
 
@@ -301,6 +306,9 @@ def _transform(self, X, y=None):
 
         Xcol = X.columns
         X = X[X.columns[0]]
+        if self.return_indices:
+            Xcol = list(Xcol)
+            Xcol.append("resampled_index")
 
         if len(X) <= self.block_length_:
             raise ValueError(
@@ -335,16 +343,17 @@ def _transform(self, X, y=None):
 
         # initialize the dataframe that will store the bootstrapped series
         if self.return_actual:
-            df_list = [
-                pd.DataFrame(
-                    X.values,
-                    index=pd.MultiIndex.from_product(
-                        iterables=[["actual"], X_index],
-                        names=["series_id", "time_index"],
-                    ),
-                    columns=[col_name],
-                )
-            ]
+            df = pd.DataFrame(
+                X.values,
+                index=pd.MultiIndex.from_product(
+                    iterables=[["actual"], X_index],
+                    names=["series_id", "time_index"],
+                ),
+                columns=[col_name],
+            )
+            if self.return_indices:
+                df["resampled_index"] = range(len(X))
+            df_list = [df]
         else:
             df_list = []
 
@@ -352,15 +361,15 @@ def _transform(self, X, y=None):
         rng = check_random_state(self.random_state)
         # create multiple series
         for i in range(self.n_series):
+            bs_ts, indices = _moving_block_bootstrap(
+                ts=resid,
+                block_length=self.block_length_,
+                replacement=self.sampling_replacement,
+                random_state=rng,
+                return_indices=True,
+            )
             new_series = self.box_cox_transformer_.inverse_transform(
-                _moving_block_bootstrap(
-                    ts=resid,
-                    block_length=self.block_length_,
-                    replacement=self.sampling_replacement,
-                    random_state=rng,
-                )
-                + seasonal
-                + trend
+                bs_ts + seasonal + trend
             )
 
             new_series_id = f"synthetic_{i}"
@@ -369,11 +378,12 @@ def _transform(self, X, y=None):
                 names=["series_id", "time_index"],
             )
 
-            df_list.append(
-                pd.DataFrame(
-                    data=new_series.values, index=new_df_index, columns=[col_name]
-                )
+            df = pd.DataFrame(
+                data=new_series.values, index=new_df_index, columns=[col_name]
             )
+            if self.return_indices:
+                df["resampled_index"] = indices
+            df_list.append(df)
 
         Xt = pd.concat(df_list)
         Xt.columns = Xcol
@@ -440,6 +450,9 @@ class MovingBlockBootstrapTransformer(BaseTransformer):
         The actual time series will be labelled as "actual"
     random_state : int, np.random.RandomState or None, by default None
         Controls the randomness of the estimator
+    return_indices : bool, optional
+        If True, the output will contain the resampled indices as extra column,
+        by default False.
 
     See Also
     --------
@@ -503,6 +516,7 @@ class MovingBlockBootstrapTransformer(BaseTransformer):
         "enforce_index_type": None,  # index type that needs to be enforced in X/y
         "fit_is_empty": True,  # is fit empty and can be skipped? Yes = True
         "transform-returns-same-time-index": False,
+        "can_return_indices": True,
     }
 
     def __init__(
@@ -512,12 +526,14 @@ def __init__(
         sampling_replacement: bool = False,
         return_actual: bool = True,
         random_state: Union[int, np.random.RandomState] = None,
+        return_indices=False,
     ):
         self.n_series = n_series
         self.block_length = block_length
         self.sampling_replacement = sampling_replacement
         self.return_actual = return_actual
         self.random_state = random_state
+        self.return_indices = return_indices
 
         super().__init__()
 
@@ -538,7 +554,9 @@ def _transform(self, X, y=None):
         -------
         transformed version of X
         """
-        Xcol = X.columns
+        Xcol = list(X.columns)
+        if self.return_indices:
+            Xcol.append("resampled_index")
         X = X[X.columns[0]]
 
         if len(X) <= self.block_length:
@@ -554,16 +572,17 @@ def _transform(self, X, y=None):
 
         # initialize the dataframe that will store the bootstrapped series
         if self.return_actual:
-            df_list = [
-                pd.DataFrame(
-                    X.values,
-                    index=pd.MultiIndex.from_product(
-                        iterables=[["actual"], X_index],
-                        names=["series_id", "time_index"],
-                    ),
-                    columns=[col_name],
-                )
-            ]
+            df = pd.DataFrame(
+                X.values,
+                index=pd.MultiIndex.from_product(
+                    iterables=[["actual"], X_index],
+                    names=["series_id", "time_index"],
+                ),
+                columns=[col_name],
+            )
+            if self.return_indices:
+                df["resampled_index"] = range(len(X))
+            df_list = [df]
         else:
             df_list = []
 
@@ -571,11 +590,12 @@ def _transform(self, X, y=None):
         rng = check_random_state(self.random_state)
         # create multiple series
         for i in range(self.n_series):
-            new_series = _moving_block_bootstrap(
+            new_series, indices = _moving_block_bootstrap(
                 ts=X,
                 block_length=self.block_length,
                 replacement=self.sampling_replacement,
                 random_state=rng,
+                return_indices=True,
             )
 
             new_series_id = f"synthetic_{i}"
@@ -583,11 +603,12 @@ def _transform(self, X, y=None):
                 iterables=[[new_series_id], new_series.index],
                 names=["series_id", "time_index"],
             )
-            df_list.append(
-                pd.DataFrame(
-                    data=new_series.values, index=new_df_index, columns=[col_name]
-                )
+            df = pd.DataFrame(
+                data=new_series.values, index=new_df_index, columns=[col_name]
             )
+            if self.return_indices:
+                df["resampled_index"] = indices
+            df_list.append(df)
 
         Xt = pd.concat(df_list)
         Xt.columns = Xcol
@@ -629,6 +650,7 @@ def _moving_block_bootstrap(
     block_length: int,
     replacement: bool = False,
     random_state: Union[int, np.random.RandomState] = None,
+    return_indices=False,
 ) -> pd.Series:
     """Create a synthetic time series using the moving block bootstrap method MBB.
 
@@ -660,26 +682,31 @@ def _moving_block_bootstrap(
         )
 
     if block_length == 1 and not replacement:
-        mbb_values = copy(ts_values)
-        rng.shuffle(mbb_values)
+        mbb_indices = np.arange(len(ts_index))
+        rng.shuffle(mbb_indices)
     elif block_length == 1:
-        mbb_values = rng.choice(ts_values, size=ts_length, replace=replacement)
+        mbb_indices = rng.choice(
+            np.arange(len(ts_index)), size=ts_length, replace=replacement
+        )
     else:
+        ts_int_indexes = np.arange(len(ts_index))
         total_num_blocks = int(ts_length / block_length) + 2
         block_origns = rng.choice(
             ts_length - block_length + 1, size=total_num_blocks, replace=replacement
         )
-        mbb_values = [
-            val for i in block_origns for val in ts_values[i : i + block_length]
+        mbb_indices = [
+            val for i in block_origns for val in ts_int_indexes[i : i + block_length]
         ]
         # remove the first few observations and ensure new series has the
         # same length as the original
         remove_first = rng.choice(block_length - 1)
-        mbb_values = mbb_values[remove_first : remove_first + ts_length]
+        mbb_indices = mbb_indices[remove_first : remove_first + ts_length]
 
-    mbb_series = pd.Series(data=mbb_values, index=ts_index)
-
-    return mbb_series
+    mbb_series = pd.Series(data=ts_values[mbb_indices], index=ts_index)
+    if not return_indices:
+        return mbb_series
+    else:
+        return mbb_series, mbb_indices
 
 
 def _get_series_name(ts: Union[pd.Series, pd.DataFrame]) -> str:
diff --git a/sktime/transformations/bootstrap/_tsbootstrap.py b/sktime/transformations/bootstrap/_tsbootstrap.py
index 93894043abb..0cecd51ab42 100644
--- a/sktime/transformations/bootstrap/_tsbootstrap.py
+++ b/sktime/transformations/bootstrap/_tsbootstrap.py
@@ -24,18 +24,19 @@ class TSBootstrapAdapter(BaseTransformer):
         The splitter used for the bootstrap splitting.
     include_actual : bool, default=False
         Whether to include the actual data in the output.
+    return_indices : bool, optional
+        If True, the output will contain the resampled indices as extra column,
+        by default False.
 
     Examples
     --------
     >>> from sktime.datasets import load_airline
     >>> from sktime.transformations.bootstrap import TSBootstrapAdapter
-    >>> from tsbootstrap import (
-    ...    MovingBlockBootstrap,
-    ...    MovingBlockBootstrapConfig
-    ... ) # doctest: +SKIP
+    >>> from tsbootstrap import MovingBlockBootstrap, # doctest: +SKIP
     >>> y = load_airline()
-    >>> config = MovingBlockBootstrapConfig(10, n_bootstraps=10)  # doctest: +SKIP
-    >>> bootstrap = TSBootstrapAdapter(MovingBlockBootstrap(config))  # doctest: +SKIP
+    >>> bootstrap = TSBootstrapAdapter(
+    ...    MovingBlockBootstrap(n_bootstrap=10, block_length=10)
+    ... )  # doctest: +SKIP
     >>> result = bootstrap.fit_transform(y)  # doctest: +SKIP
     """
 
@@ -59,15 +60,18 @@ class TSBootstrapAdapter(BaseTransformer):
         "enforce_index_type": None,  # index type that needs to be enforced in X/y
         "fit_is_empty": True,  # is fit empty and can be skipped? Yes = True
         "transform-returns-same-time-index": True,
+        "can_return_indices": True,
     }
 
     def __init__(
         self,
         bootstrap,
         include_actual=False,
+        return_indices=False,
     ):
         self.bootstrap = bootstrap
         self.include_actual = include_actual
+        self.return_indices = return_indices
         super().__init__()
 
     def _transform(self, X, y=None):
@@ -87,9 +91,20 @@ def _transform(self, X, y=None):
         -------
         transformed version of X
         """
-        bootstrapped_samples = self.bootstrap.bootstrap(X, test_ratio=0)
+        # Need to be cloned otherwise it helds a inner state and fails during update
+        bootstrapped_samples = self.bootstrap.clone().bootstrap(
+            X, test_ratio=0, return_indices=True
+        )
 
         def wrap_df(spl):
+            if isinstance(spl, tuple):
+                index = spl[1]
+                spl = spl[0]
+                df = pd.DataFrame(spl, index=X.index, columns=X.columns)
+                return pd.concat(
+                    [df, pd.Series(index, X.index, name="resampled_index")], axis=1
+                )
+
             return pd.DataFrame(spl, index=X.index, columns=X.columns)
 
         bootstrapped_samples = [wrap_df(sample) for sample in bootstrapped_samples]
@@ -100,7 +115,12 @@ def wrap_df(spl):
             bootstrapped_samples = [X_actual] + bootstrapped_samples
             spl_keys = ["actual"] + spl_keys
 
-        return pd.concat(bootstrapped_samples, keys=spl_keys, axis=0)
+        df = pd.concat(bootstrapped_samples, keys=spl_keys, axis=0)
+        if self.return_indices:
+            return df
+        else:
+            del df["resampled_index"]
+            return df
 
     @classmethod
     def get_test_params(cls, parameter_set="default"):


Code review comments:
diff_hunk: @@ -128,13 +134,41 @@ def __init__(
 
         super().__init__()
 
-        from tsbootstrap import MovingBlockBootstrap
+        from tsbootstrap import BaseTimeSeriesBootstrap, MovingBlockBootstrap
 
-        self.bootstrap_transformer_ = (
-            clone(bootstrap_transformer)
-            if bootstrap_transformer is not None
-            else MovingBlockBootstrap()
-        )
+        if isinstance(bootstrap_transformer, BaseTimeSeriesBootstrap):
+            warn(
+                "Passing directly bootstrapper from tsbootstrap is deprecated."
+                + "And will change in version 0.37.0"
+                + "Please wrap bootstrapper from tsbootstrap using tsbootstrapAdapter.",
+                DeprecationWarning,
+            )
+            self.bootstrap_transformer_ = TSBootstrapAdapter(
+                bootstrap_transformer, return_indices=True
+            )
+        else:
+            self.bootstrap_transformer_ = bootstrap_transformer
+
+        if self.bootstrap_transformer is None:
+            warn(
+                "The default value for the bootstrap_transformer will change to the"
+                + "sktime MovingBlockBootstrap in version 0.37.0."
+                + "If you want to use the current behavior"
+                + "Please bass  TSBootstrapAdapter(MovingBlockBootstrap())"
path: sktime/forecasting/enbpi.py
text: @user1:
typo, pass

diff_hunk: @@ -128,13 +134,41 @@ def __init__(
 
         super().__init__()
 
-        from tsbootstrap import MovingBlockBootstrap
+        from tsbootstrap import BaseTimeSeriesBootstrap, MovingBlockBootstrap
 
-        self.bootstrap_transformer_ = (
-            clone(bootstrap_transformer)
-            if bootstrap_transformer is not None
-            else MovingBlockBootstrap()
-        )
+        if isinstance(bootstrap_transformer, BaseTimeSeriesBootstrap):
path: sktime/forecasting/enbpi.py
text: @user1:
this inheritance check and import can be avoided by doing the type check, using a check for the tag `object_type == "bootstrap"`