Optimisation of historical forecast for regression models

CHANGELOG.md

@@ -14,6 +14,7 @@ but cannot always guarantee backwards compatibility. Changes that may **break co
 - General model improvements:
   - Added support for `PathLike` to the `save()` and `load()` functions of all non-deep learning based models. [#1754](https://github.com/unit8co/darts/pull/1754) by [Simon Sudrich](https://github.com/sudrich).
   - Improved efficiency of `historical_forecasts()` and `backtest()` for all models giving significant process time reduction for larger number of predict iterations and series. [#1801](https://github.com/unit8co/darts/pull/1801) by [Dennis Bader](https://github.com/dennisbader).
+  - Optimized `historical_forecasts()` for `RegressionModel` when `retrain=False` and `forecast_horizon <= output_chunk_length` by vectorizing the prediction. [#1885](https://github.com/unit8co/darts/pull/1885) by [Antoine Madrona](https://github.com/madtoinou).
   - Added model property `ForecastingModel.supports_multivariate` to indicate whether the model supports multivariate forecasting. [#1848](https://github.com/unit8co/darts/pull/1848) by [Felix Divo](https://github.com/felixdivo).
   - `Prophet` now supports conditional seasonalities, and properly handles all parameters passed to `Prophet.add_seasonality()` and model creation parameter `add_seasonalities` [#1829](https://github.com/unit8co/darts/pull/#1829) by [Idan Shilon](https://github.com/id5h).
   - Added support for direct prediction of the likelihood parameters to probabilistic models using a likelihood (regression and torch models). Set `predict_likelihood_parameters=True` when calling `predict()`. [#1811](https://github.com/unit8co/darts/pull/1811) by [Antoine Madrona](https://github.com/madtoinou).

darts/dataprocessing/encoders/encoder_base.py

@@ -178,7 +178,14 @@ def generate_train_inference_idx(
         inference_idx, _ = self.generate_inference_idx(
             n=n, target=target, covariates=covariates
         )
-        return train_idx.__class__.union(train_idx, inference_idx), target_end
+        # generate index end is inclusive, should not be a problem when taking union
+        gap = generate_index(
+            start=train_idx[-1], end=inference_idx[0] - target.freq, freq=target.freq
+        )
+        return (
+            train_idx.__class__.union(train_idx, gap).union(inference_idx),
+            target_end,
+        )
 
     @property
     @abstractmethod

darts/models/forecasting/regression_model.py

@@ -27,9 +27,10 @@
 if their static covariates do not have the same size, the shorter ones are padded with 0 valued features.
 """
 from collections import OrderedDict
-from typing import Any, Dict, List, Optional, Sequence, Tuple, Union
+from typing import Any, Callable, Dict, List, Optional, Sequence, Tuple, Union
 
 import numpy as np
+import pandas as pd
 from sklearn.linear_model import LinearRegression
 
 from darts.logging import get_logger, raise_if, raise_if_not, raise_log
@@ -40,6 +41,10 @@
     create_lagged_component_names,
     create_lagged_training_data,
 )
+from darts.utils.historical_forecasts import (
+    _optimized_historical_forecasts_regression_all_points,
+    _optimized_historical_forecasts_regression_last_points_only,
+)
 from darts.utils.multioutput import MultiOutputRegressor
 from darts.utils.utils import (
     _check_quantiles,
@@ -860,6 +865,121 @@ def supports_future_covariates(self) -> bool:
     def supports_static_covariates(self) -> bool:
         return True
 
+    @property
+    def supports_optimized_historical_forecasts(self) -> bool:
+        return True
+
+    def _check_optimizable_historical_forecasts(
+        self,
+        forecast_horizon: int,
+        retrain: Union[bool, int, Callable[..., bool]],
+        show_warnings=bool,
+    ) -> bool:
+        """
+        Historical forecast can be optimized only if `retrain=False` and `forecast_horizon <= self.output_chunk_length`
+        (no auto-regression required).
+        """
+
+        supported_retrain = (retrain is False) or (retrain == 0)
+        supported_forecast_horizon = forecast_horizon <= self.output_chunk_length
+        if supported_retrain and supported_forecast_horizon:
+            return True
+
+        if show_warnings:
+            if not supported_retrain:
+                logger.warning(
+                    "`enable_optimization=True` is ignored because `retrain` is not `False`"
+                    "To hide this warning, set `show_warnings=False` or `enable_optimization=False`."
+                )
+            if not supported_forecast_horizon:
+                logger.warning(
+                    "`enable_optimization=True` is ignored because "
+                    "`forecast_horizon > self.output_chunk_length`."
+                    "To hide this warning, set `show_warnings=False` or `enable_optimization=False`."
+                )
+
+        return False
+
+    def _optimized_historical_forecasts(
+        self,
+        series: Optional[Sequence[TimeSeries]],
+        past_covariates: Optional[Sequence[TimeSeries]] = None,
+        future_covariates: Optional[Sequence[TimeSeries]] = None,
+        num_samples: int = 1,
+        start: Optional[Union[pd.Timestamp, float, int]] = None,
+        forecast_horizon: int = 1,
+        stride: int = 1,
+        overlap_end: bool = False,
+        last_points_only: bool = True,
+        verbose: bool = False,
+        show_warnings: bool = True,
+        predict_likelihood_parameters: bool = False,
+    ) -> Union[
+        TimeSeries, List[TimeSeries], Sequence[TimeSeries], Sequence[List[TimeSeries]]
+    ]:
+        """
+        TODO: support forecast_horizon > output_chunk_length (auto-regression)
+        """
+        if not self._fit_called:
+            raise_log(
+                ValueError("Model has not been fit yet."),
+                logger,
+            )
+        if forecast_horizon > self.output_chunk_length:
+            raise_log(
+                ValueError(
+                    "`forecast_horizon > model.output_chunk_length` requires auto-regression which is not "
+                    "supported in this optimized routine."
+                ),
+                logger,
+            )
+
+        # manage covariates, usually handled by RegressionModel.predict()
+        if past_covariates is None and self.past_covariate_series is not None:
+            past_covariates = [self.past_covariate_series] * len(series)
+        if future_covariates is None and self.future_covariate_series is not None:
+            future_covariates = [self.future_covariate_series] * len(series)
+
+        self._verify_static_covariates(series[0].static_covariates)
+
+        if self.encoders.encoding_available:
+            past_covariates, future_covariates = self.generate_fit_predict_encodings(
+                n=forecast_horizon,
+                series=series,
+                past_covariates=past_covariates,
+                future_covariates=future_covariates,
+            )
+
+        # TODO: move the loop here instead of duplicated code in each sub-routine?
+        if last_points_only:
+            return _optimized_historical_forecasts_regression_last_points_only(
+                model=self,
+                series=series,
+                past_covariates=past_covariates,
+                future_covariates=future_covariates,
+                num_samples=num_samples,
+                start=start,
+                forecast_horizon=forecast_horizon,
+                stride=stride,
+                overlap_end=overlap_end,
+                show_warnings=show_warnings,
+                predict_likelihood_parameters=predict_likelihood_parameters,
+            )
+        else:
+            return _optimized_historical_forecasts_regression_all_points(
+                model=self,
+                series=series,
+                past_covariates=past_covariates,
+                future_covariates=future_covariates,
+                num_samples=num_samples,
+                start=start,
+                forecast_horizon=forecast_horizon,
+                stride=stride,
+                overlap_end=overlap_end,
+                show_warnings=show_warnings,
+                predict_likelihood_parameters=predict_likelihood_parameters,
+            )
+
 
 class _LikelihoodMixin:
     """

darts/tests/models/forecasting/test_historical_forecasts.py

@@ -1,4 +1,5 @@
 import unittest
+from typing import Union
 
 import numpy as np
 import pandas as pd
@@ -625,6 +626,156 @@ def test_regression_auto_start_multiple_no_cov(self):
                 f"Expected {theorical_forecast_length}, got {len(forecasts[0])} and {len(forecasts[1])}",
             )
 
+    @pytest.mark.slow
+    def test_optimized_historical_forecasts_regression(self):
+        start_ts = pd.Timestamp("2000-01-01")
+        ts_univariate = tg.linear_timeseries(
+            start_value=1, end_value=100, length=20, start=start_ts
+        )
+        ts_multivariate = ts_univariate.stack(
+            tg.sine_timeseries(length=20, start=start_ts)
+        )
+        # slightly longer to not affect the last predictable timestamp
+        ts_covs = tg.gaussian_timeseries(length=30, start=start_ts)
+        start = 14
+        model_cls = LinearRegressionModel
+        for ts in [ts_univariate, ts_multivariate]:
+            # cover several covariates combinations and several regression models
+            for _, model_kwargs, _ in (
+                models_reg_no_cov_cls_kwargs + models_reg_cov_cls_kwargs
+            ):
+                for multi_models in [True, False]:
+                    for forecast_horizon in [1, 5]:
+                        # ocl == forecast horizon
+                        model_kwargs_same = model_kwargs.copy()
+                        model_kwargs_same["output_chunk_length"] = forecast_horizon
+                        model_kwargs_same["multi_models"] = multi_models
+                        model_same = model_cls(**model_kwargs_same)
+                        model_same.fit(
+                            series=ts[:start],
+                            past_covariates=ts_covs
+                            if model_same.supports_past_covariates
+                            else None,
+                            future_covariates=ts_covs
+                            if model_same.supports_future_covariates
+                            else None,
+                        )
+                        # ocl >= forecast horizon
+                        model_kwargs_diff = model_kwargs.copy()
+                        model_kwargs_diff["output_chunk_length"] = 5
+                        model_kwargs_diff["multi_models"] = multi_models
+                        model_diff = model_cls(**model_kwargs_same)
+                        model_diff.fit(
+                            series=ts[:start],
+                            past_covariates=ts_covs
+                            if model_diff.supports_past_covariates
+                            else None,
+                            future_covariates=ts_covs
+                            if model_diff.supports_future_covariates
+                            else None,
+                        )
+                        for model in [model_same, model_diff]:
+                            for last_points_only in [True, False]:
+                                for stride in [1, 2]:
+                                    hist_fct = model.historical_forecasts(
+                                        series=ts,
+                                        past_covariates=ts_covs
+                                        if model.supports_past_covariates
+                                        else None,
+                                        future_covariates=ts_covs
+                                        if model.supports_future_covariates
+                                        else None,
+                                        start=start,
+                                        retrain=False,
+                                        last_points_only=last_points_only,
+                                        stride=stride,
+                                        forecast_horizon=forecast_horizon,
+                                        enable_optimization=False,
+                                    )
+
+                                    # manually packing the series in list to match expected inputs
+                                    opti_hist_fct = (
+                                        model._optimized_historical_forecasts(
+                                            series=[ts],
+                                            past_covariates=[ts_covs]
+                                            if model.supports_past_covariates
+                                            else None,
+                                            future_covariates=[ts_covs]
+                                            if model.supports_future_covariates
+                                            else None,
+                                            start=start,
+                                            last_points_only=last_points_only,
+                                            stride=stride,
+                                            forecast_horizon=forecast_horizon,
+                                        )
+                                    )
+                                    # pack the output to generalize the tests
+                                    if last_points_only:
+                                        hist_fct = [hist_fct]
+                                        opti_hist_fct = [opti_hist_fct]
+
+                                    for fct, opti_fct in zip(hist_fct, opti_hist_fct):
+                                        self.assertTrue(
+                                            (
+                                                fct.time_index == opti_fct.time_index
+                                            ).all()
+                                        )
+                                        np.testing.assert_array_almost_equal(
+                                            fct.all_values(), opti_fct.all_values()
+                                        )
+
+    def test_optimized_historical_forecasts_regression_with_encoders(self):
+        for use_covs in [False, True]:
+            series_train, series_val = self.ts_pass_train, self.ts_pass_val
+            model = LinearRegressionModel(
+                lags=3,
+                lags_past_covariates=2,
+                lags_future_covariates=[2, 3],
+                add_encoders={
+                    "cyclic": {"future": ["month"]},
+                    "datetime_attribute": {"past": ["dayofweek"]},
+                },
+                output_chunk_length=5,
+            )
+            if use_covs:
+                pc = tg.gaussian_timeseries(
+                    start=series_train.start_time() - 2 * series_train.freq,
+                    end=series_val.end_time(),
+                    freq=series_train.freq,
+                )
+                fc = tg.gaussian_timeseries(
+                    start=series_train.start_time() + 3 * series_train.freq,
+                    end=series_val.end_time() + 4 * series_train.freq,
+                    freq=series_train.freq,
+                )
+            else:
+                pc, fc = None, None
+
+            model.fit(self.ts_pass_train, past_covariates=pc, future_covariates=fc)
+
+            hist_fct = model.historical_forecasts(
+                series=self.ts_pass_val,
+                past_covariates=pc,
+                future_covariates=fc,
+                retrain=False,
+                last_points_only=True,
+                forecast_horizon=5,
+                enable_optimization=False,
+            )
+
+            opti_hist_fct = model._optimized_historical_forecasts(
+                series=[self.ts_pass_val],
+                past_covariates=[pc],
+                future_covariates=[fc],
+                last_points_only=True,
+                forecast_horizon=5,
+            )
+
+            self.assertTrue((hist_fct.time_index == opti_hist_fct.time_index).all())
+            np.testing.assert_array_almost_equal(
+                hist_fct.all_values(), opti_hist_fct.all_values()
+            )
+
     @pytest.mark.slow
     @unittest.skipUnless(
         TORCH_AVAILABLE,
@@ -1221,11 +1372,13 @@ def retrain_f_delayed_true(
 
         # test int
         helper_hist_forecasts(10, 0.9)
-        expected_msg = "Model has not been fit before the first predict iteration at prediction point (in time)"
+        expected_msg = "Model has not been fit yet."
         # `retrain=0` with not-trained model, encountering directly a predictable time index
         with pytest.raises(ValueError) as error_msg:
             helper_hist_forecasts(0, 0.9)
-        self.assertTrue(str(error_msg.value).startswith(expected_msg))
+        self.assertTrue(
+            str(error_msg.value).startswith(expected_msg), str(error_msg.value)
+        )
 
         # test bool
         helper_hist_forecasts(True, 0.9)
@@ -1251,7 +1404,9 @@ def test_predict_likelihood_parameters(self):
         """standard checks that historical forecasts work with direct likelihood parameter predictions
         with regression and torch models."""
 
-        def create_model(ocl, use_ll=True, model_type="regression"):
+        def create_model(
+            ocl, use_ll=True, model_type="regression"
+        ) -> Union[LinearRegressionModel, NLinearModel]:
             if model_type == "regression":
                 return LinearRegressionModel(
                     lags=3,

darts/utils/__init__.py

@@ -4,7 +4,6 @@
 """
 from .utils import (
     _build_tqdm_iterator,
-    _historical_forecasts_general_checks,
     _parallel_apply,
     _with_sanity_checks,
     retain_period_common_to_all,

darts/utils/historical_forecasts/__init__.py

@@ -0,0 +1,9 @@
+from .optimized_historical_forecasts import (
+    _optimized_historical_forecasts_regression_all_points,
+    _optimized_historical_forecasts_regression_last_points_only,
+)
+from .utils import (
+    _get_historical_forecast_boundaries,
+    _historical_forecasts_general_checks,
+    _historical_forecasts_start_warnings,
+)