diff --git a/CHANGELOG.md b/CHANGELOG.md index 3736cca7..29a2ecac 100644 --- a/CHANGELOG.md +++ b/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 diff --git a/docs/gen_modules/tslearn.forecasting.rst b/docs/gen_modules/tslearn.forecasting.rst index 0236e1ac..e62231c5 100644 --- a/docs/gen_modules/tslearn.forecasting.rst +++ b/docs/gen_modules/tslearn.forecasting.rst @@ -13,3 +13,4 @@ tslearn.forecasting VARIMA AutoVARIMA + ScaledForecastingPipeline diff --git a/tests/test_forecasting.py b/tests/test_forecasting.py index 26914136..940d9d30 100644 --- a/tests/test_forecasting.py +++ b/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 \ No newline at end of file + 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) diff --git a/tests/test_preprocessing.py b/tests/test_preprocessing.py index 4c9b2c62..eea8fa4a 100644 --- a/tests/test_preprocessing.py +++ b/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], diff --git a/tslearn/forecasting/__init__.py b/tslearn/forecasting/__init__.py index 54ff4962..0bf5ff01 100644 --- a/tslearn/forecasting/__init__.py +++ b/tslearn/forecasting/__init__.py @@ -5,5 +5,6 @@ """ from ._arima import VARIMA, AutoVARIMA +from ._pipeline import ScaledForecastingPipeline -__all__ = ["VARIMA", "AutoVARIMA"] +__all__ = ["VARIMA", "AutoVARIMA", "ScaledForecastingPipeline"] diff --git a/tslearn/forecasting/_pipeline.py b/tslearn/forecasting/_pipeline.py new file mode 100644 index 00000000..808d5c60 --- /dev/null +++ b/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) diff --git a/tslearn/preprocessing/preprocessing.py b/tslearn/preprocessing/preprocessing.py index 647da0bb..648abfd7 100644 --- a/tslearn/preprocessing/preprocessing.py +++ b/tslearn/preprocessing/preprocessing.py @@ -140,20 +140,31 @@ def __sklearn_tags__(self): class TimeSeriesScalerMinMax(TimeSeriesMixin, TransformerMixin, BaseEstimator): """Scaler for time series datasets. - When `per_timeseries` is False, scales features based on computation led on the fitted data, - so that their span in given dimensions is between ``min`` and ``max`` where ``value_range=(min, max)``. - The transformation is stateless otherwise, dealing with each timeseries individually. Parameters ---------- value_range : tuple (default: (0., 1.)) The minimum and maximum value for the output time series. - per_timeseries: bool (default: True) + per_timeseries : bool (default: True) Wether the scaling should be performed per time series. - per_feature: bool (default: True) + When `per_timeseries` is False, scales features based on computation led on the fitted data, + so that their span in given dimensions is between ``min`` and ``max`` where ``value_range=(min, max)``. + The transformation is stateless otherwise, dealing with each timeseries individually. + per_feature : bool (default: True) Wether the scaling should be performed per feature. Meaningless for univariate timeseries. + Attributes + ---------- + min_ : array-like of shape=(1, 1, d or 1) + The miminum value(s) seen in data and used for normalization. + Only available when `per_timeseries` is False. + The shape depends on the `per_feature` parameters. + max_ : array-like of shape=(1, 1, d or 1) + The maximum value(s) seen in data and used for normalization. + Only available when `per_timeseries` is False. + The shape depends on the `per_feature` parameters. + Notes ----- NaNs within a time series are ignored when calculating min and max. @@ -203,9 +214,15 @@ def fit(self, X, y=None, **kwargs): self._X_fit_dims = X_.shape self.n_features_in_ = self._X_fit_dims[-1] - if not self.per_timeseries: - self.min_, self.max_ = self._process(X_) + # Reset if needed + if hasattr(self, 'min_'): + del self.min_ + if hasattr(self, 'max_'): + del self.max_ + if not self.per_timeseries: + min_, max_ = self._process(X_) + self.min_, self.max_ = min_.reshape((-1,)), max_.reshape((-1,)) return self def fit_transform(self, X, y=None, **kwargs): @@ -248,7 +265,7 @@ def transform(self, X, y=None, **kwargs): Returns ------- - numpy.ndarray + array-like of shape (n_ts, sz, d) Rescaled time series dataset. """ if self.value_range[0] >= self.value_range[1]: @@ -260,7 +277,8 @@ def transform(self, X, y=None, **kwargs): X_ = to_time_series_dataset(X_) X_ = check_dims(X_, X_fit_dims=self._X_fit_dims, check_n_features_only=True, extend=False) - min_, max_ = self._process(X_) if self.per_timeseries else (self.min_, self.max_) + min_, max_ = (self._process(X_) + if self.per_timeseries else (self.min_.reshape(1, 1, -1), self.max_.reshape(1, 1, -1))) range_t = max_ - min_ range_t[range_t == 0.] = 1. @@ -268,6 +286,36 @@ def transform(self, X, y=None, **kwargs): X_ = nomin / range_t + self.value_range[0] return X_ + def inverse_transform(self, X): + """ + Undo the scaling of X based on fitted data. + Only available when per_timeseries is False, raises RuntimeError otherwise. + + Parameters + ---------- + X : array-like of shape (n_ts, sz, d) + Input dataset. + + Returns + ------- + array-like of shape (n_ts, sz, d) + Transformed dataset. + """ + check_is_fitted(self, '_X_fit_dims') + if self.per_timeseries: + raise RuntimeError("Cannot inverse per timeseries scaling.") + + X_ = check_array(X, allow_nd=True, force_all_finite=False) + X_ = to_time_series_dataset(X_) + X_ = check_dims(X_, X_fit_dims=self._X_fit_dims, check_n_features_only=True, extend=False) + + min_, max_ = self.min_.reshape(1, 1, -1), self.max_.reshape(1, 1, -1) + + X_ -= self.value_range[0] + X_ *= (max_ - min_) / (self.value_range[1] - self.value_range[0]) + X_ += min_ + return X_ + def _more_tags(self): tags = super()._more_tags() tags.update({'allow_nan': True, ALLOW_VARIABLE_LENGTH: True}) @@ -282,9 +330,6 @@ def __sklearn_tags__(self): class TimeSeriesScalerMeanVariance(TimeSeriesMixin, TransformerMixin, BaseEstimator): """Scaler for time series datasets. - When `per_timeseries` is False, scales features based on computation led on the fitted data, - so that their mean (resp. standard deviation) in given dimensions is mu (resp. std). - The transformation is stateless otherwise, dealing with each timeseries individually. Parameters ---------- @@ -292,12 +337,27 @@ class TimeSeriesScalerMeanVariance(TimeSeriesMixin, TransformerMixin, BaseEstima Mean of the output time series. std : float (default: 1.) Standard deviation of the output time series. - per_timeseries: bool (default: True) + per_timeseries : bool (default: True) Whether the scaling should be performed per time series. - per_feature: bool (default: True) + When `per_timeseries` is False, scales features based on computation led on the fitted data, + so that their mean (resp. standard deviation) in given dimensions is mu (resp. std). + The transformation is stateless otherwise, dealing with each timeseries individually. + per_feature : bool (default: True) Whether the scaling should be performed per feature. Meaningless for univariate timeseries. + Attributes + ---------- + mean_ : array-like of shape=(1, 1, d or 1) + The mean value(s) seen in input data and used for normalization. + Only available when `per_timeseries` is False. + The shape depends on the `per_feature` parameters. + + std_ : array-like of shape=(1, 1, d or 1) + The standard deviation values seen in input data and used for normalization. + Only available when `per_timeseries` is False. + The shape depends on the `per_feature` parameters. + Notes ----- NaNs within a time series are ignored when calculating mu and std. @@ -347,6 +407,12 @@ def fit(self, X, y=None, **kwargs): self._X_fit_dims = X_.shape self.n_features_in_ = self._X_fit_dims[-1] + # Reset if needed + if hasattr(self, 'mean_'): + del self.mean_ + if hasattr(self, 'std_'): + del self.std_ + if not self.per_timeseries: self.mean_, self.std_ = self._process(X_) @@ -402,11 +468,43 @@ def transform(self, X, y=None, **kwargs): X_ = to_time_series_dataset(X_) X_ = check_dims(X_, X_fit_dims=self._X_fit_dims, check_n_features_only=True, extend=False) - mean_, std_ = self._process(X_) if self.per_timeseries else (self.mean_, self.std_) + mean_, std_ = self._process(X_) if self.per_timeseries else (self.mean_.reshape(1, 1, -1), + self.std_.reshape(1, 1, -1)) X_ = (X_ - mean_) * self.std / std_ + self.mu return X_ + def inverse_transform(self, X): + """ + Undo the scaling of X based on fitted data. + Only available when per_timeseries is False, raises RuntimeError otherwise. + + Parameters + ---------- + X : array-like of shape (n_ts, sz, d) + Input dataset. + + Returns + ------- + array-like of shape (n_ts, sz, d) + Transformed dataset. + """ + + check_is_fitted(self, '_X_fit_dims') + if self.per_timeseries: + raise RuntimeError("Cannot inverse per timeseries scaling.") + + X_ = check_array(X, allow_nd=True, force_all_finite=False) + X_ = to_time_series_dataset(X_) + X_ = check_dims(X_, X_fit_dims=self._X_fit_dims, check_n_features_only=True, extend=False) + + mean_, std_ = self.mean_.reshape(1, 1, -1), self.std_.reshape(1, 1, -1) + + X_ += self.mu + X_ *= std_ / self.std + X_ += mean_ + return X_ + def _more_tags(self): tags = super()._more_tags() tags.update({'allow_nan': True, ALLOW_VARIABLE_LENGTH: True})