Repository: sktime/sktime

Issue:
[ENH] Let HEF fit and predict in parallel
**Is your feature request related to a problem? Please describe.**
As it is the HierarchyEnsembleForecaster fits every model one at a time which may result in long training times if the hierarchy becomes sufficiently large. 

**Describe the solution you'd like**
Instead enable the option to use the sktime backend to let the forecasters be fit in parallel (also do that for the prediction since we are already at it) 

**Describe alternatives you've considered**
I suppose instead of creating the HEF with new forecasters one could use any already existing framework that allows us to train the forecasters in parallel and then skip the HEF training by setting 

HEF.fitted_list = [already trained forecasters]

**Additional context**
The idea to solve this is to replace the for-loops in both the by=level and by=node with something like this 

    fcstr_list = parallelize(
    node_fit,
    node_dict.items(),
    meta=meta,
    backend=backend,
    backend_params=backend_params,
    )
where node_fit would then just be a function to wrap around the forecasters.fit lines (since we cant parallelize over that directly). This could already look like 

      def node_fit(params, meta):
  
      key, node = params
      z = meta["z"]
      X = meta["X"]
  
      frcstr = meta["fcstr_dict"][key]
      df = z[z.index.droplevel(-1).isin(node)]
      if X is not None:
          X = X.loc[df.index]
      frcstr.fit(df, fh=meta["fh"], X=X)
      return frcstr, df




Discussion:
LHoelper: > Looks great!
> 
> I would also be interested to hear from @VyomkeshVyas and @felipeangelimvieira.
> 
> Minor, non-blocking style point: I would suggest to move the new functions to the bottom of the file, to make it easier to read.

Will do. One question to that (since this is my first PR). Does one normally wait until the running checks are finished to push a change like this or do i just fire them as fast as possible? 
fkiraly: > or do i just fire them as fast as possible?

If you push the changes, then the jobs will be aborted - it is fine to do so.

Generally, best practice is to push the newest state when you have it - we squash PRs, so no need to minimize the number of commits.

There is no requirement for one or the other option.

Welcome to `sktime` btw, great contribution!

Patch:
diff --git a/.all-contributorsrc b/.all-contributorsrc
index ecf627af2..59b7d9938 100644
--- a/.all-contributorsrc
+++ b/.all-contributorsrc
@@ -3414,6 +3414,15 @@
         "code"
       ]
     },
+    {
+      "login": "LHoelper",
+      "name": "Luca Hölper",
+      "avatar_url": "https://avatars.githubusercontent.com/LHoelper",
+      "profile": "https://github.com/LHoelper",
+      "contributions": [
+        "code"
+      ]
+    },
     {
       "login": "sky0walker99",
       "name": "Muhammed Haroon",
diff --git a/sktime/forecasting/base/_meta.py b/sktime/forecasting/base/_meta.py
index 556b5f31a..14cf92ea1 100644
--- a/sktime/forecasting/base/_meta.py
+++ b/sktime/forecasting/base/_meta.py
@@ -92,22 +92,43 @@ class _HeterogenousEnsembleForecaster(_HeterogenousMetaEstimator, BaseForecaster
 
         return estimator_tuples
 
+    def _get_forecaster_list(self):
+        """Return list of forecasters."""
+        return [x[1] for x in self.forecasters_]
+
+    def _get_forecaster_names(self):
+        """Return list of forecaster names."""
+        return [x[0] for x in self.forecasters_]
+
     def _fit_forecasters(self, forecasters, y, X, fh):
-        """Fit all forecasters in parallel."""
+        """Fit all forecasters in parallel.
+
+        Returns
+        -------
+        list of references to fitted forecasters
+            in same order as forecasters
+        """
         from joblib import Parallel, delayed
 
         def _fit_forecaster(forecaster, y, X, fh):
             """Fit single forecaster."""
             return forecaster.fit(y, X, fh)
 
-        self.forecasters_ = Parallel(n_jobs=self.n_jobs)(
+        if forecasters is None:
+            forecasters = self._get_forecaster_list()
+
+        fitted_fcst = Parallel(n_jobs=self.n_jobs)(
             delayed(_fit_forecaster)(forecaster.clone(), y, X, fh)
             for forecaster in forecasters
         )
+        fcst_names = self._get_forecaster_names()
+        self.forecasters_ = list(zip(fcst_names, fitted_fcst))
 
-    def _predict_forecasters(self, fh=None, X=None):
+    def _predict_forecasters(self, fh=None, X=None, forecasters=None):
         """Collect results from forecaster.predict() calls."""
-        return [forecaster.predict(fh=fh, X=X) for forecaster in self.forecasters_]
+        if forecasters is None:
+            forecasters = self._get_forecaster_list()
+        return [forecaster.predict(fh=fh, X=X) for forecaster in forecasters]
 
     def _update(self, y, X=None, update_params=True):
         """Update fitted parameters.
@@ -122,6 +143,6 @@ class _HeterogenousEnsembleForecaster(_HeterogenousMetaEstimator, BaseForecaster
         -------
         self : an instance of self.
         """
-        for forecaster in self.forecasters_:
+        for forecaster in self._get_forecaster_list():
             forecaster.update(y, X, update_params=update_params)
         return self
diff --git a/sktime/forecasting/compose/_ensemble.py b/sktime/forecasting/compose/_ensemble.py
index fa7ff46ad..7b999f59b 100644
--- a/sktime/forecasting/compose/_ensemble.py
+++ b/sktime/forecasting/compose/_ensemble.py
@@ -328,18 +328,6 @@ class EnsembleForecaster(_HeterogenousEnsembleForecaster):
         "scitype:y": "both",
     }
 
-    # for default get_params/set_params from _HeterogenousMetaEstimator
-    # _steps_attr points to the attribute of self
-    # which contains the heterogeneous set of estimators
-    # this must be an iterable of (name: str, estimator, ...) tuples for the default
-    _steps_attr = "_forecasters"
-
-    # if the estimator is fittable, _HeterogenousMetaEstimator also
-    # provides an override for get_fitted_params for params from the fitted estimators
-    # the fitted estimators should be in a different attribute, _steps_fitted_attr
-    # this must be an iterable of (name: str, estimator, ...) tuples for the default
-    _steps_fitted_attr = "forecasters_"
-
     def __init__(self, forecasters, n_jobs=None, aggfunc="mean", weights=None):
         self.aggfunc = aggfunc
         self.weights = weights
@@ -391,8 +379,7 @@ class EnsembleForecaster(_HeterogenousEnsembleForecaster):
         -------
         self : returns an instance of self.
         """
-        forecasters = [f[1] for f in self._forecasters]
-        self._fit_forecasters(forecasters, y, X, fh)
+        self._fit_forecasters(None, y, X, fh)
         return self
 
     def _predict(self, fh, X):
diff --git a/sktime/forecasting/compose/_hierarchy_ensemble.py b/sktime/forecasting/compose/_hierarchy_ensemble.py
index 457c1b71f..cf7fcb622 100644
--- a/sktime/forecasting/compose/_hierarchy_ensemble.py
+++ b/sktime/forecasting/compose/_hierarchy_ensemble.py
@@ -11,6 +11,7 @@ from sktime.base._meta import flatten
 from sktime.forecasting.base._base import BaseForecaster
 from sktime.forecasting.base._meta import _HeterogenousEnsembleForecaster
 from sktime.transformations.hierarchical.aggregate import _check_index_no_total
+from sktime.utils.parallel import parallelize
 from sktime.utils.warnings import warn
 
 
@@ -54,6 +55,39 @@ class HierarchyEnsembleForecaster(_HeterogenousEnsembleForecaster):
         if passed, applies ``default`` forecaster on the nodes/levels
         not mentioned in the ``forecaster`` argument.
 
+    backend : {"dask", "loky", "multiprocessing", "threading"}, by default None.
+        Runs parallel evaluate if specified and ``strategy`` is set as "refit".
+
+        - "None": executes loop sequentally, simple list comprehension
+        - "loky", "multiprocessing" and "threading": uses ``joblib.Parallel`` loops
+        - "joblib": custom and 3rd party ``joblib`` backends, e.g., ``spark``
+        - "dask": uses ``dask``, requires ``dask`` package in environment
+        - "dask_lazy": same as "dask",
+          but changes the return to (lazy) ``dask.dataframe.DataFrame``.
+
+        Recommendation: Use "dask" or "loky" for parallel evaluate.
+        "threading" is unlikely to see speed ups due to the GIL and the serialization
+        backend (``cloudpickle``) for "dask" and "loky" is generally more robust
+        than the standard ``pickle`` library used in "multiprocessing".
+    backend_params : dict, optional
+        additional parameters passed to the backend as config.
+        Directly passed to ``utils.parallel.parallelize``.
+        Valid keys depend on the value of ``backend``:
+
+        - "None": no additional parameters, ``backend_params`` is ignored
+        - "loky", "multiprocessing" and "threading": default ``joblib`` backends
+          any valid keys for ``joblib.Parallel`` can be passed here, e.g., ``n_jobs``,
+          with the exception of ``backend`` which is directly controlled by ``backend``.
+          If ``n_jobs`` is not passed, it will default to ``-1``, other parameters
+          will default to ``joblib`` defaults.
+        - "joblib": custom and 3rd party ``joblib`` backends, e.g., ``spark``.
+          any valid keys for ``joblib.Parallel`` can be passed here, e.g., ``n_jobs``,
+          ``backend`` must be passed as a key of ``backend_params`` in this case.
+          If ``n_jobs`` is not passed, it will default to ``-1``, other parameters
+          will default to ``joblib`` defaults.
+        - "dask": any valid keys for ``dask.compute`` can be passed,
+          e.g., ``scheduler``
+
     Examples
     --------
     >>> from sktime.forecasting.compose import HierarchyEnsembleForecaster
@@ -109,11 +143,15 @@ class HierarchyEnsembleForecaster(_HeterogenousEnsembleForecaster):
 
     _steps_attr = "_forecasters"
 
-    def __init__(self, forecasters, by="level", default=None):
+    def __init__(
+        self, forecasters, by="level", default=None, backend=None, backend_params=None
+    ):
         self.forecasters = forecasters
         self.by = by
         self.default = default
 
+        self.backend = backend
+        self.backend_params = backend_params
         super().__init__(forecasters=None)
 
         if isinstance(forecasters, BaseForecaster):
@@ -220,31 +258,40 @@ class HierarchyEnsembleForecaster(_HeterogenousEnsembleForecaster):
             self.fitted_list_.append([frcstr, y.index])
             return self
 
+        meta = {"X": x, "z": z, "fh": fh}
+
         if self.by == "level":
             hier_dict = self._get_hier_dict(z)
-            for _, forecaster, level in self.forecasters_:
-                if level in hier_dict.keys():
-                    frcstr = forecaster
-                    df = z[z.index.droplevel(-1).isin(hier_dict[level])]
-                    if X is not None:
-                        x = X.loc[df.index]
-                    frcstr.fit(df, fh=fh, X=x)
-                    self.fitted_list_.append([frcstr, df.index.droplevel(-1).unique()])
-                    self.forecasters_ = [
-                        (name, frcstr if f == forecaster else f, level)
-                        for name, f, level in self.forecasters_
-                    ]
+            meta["hier_dict"] = hier_dict
+
+            fcstr_list = parallelize(
+                _level_fit,
+                iter=self.forecasters_,
+                meta=meta,
+                backend=self.backend,
+                backend_params=self.backend_params,
+            )
 
+            for i in range(len(fcstr_list)):
+                self.fitted_list_.append(
+                    [fcstr_list[i][0], fcstr_list[i][1].index.droplevel(-1).unique()]
+                )
+                name, _, level = self.forecasters_[i]
+                self.forecasters_[i] = (name, fcstr_list[i][0], level)
         else:
             node_dict, frcstr_dict = self._get_node_dict(z)
-            for key, nodes in node_dict.items():
-                frcstr = frcstr_dict[key]
-                df = z[z.index.droplevel(-1).isin(nodes)]
-                if X is not None:
-                    x = X.loc[df.index]
-                frcstr.fit(df, fh=fh, X=x)
-                self.fitted_list_.append([frcstr, df.index.droplevel(-1).unique()])
+            meta["fcstr_dict"] = frcstr_dict
+
+            fcstr_list = parallelize(
+                _node_fit,
+                node_dict.items(),
+                meta=meta,
+                backend=self.backend,
+                backend_params=self.backend_params,
+            )
 
+            for fcstr, df in fcstr_list:
+                self.fitted_list_.append([fcstr, df.index.droplevel(-1).unique()])
         return self
 
     def _get_hier_dict(self, z):
@@ -405,18 +452,20 @@ class HierarchyEnsembleForecaster(_HeterogenousEnsembleForecaster):
         y_pred : pd.Series
             Point predictions
         """
-        preds = []
-
         if X is not None:
             if _check_index_no_total(X):
                 X = self._aggregate(X)
         x = X
 
-        for forecaster, ind in self.fitted_list:
-            if X is not None and X.index.nlevels > 1:
-                x = X[X.index.droplevel(-1).isin(ind)]
-            pred = forecaster.predict(fh=fh, X=x)
-            preds.append(pred)
+        meta = {"x": x, "fh": fh}
+
+        preds = parallelize(
+            _predict_one_forecaster,
+            self.fitted_list,
+            meta,
+            backend=self.backend,
+            backend_params=self.backend_params,
+        )
 
         preds = pd.concat(preds, axis=0)
         preds.sort_index(inplace=True)
@@ -674,3 +723,98 @@ class HierarchyEnsembleForecaster(_HeterogenousEnsembleForecaster):
         params3 = {"forecasters": NaiveForecaster()}
 
         return [params1, params2, params3]
+
+
+def _level_fit(params, meta):
+    """Fit a single forecaster.
+
+    Called from _fit if by=level was set to allow for parallelization.
+
+    Parameters
+    ----------
+    params: tuples
+        (str, estimator, int or list of tuple/s)
+        If called from the _fit function this will contain the self.forecasters_ list
+    meta: dict
+        Meta dictionary for parallelization. Needs to contain (at least)
+        the following keys:
+        z: pd.Series
+            The aggregated target time series
+        X: pd.DataFrame, None
+            The aggregated input data
+        hier_dict: dict
+            The level dictionary as created by the get_hier_dict function
+        fh : int, list or np.array, optional (default=None)
+            The forecasters horizon with the steps ahead to to predict.
+    """
+    _, forecaster, level = params
+    z = meta["z"]
+    X = meta["X"]
+    hier_dict = meta["hier_dict"]
+    if level in hier_dict.keys():
+        frcstr = forecaster
+        df = z[z.index.droplevel(-1).isin(hier_dict[level])]
+        if X is not None:
+            X = X.loc[df.index]
+        frcstr.fit(df, fh=meta["fh"], X=X)
+    return frcstr, df
+
+
+def _node_fit(params, meta):
+    """Fit a single forecaster.
+
+    Called from _fit if by=node was set to allow for parallelization.
+
+    Parameters
+    ----------
+    params: tuples
+        (str, estimator, int or list of tuple/s)
+        If called from the _fit function this will contain the node_dict from the
+        get_node_dict function
+    meta: dict
+        Meta dictionary for parallelization. Needs to contain (at least) the
+        following keys:
+        z: pd.Series
+            The aggregated target time series
+        X: pd.DataFrame, None
+            The aggregated input data
+        fcstr_dict: dict
+            The forecaster dictionary as created by the get_node_dict function
+        fh : int, list or np.array, optional (default=None)
+            The forecasters horizon with the steps ahead to to predict.
+    """
+    key, node = params
+    z = meta["z"]
+    X = meta["X"]
+    frcstr = meta["fcstr_dict"][key]
+    df = z[z.index.droplevel(-1).isin(node)]
+    if X is not None:
+        X = X.loc[df.index]
+    frcstr.fit(df, fh=meta["fh"], X=X)
+    return frcstr, df
+
+
+def _predict_one_forecaster(params, meta):
+    """Predict a single forecaster.
+
+    Called from _predict to allow for parallelization.
+
+    Parameters
+    ----------
+    params: tuples
+        (estimator, str)
+        If called from the _predict function this will contain the self.fitted list
+    meta: dict
+        Meta dictionary for parallelization. Needs to contain (at least)
+        the following keys:
+        X: pd.DataFrame, None
+            The aggregated input data
+        fh : int, list or np.array, optional (default=None)
+            The forecasters horizon with the steps ahead to to predict.
+    """
+    X = meta["x"]
+    fh = meta["fh"]
+    if X is not None and X.index.nlevels > 1:
+        X = X[X.index.droplevel(-1).isin(params[1])]
+    pred = params[0].predict(fh=fh, X=X)
+    return pred
diff --git a/sktime/forecasting/compose/_stack.py b/sktime/forecasting/compose/_stack.py
index c87f990e1..e5cc68a35 100644
--- a/sktime/forecasting/compose/_stack.py
+++ b/sktime/forecasting/compose/_stack.py
@@ -143,7 +143,7 @@ class StackingForecaster(_HeterogenousEnsembleForecaster):
         """
         if update_params:
             warn("Updating `final regressor is not implemented", obj=self)
-        for forecaster in self.forecasters_:
+        for forecaster in self._get_forecaster_list():
             forecaster.update(y, X, update_params=update_params)
         return self
 
diff --git a/sktime/forecasting/online_learning/_online_ensemble.py b/sktime/forecasting/online_learning/_online_ensemble.py
index 91c8aef74..9e8721f3e 100644
--- a/sktime/forecasting/online_learning/_online_ensemble.py
+++ b/sktime/forecasting/online_learning/_online_ensemble.py
@@ -107,7 +107,7 @@ class OnlineEnsembleForecaster(EnsembleForecaster):
         if len(y) >= 1 and self.ensemble_algorithm is not None:
             self._fit_ensemble(y, X)
 
-        for forecaster in self.forecasters_:
+        for forecaster in self._get_forecaster_list():
             forecaster.update(y, X, update_params=update_params)
 
         return self


Test patch:
diff --git a/sktime/forecasting/compose/tests/test_hierarchy_ensemble.py b/sktime/forecasting/compose/tests/test_hierarchy_ensemble.py
index 760818b69..bc3d2342d 100644
--- a/sktime/forecasting/compose/tests/test_hierarchy_ensemble.py
+++ b/sktime/forecasting/compose/tests/test_hierarchy_ensemble.py
@@ -63,6 +63,109 @@ def test_hierarchy_ensemble_level_predict(forecasters):
     assert np.all(actual_pred.loc[def_pred.index] == def_pred), msg
 
 
+@pytest.mark.skipif(
+    not run_test_for_class(HierarchyEnsembleForecaster),
+    reason="run test only if softdeps are present and incrementally (if requested)",
+)
+@pytest.mark.parametrize(
+    "forecasters",
+    [
+        [("ptf", PolynomialTrendForecaster(), 0), ("naive", NaiveForecaster(), 1)],
+        [("naive", NaiveForecaster(), 0)],
+    ],
+)
+def test_hierarchy_ensemble_level_predict_parallel(forecasters):
+    """Check the level predictions."""
+    agg = Aggregator()
+
+    y = _bottom_hier_datagen(
+        no_bottom_nodes=7,
+        no_levels=2,
+        random_seed=123,
+    )
+
+    forecaster = HierarchyEnsembleForecaster(
+        forecasters,
+        default=forecasters[0][1].clone(),
+        by="level",
+        backend="joblib",
+        backend_params={"backend": "loky", "n_jobs": -1},
+    )
+
+    forecaster.fit(y, fh=[1, 2, 3])
+    actual_pred = forecaster.predict()
+
+    y = agg.fit_transform(y)
+
+    for i in range(len(forecasters)):
+        test_frcstr = forecasters[i][1].clone()
+        df = y[y.index.droplevel(-1).isin(forecaster.fitted_list[i][1])]
+        test_frcstr.fit(df, fh=[1, 2, 3])
+        test_pred = test_frcstr.predict()
+        msg = "Level predictions do not match"
+        assert np.all(actual_pred.loc[test_pred.index] == test_pred), msg
+
+    def_frcstr = forecasters[0][1].clone()
+    df = y[y.index.droplevel(-1).isin(forecaster.fitted_list[-1][1])]
+    def_frcstr.fit(df, fh=[1, 2, 3])
+    def_pred = def_frcstr.predict()
+    msg = "Level default predictions do not match"
+    assert np.all(actual_pred.loc[def_pred.index] == def_pred), msg
+
+
+@pytest.mark.skipif(
+    not run_test_for_class(HierarchyEnsembleForecaster),
+    reason="run test only if softdeps are present and incrementally (if requested)",
+)
+@pytest.mark.parametrize(
+    "forecasters",
+    [
+        [
+            ("ptf", PolynomialTrendForecaster(), [("__total", "__total")]),
+            ("naive", NaiveForecaster(), [("l2_node01", "__total")]),
+        ],
+        [("naive", NaiveForecaster(), [("l2_node01", "l1_node06")])],
+    ],
+)
+def test_hierarchy_ensemble_node_predict_parallel(forecasters):
+    """Check the node predictions."""
+    agg = Aggregator()
+
+    y = _bottom_hier_datagen(
+        no_bottom_nodes=7,
+        no_levels=2,
+        random_seed=123,
+    )
+
+    forecaster = HierarchyEnsembleForecaster(
+        forecasters,
+        by="node",
+        default=forecasters[0][1].clone(),
+        backend="joblib",
+        backend_params={"backend": "loky", "n_jobs": -1},
+    )
+
+    forecaster.fit(y, fh=[1, 2, 3])
+    actual_pred = forecaster.predict()
+
+    y = agg.fit_transform(y)
+
+    for i in range(len(forecasters)):
+        test_frcstr = forecasters[i][1].clone()
+        df = y[y.index.droplevel(-1).isin(forecaster.fitted_list[i][1])]
+        test_frcstr.fit(df, fh=[1, 2, 3])
+        test_pred = test_frcstr.predict()
+        msg = "Node predictions do not match"
+        assert np.all(actual_pred.loc[test_pred.index] == test_pred), msg
+
+    def_frcstr = forecasters[0][1].clone()
+    df = y[y.index.droplevel(-1).isin(forecaster.fitted_list[-1][1])]
+    def_frcstr.fit(df, fh=[1, 2, 3])
+    def_pred = def_frcstr.predict()
+    msg = "Node default predictions do not match"
+    assert np.all(actual_pred.loc[def_pred.index] == def_pred), msg
+
+
 @pytest.mark.skipif(
     not run_test_for_class(HierarchyEnsembleForecaster),
     reason="run test only if softdeps are present and incrementally (if requested)",