Add a new forecasting pipeline that does scale - forecast - unscale

CHANGELOG.md

@@ -15,6 +15,11 @@ Changelogs for this project are recorded in this file since v0.2.0.
 
 * soft-dtw related tools now support Itakura and Sakoe-Chiba global constraints. ([#189](https://github.com/tslearn-team/tslearn/issues/189))
 * Multithreading support added to `cdist_soft_dtw` and `softdtw_barycenter`. ([#310](https://github.com/tslearn-team/tslearn/issues/310)) 
+* `TimeSeriesScalerMinMax` and `TimeSeriesScalerMeanVariance` now offer an inverse_transform method 
+to reverse the normalization. ([#697](https://github.com/tslearn-team/tslearn/issues/697))
+* `tslearn.forecasting.ScaledForecastingPipeline` wraps a forecaster so that it is fitted on scaled data and 
+its predictions are un-scaled back, using any `tslearn.preprocessing` scaler, including per-series ones 
+(`per_timeseries=True`). ([#708](https://github.com/tslearn-team/tslearn/issues/708))
 
 ### Removed
 

docs/gen_modules/tslearn.forecasting.rst

@@ -13,3 +13,4 @@ tslearn.forecasting
 
       VARIMA
       AutoVARIMA
+      ScaledForecastingPipeline

tests/test_forecasting.py

@@ -2,8 +2,15 @@
 
 import pytest
 
+from sklearn.base import clone
+
 from tslearn.generators import random_walks
-from tslearn.forecasting import VARIMA, AutoVARIMA
+from tslearn.forecasting import VARIMA, AutoVARIMA, ScaledForecastingPipeline
+from tslearn.preprocessing import (
+    TimeSeriesScalerMeanVariance,
+    TimeSeriesScalerMinMax,
+    TimeSeriesResampler,
+)
 
 
 def test_VARIMA():
@@ -207,4 +214,100 @@ def test_verbosity(capteesys):
 
     AutoVARIMA(max_d=0, max_iter=2, verbose=1).fit(data)
     captured = capteesys.readouterr()
-    assert "Default d for non stationarity 0 is used." in captured.out
+    assert "Default d for non stationarity 0 is used." in captured.out
+
+
+@pytest.mark.parametrize(
+    "scaler",
+    [
+        TimeSeriesScalerMeanVariance(per_timeseries=False),
+        TimeSeriesScalerMinMax(per_timeseries=False),
+    ]
+)
+def test_scaled_forecaster_global_matches_manual_scaling(scaler):
+    data = random_walks(n_ts=5, sz=20, d=2, random_state=0) * 100 + 500
+
+    pipeline = ScaledForecastingPipeline(VARIMA(1, 0, 0), scaler=scaler).fit(data)
+    assert pipeline.per_timeseries_ is False
+    predicted = pipeline.predict(n=3)
+    assert predicted.shape == (5, 3, 2)
+
+    reference_scaler = clone(scaler)
+    scaled_data = reference_scaler.fit_transform(data)
+    reference_predicted = VARIMA(1, 0, 0).fit(scaled_data).predict(n=3)
+    np.testing.assert_allclose(
+        predicted,
+        reference_scaler.inverse_transform(reference_predicted)
+    )
+
+
+@pytest.mark.parametrize(
+    "scaler",
+    [None, TimeSeriesScalerMeanVariance(), TimeSeriesScalerMinMax()]
+)
+def test_scaled_forecaster_per_series_matches_manual_scaling(scaler):
+    # `per_timeseries=True` is the default of tslearn scalers, so both the
+    # default `scaler=None` and an explicit scaler exercise per-series mode.
+    data = random_walks(n_ts=5, sz=20, d=2, random_state=0)
+    # Give series wildly different levels/spreads, the case per-series
+    # scaling is meant for.
+    data = data * np.array([1, 100, 0.01, 10, 1000]).reshape(5, 1, 1)
+
+    pipeline = ScaledForecastingPipeline(VARIMA(1, 0, 0), scaler=scaler).fit(data)
+    assert pipeline.per_timeseries_ is True
+    assert len(pipeline.scalers_) == 5
+    predicted = pipeline.predict(n=3)
+    assert predicted.shape == (5, 3, 2)
+
+    template = pipeline._scaler_template_
+    scalers = [clone(template).fit(data[i:i + 1]) for i in range(5)]
+    scaled_data = np.concatenate(
+        [s.transform(data[i:i + 1]) for i, s in enumerate(scalers)], axis=0
+    )
+    reference_predicted_scaled = VARIMA(1, 0, 0).fit(scaled_data).predict(n=3)
+    reference_predicted = np.concatenate(
+        [
+            s.inverse_transform(reference_predicted_scaled[i:i + 1])
+            for i, s in enumerate(scalers)
+        ],
+        axis=0,
+    )
+    np.testing.assert_allclose(predicted, reference_predicted)
+
+
+def test_scaled_forecaster_per_series_predict_uses_fresh_statistics():
+    data = random_walks(n_ts=3, sz=20, d=1, random_state=0)
+    data = data * np.array([1, 100, 0.01]).reshape(3, 1, 1)
+    pipeline = ScaledForecastingPipeline(VARIMA(1, 0, 0)).fit(data)
+
+    new_data = random_walks(n_ts=3, sz=15, d=1, random_state=1)
+    new_data = new_data * np.array([5, 2, 50]).reshape(3, 1, 1) + 3
+    predicted_new = pipeline.predict(new_data, n=2)
+    assert predicted_new.shape == (3, 2, 1)
+
+    # Statistics used for the new prediction differ from the fit-time ones,
+    # since per-series scaling is recomputed from whatever is transformed.
+    template = pipeline._scaler_template_
+    fresh_scalers = [
+        clone(template).fit(new_data[i:i + 1]) for i in range(3)
+    ]
+    for fitted, fresh in zip(pipeline.scalers_, fresh_scalers):
+        assert not np.allclose(fitted.mean_, fresh.mean_)
+
+
+def test_scaled_forecaster_fit_predict():
+    data = random_walks(n_ts=3, sz=20, d=1, random_state=0) * 10 + 50
+    pipeline = ScaledForecastingPipeline(VARIMA(1, 0, 0))
+    np.testing.assert_allclose(
+        pipeline.fit_predict(data, n=3),
+        pipeline.fit(data).predict(n=3)
+    )
+
+
+def test_scaled_forecaster_rejects_scaler_without_inverse_transform():
+    data = random_walks(n_ts=3, sz=20, random_state=0)
+    with pytest.raises(ValueError):
+        ScaledForecastingPipeline(
+            VARIMA(1, 0, 0),
+            scaler=TimeSeriesResampler(sz=20)
+        ).fit(data)

tests/test_preprocessing.py

@@ -2,10 +2,10 @@
 
 import pytest
 
+from tslearn.generators import random_walks
 from tslearn.preprocessing import (TimeSeriesScalerMeanVariance,
                                    TimeSeriesScalerMinMax,
                                    TimeSeriesImputer)
-
 from tslearn.utils import to_time_series_dataset, to_time_series
 
 
@@ -19,6 +19,11 @@ def test_single_value_ts_no_nan():
     assert np.sum(np.isnan(minmax_scaler.fit_transform(X))) == 0
 
 
+def test_min_max_scaler_range():
+    with pytest.raises(ValueError):
+        TimeSeriesScalerMinMax((1., 0.)).fit_transform([[1, 2, 3]])
+
+
 def test_min_max_scaler_variable_length():
     X = [
         [1, np.nan],
@@ -40,6 +45,30 @@ def test_min_max_scaler_variable_length():
         np.array([[[0], [0.5], [1]]])
     )
 
+
+@pytest.mark.parametrize(
+    "scaler",
+    [TimeSeriesScalerMinMax, TimeSeriesScalerMeanVariance]
+)
+@pytest.mark.parametrize(
+    "per_timeseries, per_feature",
+    [(True, True), (True, False), (False, True), (False, False)]
+)
+def test_scaler_inverse_transform(scaler, per_timeseries, per_feature):
+    X = random_walks(10, 10, 2, mu=1, random_state=0)
+
+    estimator = scaler(per_timeseries=per_timeseries, per_feature=per_feature)
+    transformed = estimator.fit_transform(X)
+    if per_timeseries:
+        with pytest.raises(RuntimeError):
+            estimator.inverse_transform(X)
+    else:
+        np.testing.assert_array_almost_equal(
+        estimator.inverse_transform(transformed),
+        X
+    )
+
+
 def test_min_max_scaler_modes():
     univariate_dataset = [
         [1, 2, 3],

tslearn/forecasting/__init__.py

@@ -5,5 +5,6 @@
 """
 
 from ._arima import VARIMA, AutoVARIMA
+from ._pipeline import ScaledForecastingPipeline
 
-__all__ = ["VARIMA", "AutoVARIMA"]
+__all__ = ["VARIMA", "AutoVARIMA", "ScaledForecastingPipeline"]

tslearn/forecasting/_pipeline.py

@@ -0,0 +1,187 @@
+"""Scale-predict-unscale pipeline for forecasters."""
+
+from sklearn.base import BaseEstimator, clone
+from sklearn.utils.validation import check_is_fitted
+
+from tslearn.backend import instantiate_backend
+from tslearn.bases import TimeSeriesMixin
+from tslearn.preprocessing import TimeSeriesScalerMeanVariance
+from tslearn.utils import check_array, to_time_series_dataset
+
+
+class ScaledForecastingPipeline(TimeSeriesMixin, BaseEstimator):
+    """Scale a series, forecast, and un-scale the forecast.
+
+    A pattern often seen in forecasting is to scale each series before
+    handing it over to a model, and to un-scale the model's predictions
+    using the statistics computed by the scaler, so that the forecaster
+    itself only ever sees normalized data. This estimator
+    acts as a pipeline and fills that gap for :mod:`tslearn.forecasting` forecasters.
+
+    Parameters
+    ----------
+        forecaster : estimator
+          A tslearn forecaster, exposing ``fit(X, y=None)`` and
+          ``predict(X=None, n=1)``, that will be trained on the scaled data.
+        scaler : transformer or None (default: None)
+          A scaler exposing ``fit``, ``transform`` and ``inverse_transform``,
+          such as the ones in :mod:`tslearn.preprocessing`. When None, a
+          :class:`~tslearn.preprocessing.TimeSeriesScalerMeanVariance` is
+          used, scaling each series independently (its default,
+          ``per_timeseries=True``).
+
+    Attributes
+    ----------
+        per_timeseries_ : bool
+          Whether scaling is performed independently per series (mirrors
+          ``scaler.per_timeseries``, defaulting to False when the scaler
+          does not have such a parameter).
+        scalers_ : list of transformers
+          The fitted scaler(s): a single one when ``per_timeseries_`` is False,
+          one per training series otherwise.
+        forecaster_ : estimator
+          The forecaster, fitted on scaled data.
+
+    Examples
+    --------
+    >>> from tslearn.forecasting import VARIMA
+    >>> from tslearn.generators import random_walks
+    >>> X = random_walks(n_ts=2, sz=20, d=1, random_state=0)
+    >>> model = ScaledForecastingPipeline(VARIMA(1, 0, 0)).fit(X)
+    >>> model.predict(n=3).shape
+    (2, 3, 1)
+
+    See Also
+    --------
+        VARIMA, AutoVARIMA: Forecasters that can be wrapped as-is.
+    """
+
+    def __init__(self, forecaster, scaler=None):
+        self.forecaster = forecaster
+        self.scaler = scaler
+
+    def _resolve_scaler(self):
+        """Build the per-group scaler template and detect the scaling mode.
+
+        A ``per_timeseries=True`` scaler cannot serve `inverse_transform`,
+        so per-series scaling is emulated with one ``per_timeseries=False``
+        clone of the scaler per series, rather than with a single
+        ``per_timeseries=True`` instance for the whole batch.
+        """
+        scaler = (
+            TimeSeriesScalerMeanVariance()
+            if self.scaler is None
+            else clone(self.scaler)
+        )
+        if not hasattr(scaler, "inverse_transform"):
+            raise ValueError(
+                f"{scaler.__class__.__name__} does not implement "
+                "`inverse_transform`, so it cannot be used to un-scale "
+                "forecasts."
+            )
+        per_series = bool(getattr(scaler, "per_timeseries", False))
+        if per_series:
+            scaler.set_params(per_timeseries=False)
+        return scaler, per_series
+
+    def _fit_scalers(self, X):
+        if not self.per_timeseries_:
+            return [clone(self._scaler_template_).fit(X)]
+        return [
+            clone(self._scaler_template_).fit(X[i : i + 1])
+            for i in range(X.shape[0])
+        ]
+
+    def _transform(self, scalers, X):
+        if not self.per_timeseries_:
+            return scalers[0].transform(X)
+        be = instantiate_backend(X)
+        return be.vstack(
+            [scaler.transform(X[i : i + 1]) for i, scaler in enumerate(scalers)]
+        )
+
+    def _inverse_transform(self, scalers, forecast):
+        if not self.per_timeseries_:
+            return scalers[0].inverse_transform(forecast)
+        be = instantiate_backend(forecast)
+        return be.vstack(
+            [
+                scaler.inverse_transform(forecast[i : i + 1])
+                for i, scaler in enumerate(scalers)
+            ]
+        )
+
+    def fit(self, X, y=None):
+        """Fit the scaler(s), then the forecaster on the scaled data.
+
+        Parameters
+        ----------
+            X : array-like of shape=(n_ts, sz, d)
+                Time series dataset.
+            y : Ignored
+
+        Returns
+        -------
+            self
+                The fitted estimator
+        """
+        X = check_array(X, allow_nd=True, force_all_finite=False)
+        X = to_time_series_dataset(X)
+
+        self._scaler_template_, self.per_timeseries_ = self._resolve_scaler()
+        self.scalers_ = self._fit_scalers(X)
+        X_scaled = self._transform(self.scalers_, X)
+
+        self.forecaster_ = clone(self.forecaster)
+        self.forecaster_.fit(X_scaled)
+        return self
+
+    def predict(self, X=None, n=1):
+        """Forecast ``n`` timestamps ahead, on the original data scale.
+
+        Parameters
+        ----------
+            X : array-like of shape=(n_ts, sz, d) or None (default: None)
+              Time series dataset to forecast. If None, the data passed at
+              fit time is forecasted instead.
+            n : int (default: 1)
+              The number of timestamps to forecast, a.k.a. the horizon.
+
+        Returns
+        -------
+            array of shape=(n_ts, n, d)
+              Array of forecasted timestamps, on the original data scale.
+        """
+        check_is_fitted(self, "forecaster_")
+        if X is None:
+            scalers = self.scalers_
+            X_scaled = None
+        else:
+            X = check_array(X, allow_nd=True, force_all_finite=False)
+            X = to_time_series_dataset(X)
+            scalers = (
+                self.scalers_
+                if not self.per_timeseries_
+                else self._fit_scalers(X)
+            )
+            X_scaled = self._transform(scalers, X)
+        forecast_scaled = self.forecaster_.predict(X_scaled, n=n)
+        return self._inverse_transform(scalers, forecast_scaled)
+
+    def fit_predict(self, X, y=None, n=1):
+        """Fit the estimator and forecast ``n`` timestamps for the given data.
+
+        Parameters
+        ----------
+            X : array-like of shape=(n_ts, sz, d)
+                Time series dataset.
+            y : Ignored
+            n : int (default: 1)
+                The number of timestamps to forecast, a.k.a. the horizon.
+
+        Returns
+        -------
+            array of shape=(n_ts, n, d)
+              Array of forecasted timestamps, on the original data scale.
+        """
+        return self.fit(X, y).predict(X, n=n)