From aa76e2bc109a9bbcae0ed4da20541a9952a31291 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Thu, 24 Oct 2024 09:18:02 +0100 Subject: [PATCH 01/31] forecaster base and dummy --- aeon/forecasting/__init__.py | 5 ++ aeon/forecasting/_dummy.py | 25 ++++++ aeon/forecasting/base.py | 121 ++++++++++++++++++++++++++++ aeon/forecasting/tests/__init__.py | 1 + aeon/forecasting/tests/test_base.py | 16 ++++ aeon/utils/base/_register.py | 2 + 6 files changed, 170 insertions(+) create mode 100644 aeon/forecasting/__init__.py create mode 100644 aeon/forecasting/_dummy.py create mode 100644 aeon/forecasting/base.py create mode 100644 aeon/forecasting/tests/__init__.py create mode 100644 aeon/forecasting/tests/test_base.py diff --git a/aeon/forecasting/__init__.py b/aeon/forecasting/__init__.py new file mode 100644 index 0000000000..be6dc99141 --- /dev/null +++ b/aeon/forecasting/__init__.py @@ -0,0 +1,5 @@ +"""Forecasters.""" + +__all__ = ["DummyForecaster"] + +from aeon.forecasting._dummy import DummyForecaster diff --git a/aeon/forecasting/_dummy.py b/aeon/forecasting/_dummy.py new file mode 100644 index 0000000000..f708d2eb40 --- /dev/null +++ b/aeon/forecasting/_dummy.py @@ -0,0 +1,25 @@ +"""DummyForecaster always predicts the last value seen in training.""" + +from aeon.forecasting.base import BaseForecaster + + +class DummyForecaster(BaseForecaster): + """Dummy forecaster always predicts the last value seen in training.""" + + def __init__(self): + """Initialize DummyForecaster.""" + self.last_value_ = None + super().__init__() + + def _fit(self, y, exog=None): + """Fit dummy forecaster.""" + self.last_value_ = y[-1] + return self + + def _predict(self, y=None, exog=None): + """Predict using dummy forecaster.""" + return self.last_value_ + + def _forecast(self, y, X=None): + """Forecast using dummy forecaster.""" + return y[-1] diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py new file mode 100644 index 0000000000..04f05ac184 --- /dev/null +++ b/aeon/forecasting/base.py @@ -0,0 +1,121 @@ +"""BaseForecaster class. + +A simplified first base class for foreacasting models. The focus here is on a +specific form of forecasting: longer series, long winodws and single step forecasting. + +aeon enhancement proposal +https://github.com/aeon-toolkit/aeon-admin/pull/14 + +""" + +from abc import ABC, abstractmethod + +from aeon.base import BaseSeriesEstimator + + +class BaseForecaster(BaseSeriesEstimator, ABC): + """ + Abstract base class for time series forecasters. + + The base forecaster specifies the methods and method signatures that all + forecasters have to implement. Attributes with an underscore suffix are set in the + method fit. + + Parameters + ---------- + horizon : int, default =1 + The number of time steps ahead to forecast. If horizon is one, the forecaster + will learn to predict one point ahead + window : int or None + The window prior to the current time point to use in forecasting. So if + horizon is one, forecaster will train using points $i$ to $window+i-1$ to + predict value $window+i$. If horizon is 4, forecaster will used points $i$ + to $window+i-1$ to predict value $window+i+3$. If None, the algorithm will + internally determine what data to use to predict `horizon` steps ahead. + """ + + # TODO: add any forecasting specific tags + _tags = { + "capability:univariate": True, + "capability:multivariate": False, + "capability:missing_values": False, + "X_inner_type": "np.ndarray", + } + + def __init__(self, horizon=1, window=None, axis=1): + self.horizon = horizon + self.window = window + self._is_fitted = False + super().__init__(axis) + + def fit(self, y, exog=None): + """Fit forecaster to series y. + + Fit a forecaster to predict self.horizon steps ahead using y. + + Parameters + ---------- + y : np.ndarray + A time series on which to learn a forecaster to predict horizon ahead + exog : np.ndarray, default =None + Optional exogenous time series data assumed to be aligned with y + + Returns + ------- + self + Fitted BaseForecaster. + """ + # Validate y + + # Convert if necessary + + if exog is not None: + raise NotImplementedError("Exogenous variables not yet supported") + # Validate exog + self._is_fitted = True + return self._fit(y, exog) + + @abstractmethod + def _fit(self, y, exog=None): ... + + def predict(self, y=None, exog=None): + """Predict the next horizon steps ahead. + + Parameters + ---------- + y : np.ndarray, default = None + A time series to predict the next horizon value for. If None, + predict the next horizon value after series seen in fit. + exog : np.ndarray, default =None + Optional exogenous time series data assumed to be aligned with y + + Returns + ------- + float + single prediction self.horizon steps ahead of y. + """ + if not self._is_fitted: + raise ValueError("Forecaster must be fitted before predicting") + if exog is not None: + raise NotImplementedError("Exogenous variables not yet supported") + # Validate exog + self._is_fitted = True + return self._predict(y, exog) + + @abstractmethod + def _predict(self, y=None, exog=None): ... + + def forecast(self, y, X=None): + """ + + Forecast basically fit_predict. + + Returns + ------- + np.ndarray + single prediction directly after the last point in X. + """ + return self._forecast(y, X) + + @abstractmethod + def _forecast(self, y=None, exog=None): ... diff --git a/aeon/forecasting/tests/__init__.py b/aeon/forecasting/tests/__init__.py new file mode 100644 index 0000000000..90b32266a4 --- /dev/null +++ b/aeon/forecasting/tests/__init__.py @@ -0,0 +1 @@ +"""Forecaster tests.""" diff --git a/aeon/forecasting/tests/test_base.py b/aeon/forecasting/tests/test_base.py new file mode 100644 index 0000000000..e1a634eba3 --- /dev/null +++ b/aeon/forecasting/tests/test_base.py @@ -0,0 +1,16 @@ +"""Test base forecaster.""" + +import numpy as np + +from aeon.forecasting import DummyForecaster + + +def test_base_forecaster(): + """Test base forecaster functionality.""" + f = DummyForecaster() + y = np.random.rand(50) + f.fit(y) + p1 = f.predict() + assert p1 == y[-1] + p2 = f.forecast(y) + assert p2 == p1 diff --git a/aeon/utils/base/_register.py b/aeon/utils/base/_register.py index 1327d626ef..1d81c2512c 100644 --- a/aeon/utils/base/_register.py +++ b/aeon/utils/base/_register.py @@ -21,6 +21,7 @@ from aeon.classification.base import BaseClassifier from aeon.classification.early_classification import BaseEarlyClassifier from aeon.clustering.base import BaseClusterer +from aeon.forecasting.base import BaseForecaster from aeon.regression.base import BaseRegressor from aeon.segmentation.base import BaseSegmenter from aeon.similarity_search.base import BaseSimilaritySearch @@ -45,6 +46,7 @@ "segmenter": BaseSegmenter, "similarity_searcher": BaseSimilaritySearch, "series-transformer": BaseSeriesTransformer, + "forecaster": BaseForecaster, } # base classes which are valid for estimator to directly inherit from From a5dbc282c6951713285a28dbe9ce17eb93b379c5 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Fri, 25 Oct 2024 10:20:06 +0100 Subject: [PATCH 02/31] forecasting tests --- aeon/forecasting/base.py | 4 +- .../_yield_forecasting_checks.py | 52 +++++++++++++++++++ .../mock_estimators/_mock_forecasters.py | 22 ++++++++ aeon/testing/testing_data.py | 2 + 4 files changed, 78 insertions(+), 2 deletions(-) create mode 100644 aeon/testing/estimator_checking/_yield_forecasting_checks.py create mode 100644 aeon/testing/mock_estimators/_mock_forecasters.py diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index 04f05ac184..84fd33d9fa 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -8,12 +8,12 @@ """ -from abc import ABC, abstractmethod +from abc import abstractmethod from aeon.base import BaseSeriesEstimator -class BaseForecaster(BaseSeriesEstimator, ABC): +class BaseForecaster(BaseSeriesEstimator): """ Abstract base class for time series forecasters. diff --git a/aeon/testing/estimator_checking/_yield_forecasting_checks.py b/aeon/testing/estimator_checking/_yield_forecasting_checks.py new file mode 100644 index 0000000000..50e19d11aa --- /dev/null +++ b/aeon/testing/estimator_checking/_yield_forecasting_checks.py @@ -0,0 +1,52 @@ +"""Tests for all forecasters.""" + +from functools import partial + +import numpy as np + +from aeon.base._base import _clone_estimator +from aeon.base._base_series import VALID_INNER_TYPES + + +def _yield_forecasting_checks(estimator_class, estimator_instances, datatypes): + """Yield all forecasting checks for an aeon forecaster.""" + # only class required + yield partial(check_forecasting_base_functionality, estimator_class=estimator_class) + + # test class instances + for _, estimator in enumerate(estimator_instances): + # no data needed + yield partial(check_forecaster_instance, estimator=estimator) + + +def check_forecasting_base_functionality(estimator_class): + """Test compliance with the base class contract.""" + # Test they dont override final methods, because python does not enforce this + assert "fit" not in estimator_class.__dict__ + assert "predict" not in estimator_class.__dict__ + assert "forecast" not in estimator_class.__dict__ + fit_is_empty = estimator_class.get_class_tag(tag_name="fit_is_empty") + assert not fit_is_empty == "_fit" not in estimator_class.__dict__ + # Test valid tag for X_inner_type + X_inner_type = estimator_class.get_class_tag(tag_name="X_inner_type") + assert X_inner_type in VALID_INNER_TYPES + # Must have at least one set to True + multi = estimator_class.get_class_tag(tag_name="capability:multivariate") + uni = estimator_class.get_class_tag(tag_name="capability:univariate") + assert multi or uni + + +def check_forecaster_instance(estimator): + """Test forecasters.""" + estimator = _clone_estimator(estimator) + pass + # Sort + # Check output correct: predict should return a float + y = np.array([0.5, 0.7, 0.8, 0.9, 1.0]) + estimator.fit(y) + p = estimator.predict() + assert isinstance(p, float) + # forecast should return a float equal to fit/predict + p2 = estimator.forecast(y) + assert p == p2 + # Add other tests when diff --git a/aeon/testing/mock_estimators/_mock_forecasters.py b/aeon/testing/mock_estimators/_mock_forecasters.py new file mode 100644 index 0000000000..e5b5403b08 --- /dev/null +++ b/aeon/testing/mock_estimators/_mock_forecasters.py @@ -0,0 +1,22 @@ +"""Mock forecasters useful for testing and debugging. + +Used in tests for the forecasting base class. +""" + +from aeon.forecasting.base import BaseForecaster + + +class MockForecaster(BaseForecaster): + """Mock segmenter for testing.""" + + def __init__(self): + super().__init__() + + def _fit(self, y, X=None): + return self + + def _predict(self, y): + return 1.0 + + def _forecast(self, y, X=None): + return 1.0 diff --git a/aeon/testing/testing_data.py b/aeon/testing/testing_data.py index 7ab34d2be3..a5dcf7acf6 100644 --- a/aeon/testing/testing_data.py +++ b/aeon/testing/testing_data.py @@ -7,6 +7,7 @@ from aeon.classification import BaseClassifier from aeon.classification.early_classification import BaseEarlyClassifier from aeon.clustering import BaseClusterer +from aeon.forecasting import BaseForecaster from aeon.regression import BaseRegressor from aeon.segmentation import BaseSegmenter from aeon.similarity_search import BaseSimilaritySearch @@ -829,6 +830,7 @@ def _get_label_type_for_estimator(estimator): isinstance(estimator, BaseAnomalyDetector) or isinstance(estimator, BaseSegmenter) or isinstance(estimator, BaseSeriesTransformer) + or isinstance(estimator, BaseForecaster) ): label_type = "NoLabel" else: From 7fee274730f2aa9ced1acbb9e8c0ce5b174aecec Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Fri, 25 Oct 2024 10:47:33 +0100 Subject: [PATCH 03/31] forecasting tests --- aeon/testing/testing_data.py | 10 +++++++++- 1 file changed, 9 insertions(+), 1 deletion(-) diff --git a/aeon/testing/testing_data.py b/aeon/testing/testing_data.py index a5dcf7acf6..c5b1dfe8ff 100644 --- a/aeon/testing/testing_data.py +++ b/aeon/testing/testing_data.py @@ -814,7 +814,15 @@ def _get_label_type_for_estimator(estimator): Returns ------- label_type : str - Label type key for the estimator for use in TEST_LABEL_DICT. + Label type key for the estimator for use in FULL_TEST_DATA_DICT. Indicates + whether estimator can take labels for training data, and if so what kind. + Classification indicates the estimator takes a discrete target variable, + modelled by an np.array of integers. + Regression indicates the estimator takes a continuous target variable, + modelled by an np.array of floats. + NoLabel indicates the estimator does not take a target variable, or that no + estimator is yet implemented to take labels. + """ if ( isinstance(estimator, BaseClassifier) From c4628ca757f801de015c1b3b713c810e4979c28b Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Fri, 25 Oct 2024 10:49:20 +0100 Subject: [PATCH 04/31] forecasting tests --- aeon/testing/testing_data.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/aeon/testing/testing_data.py b/aeon/testing/testing_data.py index c5b1dfe8ff..13105025f1 100644 --- a/aeon/testing/testing_data.py +++ b/aeon/testing/testing_data.py @@ -722,8 +722,8 @@ def _get_datatypes_for_estimator(estimator): Returns ------- datatypes : list of tuple - List of valid data types keys for the estimator usable in FULL_TEST_DATA_DICT - and TEST_LABEL_DICT. Each tuple is formatted (data_key, label_key). + List of valid data types keys for the estimator usable in + FULL_TEST_DATA_DICT. Each tuple is formatted (data_key, label_key). """ datatypes = [] univariate, multivariate, unequal_length, missing_values = ( From e1540788975220c9e10582c10bc683165dbfe004 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Fri, 25 Oct 2024 11:27:57 +0100 Subject: [PATCH 05/31] forecasting tests --- aeon/forecasting/__init__.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/aeon/forecasting/__init__.py b/aeon/forecasting/__init__.py index be6dc99141..3f2263f987 100644 --- a/aeon/forecasting/__init__.py +++ b/aeon/forecasting/__init__.py @@ -1,5 +1,6 @@ """Forecasters.""" -__all__ = ["DummyForecaster"] +__all__ = ["DummyForecaster", "BaseForecaster"] from aeon.forecasting._dummy import DummyForecaster +from aeon.forecasting.base import BaseForecaster From a6043d63108827f257a4581118caeaf7117f90ca Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Fri, 25 Oct 2024 17:24:59 +0100 Subject: [PATCH 06/31] regression --- aeon/forecasting/_regression.py | 61 +++++++++++++++++++++++++++++++++ 1 file changed, 61 insertions(+) create mode 100644 aeon/forecasting/_regression.py diff --git a/aeon/forecasting/_regression.py b/aeon/forecasting/_regression.py new file mode 100644 index 0000000000..3cce7b0704 --- /dev/null +++ b/aeon/forecasting/_regression.py @@ -0,0 +1,61 @@ +"""Window-based regression forecaster. + +General purpose forecaster to use with any scikit learn or aeon compatible +regressor. Simply forms a collection of windows from the time series and trains to +predict the next +""" + +import numpy as np +from sklearn.linear_model import LinearRegression + +from aeon.forecasting.base import BaseForecaster + + +class RegressionForecaster(BaseForecaster): + def __init__(self, window, horizon=1, regressor=None): + self.regressor = regressor + super().__init__(horizon, window) + + def _fit(self, y, exog=None): + """Fit forecaster to time series. + + Split X into windows of length window and train the forecaster on each window + to predict the horizon ahead. + + Parameters + ---------- + X : Time series on which to learn a forecaster + + Returns + ------- + self + Fitted estimator + """ + # Window data + if self.regressor is None: + self.regressor_ = LinearRegression() + else: + self.regressor_ = self.regressor + X = np.lib.stride_tricks.sliding_window_view(y, window_shape=self.window) + # Ignore the final horizon values: need to store these for pred with empty y + X = X[: -self.horizon] + # Extract y + y = y[self.window + self.horizon - 1 :] + self.last_ = y[-self.window :] + self.regressor_.fit(y, exog) + return self + + def _predict(self, y=None, exog=None): + """Predict values for time series X.""" + if y is None: + return self.regressor_.predict(self.last_) + + return self.regressor_.predict(y[-self.window :]) + + def _forecast(self, y, exog=None): + """Forecast values for time series X. + + NOTE: deal with horizons + """ + self.fit(y, exog) + return self.predict(y) From 5db044b129f479a60ffc349f9905ede5d867b674 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Fri, 25 Oct 2024 20:25:36 +0100 Subject: [PATCH 07/31] notebook --- aeon/forecasting/__init__.py | 3 +- aeon/forecasting/_dummy.py | 4 +- aeon/forecasting/_regression.py | 5 +- aeon/forecasting/base.py | 7 +- aeon/utils/tags/_tags.py | 8 + examples/forecasting/forecasting.ipynb | 312 +++++++++++++++++++++++++ 6 files changed, 333 insertions(+), 6 deletions(-) create mode 100644 examples/forecasting/forecasting.ipynb diff --git a/aeon/forecasting/__init__.py b/aeon/forecasting/__init__.py index 3f2263f987..e114736cbc 100644 --- a/aeon/forecasting/__init__.py +++ b/aeon/forecasting/__init__.py @@ -1,6 +1,7 @@ """Forecasters.""" -__all__ = ["DummyForecaster", "BaseForecaster"] +__all__ = ["DummyForecaster", "BaseForecaster", "RegressionForecaster"] from aeon.forecasting._dummy import DummyForecaster +from aeon.forecasting._regression import RegressionForecaster from aeon.forecasting.base import BaseForecaster diff --git a/aeon/forecasting/_dummy.py b/aeon/forecasting/_dummy.py index f708d2eb40..ab2ac4b8ee 100644 --- a/aeon/forecasting/_dummy.py +++ b/aeon/forecasting/_dummy.py @@ -13,6 +13,7 @@ def __init__(self): def _fit(self, y, exog=None): """Fit dummy forecaster.""" + y = y.squeeze() self.last_value_ = y[-1] return self @@ -20,6 +21,7 @@ def _predict(self, y=None, exog=None): """Predict using dummy forecaster.""" return self.last_value_ - def _forecast(self, y, X=None): + def _forecast(self, y, exog=None): """Forecast using dummy forecaster.""" + y = y.squeeze() return y[-1] diff --git a/aeon/forecasting/_regression.py b/aeon/forecasting/_regression.py index 3cce7b0704..bc01213d17 100644 --- a/aeon/forecasting/_regression.py +++ b/aeon/forecasting/_regression.py @@ -36,20 +36,21 @@ def _fit(self, y, exog=None): self.regressor_ = LinearRegression() else: self.regressor_ = self.regressor + y = y.squeeze() X = np.lib.stride_tricks.sliding_window_view(y, window_shape=self.window) # Ignore the final horizon values: need to store these for pred with empty y X = X[: -self.horizon] # Extract y y = y[self.window + self.horizon - 1 :] self.last_ = y[-self.window :] - self.regressor_.fit(y, exog) + self.regressor_.fit(X=X, y=y) return self def _predict(self, y=None, exog=None): """Predict values for time series X.""" if y is None: return self.regressor_.predict(self.last_) - + y = y.squeeze() return self.regressor_.predict(y[-self.window :]) def _forecast(self, y, exog=None): diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index 84fd33d9fa..9c213a2c41 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -39,7 +39,7 @@ class BaseForecaster(BaseSeriesEstimator): "capability:univariate": True, "capability:multivariate": False, "capability:missing_values": False, - "X_inner_type": "np.ndarray", + "y_inner_type": "np.ndarray", } def __init__(self, horizon=1, window=None, axis=1): @@ -68,7 +68,7 @@ def fit(self, y, exog=None): # Validate y # Convert if necessary - + y = self._preprocess_series(y, axis=self.axis, store_metadata=False) if exog is not None: raise NotImplementedError("Exogenous variables not yet supported") # Validate exog @@ -94,6 +94,8 @@ def predict(self, y=None, exog=None): float single prediction self.horizon steps ahead of y. """ + if y is not None: + y = self._preprocess_series(y, axis=self.axis, store_metadata=False) if not self._is_fitted: raise ValueError("Forecaster must be fitted before predicting") if exog is not None: @@ -115,6 +117,7 @@ def forecast(self, y, X=None): np.ndarray single prediction directly after the last point in X. """ + y = self._preprocess_series(y, axis=self.axis, store_metadata=False) return self._forecast(y, X) @abstractmethod diff --git a/aeon/utils/tags/_tags.py b/aeon/utils/tags/_tags.py index 650a7c4dc4..8be4d787ee 100644 --- a/aeon/utils/tags/_tags.py +++ b/aeon/utils/tags/_tags.py @@ -53,6 +53,14 @@ class : identifier for the base class of objects this tag applies to "description": "What data structure(s) the estimator uses internally for " "fit/predict.", }, + "y_inner_type": { + "class": "estimator", + "type": [ + ("list||str", SERIES_DATA_TYPES), + ], + "description": "What data structure(s) the estimator uses internally for " + "fit/predict.", + }, "algorithm_type": { "class": "estimator", "type": [ diff --git a/examples/forecasting/forecasting.ipynb b/examples/forecasting/forecasting.ipynb new file mode 100644 index 0000000000..8db44440a5 --- /dev/null +++ b/examples/forecasting/forecasting.ipynb @@ -0,0 +1,312 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "source": [ + "# Time series forecasting with aeon\n", + "\n", + "This notebook describes the new, experimental, forecasting module in aeon. We have\n", + "recently removed a lot of legacy code that was almost entirely wrappers around other\n", + "projects, mostly statsmodels. Most of\n", + "the contributors to aeon are from a computer science/machine learning background\n", + "rather than stats and forecasting, and our objectives for forecasting have changed to\n", + " reflect this. Our focus is on:\n", + "\n", + "1. not attempting to be a comprehensive forecasting package.\n", + "\n", + "\n", + "We are not trying to do all things in forecasting. There are many packages that do\n", + "that already. Our previous experience taught us that trying to be all things to all\n", + "people resulted in over engineered spaghetti that was hard to maintain, understand or\n", + " use.\n", + "\n", + "2. fast forecasting with numpy arrays.\n", + "\n", + "Whilst our forecasters will work with data frames, our design principle is to write\n", + "code optimised with numba and numpy. We found that extensive use of data frames in\n", + "the internal calculations of forecasters makes them much slower and harder to\n", + "understand for those not used to using dataframes daily.\n", + "\n", + "3. forecasting using machine learning and deep learning.\n", + "\n", + "we want to implement and assess the latest machine learning and deep learning\n", + "forecasting for scenarios where it makes sense to use them. Our initial experimental\n", + "focus will be on forecasting with long series for a single forecasting horizon.\n" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "## Base Class\n", + "\n", + "Our first design choice for forecasting is to pass the forecasting horizon in the\n", + "constructor (default is 1). This is because we want a simpler use case: a forecaster\n", + "trains to predict so many places in the future, then for unseen data, it predicts the\n", + " same number of steps ahead. We recognise there are other scenarios, but this is the\n", + " cleanest way to start.\n", + "\n", + " The base class for all forecasters is `BaseForecaster`. It inherits from\n", + " `BaseSeriesEstimator`, which is also the base class for the other series estimators\n", + " in aeon: `BaseSegmenter`, `BaseAnomalyDetector` and `BaseSeriesTransformer`. The\n", + " base class `BaseSeriesEstimator` contains a method to validate and possibly convert\n", + " an input series.\n", + "The `BaseForecaster` has three core methods: `fit`, `predict` and `forecast`. It is\n", + "an abstract class, and each of these methods calls a protected method `_fit`,\n", + "`_predict` and `_forecast`.\n" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 1, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "['clone', 'fit', 'forecast', 'get_fitted_params', 'get_metadata_routing', 'get_params', 'get_tag', 'get_tags', 'predict', 'reset', 'set_params', 'set_tags']\n" + ] + } + ], + "source": [ + "import inspect\n", + "\n", + "from aeon.forecasting import BaseForecaster\n", + "\n", + "# List methods\n", + "public_methods = [\n", + " func[0]\n", + " for func in inspect.getmembers(BaseForecaster, predicate=inspect.isfunction)\n", + " if not func[0].startswith(\"_\")\n", + "]\n", + "print(public_methods)" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + " All estimators in `aeon` have tags. One specific to\n", + "forecasting is `y_inner_type`. This specifies the inner type the sub class of\n", + "BaseForecaster needs to input the method `_fit` and `_predict`. The default is `np\n", + ".ndarray` but it can also be `pd.DataFrame` or `pd.Series`. You can pass\n", + "forecaster and of `SERIES_DATA_TYPES` and it will be converted to `y_inner_type` in\n", + "`fit`, `predict` and `forecast`." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 2, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + " Possible data structures for input to forecaster ['pd.Series', 'pd.DataFrame', 'np.ndarray']\n", + "\n", + " Tags for BaseForecaster: {'python_version': None, 'python_dependencies': None, 'cant_pickle': False, 'non_deterministic': False, 'algorithm_type': None, 'capability:missing_values': False, 'capability:multithreading': False, 'capability:univariate': True, 'capability:multivariate': False, 'X_inner_type': 'np.ndarray', 'y_inner_type': 'np.ndarray'}\n" + ] + } + ], + "source": [ + "from aeon.utils import SERIES_DATA_TYPES\n", + "\n", + "print(\" Possible data structures for input to forecaster \", SERIES_DATA_TYPES)\n", + "print(\"\\n Tags for BaseForecaster: \", BaseForecaster.get_class_tags())" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "We use the standard airline dataset for examples. This can be stored as a pd.Series,\n", + "pd.DataFrame or np.ndarray." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 3, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "\n" + ] + } + ], + "source": [ + "import pandas as pd\n", + "\n", + "from aeon.datasets import load_airline\n", + "\n", + "y = load_airline()\n", + "print(type(y))\n", + "y2 = y.to_numpy()\n", + "y3 = pd.DataFrame(y)" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "## DummyForecaster\n", + "\n", + "A dummy forecaster can illustrate the use cases for forecasting. This\n", + "forecaster simply returns the last value of the train data for the forecast. By\n", + "default the horizon is 1. It makes no difference for this forecaster. It's inner type\n", + " is `np.ndarray` so all three allowable input types are internally converted to numpy\n", + " arrays." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 4, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "np.ndarray\n", + "432.0\n" + ] + } + ], + "source": [ + "# Fit then predict\n", + "from aeon.forecasting import DummyForecaster\n", + "\n", + "d = DummyForecaster()\n", + "print(d.get_tag(\"y_inner_type\"))\n", + "d.fit(y)\n", + "p = d.predict()\n", + "print(p)" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 5, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "432.0\n" + ] + } + ], + "source": [ + "# forecast is equivalent to fit_predict in other estimators\n", + "p2 = d.forecast(y)\n", + "print(p2)" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "## Regression based forecasting\n", + "\n", + "Our main focus will be forecasting through a sliding window and a regressor. We\n", + "provide a basic implementation of this in `RegressionForecaster`. This class can take\n", + " a regressor as a constructor parameter. It will train the regressor on the windowed\n", + " series, then apply the data to new series. There will be a notebook for more details\n", + " of the use of RegressionForecaster. We just use it to demonstrate different horizons" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 6, + "outputs": [ + { + "ename": "ValueError", + "evalue": "Expected 2D array, got 1D array instead:\narray=[420. 472. 548. 559. 463. 407. 362. 405. 417. 391. 419. 461. 472. 535.\n 622. 606. 508. 461. 390. 432.].\nReshape your data either using array.reshape(-1, 1) if your data has a single feature or array.reshape(1, -1) if it contains a single sample.", + "output_type": "error", + "traceback": [ + "\u001B[1;31m---------------------------------------------------------------------------\u001B[0m", + "\u001B[1;31mValueError\u001B[0m Traceback (most recent call last)", + "Cell \u001B[1;32mIn[6], line 4\u001B[0m\n\u001B[0;32m 2\u001B[0m r \u001B[38;5;241m=\u001B[39m RegressionForecaster(window \u001B[38;5;241m=\u001B[39m \u001B[38;5;241m20\u001B[39m)\n\u001B[0;32m 3\u001B[0m r\u001B[38;5;241m.\u001B[39mfit(y)\n\u001B[1;32m----> 4\u001B[0m p \u001B[38;5;241m=\u001B[39m \u001B[43mr\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mpredict\u001B[49m\u001B[43m(\u001B[49m\u001B[43m)\u001B[49m\n\u001B[0;32m 5\u001B[0m \u001B[38;5;28mprint\u001B[39m(p)\n\u001B[0;32m 6\u001B[0m r2 \u001B[38;5;241m=\u001B[39m RegressionForecaster(window \u001B[38;5;241m=\u001B[39m \u001B[38;5;241m20\u001B[39m, horizon \u001B[38;5;241m=\u001B[39m \u001B[38;5;241m5\u001B[39m)\n", + "File \u001B[1;32mC:\\Code\\aeon\\aeon\\forecasting\\base.py:105\u001B[0m, in \u001B[0;36mBaseForecaster.predict\u001B[1;34m(self, y, exog)\u001B[0m\n\u001B[0;32m 103\u001B[0m \u001B[38;5;66;03m# Validate exog\u001B[39;00m\n\u001B[0;32m 104\u001B[0m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39m_is_fitted \u001B[38;5;241m=\u001B[39m \u001B[38;5;28;01mTrue\u001B[39;00m\n\u001B[1;32m--> 105\u001B[0m \u001B[38;5;28;01mreturn\u001B[39;00m \u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43m_predict\u001B[49m\u001B[43m(\u001B[49m\u001B[43my\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mexog\u001B[49m\u001B[43m)\u001B[49m\n", + "File \u001B[1;32mC:\\Code\\aeon\\aeon\\forecasting\\_regression.py:52\u001B[0m, in \u001B[0;36mRegressionForecaster._predict\u001B[1;34m(self, y, exog)\u001B[0m\n\u001B[0;32m 50\u001B[0m \u001B[38;5;250m\u001B[39m\u001B[38;5;124;03m\"\"\"Predict values for time series X.\"\"\"\u001B[39;00m\n\u001B[0;32m 51\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m y \u001B[38;5;129;01mis\u001B[39;00m \u001B[38;5;28;01mNone\u001B[39;00m:\n\u001B[1;32m---> 52\u001B[0m \u001B[38;5;28;01mreturn\u001B[39;00m \u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mregressor_\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mpredict\u001B[49m\u001B[43m(\u001B[49m\u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mlast_\u001B[49m\u001B[43m)\u001B[49m\n\u001B[0;32m 53\u001B[0m y \u001B[38;5;241m=\u001B[39m y\u001B[38;5;241m.\u001B[39msqueeze()\n\u001B[0;32m 54\u001B[0m \u001B[38;5;28;01mreturn\u001B[39;00m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mregressor_\u001B[38;5;241m.\u001B[39mpredict(y[\u001B[38;5;241m-\u001B[39m\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mwindow :])\n", + "File \u001B[1;32mC:\\Code\\aeon\\venv\\lib\\site-packages\\sklearn\\linear_model\\_base.py:306\u001B[0m, in \u001B[0;36mLinearModel.predict\u001B[1;34m(self, X)\u001B[0m\n\u001B[0;32m 292\u001B[0m \u001B[38;5;28;01mdef\u001B[39;00m \u001B[38;5;21mpredict\u001B[39m(\u001B[38;5;28mself\u001B[39m, X):\n\u001B[0;32m 293\u001B[0m \u001B[38;5;250m \u001B[39m\u001B[38;5;124;03m\"\"\"\u001B[39;00m\n\u001B[0;32m 294\u001B[0m \u001B[38;5;124;03m Predict using the linear model.\u001B[39;00m\n\u001B[0;32m 295\u001B[0m \n\u001B[1;32m (...)\u001B[0m\n\u001B[0;32m 304\u001B[0m \u001B[38;5;124;03m Returns predicted values.\u001B[39;00m\n\u001B[0;32m 305\u001B[0m \u001B[38;5;124;03m \"\"\"\u001B[39;00m\n\u001B[1;32m--> 306\u001B[0m \u001B[38;5;28;01mreturn\u001B[39;00m \u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43m_decision_function\u001B[49m\u001B[43m(\u001B[49m\u001B[43mX\u001B[49m\u001B[43m)\u001B[49m\n", + "File \u001B[1;32mC:\\Code\\aeon\\venv\\lib\\site-packages\\sklearn\\linear_model\\_base.py:285\u001B[0m, in \u001B[0;36mLinearModel._decision_function\u001B[1;34m(self, X)\u001B[0m\n\u001B[0;32m 282\u001B[0m \u001B[38;5;28;01mdef\u001B[39;00m \u001B[38;5;21m_decision_function\u001B[39m(\u001B[38;5;28mself\u001B[39m, X):\n\u001B[0;32m 283\u001B[0m check_is_fitted(\u001B[38;5;28mself\u001B[39m)\n\u001B[1;32m--> 285\u001B[0m X \u001B[38;5;241m=\u001B[39m \u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43m_validate_data\u001B[49m\u001B[43m(\u001B[49m\u001B[43mX\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43maccept_sparse\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43m[\u001B[49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[38;5;124;43mcsr\u001B[39;49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[38;5;124;43mcsc\u001B[39;49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[38;5;124;43mcoo\u001B[39;49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[43m]\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mreset\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[38;5;28;43;01mFalse\u001B[39;49;00m\u001B[43m)\u001B[49m\n\u001B[0;32m 286\u001B[0m coef_ \u001B[38;5;241m=\u001B[39m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mcoef_\n\u001B[0;32m 287\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m coef_\u001B[38;5;241m.\u001B[39mndim \u001B[38;5;241m==\u001B[39m \u001B[38;5;241m1\u001B[39m:\n", + "File \u001B[1;32mC:\\Code\\aeon\\venv\\lib\\site-packages\\sklearn\\base.py:633\u001B[0m, in \u001B[0;36mBaseEstimator._validate_data\u001B[1;34m(self, X, y, reset, validate_separately, cast_to_ndarray, **check_params)\u001B[0m\n\u001B[0;32m 631\u001B[0m out \u001B[38;5;241m=\u001B[39m X, y\n\u001B[0;32m 632\u001B[0m \u001B[38;5;28;01melif\u001B[39;00m \u001B[38;5;129;01mnot\u001B[39;00m no_val_X \u001B[38;5;129;01mand\u001B[39;00m no_val_y:\n\u001B[1;32m--> 633\u001B[0m out \u001B[38;5;241m=\u001B[39m check_array(X, input_name\u001B[38;5;241m=\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mX\u001B[39m\u001B[38;5;124m\"\u001B[39m, \u001B[38;5;241m*\u001B[39m\u001B[38;5;241m*\u001B[39mcheck_params)\n\u001B[0;32m 634\u001B[0m \u001B[38;5;28;01melif\u001B[39;00m no_val_X \u001B[38;5;129;01mand\u001B[39;00m \u001B[38;5;129;01mnot\u001B[39;00m no_val_y:\n\u001B[0;32m 635\u001B[0m out \u001B[38;5;241m=\u001B[39m _check_y(y, \u001B[38;5;241m*\u001B[39m\u001B[38;5;241m*\u001B[39mcheck_params)\n", + "File \u001B[1;32mC:\\Code\\aeon\\venv\\lib\\site-packages\\sklearn\\utils\\validation.py:1050\u001B[0m, in \u001B[0;36mcheck_array\u001B[1;34m(array, accept_sparse, accept_large_sparse, dtype, order, copy, force_writeable, force_all_finite, ensure_2d, allow_nd, ensure_min_samples, ensure_min_features, estimator, input_name)\u001B[0m\n\u001B[0;32m 1043\u001B[0m \u001B[38;5;28;01melse\u001B[39;00m:\n\u001B[0;32m 1044\u001B[0m msg \u001B[38;5;241m=\u001B[39m (\n\u001B[0;32m 1045\u001B[0m \u001B[38;5;124mf\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mExpected 2D array, got 1D array instead:\u001B[39m\u001B[38;5;130;01m\\n\u001B[39;00m\u001B[38;5;124marray=\u001B[39m\u001B[38;5;132;01m{\u001B[39;00marray\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m.\u001B[39m\u001B[38;5;130;01m\\n\u001B[39;00m\u001B[38;5;124m\"\u001B[39m\n\u001B[0;32m 1046\u001B[0m \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mReshape your data either using array.reshape(-1, 1) if \u001B[39m\u001B[38;5;124m\"\u001B[39m\n\u001B[0;32m 1047\u001B[0m \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124myour data has a single feature or array.reshape(1, -1) \u001B[39m\u001B[38;5;124m\"\u001B[39m\n\u001B[0;32m 1048\u001B[0m \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mif it contains a single sample.\u001B[39m\u001B[38;5;124m\"\u001B[39m\n\u001B[0;32m 1049\u001B[0m )\n\u001B[1;32m-> 1050\u001B[0m \u001B[38;5;28;01mraise\u001B[39;00m \u001B[38;5;167;01mValueError\u001B[39;00m(msg)\n\u001B[0;32m 1052\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m dtype_numeric \u001B[38;5;129;01mand\u001B[39;00m \u001B[38;5;28mhasattr\u001B[39m(array\u001B[38;5;241m.\u001B[39mdtype, \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mkind\u001B[39m\u001B[38;5;124m\"\u001B[39m) \u001B[38;5;129;01mand\u001B[39;00m array\u001B[38;5;241m.\u001B[39mdtype\u001B[38;5;241m.\u001B[39mkind \u001B[38;5;129;01min\u001B[39;00m \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mUSV\u001B[39m\u001B[38;5;124m\"\u001B[39m:\n\u001B[0;32m 1053\u001B[0m \u001B[38;5;28;01mraise\u001B[39;00m \u001B[38;5;167;01mValueError\u001B[39;00m(\n\u001B[0;32m 1054\u001B[0m \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mdtype=\u001B[39m\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mnumeric\u001B[39m\u001B[38;5;124m'\u001B[39m\u001B[38;5;124m is not compatible with arrays of bytes/strings.\u001B[39m\u001B[38;5;124m\"\u001B[39m\n\u001B[0;32m 1055\u001B[0m \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mConvert your data to numeric values explicitly instead.\u001B[39m\u001B[38;5;124m\"\u001B[39m\n\u001B[0;32m 1056\u001B[0m )\n", + "\u001B[1;31mValueError\u001B[0m: Expected 2D array, got 1D array instead:\narray=[420. 472. 548. 559. 463. 407. 362. 405. 417. 391. 419. 461. 472. 535.\n 622. 606. 508. 461. 390. 432.].\nReshape your data either using array.reshape(-1, 1) if your data has a single feature or array.reshape(1, -1) if it contains a single sample." + ] + } + ], + "source": [ + "from aeon.forecasting import RegressionForecaster\n", + "\n", + "r = RegressionForecaster(window=20)\n", + "r.fit(y)\n", + "p = r.predict()\n", + "print(p)\n", + "r2 = RegressionForecaster(window=20, horizon=5)\n", + "r.fit(y)\n", + "p = r.predict()\n", + "print(p)" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [], + "metadata": { + "collapsed": false + } + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3", + "language": "python", + "name": "python3" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 2 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython2", + "version": "2.7.6" + } + }, + "nbformat": 4, + "nbformat_minor": 0 +} From 7d932a0da7506cac6439837ed6d4c2c2bddb77b5 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Thu, 31 Oct 2024 17:35:59 +0000 Subject: [PATCH 08/31] regressor --- aeon/forecasting/_regression.py | 20 ++++++++++++++++++-- 1 file changed, 18 insertions(+), 2 deletions(-) diff --git a/aeon/forecasting/_regression.py b/aeon/forecasting/_regression.py index bc01213d17..921591133e 100644 --- a/aeon/forecasting/_regression.py +++ b/aeon/forecasting/_regression.py @@ -43,6 +43,7 @@ def _fit(self, y, exog=None): # Extract y y = y[self.window + self.horizon - 1 :] self.last_ = y[-self.window :] + self.last_ = self.last_.reshape(1, -1) self.regressor_.fit(X=X, y=y) return self @@ -50,8 +51,8 @@ def _predict(self, y=None, exog=None): """Predict values for time series X.""" if y is None: return self.regressor_.predict(self.last_) - y = y.squeeze() - return self.regressor_.predict(y[-self.window :]) + last = y[:, -self.window :] + return self.regressor_.predict(last) def _forecast(self, y, exog=None): """Forecast values for time series X. @@ -60,3 +61,18 @@ def _forecast(self, y, exog=None): """ self.fit(y, exog) return self.predict(y) + + def _get_test_params(cls, parameter_set="default"): + """Return testing parameter settings for the estimator. + + Parameters + ---------- + parameter_set : str, default='default' + Name of the parameter set to return. + + Returns + ------- + dict + Dictionary of testing parameter settings. + """ + return {"window": 4} From b89f165fa9fc14f8dde1c1b8245acfa1290c7c90 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Thu, 31 Oct 2024 17:42:52 +0000 Subject: [PATCH 09/31] regressor --- aeon/forecasting/_regression.py | 27 ++++++++++++++++++++++++++- aeon/forecasting/base.py | 9 +-------- 2 files changed, 27 insertions(+), 9 deletions(-) diff --git a/aeon/forecasting/_regression.py b/aeon/forecasting/_regression.py index 921591133e..992013f181 100644 --- a/aeon/forecasting/_regression.py +++ b/aeon/forecasting/_regression.py @@ -12,9 +12,34 @@ class RegressionForecaster(BaseForecaster): + """ + Regression based forecasting. + + Container for forecaster that reduces forecasting to regression through a + window. Form a collection of sub series of length `window` through a sliding + winodw to form X, take `horizon` points ahead to form `y`, then apply an aeon or + sklearn regressor. + + + Parameters + ---------- + window : int + The window prior to the current time point to use in forecasting. So if + horizon is one, forecaster will train using points $i$ to $window+i-1$ to + predict value $window+i$. If horizon is 4, forecaster will used points $i$ + to $window+i-1$ to predict value $window+i+3$. If None, the algorithm will + internally determine what data to use to predict `horizon` steps ahead. + horizon : int, default =1 + The number of time steps ahead to forecast. If horizon is one, the forecaster + will learn to predict one point ahead + regressor : object, default =None + Regression estimator that implements BaseRegressor or is otherwise compatible + with sklearn regressors. + """ + def __init__(self, window, horizon=1, regressor=None): self.regressor = regressor - super().__init__(horizon, window) + super().__init__(horizon, axis=1) def _fit(self, y, exog=None): """Fit forecaster to time series. diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index 9c213a2c41..70cabc7c9a 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -26,12 +26,6 @@ class BaseForecaster(BaseSeriesEstimator): horizon : int, default =1 The number of time steps ahead to forecast. If horizon is one, the forecaster will learn to predict one point ahead - window : int or None - The window prior to the current time point to use in forecasting. So if - horizon is one, forecaster will train using points $i$ to $window+i-1$ to - predict value $window+i$. If horizon is 4, forecaster will used points $i$ - to $window+i-1$ to predict value $window+i+3$. If None, the algorithm will - internally determine what data to use to predict `horizon` steps ahead. """ # TODO: add any forecasting specific tags @@ -42,9 +36,8 @@ class BaseForecaster(BaseSeriesEstimator): "y_inner_type": "np.ndarray", } - def __init__(self, horizon=1, window=None, axis=1): + def __init__(self, horizon=1, axis=1): self.horizon = horizon - self.window = window self._is_fitted = False super().__init__(axis) From 72fe6313b82a4d1990533f584a1a46d22ddaaff7 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Thu, 31 Oct 2024 17:50:32 +0000 Subject: [PATCH 10/31] regressor --- aeon/forecasting/_regression.py | 1 + 1 file changed, 1 insertion(+) diff --git a/aeon/forecasting/_regression.py b/aeon/forecasting/_regression.py index 992013f181..7f576d2990 100644 --- a/aeon/forecasting/_regression.py +++ b/aeon/forecasting/_regression.py @@ -87,6 +87,7 @@ def _forecast(self, y, exog=None): self.fit(y, exog) return self.predict(y) + @classmethod def _get_test_params(cls, parameter_set="default"): """Return testing parameter settings for the estimator. From 80f8bd521f58f0304a9f80e8fe10a9fae6e96f30 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Thu, 31 Oct 2024 18:00:08 +0000 Subject: [PATCH 11/31] tags --- aeon/utils/tags/_tags.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/aeon/utils/tags/_tags.py b/aeon/utils/tags/_tags.py index 8be4d787ee..b5d4bc4d7f 100644 --- a/aeon/utils/tags/_tags.py +++ b/aeon/utils/tags/_tags.py @@ -10,6 +10,7 @@ class : identifier for the base class of objects this tag applies to ("str", list_of_string) - any string in list_of_string is valid ("list", list_of_string) - any sub-list is valid ("list||str", list_of_string) - combination of the above + ("list||str", list_of_string) - combination of the above None - no value for the tag description : plain English description of the tag """ @@ -54,7 +55,7 @@ class : identifier for the base class of objects this tag applies to "fit/predict.", }, "y_inner_type": { - "class": "estimator", + "class": "forecaster", "type": [ ("list||str", SERIES_DATA_TYPES), ], From f4f828ff50aa80b2c89f364da5aa31bf9b636a59 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Thu, 31 Oct 2024 18:11:31 +0000 Subject: [PATCH 12/31] tags --- aeon/forecasting/base.py | 1 + 1 file changed, 1 insertion(+) diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index 70cabc7c9a..a18a58dfd6 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -34,6 +34,7 @@ class BaseForecaster(BaseSeriesEstimator): "capability:multivariate": False, "capability:missing_values": False, "y_inner_type": "np.ndarray", + "fit_is_empty": False, } def __init__(self, horizon=1, axis=1): From aeb8bbbaf9682a4aa9e8de879a4e5d44dbefe58b Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Thu, 31 Oct 2024 18:16:42 +0000 Subject: [PATCH 13/31] requires_y --- aeon/forecasting/base.py | 2 +- aeon/utils/tags/_tags.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index a18a58dfd6..b109d915e2 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -28,13 +28,13 @@ class BaseForecaster(BaseSeriesEstimator): will learn to predict one point ahead """ - # TODO: add any forecasting specific tags _tags = { "capability:univariate": True, "capability:multivariate": False, "capability:missing_values": False, "y_inner_type": "np.ndarray", "fit_is_empty": False, + "requires_y": True, } def __init__(self, horizon=1, axis=1): diff --git a/aeon/utils/tags/_tags.py b/aeon/utils/tags/_tags.py index b5d4bc4d7f..55a5b99897 100644 --- a/aeon/utils/tags/_tags.py +++ b/aeon/utils/tags/_tags.py @@ -139,7 +139,7 @@ class : identifier for the base class of objects this tag applies to "point belongs to.", }, "requires_y": { - "class": ["transformer", "anomaly-detector", "segmenter"], + "class": ["transformer", "anomaly-detector", "segmenter", "forecaster"], "type": "bool", "description": "Does this estimator require y to be passed in its methods?", }, From f8209842e5381c885f40eb8ecb4ec53b43e63479 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Fri, 1 Nov 2024 18:15:55 +0000 Subject: [PATCH 14/31] forecasting notebook --- aeon/forecasting/_regression.py | 1 + examples/forecasting/forecasting.ipynb | 81 ++++++++++++++++++-------- 2 files changed, 58 insertions(+), 24 deletions(-) diff --git a/aeon/forecasting/_regression.py b/aeon/forecasting/_regression.py index 7f576d2990..d113775eef 100644 --- a/aeon/forecasting/_regression.py +++ b/aeon/forecasting/_regression.py @@ -38,6 +38,7 @@ class RegressionForecaster(BaseForecaster): """ def __init__(self, window, horizon=1, regressor=None): + self.window = window self.regressor = regressor super().__init__(horizon, axis=1) diff --git a/examples/forecasting/forecasting.ipynb b/examples/forecasting/forecasting.ipynb index 8db44440a5..d53ec82f78 100644 --- a/examples/forecasting/forecasting.ipynb +++ b/examples/forecasting/forecasting.ipynb @@ -63,7 +63,7 @@ }, { "cell_type": "code", - "execution_count": 1, + "execution_count": 13, "outputs": [ { "name": "stdout", @@ -106,7 +106,7 @@ }, { "cell_type": "code", - "execution_count": 2, + "execution_count": 14, "outputs": [ { "name": "stdout", @@ -114,7 +114,7 @@ "text": [ " Possible data structures for input to forecaster ['pd.Series', 'pd.DataFrame', 'np.ndarray']\n", "\n", - " Tags for BaseForecaster: {'python_version': None, 'python_dependencies': None, 'cant_pickle': False, 'non_deterministic': False, 'algorithm_type': None, 'capability:missing_values': False, 'capability:multithreading': False, 'capability:univariate': True, 'capability:multivariate': False, 'X_inner_type': 'np.ndarray', 'y_inner_type': 'np.ndarray'}\n" + " Tags for BaseForecaster: {'python_version': None, 'python_dependencies': None, 'cant_pickle': False, 'non_deterministic': False, 'algorithm_type': None, 'capability:missing_values': False, 'capability:multithreading': False, 'capability:univariate': True, 'capability:multivariate': False, 'X_inner_type': 'np.ndarray', 'y_inner_type': 'np.ndarray', 'fit_is_empty': False, 'requires_y': True}\n" ] } ], @@ -140,7 +140,7 @@ }, { "cell_type": "code", - "execution_count": 3, + "execution_count": 15, "outputs": [ { "name": "stdout", @@ -181,7 +181,7 @@ }, { "cell_type": "code", - "execution_count": 4, + "execution_count": 16, "outputs": [ { "name": "stdout", @@ -208,7 +208,7 @@ }, { "cell_type": "code", - "execution_count": 5, + "execution_count": 17, "outputs": [ { "name": "stdout", @@ -244,23 +244,15 @@ }, { "cell_type": "code", - "execution_count": 6, + "execution_count": 19, "outputs": [ { - "ename": "ValueError", - "evalue": "Expected 2D array, got 1D array instead:\narray=[420. 472. 548. 559. 463. 407. 362. 405. 417. 391. 419. 461. 472. 535.\n 622. 606. 508. 461. 390. 432.].\nReshape your data either using array.reshape(-1, 1) if your data has a single feature or array.reshape(1, -1) if it contains a single sample.", - "output_type": "error", - "traceback": [ - "\u001B[1;31m---------------------------------------------------------------------------\u001B[0m", - "\u001B[1;31mValueError\u001B[0m Traceback (most recent call last)", - "Cell \u001B[1;32mIn[6], line 4\u001B[0m\n\u001B[0;32m 2\u001B[0m r \u001B[38;5;241m=\u001B[39m RegressionForecaster(window \u001B[38;5;241m=\u001B[39m \u001B[38;5;241m20\u001B[39m)\n\u001B[0;32m 3\u001B[0m r\u001B[38;5;241m.\u001B[39mfit(y)\n\u001B[1;32m----> 4\u001B[0m p \u001B[38;5;241m=\u001B[39m \u001B[43mr\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mpredict\u001B[49m\u001B[43m(\u001B[49m\u001B[43m)\u001B[49m\n\u001B[0;32m 5\u001B[0m \u001B[38;5;28mprint\u001B[39m(p)\n\u001B[0;32m 6\u001B[0m r2 \u001B[38;5;241m=\u001B[39m RegressionForecaster(window \u001B[38;5;241m=\u001B[39m \u001B[38;5;241m20\u001B[39m, horizon \u001B[38;5;241m=\u001B[39m \u001B[38;5;241m5\u001B[39m)\n", - "File \u001B[1;32mC:\\Code\\aeon\\aeon\\forecasting\\base.py:105\u001B[0m, in \u001B[0;36mBaseForecaster.predict\u001B[1;34m(self, y, exog)\u001B[0m\n\u001B[0;32m 103\u001B[0m \u001B[38;5;66;03m# Validate exog\u001B[39;00m\n\u001B[0;32m 104\u001B[0m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39m_is_fitted \u001B[38;5;241m=\u001B[39m \u001B[38;5;28;01mTrue\u001B[39;00m\n\u001B[1;32m--> 105\u001B[0m \u001B[38;5;28;01mreturn\u001B[39;00m \u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43m_predict\u001B[49m\u001B[43m(\u001B[49m\u001B[43my\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mexog\u001B[49m\u001B[43m)\u001B[49m\n", - "File \u001B[1;32mC:\\Code\\aeon\\aeon\\forecasting\\_regression.py:52\u001B[0m, in \u001B[0;36mRegressionForecaster._predict\u001B[1;34m(self, y, exog)\u001B[0m\n\u001B[0;32m 50\u001B[0m \u001B[38;5;250m\u001B[39m\u001B[38;5;124;03m\"\"\"Predict values for time series X.\"\"\"\u001B[39;00m\n\u001B[0;32m 51\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m y \u001B[38;5;129;01mis\u001B[39;00m \u001B[38;5;28;01mNone\u001B[39;00m:\n\u001B[1;32m---> 52\u001B[0m \u001B[38;5;28;01mreturn\u001B[39;00m \u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mregressor_\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mpredict\u001B[49m\u001B[43m(\u001B[49m\u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mlast_\u001B[49m\u001B[43m)\u001B[49m\n\u001B[0;32m 53\u001B[0m y \u001B[38;5;241m=\u001B[39m y\u001B[38;5;241m.\u001B[39msqueeze()\n\u001B[0;32m 54\u001B[0m \u001B[38;5;28;01mreturn\u001B[39;00m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mregressor_\u001B[38;5;241m.\u001B[39mpredict(y[\u001B[38;5;241m-\u001B[39m\u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mwindow :])\n", - "File \u001B[1;32mC:\\Code\\aeon\\venv\\lib\\site-packages\\sklearn\\linear_model\\_base.py:306\u001B[0m, in \u001B[0;36mLinearModel.predict\u001B[1;34m(self, X)\u001B[0m\n\u001B[0;32m 292\u001B[0m \u001B[38;5;28;01mdef\u001B[39;00m \u001B[38;5;21mpredict\u001B[39m(\u001B[38;5;28mself\u001B[39m, X):\n\u001B[0;32m 293\u001B[0m \u001B[38;5;250m \u001B[39m\u001B[38;5;124;03m\"\"\"\u001B[39;00m\n\u001B[0;32m 294\u001B[0m \u001B[38;5;124;03m Predict using the linear model.\u001B[39;00m\n\u001B[0;32m 295\u001B[0m \n\u001B[1;32m (...)\u001B[0m\n\u001B[0;32m 304\u001B[0m \u001B[38;5;124;03m Returns predicted values.\u001B[39;00m\n\u001B[0;32m 305\u001B[0m \u001B[38;5;124;03m \"\"\"\u001B[39;00m\n\u001B[1;32m--> 306\u001B[0m \u001B[38;5;28;01mreturn\u001B[39;00m \u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43m_decision_function\u001B[49m\u001B[43m(\u001B[49m\u001B[43mX\u001B[49m\u001B[43m)\u001B[49m\n", - "File \u001B[1;32mC:\\Code\\aeon\\venv\\lib\\site-packages\\sklearn\\linear_model\\_base.py:285\u001B[0m, in \u001B[0;36mLinearModel._decision_function\u001B[1;34m(self, X)\u001B[0m\n\u001B[0;32m 282\u001B[0m \u001B[38;5;28;01mdef\u001B[39;00m \u001B[38;5;21m_decision_function\u001B[39m(\u001B[38;5;28mself\u001B[39m, X):\n\u001B[0;32m 283\u001B[0m check_is_fitted(\u001B[38;5;28mself\u001B[39m)\n\u001B[1;32m--> 285\u001B[0m X \u001B[38;5;241m=\u001B[39m \u001B[38;5;28;43mself\u001B[39;49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43m_validate_data\u001B[49m\u001B[43m(\u001B[49m\u001B[43mX\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43maccept_sparse\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[43m[\u001B[49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[38;5;124;43mcsr\u001B[39;49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[38;5;124;43mcsc\u001B[39;49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[38;5;124;43mcoo\u001B[39;49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[43m]\u001B[49m\u001B[43m,\u001B[49m\u001B[43m \u001B[49m\u001B[43mreset\u001B[49m\u001B[38;5;241;43m=\u001B[39;49m\u001B[38;5;28;43;01mFalse\u001B[39;49;00m\u001B[43m)\u001B[49m\n\u001B[0;32m 286\u001B[0m coef_ \u001B[38;5;241m=\u001B[39m \u001B[38;5;28mself\u001B[39m\u001B[38;5;241m.\u001B[39mcoef_\n\u001B[0;32m 287\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m coef_\u001B[38;5;241m.\u001B[39mndim \u001B[38;5;241m==\u001B[39m \u001B[38;5;241m1\u001B[39m:\n", - "File \u001B[1;32mC:\\Code\\aeon\\venv\\lib\\site-packages\\sklearn\\base.py:633\u001B[0m, in \u001B[0;36mBaseEstimator._validate_data\u001B[1;34m(self, X, y, reset, validate_separately, cast_to_ndarray, **check_params)\u001B[0m\n\u001B[0;32m 631\u001B[0m out \u001B[38;5;241m=\u001B[39m X, y\n\u001B[0;32m 632\u001B[0m \u001B[38;5;28;01melif\u001B[39;00m \u001B[38;5;129;01mnot\u001B[39;00m no_val_X \u001B[38;5;129;01mand\u001B[39;00m no_val_y:\n\u001B[1;32m--> 633\u001B[0m out \u001B[38;5;241m=\u001B[39m check_array(X, input_name\u001B[38;5;241m=\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mX\u001B[39m\u001B[38;5;124m\"\u001B[39m, \u001B[38;5;241m*\u001B[39m\u001B[38;5;241m*\u001B[39mcheck_params)\n\u001B[0;32m 634\u001B[0m \u001B[38;5;28;01melif\u001B[39;00m no_val_X \u001B[38;5;129;01mand\u001B[39;00m \u001B[38;5;129;01mnot\u001B[39;00m no_val_y:\n\u001B[0;32m 635\u001B[0m out \u001B[38;5;241m=\u001B[39m _check_y(y, \u001B[38;5;241m*\u001B[39m\u001B[38;5;241m*\u001B[39mcheck_params)\n", - "File \u001B[1;32mC:\\Code\\aeon\\venv\\lib\\site-packages\\sklearn\\utils\\validation.py:1050\u001B[0m, in \u001B[0;36mcheck_array\u001B[1;34m(array, accept_sparse, accept_large_sparse, dtype, order, copy, force_writeable, force_all_finite, ensure_2d, allow_nd, ensure_min_samples, ensure_min_features, estimator, input_name)\u001B[0m\n\u001B[0;32m 1043\u001B[0m \u001B[38;5;28;01melse\u001B[39;00m:\n\u001B[0;32m 1044\u001B[0m msg \u001B[38;5;241m=\u001B[39m (\n\u001B[0;32m 1045\u001B[0m \u001B[38;5;124mf\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mExpected 2D array, got 1D array instead:\u001B[39m\u001B[38;5;130;01m\\n\u001B[39;00m\u001B[38;5;124marray=\u001B[39m\u001B[38;5;132;01m{\u001B[39;00marray\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m.\u001B[39m\u001B[38;5;130;01m\\n\u001B[39;00m\u001B[38;5;124m\"\u001B[39m\n\u001B[0;32m 1046\u001B[0m \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mReshape your data either using array.reshape(-1, 1) if \u001B[39m\u001B[38;5;124m\"\u001B[39m\n\u001B[0;32m 1047\u001B[0m \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124myour data has a single feature or array.reshape(1, -1) \u001B[39m\u001B[38;5;124m\"\u001B[39m\n\u001B[0;32m 1048\u001B[0m \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mif it contains a single sample.\u001B[39m\u001B[38;5;124m\"\u001B[39m\n\u001B[0;32m 1049\u001B[0m )\n\u001B[1;32m-> 1050\u001B[0m \u001B[38;5;28;01mraise\u001B[39;00m \u001B[38;5;167;01mValueError\u001B[39;00m(msg)\n\u001B[0;32m 1052\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m dtype_numeric \u001B[38;5;129;01mand\u001B[39;00m \u001B[38;5;28mhasattr\u001B[39m(array\u001B[38;5;241m.\u001B[39mdtype, \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mkind\u001B[39m\u001B[38;5;124m\"\u001B[39m) \u001B[38;5;129;01mand\u001B[39;00m array\u001B[38;5;241m.\u001B[39mdtype\u001B[38;5;241m.\u001B[39mkind \u001B[38;5;129;01min\u001B[39;00m \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mUSV\u001B[39m\u001B[38;5;124m\"\u001B[39m:\n\u001B[0;32m 1053\u001B[0m \u001B[38;5;28;01mraise\u001B[39;00m \u001B[38;5;167;01mValueError\u001B[39;00m(\n\u001B[0;32m 1054\u001B[0m \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mdtype=\u001B[39m\u001B[38;5;124m'\u001B[39m\u001B[38;5;124mnumeric\u001B[39m\u001B[38;5;124m'\u001B[39m\u001B[38;5;124m is not compatible with arrays of bytes/strings.\u001B[39m\u001B[38;5;124m\"\u001B[39m\n\u001B[0;32m 1055\u001B[0m \u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mConvert your data to numeric values explicitly instead.\u001B[39m\u001B[38;5;124m\"\u001B[39m\n\u001B[0;32m 1056\u001B[0m )\n", - "\u001B[1;31mValueError\u001B[0m: Expected 2D array, got 1D array instead:\narray=[420. 472. 548. 559. 463. 407. 362. 405. 417. 391. 419. 461. 472. 535.\n 622. 606. 508. 461. 390. 432.].\nReshape your data either using array.reshape(-1, 1) if your data has a single feature or array.reshape(1, -1) if it contains a single sample." + "name": "stdout", + "output_type": "stream", + "text": [ + "[451.67541971]\n", + "[505.91690237]\n", + "[505.91690237]\n" ] } ], @@ -272,8 +264,10 @@ "p = r.predict()\n", "print(p)\n", "r2 = RegressionForecaster(window=20, horizon=5)\n", - "r.fit(y)\n", - "p = r.predict()\n", + "r2.fit(y)\n", + "p = r2.predict()\n", + "print(p)\n", + "p = r2.forecast(y)\n", "print(p)" ], "metadata": { @@ -282,7 +276,46 @@ }, { "cell_type": "markdown", - "source": [], + "source": [ + "With our set up, we can make predictions with previously unseen data, thus more\n", + "closely modelling machine learning approaches. Or we can use the forecast method to\n", + "fit/predict at the same time." + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "code", + "execution_count": 20, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "[5.38765716] ,\n", + " [17.05612701]\n" + ] + } + ], + "source": [ + "import numpy as np\n", + "\n", + "y = np.random.rand(20)\n", + "p1 = r.predict(y)\n", + "p2 = r2.predict(y)\n", + "print(p1, \",\\n\", p2)" + ], + "metadata": { + "collapsed": false + } + }, + { + "cell_type": "markdown", + "source": [ + "This is our first step in the forecasting module. The next will be to add traditional\n", + " stats models in the same context." + ], "metadata": { "collapsed": false } From 5f3775798356a3efca9bbd6eb82e19264727be7f Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Fri, 1 Nov 2024 18:22:16 +0000 Subject: [PATCH 15/31] forecasting notebook --- aeon/forecasting/base.py | 1 - aeon/testing/mock_estimators/_mock_forecasters.py | 2 +- aeon/utils/tags/_tags.py | 1 - 3 files changed, 1 insertion(+), 3 deletions(-) diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index b109d915e2..0b1bd14812 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -39,7 +39,6 @@ class BaseForecaster(BaseSeriesEstimator): def __init__(self, horizon=1, axis=1): self.horizon = horizon - self._is_fitted = False super().__init__(axis) def fit(self, y, exog=None): diff --git a/aeon/testing/mock_estimators/_mock_forecasters.py b/aeon/testing/mock_estimators/_mock_forecasters.py index e5b5403b08..f5bb86d249 100644 --- a/aeon/testing/mock_estimators/_mock_forecasters.py +++ b/aeon/testing/mock_estimators/_mock_forecasters.py @@ -7,7 +7,7 @@ class MockForecaster(BaseForecaster): - """Mock segmenter for testing.""" + """Mock forecaster for testing.""" def __init__(self): super().__init__() diff --git a/aeon/utils/tags/_tags.py b/aeon/utils/tags/_tags.py index 55a5b99897..5fe18d5162 100644 --- a/aeon/utils/tags/_tags.py +++ b/aeon/utils/tags/_tags.py @@ -10,7 +10,6 @@ class : identifier for the base class of objects this tag applies to ("str", list_of_string) - any string in list_of_string is valid ("list", list_of_string) - any sub-list is valid ("list||str", list_of_string) - combination of the above - ("list||str", list_of_string) - combination of the above None - no value for the tag description : plain English description of the tag """ From 560f664b9d3031f2a4ae265483090b9dfa8b0e87 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Sat, 2 Nov 2024 11:55:49 +0000 Subject: [PATCH 16/31] remove tags --- aeon/forecasting/base.py | 2 -- aeon/utils/tags/_tags.py | 10 +--------- 2 files changed, 1 insertion(+), 11 deletions(-) diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index 0b1bd14812..6837af51ea 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -32,9 +32,7 @@ class BaseForecaster(BaseSeriesEstimator): "capability:univariate": True, "capability:multivariate": False, "capability:missing_values": False, - "y_inner_type": "np.ndarray", "fit_is_empty": False, - "requires_y": True, } def __init__(self, horizon=1, axis=1): diff --git a/aeon/utils/tags/_tags.py b/aeon/utils/tags/_tags.py index 5fe18d5162..650a7c4dc4 100644 --- a/aeon/utils/tags/_tags.py +++ b/aeon/utils/tags/_tags.py @@ -53,14 +53,6 @@ class : identifier for the base class of objects this tag applies to "description": "What data structure(s) the estimator uses internally for " "fit/predict.", }, - "y_inner_type": { - "class": "forecaster", - "type": [ - ("list||str", SERIES_DATA_TYPES), - ], - "description": "What data structure(s) the estimator uses internally for " - "fit/predict.", - }, "algorithm_type": { "class": "estimator", "type": [ @@ -138,7 +130,7 @@ class : identifier for the base class of objects this tag applies to "point belongs to.", }, "requires_y": { - "class": ["transformer", "anomaly-detector", "segmenter", "forecaster"], + "class": ["transformer", "anomaly-detector", "segmenter"], "type": "bool", "description": "Does this estimator require y to be passed in its methods?", }, From 26aec5351004b7b5d72250f3a1185301389f6a74 Mon Sep 17 00:00:00 2001 From: MatthewMiddlehurst Date: Tue, 5 Nov 2024 00:01:18 +0000 Subject: [PATCH 17/31] fix forecasting testing (they still fail though) --- aeon/clustering/deep_learning/_ae_bgru.py | 2 +- aeon/forecasting/base.py | 4 +- aeon/testing/testing_data.py | 105 +++++++++++++--------- aeon/testing/tests/test_testing_data.py | 3 + aeon/testing/utils/estimator_checks.py | 7 +- 5 files changed, 77 insertions(+), 44 deletions(-) diff --git a/aeon/clustering/deep_learning/_ae_bgru.py b/aeon/clustering/deep_learning/_ae_bgru.py index 9b7df32716..269eb90176 100644 --- a/aeon/clustering/deep_learning/_ae_bgru.py +++ b/aeon/clustering/deep_learning/_ae_bgru.py @@ -290,7 +290,7 @@ def _score(self, X, y=None): return self._estimator.score(latent_space) @classmethod - def get_test_params(cls, parameter_set="default"): + def _get_test_params(cls, parameter_set="default"): """Return testing parameter settings for the estimator. Parameters diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index 6837af51ea..48228a8226 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -87,12 +87,12 @@ def predict(self, y=None, exog=None): """ if y is not None: y = self._preprocess_series(y, axis=self.axis, store_metadata=False) - if not self._is_fitted: + if not self.is_fitted: raise ValueError("Forecaster must be fitted before predicting") if exog is not None: raise NotImplementedError("Exogenous variables not yet supported") # Validate exog - self._is_fitted = True + self.is_fitted = True return self._predict(y, exog) @abstractmethod diff --git a/aeon/testing/testing_data.py b/aeon/testing/testing_data.py index 7818527b0c..9318a37193 100644 --- a/aeon/testing/testing_data.py +++ b/aeon/testing/testing_data.py @@ -23,9 +23,11 @@ ) from aeon.transformations.collection import BaseCollectionTransformer from aeon.transformations.series import BaseSeriesTransformer +from aeon.utils.conversion import convert_collection data_rng = np.random.RandomState(42) +# Collection testing data EQUAL_LENGTH_UNIVARIATE_CLASSIFICATION = { "numpy3D": { @@ -590,39 +592,61 @@ "numpy3D": { "train": (X_classification_missing_train, y_classification_missing_train), "test": (X_classification_missing_test, y_classification_missing_test), - } + }, + "np-list": { + "train": ( + convert_collection(X_classification_missing_train, "np-list"), + y_classification_missing_train, + ), + "test": ( + convert_collection(X_classification_missing_test, "np-list"), + y_classification_missing_test, + ), + }, } -X_classification_missing_train, y_classification_missing_train = make_example_3d_numpy( +X_regression_missing_train, y_regression_missing_train = make_example_3d_numpy( n_cases=10, n_channels=1, n_timepoints=20, random_state=data_rng.randint(np.iinfo(np.int32).max), regression_target=True, ) -X_classification_missing_test, y_classification_missing_test = make_example_3d_numpy( +X_regression_missing_test, y_regression_missing_test = make_example_3d_numpy( n_cases=5, n_channels=1, n_timepoints=20, random_state=data_rng.randint(np.iinfo(np.int32).max), regression_target=True, ) -X_classification_missing_train[:, :, data_rng.choice(20, 2)] = np.nan -X_classification_missing_test[:, :, data_rng.choice(20, 2)] = np.nan +X_regression_missing_train[:, :, data_rng.choice(20, 2)] = np.nan +X_regression_missing_test[:, :, data_rng.choice(20, 2)] = np.nan MISSING_VALUES_REGRESSION = { "numpy3D": { - "train": (X_classification_missing_train, y_classification_missing_train), - "test": (X_classification_missing_test, y_classification_missing_test), - } + "train": (X_regression_missing_train, y_regression_missing_train), + "test": (X_regression_missing_test, y_regression_missing_test), + }, + "np-list": { + "train": ( + convert_collection(X_regression_missing_train, "np-list"), + y_regression_missing_train, + ), + "test": ( + convert_collection(X_regression_missing_test, "np-list"), + y_regression_missing_test, + ), + }, } +# Series testing data + X_series = make_example_1d_numpy( n_timepoints=40, random_state=data_rng.randint(np.iinfo(np.int32).max) ) X_series2 = X_series[20:40] X_series = X_series[:20] -UNIVARIATE_SERIES_NOLABEL = {"train": (X_series, None), "test": (X_series2, None)} +UNIVARIATE_SERIES_NONE = {"train": (X_series, None), "test": (X_series2, None)} X_series_mv = make_example_2d_numpy_series( n_timepoints=40, @@ -632,7 +656,7 @@ ) X_series_mv2 = X_series_mv[:, 20:40] X_series_mv = X_series_mv[:, :20] -MULTIVARIATE_SERIES_NOLABEL = { +MULTIVARIATE_SERIES_NONE = { "train": (X_series_mv, None), "test": (X_series_mv2, None), } @@ -644,7 +668,12 @@ X_series_mi2[data_rng.choice(20, 1)] = np.nan X_series_mi = X_series_mi[:20] X_series_mi[data_rng.choice(20, 2)] = np.nan -MISSING_VALUES_NOLABEL = {"train": (X_series_mi, None), "test": (X_series_mi2, None)} +MISSING_VALUES_SERIES_NONE = { + "train": (X_series_mi, None), + "test": (X_series_mi2, None), +} + +# All testing data FULL_TEST_DATA_DICT = {} # Collection @@ -705,10 +734,11 @@ FULL_TEST_DATA_DICT.update( {f"MissingValues-Regression-{k}": v for k, v in MISSING_VALUES_REGRESSION.items()} ) + # Series -FULL_TEST_DATA_DICT.update({"UnivariateSeries-NoLabel": UNIVARIATE_SERIES_NOLABEL}) -FULL_TEST_DATA_DICT.update({"MultivariateSeries-NoLabel": MULTIVARIATE_SERIES_NOLABEL}) -FULL_TEST_DATA_DICT.update({"MissingValues-NoLabel": MISSING_VALUES_NOLABEL}) +FULL_TEST_DATA_DICT.update({"UnivariateSeries-None": UNIVARIATE_SERIES_NONE}) +FULL_TEST_DATA_DICT.update({"MultivariateSeries-None": MULTIVARIATE_SERIES_NONE}) +FULL_TEST_DATA_DICT.update({"MissingValues-None": MISSING_VALUES_SERIES_NONE}) def _get_datatypes_for_estimator(estimator): @@ -729,7 +759,7 @@ def _get_datatypes_for_estimator(estimator): univariate, multivariate, unequal_length, missing_values = ( _get_capabilities_for_estimator(estimator) ) - label_type = _get_label_type_for_estimator(estimator) + task = _get_task_for_estimator(estimator) inner_types = estimator.get_tag("X_inner_type") if not isinstance(inner_types, list): @@ -738,34 +768,34 @@ def _get_datatypes_for_estimator(estimator): if isinstance(estimator, BaseCollectionEstimator): for inner_type in inner_types: if univariate: - s = f"EqualLengthUnivariate-{label_type}-{inner_type}" + s = f"EqualLengthUnivariate-{task}-{inner_type}" if s in FULL_TEST_DATA_DICT: datatypes.append(s) if unequal_length: - s = f"UnequalLengthUnivariate-{label_type}-{inner_type}" + s = f"UnequalLengthUnivariate-{task}-{inner_type}" if s in FULL_TEST_DATA_DICT: datatypes.append(s) if multivariate: - s = f"EqualLengthMultivariate-{label_type}-{inner_type}" + s = f"EqualLengthMultivariate-{task}-{inner_type}" if s in FULL_TEST_DATA_DICT: datatypes.append(s) if unequal_length: - s = f"UnequalLengthMultivariate-{label_type}-{inner_type}" + s = f"UnequalLengthMultivariate-{task}-{inner_type}" if s in FULL_TEST_DATA_DICT: datatypes.append(s) if missing_values: - datatypes.append(f"MissingValues-{label_type}-numpy3D") + datatypes.append(f"MissingValues-{task}-numpy3D") elif isinstance(estimator, BaseSeriesEstimator): if univariate: - datatypes.append("UnivariateSeries-NoLabel") + datatypes.append(f"UnivariateSeries-{task}") if multivariate: - datatypes.append("MultivariateSeries-NoLabel") + datatypes.append(f"MultivariateSeries-{task}") if missing_values: - datatypes.append("MissingValues-NoLabel") + datatypes.append(f"MissingValues-{task}") else: raise ValueError(f"Unknown estimator type: {type(estimator)}") @@ -803,27 +833,20 @@ def _get_capabilities_for_estimator(estimator): return univariate, multivariate, unequal_length, missing_values -def _get_label_type_for_estimator(estimator): - """Get label type for estimator. +def _get_task_for_estimator(estimator): + """Get task string used to select the correct test data for the estimator. Parameters ---------- estimator : BaseAeonEstimator instance or class - Estimator instance or class to check for valid input data types. + Estimator instance or class to find the task string for. Returns ------- - label_type : str - Label type key for the estimator for use in FULL_TEST_DATA_DICT. Indicates - whether estimator can take labels for training data, and if so what kind. - Classification indicates the estimator takes a discrete target variable, - modelled by an np.array of integers. - Regression indicates the estimator takes a continuous target variable, - modelled by an np.array of floats. - NoLabel indicates the estimator does not take a target variable, or that no - estimator is yet implemented to take labels. - + data_label : str + Task string for the estimator used in forming a key from FULL_TEST_DATA_DICT. """ + # collection data with class labels if ( isinstance(estimator, BaseClassifier) or isinstance(estimator, BaseEarlyClassifier) @@ -831,17 +854,19 @@ def _get_label_type_for_estimator(estimator): or isinstance(estimator, BaseCollectionTransformer) or isinstance(estimator, BaseSimilaritySearch) ): - label_type = "Classification" + data_label = "Classification" + # collection data with continuous target labels elif isinstance(estimator, BaseRegressor): - label_type = "Regression" + data_label = "Regression" + # series data with no secondary input elif ( isinstance(estimator, BaseAnomalyDetector) or isinstance(estimator, BaseSegmenter) or isinstance(estimator, BaseSeriesTransformer) or isinstance(estimator, BaseForecaster) ): - label_type = "NoLabel" + data_label = "None" else: raise ValueError(f"Unknown estimator type: {type(estimator)}") - return label_type + return data_label diff --git a/aeon/testing/tests/test_testing_data.py b/aeon/testing/tests/test_testing_data.py index 97e1ad8572..f7b918a4b7 100644 --- a/aeon/testing/tests/test_testing_data.py +++ b/aeon/testing/tests/test_testing_data.py @@ -298,3 +298,6 @@ def test_missing_values_collection(): assert np.issubdtype( MISSING_VALUES_REGRESSION[key]["test"][1].dtype, np.integer ) or np.issubdtype(MISSING_VALUES_REGRESSION[key]["test"][1].dtype, np.floating) + + +# todo series testing data diff --git a/aeon/testing/utils/estimator_checks.py b/aeon/testing/utils/estimator_checks.py index 1c9e8f8cb3..cd936cb40a 100644 --- a/aeon/testing/utils/estimator_checks.py +++ b/aeon/testing/utils/estimator_checks.py @@ -18,7 +18,12 @@ def _run_estimator_method(estimator, method_name, datatype, split): method = getattr(estimator, method_name) args = inspect.getfullargspec(method)[0] try: - if "X" in args and "y" in args: + if "y" in args and "exog" in args: + return method( + y=FULL_TEST_DATA_DICT[datatype][split][0], + exog=FULL_TEST_DATA_DICT[datatype][split][1], + ) + elif "X" in args and "y" in args: return method( X=FULL_TEST_DATA_DICT[datatype][split][0], y=FULL_TEST_DATA_DICT[datatype][split][1], From 54743a0eab82e46d38c77be6a589524b05d98d35 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Tue, 5 Nov 2024 19:10:16 +0000 Subject: [PATCH 18/31] _is_fitted -> is_fitted --- aeon/forecasting/base.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index 48228a8226..fcb9728ba8 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -63,7 +63,7 @@ def fit(self, y, exog=None): if exog is not None: raise NotImplementedError("Exogenous variables not yet supported") # Validate exog - self._is_fitted = True + self.is_fitted = True return self._fit(y, exog) @abstractmethod From b9bba0fb673f49d33223d127c6a5428a7ae714b0 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Tue, 5 Nov 2024 19:11:16 +0000 Subject: [PATCH 19/31] _is_fitted -> is_fitted --- aeon/forecasting/base.py | 2 +- examples/forecasting/forecasting.ipynb | 16 +++++++--------- 2 files changed, 8 insertions(+), 10 deletions(-) diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index fcb9728ba8..b7cd791613 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -25,7 +25,7 @@ class BaseForecaster(BaseSeriesEstimator): ---------- horizon : int, default =1 The number of time steps ahead to forecast. If horizon is one, the forecaster - will learn to predict one point ahead + will learn to predict one point ahead. """ _tags = { diff --git a/examples/forecasting/forecasting.ipynb b/examples/forecasting/forecasting.ipynb index d53ec82f78..402584d09e 100644 --- a/examples/forecasting/forecasting.ipynb +++ b/examples/forecasting/forecasting.ipynb @@ -7,18 +7,16 @@ "\n", "This notebook describes the new, experimental, forecasting module in aeon. We have\n", "recently removed a lot of legacy code that was almost entirely wrappers around other\n", - "projects, mostly statsmodels. Most of\n", - "the contributors to aeon are from a computer science/machine learning background\n", - "rather than stats and forecasting, and our objectives for forecasting have changed to\n", - " reflect this. Our focus is on:\n", + "projects, mostly statsmodels. Most of the contributors to aeon are from a computer\n", + "science/machine learning background rather than stats and forecasting, and our\n", + "objectives for forecasting have changed to reflect this. Our focus is on:\n", "\n", "1. not attempting to be a comprehensive forecasting package.\n", "\n", - "\n", - "We are not trying to do all things in forecasting. There are many packages that do\n", - "that already. Our previous experience taught us that trying to be all things to all\n", - "people resulted in over engineered spaghetti that was hard to maintain, understand or\n", - " use.\n", + "Forecasting is a wide field with lots of specific variants and use cases. The open\n", + "source landscape is crowded with packages that focus primarily or exclusively on\n", + "forecasting. We are not trying to do all things in forecasting. We want to focus on a\n", + " few key use cases that reflect our research interests.\n", "\n", "2. fast forecasting with numpy arrays.\n", "\n", From 9f8bed64bc1bc5a1bf1531e4868997b8c632835c Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Tue, 5 Nov 2024 20:19:03 +0000 Subject: [PATCH 20/31] _forecast --- aeon/forecasting/base.py | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index b7cd791613..023e9d6db7 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -109,7 +109,13 @@ def forecast(self, y, X=None): single prediction directly after the last point in X. """ y = self._preprocess_series(y, axis=self.axis, store_metadata=False) + self.is_fitted = True return self._forecast(y, X) - @abstractmethod - def _forecast(self, y=None, exog=None): ... + def _forecast(self, y=None, exog=None): + """Forecast values for time series X. + + NOTE: deal with horizons + """ + self.fit(y, exog) + return self.predict(y, exog) From 63d22ff955e5cf518fd76348ac0869e8a72c9708 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Tue, 5 Nov 2024 20:32:51 +0000 Subject: [PATCH 21/31] notebook --- examples/forecasting/forecasting.ipynb | 76 +++++++++++--------------- 1 file changed, 31 insertions(+), 45 deletions(-) diff --git a/examples/forecasting/forecasting.ipynb b/examples/forecasting/forecasting.ipynb index 402584d09e..993bf0e5a4 100644 --- a/examples/forecasting/forecasting.ipynb +++ b/examples/forecasting/forecasting.ipynb @@ -61,7 +61,7 @@ }, { "cell_type": "code", - "execution_count": 13, + "execution_count": 1, "outputs": [ { "name": "stdout", @@ -104,7 +104,7 @@ }, { "cell_type": "code", - "execution_count": 14, + "execution_count": 2, "outputs": [ { "name": "stdout", @@ -112,7 +112,7 @@ "text": [ " Possible data structures for input to forecaster ['pd.Series', 'pd.DataFrame', 'np.ndarray']\n", "\n", - " Tags for BaseForecaster: {'python_version': None, 'python_dependencies': None, 'cant_pickle': False, 'non_deterministic': False, 'algorithm_type': None, 'capability:missing_values': False, 'capability:multithreading': False, 'capability:univariate': True, 'capability:multivariate': False, 'X_inner_type': 'np.ndarray', 'y_inner_type': 'np.ndarray', 'fit_is_empty': False, 'requires_y': True}\n" + " Tags for BaseForecaster: {'python_version': None, 'python_dependencies': None, 'cant_pickle': False, 'non_deterministic': False, 'algorithm_type': None, 'capability:missing_values': False, 'capability:multithreading': False, 'capability:univariate': True, 'capability:multivariate': False, 'X_inner_type': 'np.ndarray', 'fit_is_empty': False}\n" ] } ], @@ -138,13 +138,13 @@ }, { "cell_type": "code", - "execution_count": 15, + "execution_count": 3, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ - "\n" + "\n" ] } ], @@ -155,7 +155,7 @@ "\n", "y = load_airline()\n", "print(type(y))\n", - "y2 = y.to_numpy()\n", + "y2 = pd.Series(y)\n", "y3 = pd.DataFrame(y)" ], "metadata": { @@ -179,14 +179,18 @@ }, { "cell_type": "code", - "execution_count": 16, + "execution_count": 4, "outputs": [ { - "name": "stdout", - "output_type": "stream", - "text": [ - "np.ndarray\n", - "432.0\n" + "ename": "ValueError", + "evalue": "Tag with name y_inner_type could not be found.", + "output_type": "error", + "traceback": [ + "\u001B[1;31m---------------------------------------------------------------------------\u001B[0m", + "\u001B[1;31mValueError\u001B[0m Traceback (most recent call last)", + "Cell \u001B[1;32mIn[4], line 5\u001B[0m\n\u001B[0;32m 2\u001B[0m \u001B[38;5;28;01mfrom\u001B[39;00m \u001B[38;5;21;01maeon\u001B[39;00m\u001B[38;5;21;01m.\u001B[39;00m\u001B[38;5;21;01mforecasting\u001B[39;00m \u001B[38;5;28;01mimport\u001B[39;00m DummyForecaster\n\u001B[0;32m 4\u001B[0m d \u001B[38;5;241m=\u001B[39m DummyForecaster()\n\u001B[1;32m----> 5\u001B[0m \u001B[38;5;28mprint\u001B[39m(\u001B[43md\u001B[49m\u001B[38;5;241;43m.\u001B[39;49m\u001B[43mget_tag\u001B[49m\u001B[43m(\u001B[49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[38;5;124;43my_inner_type\u001B[39;49m\u001B[38;5;124;43m\"\u001B[39;49m\u001B[43m)\u001B[49m)\n\u001B[0;32m 6\u001B[0m d\u001B[38;5;241m.\u001B[39mfit(y)\n\u001B[0;32m 7\u001B[0m p \u001B[38;5;241m=\u001B[39m d\u001B[38;5;241m.\u001B[39mpredict()\n", + "File \u001B[1;32mC:\\Code\\aeon\\aeon\\base\\_base.py:269\u001B[0m, in \u001B[0;36mBaseAeonEstimator.get_tag\u001B[1;34m(self, tag_name, raise_error, tag_value_default)\u001B[0m\n\u001B[0;32m 266\u001B[0m tag_value \u001B[38;5;241m=\u001B[39m collected_tags\u001B[38;5;241m.\u001B[39mget(tag_name, tag_value_default)\n\u001B[0;32m 268\u001B[0m \u001B[38;5;28;01mif\u001B[39;00m raise_error \u001B[38;5;129;01mand\u001B[39;00m tag_name \u001B[38;5;129;01mnot\u001B[39;00m \u001B[38;5;129;01min\u001B[39;00m collected_tags\u001B[38;5;241m.\u001B[39mkeys():\n\u001B[1;32m--> 269\u001B[0m \u001B[38;5;28;01mraise\u001B[39;00m \u001B[38;5;167;01mValueError\u001B[39;00m(\u001B[38;5;124mf\u001B[39m\u001B[38;5;124m\"\u001B[39m\u001B[38;5;124mTag with name \u001B[39m\u001B[38;5;132;01m{\u001B[39;00mtag_name\u001B[38;5;132;01m}\u001B[39;00m\u001B[38;5;124m could not be found.\u001B[39m\u001B[38;5;124m\"\u001B[39m)\n\u001B[0;32m 271\u001B[0m \u001B[38;5;28;01mreturn\u001B[39;00m tag_value\n", + "\u001B[1;31mValueError\u001B[0m: Tag with name y_inner_type could not be found." ] } ], @@ -206,16 +210,8 @@ }, { "cell_type": "code", - "execution_count": 17, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "432.0\n" - ] - } - ], + "execution_count": null, + "outputs": [], "source": [ "# forecast is equivalent to fit_predict in other estimators\n", "p2 = d.forecast(y)\n", @@ -242,18 +238,8 @@ }, { "cell_type": "code", - "execution_count": 19, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[451.67541971]\n", - "[505.91690237]\n", - "[505.91690237]\n" - ] - } - ], + "execution_count": null, + "outputs": [], "source": [ "from aeon.forecasting import RegressionForecaster\n", "\n", @@ -285,17 +271,8 @@ }, { "cell_type": "code", - "execution_count": 20, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[5.38765716] ,\n", - " [17.05612701]\n" - ] - } - ], + "execution_count": null, + "outputs": [], "source": [ "import numpy as np\n", "\n", @@ -317,6 +294,15 @@ "metadata": { "collapsed": false } + }, + { + "cell_type": "code", + "execution_count": null, + "outputs": [], + "source": [], + "metadata": { + "collapsed": false + } } ], "metadata": { From 107cab73736bba7ea6ebfbb11f03a8c78ba48002 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Tue, 5 Nov 2024 20:34:11 +0000 Subject: [PATCH 22/31] is_fitted --- aeon/forecasting/base.py | 3 --- 1 file changed, 3 deletions(-) diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index 023e9d6db7..884119abad 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -91,8 +91,6 @@ def predict(self, y=None, exog=None): raise ValueError("Forecaster must be fitted before predicting") if exog is not None: raise NotImplementedError("Exogenous variables not yet supported") - # Validate exog - self.is_fitted = True return self._predict(y, exog) @abstractmethod @@ -109,7 +107,6 @@ def forecast(self, y, X=None): single prediction directly after the last point in X. """ y = self._preprocess_series(y, axis=self.axis, store_metadata=False) - self.is_fitted = True return self._forecast(y, X) def _forecast(self, y=None, exog=None): From 4474593995ac8cc8c0052b2e0af2163778fbdbba Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Wed, 6 Nov 2024 08:48:38 +0000 Subject: [PATCH 23/31] y_fitted --- aeon/forecasting/base.py | 1 + 1 file changed, 1 insertion(+) diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index 884119abad..b1b54f8b36 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -33,6 +33,7 @@ class BaseForecaster(BaseSeriesEstimator): "capability:multivariate": False, "capability:missing_values": False, "fit_is_empty": False, + "y_inner_type": "np.ndarray", } def __init__(self, horizon=1, axis=1): From 2d28b3ac2053b5c25db2dd72560e096dd21f3b70 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Wed, 6 Nov 2024 10:31:51 +0000 Subject: [PATCH 24/31] ETS forecaster --- aeon/forecasting/__init__.py | 3 +- aeon/forecasting/_ets.py | 83 ++++++++++++++++++++++++++++++ aeon/forecasting/tests/test_ets.py | 16 ++++++ 3 files changed, 101 insertions(+), 1 deletion(-) create mode 100644 aeon/forecasting/_ets.py create mode 100644 aeon/forecasting/tests/test_ets.py diff --git a/aeon/forecasting/__init__.py b/aeon/forecasting/__init__.py index e114736cbc..c6ca4299ca 100644 --- a/aeon/forecasting/__init__.py +++ b/aeon/forecasting/__init__.py @@ -1,7 +1,8 @@ """Forecasters.""" -__all__ = ["DummyForecaster", "BaseForecaster", "RegressionForecaster"] +__all__ = ["DummyForecaster", "BaseForecaster", "RegressionForecaster", "ETSForecaster"] from aeon.forecasting._dummy import DummyForecaster +from aeon.forecasting._ets import ETSForecaster from aeon.forecasting._regression import RegressionForecaster from aeon.forecasting.base import BaseForecaster diff --git a/aeon/forecasting/_ets.py b/aeon/forecasting/_ets.py new file mode 100644 index 0000000000..57eb9aea4d --- /dev/null +++ b/aeon/forecasting/_ets.py @@ -0,0 +1,83 @@ +import numpy as np + +from aeon.forecasting.base import BaseForecaster + + +class ETSForecaster(BaseForecaster): + """Exponential Smoothing forecaster. + + Simple first implementation with Holt-Winters method + and no seasonality. + + Parameters + ---------- + alpha : float, default = 0.2 + Level smoothing parameter. + beta : float, default = 0.2 + Trend smoothing parameter. + gamma : float, default = 0.2 + Seasonal smoothing parameter. + season_length : int, default = 1 + The length of the seasonality period. + """ + + def __init__(self, alpha=0.2, beta=0.2, gamma=0.2, season_length=1, horizon=1): + self.alpha = alpha + self.beta = beta + self.gamma = gamma + self.season_length = season_length + self.forecast_val_ = 0.0 + self.level_ = 0.0 + self.trend_ = 0.0 + self.seasonals_ = None + super().__init__(horizon=horizon, axis=1) + + def _fit(self, y, exog=None): + """Fit Exponential Smoothing forecaster to series y. + + Fit a forecaster to predict self.horizon steps ahead using y. + + Parameters + ---------- + y : np.ndarray + A time series on which to learn a forecaster to predict horizon ahead + exog : np.ndarray, default =None + Optional exogenous time series data assumed to be aligned with y + + Returns + ------- + self + Fitted BaseForecaster. + """ + data = y.squeeze() + self.n_timepoints = len(data) + sl = self.season_length + # Initialize components + self.level_ = data[0] + self.trend_ = np.mean(data[sl : 2 * sl]) - np.mean(data[:sl]) + self.seasonals_ = [data[i] / data[0] for i in range(sl)] + for t in range(sl, self.n_timepoints): + # Calculate level, trend, and seasonal components + level_prev = self.level_ + self.level_ = self.alpha * ( + data[t] / self.seasonals_[t % self.season_length] + ) + (1 - self.alpha) * (self.level_ + self.trend_) + self.trend_ = ( + self.beta * (self.level_ - level_prev) + (1 - self.beta) * self.trend_ + ) + self.seasonals_[t % self.season_length] = ( + self.gamma * (data[t] / self.level_) + + (1 - self.gamma) * self.seasonals_[t % sl] + ) + return self + + def _predict(self, y=None, exog=None): + # Generate forecasts based on the final values of level, trend, and seasonals + forecast = (self.level_ + (self.horizon + 1) * self.trend_) * self.seasonals_[ + (self.n_timepoints + self.horizon) % self.season_length + ] + return forecast + + def _forecast(self, y): + self.fit(y) + return self.predict() diff --git a/aeon/forecasting/tests/test_ets.py b/aeon/forecasting/tests/test_ets.py new file mode 100644 index 0000000000..5ee7d74329 --- /dev/null +++ b/aeon/forecasting/tests/test_ets.py @@ -0,0 +1,16 @@ +"""Test ETS.""" + +import numpy as np + +from aeon.forecasting import ETSForecaster + + +def test_ets_forecaster(): + """TestETSForecaster.""" + data = np.array( + [3, 10, 12, 13, 12, 10, 12, 3, 10, 12, 13, 12, 10, 12] + ) # Sample seasonal data + forecaster = ETSForecaster(alpha=0.5, beta=0.3, gamma=0.4, season_length=4) + forecaster.fit(data) + p = forecaster.predict() + assert np.isclose(p, 15.85174501127) From 87c5969aa820c549b6fc5cdd63d2f2a64dcc5564 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Wed, 6 Nov 2024 11:05:41 +0000 Subject: [PATCH 25/31] add y checks and conversion --- aeon/forecasting/base.py | 50 +++++++++++++++++++++++++++++++++++----- 1 file changed, 44 insertions(+), 6 deletions(-) diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index b1b54f8b36..9b4aef205d 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -10,7 +10,11 @@ from abc import abstractmethod +import numpy as np +import pandas as pd + from aeon.base import BaseSeriesEstimator +from aeon.base._base_series import VALID_INPUT_TYPES class BaseForecaster(BaseSeriesEstimator): @@ -38,6 +42,7 @@ class BaseForecaster(BaseSeriesEstimator): def __init__(self, horizon=1, axis=1): self.horizon = horizon + self.meta_ = None # Meta data related to y on the last fit super().__init__(axis) def fit(self, y, exog=None): @@ -60,7 +65,8 @@ def fit(self, y, exog=None): # Validate y # Convert if necessary - y = self._preprocess_series(y, axis=self.axis, store_metadata=False) + self._check_X(y, self.axis) + y = self._convert_y(y, self.axis) if exog is not None: raise NotImplementedError("Exogenous variables not yet supported") # Validate exog @@ -87,7 +93,8 @@ def predict(self, y=None, exog=None): single prediction self.horizon steps ahead of y. """ if y is not None: - y = self._preprocess_series(y, axis=self.axis, store_metadata=False) + self._check_X(y, self.axis) + y = self._convert_y(y, self.axis) if not self.is_fitted: raise ValueError("Forecaster must be fitted before predicting") if exog is not None: @@ -107,13 +114,44 @@ def forecast(self, y, X=None): np.ndarray single prediction directly after the last point in X. """ + self._check_X(y, self.axis) + y = self._convert_y(y, self.axis) y = self._preprocess_series(y, axis=self.axis, store_metadata=False) return self._forecast(y, X) def _forecast(self, y=None, exog=None): - """Forecast values for time series X. - - NOTE: deal with horizons - """ + """Forecast values for time series X.""" self.fit(y, exog) return self.predict(y, exog) + + def _convert_y(self, y: VALID_INPUT_TYPES, axis: int): + """Convert y to self.get_tag("y_inner_type").""" + if axis > 1 or axis < 0: + raise ValueError(f"Input axis should be 0 or 1, saw {axis}") + + inner_type = self.get_tag("y_inner_type") + if not isinstance(inner_type, list): + inner_type = [inner_type] + inner_names = [i.split(".")[-1] for i in inner_type] + + input = type(y).__name__ + if input not in inner_names: + if inner_names[0] == "ndarray": + y = y.to_numpy() + elif inner_names[0] == "DataFrame": + # converting a 1d array will create a 2d array in axis 0 format + transpose = False + if y.ndim == 1 and axis == 1: + transpose = True + y = pd.DataFrame(y) + if transpose: + y = y.T + else: + raise ValueError( + f"Unsupported inner type {inner_names[0]} derived from {inner_type}" + ) + if y.ndim > 1 and self.axis != axis: + y = y.T + elif y.ndim == 1 and isinstance(y, np.ndarray): + y = y[np.newaxis, :] if self.axis == 1 else y[:, np.newaxis] + return y From afde1f58a92c4bdf8c255d098ed2a33c9c72e125 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Wed, 6 Nov 2024 21:39:58 +0000 Subject: [PATCH 26/31] add tag --- aeon/utils/tags/_tags.py | 8 ++++++++ 1 file changed, 8 insertions(+) diff --git a/aeon/utils/tags/_tags.py b/aeon/utils/tags/_tags.py index 650a7c4dc4..554584115e 100644 --- a/aeon/utils/tags/_tags.py +++ b/aeon/utils/tags/_tags.py @@ -53,6 +53,14 @@ class : identifier for the base class of objects this tag applies to "description": "What data structure(s) the estimator uses internally for " "fit/predict.", }, + "y_inner_type": { + "class": "forecaster", + "type": [ + ("list||str", SERIES_DATA_TYPES), + ], + "description": "What data structure(s) the estimator uses internally for " + "fit/predict.", + }, "algorithm_type": { "class": "estimator", "type": [ From 5bf2828284bf7f8c80ca71907c71241e78ed2782 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Thu, 7 Nov 2024 15:04:42 +0000 Subject: [PATCH 27/31] tidy --- aeon/forecasting/_ets.py | 34 ++++++++++++++++------------------ aeon/forecasting/base.py | 7 ++++--- 2 files changed, 20 insertions(+), 21 deletions(-) diff --git a/aeon/forecasting/_ets.py b/aeon/forecasting/_ets.py index 57eb9aea4d..9f20bf9244 100644 --- a/aeon/forecasting/_ets.py +++ b/aeon/forecasting/_ets.py @@ -21,15 +21,15 @@ class ETSForecaster(BaseForecaster): The length of the seasonality period. """ - def __init__(self, alpha=0.2, beta=0.2, gamma=0.2, season_length=1, horizon=1): + def __init__(self, alpha=0.2, beta=0.2, gamma=0.2, season_len=1, horizon=1): self.alpha = alpha self.beta = beta self.gamma = gamma - self.season_length = season_length + self.season_len = season_len self.forecast_val_ = 0.0 self.level_ = 0.0 self.trend_ = 0.0 - self.seasonals_ = None + self.season_ = None super().__init__(horizon=horizon, axis=1) def _fit(self, y, exog=None): @@ -51,31 +51,29 @@ def _fit(self, y, exog=None): """ data = y.squeeze() self.n_timepoints = len(data) - sl = self.season_length + sl = self.season_len # Initialize components self.level_ = data[0] self.trend_ = np.mean(data[sl : 2 * sl]) - np.mean(data[:sl]) - self.seasonals_ = [data[i] / data[0] for i in range(sl)] + self.season_ = [data[i] / data[0] for i in range(sl)] for t in range(sl, self.n_timepoints): # Calculate level, trend, and seasonal components level_prev = self.level_ - self.level_ = self.alpha * ( - data[t] / self.seasonals_[t % self.season_length] - ) + (1 - self.alpha) * (self.level_ + self.trend_) - self.trend_ = ( - self.beta * (self.level_ - level_prev) + (1 - self.beta) * self.trend_ - ) - self.seasonals_[t % self.season_length] = ( - self.gamma * (data[t] / self.level_) - + (1 - self.gamma) * self.seasonals_[t % sl] - ) + l1 = data[t] / self.season_[t % self.season_len] + l2 = self.level_ + self.trend_ + self.level_ = self.alpha * l1 + (1 - self.alpha) * l2 + trend = self.level_ - level_prev + self.trend_ = self.beta * trend + (1 - self.beta) * self.trend_ + s1 = data[t] / self.level_ + s2 = self.season_[t % sl] + self.season_[t % self.season_len] = self.gamma * s1 + (1 - self.gamma) * s2 return self def _predict(self, y=None, exog=None): # Generate forecasts based on the final values of level, trend, and seasonals - forecast = (self.level_ + (self.horizon + 1) * self.trend_) * self.seasonals_[ - (self.n_timepoints + self.horizon) % self.season_length - ] + trend = (self.horizon + 1) * self.trend_ + seasonal = self.season_[(self.n_timepoints + self.horizon) % self.season_len] + forecast = (self.level_ + trend) * seasonal return forecast def _forecast(self, y): diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index 9b4aef205d..aff6939b2f 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -104,7 +104,7 @@ def predict(self, y=None, exog=None): @abstractmethod def _predict(self, y=None, exog=None): ... - def forecast(self, y, X=None): + def forecast(self, y, exog=None): """ Forecast basically fit_predict. @@ -116,8 +116,9 @@ def forecast(self, y, X=None): """ self._check_X(y, self.axis) y = self._convert_y(y, self.axis) - y = self._preprocess_series(y, axis=self.axis, store_metadata=False) - return self._forecast(y, X) + self._check_X(y, self.axis) + y = self._convert_y(y, self.axis) + return self._forecast(y, exog) def _forecast(self, y=None, exog=None): """Forecast values for time series X.""" From c292f5dfb1d578d94a456bd8080f85be87484d76 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Thu, 7 Nov 2024 18:15:46 +0000 Subject: [PATCH 28/31] _check_is_fitted() --- aeon/forecasting/base.py | 3 +-- 1 file changed, 1 insertion(+), 2 deletions(-) diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index aff6939b2f..814a7531a9 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -92,11 +92,10 @@ def predict(self, y=None, exog=None): float single prediction self.horizon steps ahead of y. """ + self._check_is_fitted() if y is not None: self._check_X(y, self.axis) y = self._convert_y(y, self.axis) - if not self.is_fitted: - raise ValueError("Forecaster must be fitted before predicting") if exog is not None: raise NotImplementedError("Exogenous variables not yet supported") return self._predict(y, exog) From 608558ddf3c714e7c4784b37135eb31d1521b3a5 Mon Sep 17 00:00:00 2001 From: Tony Bagnall Date: Thu, 7 Nov 2024 19:08:01 +0000 Subject: [PATCH 29/31] _check_is_fitted() --- aeon/forecasting/_ets.py | 2 +- aeon/forecasting/tests/test_ets.py | 2 +- 2 files changed, 2 insertions(+), 2 deletions(-) diff --git a/aeon/forecasting/_ets.py b/aeon/forecasting/_ets.py index 9f20bf9244..0b4172ca3c 100644 --- a/aeon/forecasting/_ets.py +++ b/aeon/forecasting/_ets.py @@ -17,7 +17,7 @@ class ETSForecaster(BaseForecaster): Trend smoothing parameter. gamma : float, default = 0.2 Seasonal smoothing parameter. - season_length : int, default = 1 + season_len : int, default = 1 The length of the seasonality period. """ diff --git a/aeon/forecasting/tests/test_ets.py b/aeon/forecasting/tests/test_ets.py index 5ee7d74329..ba92d76cb1 100644 --- a/aeon/forecasting/tests/test_ets.py +++ b/aeon/forecasting/tests/test_ets.py @@ -10,7 +10,7 @@ def test_ets_forecaster(): data = np.array( [3, 10, 12, 13, 12, 10, 12, 3, 10, 12, 13, 12, 10, 12] ) # Sample seasonal data - forecaster = ETSForecaster(alpha=0.5, beta=0.3, gamma=0.4, season_length=4) + forecaster = ETSForecaster(alpha=0.5, beta=0.3, gamma=0.4, season_len=4) forecaster.fit(data) p = forecaster.predict() assert np.isclose(p, 15.85174501127) From 5f0a3f55b777d503ce7fe8cba0a2db80375a1cc9 Mon Sep 17 00:00:00 2001 From: alexbanwell1 <31886108+alexbanwell1@users.noreply.github.com> Date: Fri, 8 Nov 2024 09:23:10 +0000 Subject: [PATCH 30/31] Add fully functional ETS Forecaster. Modify base to not set default y in forecast. Update tests for ETS Forecaster. Add script to verify ETS Forecaster against statsforecast module using a large number of random parameter inputs. (#2318) Co-authored-by: Alex Banwell --- aeon/forecasting/__init__.py | 10 +- aeon/forecasting/_ets.py | 338 ++++++++++++++++++++++++++--- aeon/forecasting/_verify_ets.py | 168 ++++++++++++++ aeon/forecasting/base.py | 2 +- aeon/forecasting/tests/test_ets.py | 72 +++++- 5 files changed, 551 insertions(+), 39 deletions(-) create mode 100644 aeon/forecasting/_verify_ets.py diff --git a/aeon/forecasting/__init__.py b/aeon/forecasting/__init__.py index c6ca4299ca..45eaa41715 100644 --- a/aeon/forecasting/__init__.py +++ b/aeon/forecasting/__init__.py @@ -1,8 +1,14 @@ """Forecasters.""" -__all__ = ["DummyForecaster", "BaseForecaster", "RegressionForecaster", "ETSForecaster"] +__all__ = [ + "DummyForecaster", + "BaseForecaster", + "RegressionForecaster", + "ETSForecaster", + "ModelType", +] from aeon.forecasting._dummy import DummyForecaster -from aeon.forecasting._ets import ETSForecaster +from aeon.forecasting._ets import ETSForecaster, ModelType from aeon.forecasting._regression import RegressionForecaster from aeon.forecasting.base import BaseForecaster diff --git a/aeon/forecasting/_ets.py b/aeon/forecasting/_ets.py index 0b4172ca3c..7c97378002 100644 --- a/aeon/forecasting/_ets.py +++ b/aeon/forecasting/_ets.py @@ -1,35 +1,129 @@ +"""ETSForecaster class. + +An implementation of the exponential smoothing statistics forecasting algorithm. +Implements additive and multiplicative error models, +None, additive and multiplicative (including damped) trend and +None, additive and mutliplicative seasonality + +aeon enhancement proposal +https://github.com/aeon-toolkit/aeon/pull/2244/ + +""" + +__maintainer__ = [] +__all__ = ["ETSForecaster", "ModelType"] + import numpy as np from aeon.forecasting.base import BaseForecaster +NONE = 0 +ADDITIVE = 1 +MULTIPLICATIVE = 2 + + +class ModelType: + """ + Class describing the error, trend and seasonality model of an ETS forecaster. + + Attributes + ---------- + error_type : int + The type of error model; either Additive(1) or Multiplicative(2) + trend_type : int + The type of trend model; one of None(0), additive(1) or multiplicative(2). + seasonality_type : int + The type of seasonality model; one of None(0), additive(1) or multiplicative(2). + seasonal_period : int + The period of the seasonality (m) (e.g., for quaterly data seasonal_period = 4). + """ + + error_type: int + trend_type: int + seasonality_type: int + seasonal_period: int + + def __init__( + self, + error_type=ADDITIVE, + trend_type=NONE, + seasonality_type=NONE, + seasonal_period=1, + ): + assert error_type != NONE, "Error must be either additive or multiplicative" + if seasonal_period < 1 or seasonality_type == NONE: + seasonal_period = 1 + self.error_type = error_type + self.trend_type = trend_type + self.seasonality_type = seasonality_type + self.seasonal_period = seasonal_period + class ETSForecaster(BaseForecaster): """Exponential Smoothing forecaster. - Simple first implementation with Holt-Winters method - and no seasonality. + An implementation of the exponential smoothing statistics forecasting algorithm. + Implements additive and multiplicative error models, + None, additive and multiplicative (including damped) trend and + None, additive and mutliplicative seasonality[1]_. Parameters ---------- - alpha : float, default = 0.2 + alpha : float, default = 0.1 Level smoothing parameter. - beta : float, default = 0.2 + beta : float, default = 0.01 Trend smoothing parameter. - gamma : float, default = 0.2 + gamma : float, default = 0.01 Seasonal smoothing parameter. - season_len : int, default = 1 - The length of the seasonality period. + phi : float, default = 0.99 + Trend damping smoothing parameters + horizon : int, default = 1 + The horizon to forecast to. + model_type : ModelType, default = ModelType() + A object of type ModelType, describing the error, + trend and seasonality type of this ETS model. + + References + ---------- + .. [1] R. J. Hyndman and G. Athanasopoulos, + Forecasting: Principles and Practice. Melbourne, Australia: OTexts, 2014. + + Examples + -------- + >>> from aeon.forecasting import ETSForecaster, ModelType + >>> from aeon.datasets import load_airline + >>> y = load_airline() + >>> forecaster = ETSForecaster(alpha=0.4, beta=0.2, gamma=0.5, phi=0.8, horizon=1, + model_type=ModelType(1,2,2,4)) + >>> forecaster.fit(y) + >>> forecaster.predict() + 366.90200486015596 """ - def __init__(self, alpha=0.2, beta=0.2, gamma=0.2, season_len=1, horizon=1): + default_model_type = ModelType() + + def __init__( + self, + alpha=0.1, + beta=0.01, + gamma=0.01, + phi=0.99, + horizon=1, + model_type=default_model_type, + ): self.alpha = alpha self.beta = beta self.gamma = gamma - self.season_len = season_len + self.phi = phi self.forecast_val_ = 0.0 self.level_ = 0.0 self.trend_ = 0.0 self.season_ = None + self.n_timepoints = 0 + self.avg_mean_sq_err_ = 0 + self.liklihood_ = 0 + self.residuals_ = [] + self.model_type = model_type super().__init__(horizon=horizon, axis=1) def _fit(self, y, exog=None): @@ -51,31 +145,211 @@ def _fit(self, y, exog=None): """ data = y.squeeze() self.n_timepoints = len(data) - sl = self.season_len - # Initialize components - self.level_ = data[0] - self.trend_ = np.mean(data[sl : 2 * sl]) - np.mean(data[:sl]) - self.season_ = [data[i] / data[0] for i in range(sl)] - for t in range(sl, self.n_timepoints): + self._initialise(data) + self.avg_mean_sq_err_ = 0 + self.liklihood_ = 0 + mul_liklihood_pt2 = 0 + self.residuals_ = np.zeros( + self.n_timepoints + ) # 1 Less residual than data points + for t, data_item in enumerate(data[self.model_type.seasonal_period :]): # Calculate level, trend, and seasonal components - level_prev = self.level_ - l1 = data[t] / self.season_[t % self.season_len] - l2 = self.level_ + self.trend_ - self.level_ = self.alpha * l1 + (1 - self.alpha) * l2 - trend = self.level_ - level_prev - self.trend_ = self.beta * trend + (1 - self.beta) * self.trend_ - s1 = data[t] / self.level_ - s2 = self.season_[t % sl] - self.season_[t % self.season_len] = self.gamma * s1 + (1 - self.gamma) * s2 + fitted_value, error = self._update_states( + data_item, t % self.model_type.seasonal_period + ) + self.residuals_[t] = error + self.avg_mean_sq_err_ += (data_item - fitted_value) ** 2 + self.liklihood_ += error * error + mul_liklihood_pt2 += np.log(np.fabs(fitted_value)) + self.avg_mean_sq_err_ /= self.n_timepoints - self.model_type.seasonal_period + self.liklihood_ = ( + self.n_timepoints - self.model_type.seasonal_period + ) * np.log(self.liklihood_) + if self.model_type.error_type == MULTIPLICATIVE: + self.liklihood_ += 2 * mul_liklihood_pt2 return self + def _update_states(self, data_item, seasonal_index): + """ + Update level, trend, and seasonality components. + + Using state space equations for an ETS model. + + Parameters + ---------- + data_item: float + The current value of the time series. + seasonal_index: int + The index to update the seasonal component. + """ + model = self.model_type + # Retrieve the current state values + level = self.level_ + trend = self.trend_ + seasonality = self.season_[seasonal_index] + fitted_value, damped_trend, trend_level_combination = self._predict_value( + trend, level, seasonality, self.phi + ) + # Calculate the error term (observed value - fitted value) + if model.error_type == MULTIPLICATIVE: + error = data_item / fitted_value - 1 # Multiplicative error + else: + error = data_item - fitted_value # Additive error + # Update level + if model.error_type == MULTIPLICATIVE: + self.level_ = trend_level_combination * (1 + self.alpha * error) + self.trend_ = damped_trend * (1 + self.beta * error) + self.season_[seasonal_index] = seasonality * (1 + self.gamma * error) + if model.seasonality_type == ADDITIVE: + self.level_ += ( + self.alpha * error * seasonality + ) # Add seasonality correction + self.season_[seasonal_index] += ( + self.gamma * error * trend_level_combination + ) + if model.trend_type == ADDITIVE: + self.trend_ += (level + seasonality) * self.beta * error + else: + self.trend_ += seasonality / level * self.beta * error + elif model.trend_type == ADDITIVE: + self.trend_ += level * self.beta * error + else: + level_correction = 1 + trend_correction = 1 + seasonality_correction = 1 + if model.seasonality_type == MULTIPLICATIVE: + # Add seasonality correction + level_correction *= seasonality + trend_correction *= seasonality + seasonality_correction *= trend_level_combination + if model.trend_type == MULTIPLICATIVE: + trend_correction *= level + self.level_ = ( + trend_level_combination + self.alpha * error / level_correction + ) + self.trend_ = damped_trend + self.beta * error / trend_correction + self.season_[seasonal_index] = ( + seasonality + self.gamma * error / seasonality_correction + ) + return (fitted_value, error) + + def _initialise(self, data): + """ + Initialize level, trend, and seasonality values for the ETS model. + + Parameters + ---------- + data : array-like + The time series data + (should contain at least two full seasons if seasonality is specified) + """ + model = self.model_type + # Initial Level: Mean of the first season + self.level_ = np.mean(data[: model.seasonal_period]) + # Initial Trend + if model.trend_type == ADDITIVE: + # Average difference between corresponding points in the first two seasons + self.trend_ = np.mean( + data[model.seasonal_period : 2 * model.seasonal_period] + - data[: model.seasonal_period] + ) + elif model.trend_type == MULTIPLICATIVE: + # Average ratio between corresponding points in the first two seasons + self.trend_ = np.mean( + data[model.seasonal_period : 2 * model.seasonal_period] + / data[: model.seasonal_period] + ) + else: + # No trend + self.trend_ = 0 + self.beta = ( + 0 # Required for the equations in _update_states to work correctly + ) + # Initial Seasonality + if model.seasonality_type == ADDITIVE: + # Seasonal component is the difference + # from the initial level for each point in the first season + self.season_ = data[: model.seasonal_period] - self.level_ + elif model.seasonality_type == MULTIPLICATIVE: + # Seasonal component is the ratio of each point in the first season + # to the initial level + self.season_ = data[: model.seasonal_period] / self.level_ + else: + # No seasonality + self.season_ = [0] + self.gamma = ( + 0 # Required for the equations in _update_states to work correctly + ) + def _predict(self, y=None, exog=None): + """ + Predict the next horizon steps ahead. + + Parameters + ---------- + y : np.ndarray, default = None + A time series to predict the next horizon value for. If None, + predict the next horizon value after series seen in fit. + exog : np.ndarray, default =None + Optional exogenous time series data assumed to be aligned with y + + Returns + ------- + float + single prediction self.horizon steps ahead of y. + """ # Generate forecasts based on the final values of level, trend, and seasonals - trend = (self.horizon + 1) * self.trend_ - seasonal = self.season_[(self.n_timepoints + self.horizon) % self.season_len] - forecast = (self.level_ + trend) * seasonal - return forecast - - def _forecast(self, y): - self.fit(y) - return self.predict() + if self.phi == 1: # No damping case + phi_h = float(self.horizon) + else: + # Geometric series formula for calculating phi + phi^2 + ... + phi^h + phi_h = self.phi * (1 - self.phi**self.horizon) / (1 - self.phi) + seasonality = self.season_[ + (self.n_timepoints + self.horizon) % self.model_type.seasonal_period + ] + fitted_value = self._predict_value( + self.trend_, self.level_, seasonality, phi_h + )[0] + return fitted_value + + def _predict_value(self, trend, level, seasonality, phi): + """ + + Generate various useful values, including the next fitted value. + + Parameters + ---------- + trend : float + The current trend value for the model + level : float + The current level value for the model + seasonality : float + The current seasonality value for the model + phi : float + The damping parameter for the model + + Returns + ------- + fitted_value : float + single prediction based on the current state variables. + damped_trend : float + The damping parameter combined with the trend dependant on the model type + trend_level_combination : float + Combination of the trend and level based on the model type. + """ + model = self.model_type + # Apply damping parameter and + # calculate commonly used combination of trend and level components + if model.trend_type == MULTIPLICATIVE: + damped_trend = trend**phi + trend_level_combination = level * damped_trend + else: # Additive trend, if no trend, then trend = 0 + damped_trend = trend * phi + trend_level_combination = level + damped_trend + + # Calculate forecast (fitted value) based on the current components + if model.seasonality_type == MULTIPLICATIVE: + fitted_value = trend_level_combination * seasonality + else: # Additive seasonality, if no seasonality, then seasonality = 0 + fitted_value = trend_level_combination + seasonality + return fitted_value, damped_trend, trend_level_combination diff --git a/aeon/forecasting/_verify_ets.py b/aeon/forecasting/_verify_ets.py new file mode 100644 index 0000000000..d93b64631f --- /dev/null +++ b/aeon/forecasting/_verify_ets.py @@ -0,0 +1,168 @@ +import random +import time + +import numpy as np +from statsforecast.ets import etscalc +from statsforecast.utils import AirPassengers as ap + +from aeon.forecasting import ETSForecaster, ModelType + +NA = -99999.0 +MAX_NMSE = 30 +MAX_SEASONAL_PERIOD = 24 + + +def setup(): + """Generate parameters required for ETS algorithms.""" + y = ap + n = len(ap) + m = random.randint(1, 24) + error = random.randint(1, 2) + trend = random.randint(0, 2) + season = random.randint(0, 2) + alpha = round(random.random(), 4) + if alpha == 0: + alpha = round(random.random(), 4) + beta = round(random.random() * alpha, 4) # 0 < beta < alpha + if beta == 0: + beta = round(random.random() * alpha, 4) + gamma = round(random.random() * (1 - alpha), 4) # 0 < beta < alpha + if gamma == 0: + gamma = round(random.random() * (1 - alpha), 4) + phi = round( + random.random() * 0.18 + 0.8, 4 + ) # Common constraint for phi is 0.8 < phi < 0.98 + e = np.zeros(n) + lik_fitets = np.zeros(1) + amse = np.zeros(MAX_NMSE) + nmse = 3 + return ( + y, + n, + m, + error, + trend, + season, + alpha, + beta, + gamma, + phi, + e, + lik_fitets, + amse, + nmse, + ) + + +def test_ets_comparison(setup_func, random_seed, catch_errors): + """Run both our statsforecast and our implementation and crosschecks.""" + random.seed(random_seed) + ( + y, + n, + m, + error, + trend, + season, + alpha, + beta, + gamma, + phi, + e, + lik_fitets, + amse, + nmse, + ) = setup_func() + # tsml-eval implementation + start = time.time() + f1 = ETSForecaster(alpha, beta, gamma, phi, 1, ModelType(error, trend, season, m)) + f1.fit(y) + end = time.time() + time_fitets = end - start + e_fitets = f1.residuals_ + amse_fitets = f1.avg_mean_sq_err_ + lik_fitets = f1.liklihood_ + # Reinitialise arrays + e.fill(0) + amse.fill(0) + f1 = ETSForecaster(alpha, beta, gamma, phi, 1, ModelType(error, trend, season, m)) + f1._initialise(y) + init_states_etscalc = np.zeros(n * (1 + (trend > 0) + m * (season > 0) + 1)) + init_states_etscalc[0] = f1.level_ + init_states_etscalc[1] = f1.trend_ + init_states_etscalc[1 + (trend != 0) : m + 1 + (trend != 0)] = f1.season_[::-1] + if season == 0: + m = 1 + # Nixtla/statsforcast implementation + start = time.time() + lik_etscalc = etscalc( + y[m:], + n - m, + init_states_etscalc, + m, + error, + trend, + season, + alpha, + beta, + gamma, + phi, + e, + amse, + nmse, + ) + end = time.time() + time_etscalc = end - start + e_etscalc = e.copy() + amse_etscalc = amse.copy()[0] + + if catch_errors: + try: + # Comparing outputs and runtime + assert np.allclose(e_fitets, e_etscalc), "Residuals Compare failed" + assert np.allclose(amse_fitets, amse_etscalc), "AMSE Compare failed" + assert np.isclose(lik_fitets, lik_etscalc), "Liklihood Compare failed" + return True + except AssertionError as e: + print(e) # noqa + print( # noqa + f"Seed: {random_seed}, Model: Error={error}, Trend={trend},\ + Seasonality={season}, seasonal period={m},\ + alpha={alpha}, beta={beta}, gamma={gamma}, phi={phi}" + ) + return False + else: + print( # noqa + f"Seed: {random_seed}, Model: Error={error}, Trend={trend},\ + Seasonality={season}, seasonal period={m}, alpha={alpha},\ + beta={beta}, gamma={gamma}, phi={phi}" + ) + diff_indices = np.where( + np.abs(e_fitets - e_etscalc) > 1e-5 * np.abs(e_etscalc) + 1e-8 + )[0] + for index in diff_indices: + print( # noqa + f"Index {index}: e_fitets = {e_fitets[index]},\ + e_etscalc = {e_etscalc[index]}" + ) + print(amse_fitets) # noqa + print(amse_etscalc) # noqa + print(lik_fitets) # noqa + print(lik_etscalc) # noqa + # assert np.allclose(init_states_fitets, init_states_etscalc) + assert np.allclose(e_fitets, e_etscalc) + assert np.allclose(amse_fitets, amse_etscalc) + assert np.isclose(lik_fitets, lik_etscalc) + print(time_fitets) # noqa + print(time_etscalc) # noqa + return True + + +if __name__ == "__main__": + # np.set_printoptions(threshold=np.inf) + # test_ets_comparison(setup, 241, False) + SUCCESSES = True + for i in range(0, 30000): + SUCCESSES &= test_ets_comparison(setup, i, True) + if SUCCESSES: + print("Test Completed Successfully with no errors") # noqa diff --git a/aeon/forecasting/base.py b/aeon/forecasting/base.py index 814a7531a9..03ac18ee18 100644 --- a/aeon/forecasting/base.py +++ b/aeon/forecasting/base.py @@ -119,7 +119,7 @@ def forecast(self, y, exog=None): y = self._convert_y(y, self.axis) return self._forecast(y, exog) - def _forecast(self, y=None, exog=None): + def _forecast(self, y, exog=None): """Forecast values for time series X.""" self.fit(y, exog) return self.predict(y, exog) diff --git a/aeon/forecasting/tests/test_ets.py b/aeon/forecasting/tests/test_ets.py index ba92d76cb1..8be3e9f799 100644 --- a/aeon/forecasting/tests/test_ets.py +++ b/aeon/forecasting/tests/test_ets.py @@ -1,16 +1,80 @@ """Test ETS.""" +__maintainer__ = [] +__all__ = [] + import numpy as np -from aeon.forecasting import ETSForecaster +from aeon.forecasting import ETSForecaster, ModelType + + +def test_ets_forecaster_additive(): + """TestETSForecaster.""" + data = np.array( + [3, 10, 12, 13, 12, 10, 12, 3, 10, 12, 13, 12, 10, 12] + ) # Sample seasonal data + forecaster = ETSForecaster( + alpha=0.5, + beta=0.3, + gamma=0.4, + phi=1, + horizon=1, + model_type=ModelType(1, 1, 1, 4), + ) + forecaster.fit(data) + p = forecaster.predict() + assert np.isclose(p, 9.191190608800001) + + +def test_ets_forecaster_mult_error(): + """TestETSForecaster.""" + data = np.array( + [3, 10, 12, 13, 12, 10, 12, 3, 10, 12, 13, 12, 10, 12] + ) # Sample seasonal data + forecaster = ETSForecaster( + alpha=0.7, + beta=0.6, + gamma=0.1, + phi=0.97, + horizon=1, + model_type=ModelType(2, 1, 1, 4), + ) + forecaster.fit(data) + p = forecaster.predict() + assert np.isclose(p, 16.20176819429869) + + +def test_ets_forecaster_mult_compnents(): + """TestETSForecaster.""" + data = np.array( + [3, 10, 12, 13, 12, 10, 12, 3, 10, 12, 13, 12, 10, 12] + ) # Sample seasonal data + forecaster = ETSForecaster( + alpha=0.4, + beta=0.2, + gamma=0.5, + phi=0.8, + horizon=1, + model_type=ModelType(1, 2, 2, 4), + ) + forecaster.fit(data) + p = forecaster.predict() + assert np.isclose(p, 12.301259229712382) -def test_ets_forecaster(): +def test_ets_forecaster_multiplicative(): """TestETSForecaster.""" data = np.array( [3, 10, 12, 13, 12, 10, 12, 3, 10, 12, 13, 12, 10, 12] ) # Sample seasonal data - forecaster = ETSForecaster(alpha=0.5, beta=0.3, gamma=0.4, season_len=4) + forecaster = ETSForecaster( + alpha=0.7, + beta=0.5, + gamma=0.2, + phi=0.85, + horizon=1, + model_type=ModelType(2, 2, 2, 4), + ) forecaster.fit(data) p = forecaster.predict() - assert np.isclose(p, 15.85174501127) + assert np.isclose(p, 16.811888294476528) From a8c00e22b581e45e7e159c235238848213a4e38e Mon Sep 17 00:00:00 2001 From: alexbanwell1 <31886108+alexbanwell1@users.noreply.github.com> Date: Fri, 15 Nov 2024 08:38:00 +0000 Subject: [PATCH 31/31] Ajb/forecasting (#2357) * Add fully functional ETS Forecaster. Modify base to not set default y in forecast. Update tests for ETS Forecaster. Add script to verify ETS Forecaster against statsforecast module using a large number of random parameter inputs. * Add faster numba version of ETS forecaster * Seperate out predict code, and add test to test without creating a class - significantly faster! * Modify _verify_ets.py to allow easy switching between statsforecast versions. This confirms that my algorithms without class overheads is significantly faster than nixtla statsforecast, and with class overheads, it is faster than their current algorithm * Add basic gradient decent optimization algorithm for smoothing parameters --------- Co-authored-by: Alex Banwell --- aeon/forecasting/_autoets_gradient_params.py | 391 ++++++++++++++++ aeon/forecasting/_ets.py | 2 +- aeon/forecasting/_ets_fast.py | 408 +++++++++++++++++ aeon/forecasting/_ets_fast_structtest.py | 431 ++++++++++++++++++ aeon/forecasting/_verify_ets.py | 244 +++++++--- .../tests/test_autoets_gradient_params.py | 21 + 6 files changed, 1441 insertions(+), 56 deletions(-) create mode 100644 aeon/forecasting/_autoets_gradient_params.py create mode 100644 aeon/forecasting/_ets_fast.py create mode 100644 aeon/forecasting/_ets_fast_structtest.py create mode 100644 aeon/forecasting/tests/test_autoets_gradient_params.py diff --git a/aeon/forecasting/_autoets_gradient_params.py b/aeon/forecasting/_autoets_gradient_params.py new file mode 100644 index 0000000000..b49e5f29f5 --- /dev/null +++ b/aeon/forecasting/_autoets_gradient_params.py @@ -0,0 +1,391 @@ +"""AutoETSForecaster class. + +Extends the ETSForecaster to automatically calculate the smoothing parameters + +aeon enhancement proposal +https://github.com/aeon-toolkit/aeon/pull/2244/ + +""" + +__maintainer__ = [] +__all__ = ["AutoETSForecaster"] + +import numpy as np +import torch + +from aeon.forecasting._ets_fast import ADDITIVE, MULTIPLICATIVE, NONE +from aeon.forecasting.base import BaseForecaster + + +class AutoETSForecaster(BaseForecaster): + """Exponential Smoothing forecaster. + + An implementation of the exponential smoothing statistics forecasting algorithm. + Implements additive and multiplicative error models, + None, additive and multiplicative (including damped) trend and + None, additive and mutliplicative seasonality[1]_. + + Parameters + ---------- + alpha : float, default = 0.1 + Level smoothing parameter. + beta : float, default = 0.01 + Trend smoothing parameter. + gamma : float, default = 0.01 + Seasonal smoothing parameter. + phi : float, default = 0.99 + Trend damping smoothing parameters + horizon : int, default = 1 + The horizon to forecast to. + model_type : ModelType, default = ModelType() + A object of type ModelType, describing the error, + trend and seasonality type of this ETS model. + + References + ---------- + .. [1] R. J. Hyndman and G. Athanasopoulos, + Forecasting: Principles and Practice. Melbourne, Australia: OTexts, 2014. + + Examples + -------- + >>> from aeon.forecasting import ETSForecaster, ModelType + >>> from aeon.datasets import load_airline + >>> y = load_airline() + >>> forecaster = ETSForecaster(alpha=0.4, beta=0.2, gamma=0.5, phi=0.8, horizon=1, + model_type=ModelType(1,2,2,4)) + >>> forecaster.fit(y) + >>> forecaster.predict() + 366.90200486015596 + """ + + def __init__( + self, + error_type=ADDITIVE, + trend_type=NONE, + seasonality_type=NONE, + seasonal_period=1, + horizon=1, + ): + assert error_type != NONE, "Error must be either additive or multiplicative" + if seasonal_period < 1 or seasonality_type == NONE: + seasonal_period = 1 + self.alpha = torch.tensor(0.1, requires_grad=True) # Level smoothing + self.beta = torch.tensor(0.05, requires_grad=True) # Trend smoothing + self.gamma = torch.tensor(0.05, requires_grad=True) # Seasonality smoothing + self.phi = torch.tensor(0.98, requires_grad=True) # Damping factor + if trend_type == NONE: + self.beta = 0 + if seasonality_type == NONE: + self.gamma = 0 + self.forecast_val_ = 0.0 + self.level = (0,) + self.trend = (0,) + self.seasonality = np.zeros(1, dtype=np.float64) + self.n_timepoints = 0 + self.avg_mean_sq_err_ = 0 + self.liklihood_ = 0 + self.residuals_ = [] + self.error_type = error_type + self.trend_type = trend_type + self.seasonality_type = seasonality_type + self.seasonal_period = seasonal_period + super().__init__(horizon=horizon, axis=1) + + def _fit(self, y, exog=None): + """Fit Exponential Smoothing forecaster to series y. + + Fit a forecaster to predict self.horizon steps ahead using y. + + Parameters + ---------- + y : np.ndarray + A time series on which to learn a forecaster to predict horizon ahead + exog : np.ndarray, default =None + Optional exogenous time series data assumed to be aligned with y + + Returns + ------- + self + Fitted BaseForecaster. + """ + data = np.array(y.squeeze(), dtype=np.float64) + ( + self.level, + self.trend, + self.seasonality, + self.residuals_, + self.avg_mean_sq_err_, + self.liklihood_, + ) = _fit( + data, + self.error_type, + self.trend_type, + self.seasonality_type, + self.seasonal_period, + ) + return self + + def _predict(self, y=None, exog=None): + """ + Predict the next horizon steps ahead. + + Parameters + ---------- + y : np.ndarray, default = None + A time series to predict the next horizon value for. If None, + predict the next horizon value after series seen in fit. + exog : np.ndarray, default =None + Optional exogenous time series data assumed to be aligned with y + + Returns + ------- + float + single prediction self.horizon steps ahead of y. + """ + y = np.array(y, dtype=np.float64) + + return _predict( + self.trend_type, + self.seasonality_type, + self.level, + self.trend, + self.seasonality, + self.phi, + self.horizon, + self.n_timepoints, + self.seasonal_period, + ) + + +def _fit(data, error_type, trend_type, seasonality_type, seasonal_period): + torch.autograd.set_detect_anomaly(True) + data = torch.tensor(data) + n_timepoints = len(data) + # print(typeof(self.states.level)) + # print(typeof(data)) + # print(typeof(self.states.seasonality)) + # print(typeof(np.full(self.model_type.seasonal_period, self.states.level))) + # print(typeof(data[: self.model_type.seasonal_period])) + level, trend, seasonality = _initialise( + trend_type, seasonality_type, seasonal_period, data + ) + alpha = torch.tensor(0.1, requires_grad=True) # Level smoothing + beta = torch.tensor(0.05, requires_grad=True) # Trend smoothing + gamma = torch.tensor(0.05, requires_grad=True) # Seasonality smoothing + phi = torch.tensor(0.98, requires_grad=True) # Damping factor + batch_size = seasonal_period * 2 + num_batches = len(data) // batch_size + # residuals_ = torch.zeros(n_timepoints) # 1 Less residual than data points + optimizer = torch.optim.SGD([alpha, beta, gamma, phi], lr=0.001) + for _epoch in range(100): # number of epochs + for i in range(1, num_batches): + batch_of_data = data[i * batch_size : (i + 1) * batch_size] + liklihood_ = torch.tensor(0, dtype=torch.float64) + mul_liklihood_pt2 = torch.tensor(0, dtype=torch.float64) + for t, data_item in enumerate(batch_of_data): + # Calculate level, trend, and seasonal components + fitted_value, error, level, trend, seasonality[t % seasonal_period] = ( + _update_states( + error_type, + trend_type, + seasonality_type, + level, + trend, + seasonality[t % seasonal_period], + data_item, + alpha, + beta, + gamma, + phi, + ) + ) + # residuals_[t] = error + liklihood_ += error * error + mul_liklihood_pt2 += torch.log(torch.abs(fitted_value)) + liklihood_ = (n_timepoints - seasonal_period) * torch.log(liklihood_) + if error_type == MULTIPLICATIVE: + liklihood_ += 2 * mul_liklihood_pt2 + liklihood_.backward() + optimizer.step() + optimizer.zero_grad() + level = level.clone().detach() + trend = trend.clone().detach() + seasonality = seasonality.clone().detach() + return alpha, beta, gamma, phi + + +def _predict( + trend_type, + seasonality_type, + level, + trend, + seasonality, + phi, + horizon, + n_timepoints, + seasonal_period, +): + # Generate forecasts based on the final values of level, trend, and seasonals + if phi == 1: # No damping case + phi_h = float(horizon) + else: + # Geometric series formula for calculating phi + phi^2 + ... + phi^h + phi_h = phi * (1 - phi**horizon) / (1 - phi) + seasonal_index = (n_timepoints + horizon) % seasonal_period + return _predict_value( + trend_type, seasonality_type, level, trend, seasonality[seasonal_index], phi_h + )[0] + + +def _initialise(trend_type, seasonality_type, seasonal_period, data): + """ + Initialize level, trend, and seasonality values for the ETS model. + + Parameters + ---------- + data : array-like + The time series data + (should contain at least two full seasons if seasonality is specified) + """ + # Initial Level: Mean of the first season + level = torch.mean(data[:seasonal_period]) + # Initial Trend + if trend_type == ADDITIVE: + # Average difference between corresponding points in the first two seasons + trend = torch.mean( + data[seasonal_period : 2 * seasonal_period] - data[:seasonal_period] + ) + elif trend_type == MULTIPLICATIVE: + # Average ratio between corresponding points in the first two seasons + trend = torch.mean( + data[seasonal_period : 2 * seasonal_period] / data[:seasonal_period] + ) + else: + # No trend + trend = 0 + # Initial Seasonality + if seasonality_type == ADDITIVE: + # Seasonal component is the difference + # from the initial level for each point in the first season + seasonality = data[:seasonal_period] - level + elif seasonality_type == MULTIPLICATIVE: + # Seasonal component is the ratio of each point in the first season + # to the initial level + seasonality = data[:seasonal_period] / level + else: + # No seasonality + seasonality = np.zeros(1) + return level, trend, seasonality + + +def _update_states( + error_type, + trend_type, + seasonality_type, + curr_level, + curr_trend, + curr_seasonality, + data_item: int, + alpha, + beta, + gamma, + phi, +): + """ + Update level, trend, and seasonality components. + + Using state space equations for an ETS model. + + Parameters + ---------- + data_item: float + The current value of the time series. + seasonal_index: int + The index to update the seasonal component. + """ + # Retrieve the current state values + fitted_value, damped_trend, trend_level_combination = _predict_value( + trend_type, seasonality_type, curr_level, curr_trend, curr_seasonality, phi + ) + # Calculate the error term (observed value - fitted value) + if error_type == MULTIPLICATIVE: + error = data_item / fitted_value - 1 # Multiplicative error + else: + error = data_item - fitted_value # Additive error + # Update level + if error_type == MULTIPLICATIVE: + level = trend_level_combination * (1 + alpha * error) + trend = damped_trend * (1 + beta * error) + seasonality = curr_seasonality * (1 + gamma * error) + if seasonality_type == ADDITIVE: + level += alpha * error * curr_seasonality # Add seasonality correction + seasonality += gamma * error * trend_level_combination + if trend_type == ADDITIVE: + trend += (curr_level + curr_seasonality) * beta * error + else: + trend += curr_seasonality / curr_level * beta * error + elif trend_type == ADDITIVE: + trend += curr_level * beta * error + else: + level_correction = 1 + trend_correction = 1 + seasonality_correction = 1 + if seasonality_type == MULTIPLICATIVE: + # Add seasonality correction + level_correction *= curr_seasonality.clone() + trend_correction *= curr_seasonality.clone() + seasonality_correction *= trend_level_combination.clone() + if trend_type == MULTIPLICATIVE: + trend_correction *= curr_level.clone() + level = ( + trend_level_combination.clone() + + alpha.clone() * error.clone() / level_correction + ) + trend = damped_trend.clone() + beta.clone() * error.clone() / trend_correction + seasonality = ( + curr_seasonality.clone() + + gamma.clone() * error.clone() / seasonality_correction + ) + return (fitted_value, error, level, trend, seasonality) + + +def _predict_value(trend_type, seasonality_type, level, trend, seasonality, phi): + """ + + Generate various useful values, including the next fitted value. + + Parameters + ---------- + trend : float + The current trend value for the model + level : float + The current level value for the model + seasonality : float + The current seasonality value for the model + phi : float + The damping parameter for the model + + Returns + ------- + fitted_value : float + single prediction based on the current state variables. + damped_trend : float + The damping parameter combined with the trend dependant on the model type + trend_level_combination : float + Combination of the trend and level based on the model type. + """ + # Apply damping parameter and + # calculate commonly used combination of trend and level components + if trend_type == MULTIPLICATIVE: + damped_trend = trend.clone() ** phi.clone() + trend_level_combination = level.clone() * damped_trend.clone() + else: # Additive trend, if no trend, then trend = 0 + damped_trend = trend.clone() * phi.clone() + trend_level_combination = level.clone() + damped_trend.clone() + + # Calculate forecast (fitted value) based on the current components + if seasonality_type == MULTIPLICATIVE: + fitted_value = trend_level_combination.clone() * seasonality.clone() + else: # Additive seasonality, if no seasonality, then seasonality = 0 + fitted_value = trend_level_combination.clone() + seasonality.clone() + return fitted_value, damped_trend, trend_level_combination diff --git a/aeon/forecasting/_ets.py b/aeon/forecasting/_ets.py index 7c97378002..4d775821f3 100644 --- a/aeon/forecasting/_ets.py +++ b/aeon/forecasting/_ets.py @@ -104,12 +104,12 @@ class ETSForecaster(BaseForecaster): def __init__( self, + model_type=default_model_type, alpha=0.1, beta=0.01, gamma=0.01, phi=0.99, horizon=1, - model_type=default_model_type, ): self.alpha = alpha self.beta = beta diff --git a/aeon/forecasting/_ets_fast.py b/aeon/forecasting/_ets_fast.py new file mode 100644 index 0000000000..33397f7637 --- /dev/null +++ b/aeon/forecasting/_ets_fast.py @@ -0,0 +1,408 @@ +"""ETSForecaster class. + +An implementation of the exponential smoothing statistics forecasting algorithm. +Implements additive and multiplicative error models, +None, additive and multiplicative (including damped) trend and +None, additive and mutliplicative seasonality + +aeon enhancement proposal +https://github.com/aeon-toolkit/aeon/pull/2244/ + +""" + +__maintainer__ = [] +__all__ = ["ETSForecaster"] + +import numpy as np +from numba import njit + +from aeon.forecasting.base import BaseForecaster + +NOGIL = False +CACHE = True + +NONE = 0 +ADDITIVE = 1 +MULTIPLICATIVE = 2 + + +class ETSForecaster(BaseForecaster): + """Exponential Smoothing forecaster. + + An implementation of the exponential smoothing statistics forecasting algorithm. + Implements additive and multiplicative error models, + None, additive and multiplicative (including damped) trend and + None, additive and mutliplicative seasonality[1]_. + + Parameters + ---------- + alpha : float, default = 0.1 + Level smoothing parameter. + beta : float, default = 0.01 + Trend smoothing parameter. + gamma : float, default = 0.01 + Seasonal smoothing parameter. + phi : float, default = 0.99 + Trend damping smoothing parameters + horizon : int, default = 1 + The horizon to forecast to. + model_type : ModelType, default = ModelType() + A object of type ModelType, describing the error, + trend and seasonality type of this ETS model. + + References + ---------- + .. [1] R. J. Hyndman and G. Athanasopoulos, + Forecasting: Principles and Practice. Melbourne, Australia: OTexts, 2014. + + Examples + -------- + >>> from aeon.forecasting import ETSForecaster, ModelType + >>> from aeon.datasets import load_airline + >>> y = load_airline() + >>> forecaster = ETSForecaster(alpha=0.4, beta=0.2, gamma=0.5, phi=0.8, horizon=1, + model_type=ModelType(1,2,2,4)) + >>> forecaster.fit(y) + >>> forecaster.predict() + 366.90200486015596 + """ + + def __init__( + self, + error_type=ADDITIVE, + trend_type=NONE, + seasonality_type=NONE, + seasonal_period=1, + alpha=0.1, + beta=0.01, + gamma=0.01, + phi=0.99, + horizon=1, + ): + assert error_type != NONE, "Error must be either additive or multiplicative" + if seasonal_period < 1 or seasonality_type == NONE: + seasonal_period = 1 + self.alpha = alpha + self.beta = beta + self.gamma = gamma + self.phi = phi + if trend_type == NONE: + self.beta = 0 + if seasonality_type == NONE: + self.gamma = 0 + self.forecast_val_ = 0.0 + self.level = (0,) + self.trend = (0,) + self.seasonality = np.zeros(1, dtype=np.float64) + self.n_timepoints = 0 + self.avg_mean_sq_err_ = 0 + self.liklihood_ = 0 + self.residuals_ = [] + self.error_type = error_type + self.trend_type = trend_type + self.seasonality_type = seasonality_type + self.seasonal_period = seasonal_period + super().__init__(horizon=horizon, axis=1) + + def _fit(self, y, exog=None): + """Fit Exponential Smoothing forecaster to series y. + + Fit a forecaster to predict self.horizon steps ahead using y. + + Parameters + ---------- + y : np.ndarray + A time series on which to learn a forecaster to predict horizon ahead + exog : np.ndarray, default =None + Optional exogenous time series data assumed to be aligned with y + + Returns + ------- + self + Fitted BaseForecaster. + """ + data = np.array(y.squeeze(), dtype=np.float64) + ( + self.level, + self.trend, + self.seasonality, + self.residuals_, + self.avg_mean_sq_err_, + self.liklihood_, + ) = _fit( + data, + self.error_type, + self.trend_type, + self.seasonality_type, + self.seasonal_period, + self.alpha, + self.beta, + self.gamma, + self.phi, + ) + return self + + def _predict(self, y=None, exog=None): + """ + Predict the next horizon steps ahead. + + Parameters + ---------- + y : np.ndarray, default = None + A time series to predict the next horizon value for. If None, + predict the next horizon value after series seen in fit. + exog : np.ndarray, default =None + Optional exogenous time series data assumed to be aligned with y + + Returns + ------- + float + single prediction self.horizon steps ahead of y. + """ + y = np.array(y, dtype=np.float64) + + return _predict( + self.trend_type, + self.seasonality_type, + self.level, + self.trend, + self.seasonality, + self.phi, + self.horizon, + self.n_timepoints, + self.seasonal_period, + ) + + +@njit(nogil=NOGIL, cache=CACHE) +def _fit( + data, + error_type, + trend_type, + seasonality_type, + seasonal_period, + alpha, + beta, + gamma, + phi, +): + n_timepoints = len(data) + # print(typeof(self.states.level)) + # print(typeof(data)) + # print(typeof(self.states.seasonality)) + # print(typeof(np.full(self.model_type.seasonal_period, self.states.level))) + # print(typeof(data[: self.model_type.seasonal_period])) + level, trend, seasonality = _initialise( + trend_type, seasonality_type, seasonal_period, data + ) + avg_mean_sq_err_ = 0 + liklihood_ = 0 + mul_liklihood_pt2 = 0 + residuals_ = np.zeros(n_timepoints) # 1 Less residual than data points + for t, data_item in enumerate(data[seasonal_period:]): + # Calculate level, trend, and seasonal components + fitted_value, error, level, trend, seasonality[t % seasonal_period] = ( + _update_states( + error_type, + trend_type, + seasonality_type, + level, + trend, + seasonality[t % seasonal_period], + data_item, + alpha, + beta, + gamma, + phi, + ) + ) + residuals_[t] = error + avg_mean_sq_err_ += (data_item - fitted_value) ** 2 + liklihood_ += error * error + mul_liklihood_pt2 += np.log(np.fabs(fitted_value)) + avg_mean_sq_err_ /= n_timepoints - seasonal_period + liklihood_ = (n_timepoints - seasonal_period) * np.log(liklihood_) + if error_type == MULTIPLICATIVE: + liklihood_ += 2 * mul_liklihood_pt2 + return level, trend, seasonality, residuals_, avg_mean_sq_err_, liklihood_ + + +def _predict( + trend_type, + seasonality_type, + level, + trend, + seasonality, + phi, + horizon, + n_timepoints, + seasonal_period, +): + # Generate forecasts based on the final values of level, trend, and seasonals + if phi == 1: # No damping case + phi_h = float(horizon) + else: + # Geometric series formula for calculating phi + phi^2 + ... + phi^h + phi_h = phi * (1 - phi**horizon) / (1 - phi) + seasonal_index = (n_timepoints + horizon) % seasonal_period + return _predict_value( + trend_type, + seasonality_type, + level, + trend, + seasonality[seasonal_index], + phi_h, + )[0] + + +@njit(nogil=NOGIL, cache=CACHE) +def _initialise(trend_type, seasonality_type, seasonal_period, data): + """ + Initialize level, trend, and seasonality values for the ETS model. + + Parameters + ---------- + data : array-like + The time series data + (should contain at least two full seasons if seasonality is specified) + """ + # Initial Level: Mean of the first season + level = np.mean(data[:seasonal_period]) + # Initial Trend + if trend_type == ADDITIVE: + # Average difference between corresponding points in the first two seasons + trend = np.mean( + data[seasonal_period : 2 * seasonal_period] - data[:seasonal_period] + ) + elif trend_type == MULTIPLICATIVE: + # Average ratio between corresponding points in the first two seasons + trend = np.mean( + data[seasonal_period : 2 * seasonal_period] / data[:seasonal_period] + ) + else: + # No trend + trend = 0 + # Initial Seasonality + if seasonality_type == ADDITIVE: + # Seasonal component is the difference + # from the initial level for each point in the first season + seasonality = data[:seasonal_period] - level + elif seasonality_type == MULTIPLICATIVE: + # Seasonal component is the ratio of each point in the first season + # to the initial level + seasonality = data[:seasonal_period] / level + else: + # No seasonality + seasonality = np.zeros(1) + return level, trend, seasonality + + +@njit(nogil=NOGIL, cache=CACHE) +def _update_states( + error_type, + trend_type, + seasonality_type, + level, + trend, + seasonality, + data_item: int, + alpha, + beta, + gamma, + phi, +): + """ + Update level, trend, and seasonality components. + + Using state space equations for an ETS model. + + Parameters + ---------- + data_item: float + The current value of the time series. + seasonal_index: int + The index to update the seasonal component. + """ + # Retrieve the current state values + curr_level = level + curr_seasonality = seasonality + fitted_value, damped_trend, trend_level_combination = _predict_value( + trend_type, seasonality_type, level, trend, seasonality, phi + ) + # Calculate the error term (observed value - fitted value) + if error_type == MULTIPLICATIVE: + error = data_item / fitted_value - 1 # Multiplicative error + else: + error = data_item - fitted_value # Additive error + # Update level + if error_type == MULTIPLICATIVE: + level = trend_level_combination * (1 + alpha * error) + trend = damped_trend * (1 + beta * error) + seasonality = curr_seasonality * (1 + gamma * error) + if seasonality_type == ADDITIVE: + level += alpha * error * curr_seasonality # Add seasonality correction + seasonality += gamma * error * trend_level_combination + if trend_type == ADDITIVE: + trend += (curr_level + curr_seasonality) * beta * error + else: + trend += curr_seasonality / curr_level * beta * error + elif trend_type == ADDITIVE: + trend += curr_level * beta * error + else: + level_correction = 1 + trend_correction = 1 + seasonality_correction = 1 + if seasonality_type == MULTIPLICATIVE: + # Add seasonality correction + level_correction *= curr_seasonality + trend_correction *= curr_seasonality + seasonality_correction *= trend_level_combination + if trend_type == MULTIPLICATIVE: + trend_correction *= curr_level + level = trend_level_combination + alpha * error / level_correction + trend = damped_trend + beta * error / trend_correction + seasonality = curr_seasonality + gamma * error / seasonality_correction + return (fitted_value, error, level, trend, seasonality) + + +@njit(nogil=NOGIL, cache=CACHE) +def _predict_value(trend_type, seasonality_type, level, trend, seasonality, phi): + """ + + Generate various useful values, including the next fitted value. + + Parameters + ---------- + trend : float + The current trend value for the model + level : float + The current level value for the model + seasonality : float + The current seasonality value for the model + phi : float + The damping parameter for the model + + Returns + ------- + fitted_value : float + single prediction based on the current state variables. + damped_trend : float + The damping parameter combined with the trend dependant on the model type + trend_level_combination : float + Combination of the trend and level based on the model type. + """ + # Apply damping parameter and + # calculate commonly used combination of trend and level components + if trend_type == MULTIPLICATIVE: + damped_trend = trend**phi + trend_level_combination = level * damped_trend + else: # Additive trend, if no trend, then trend = 0 + damped_trend = trend * phi + trend_level_combination = level + damped_trend + + # Calculate forecast (fitted value) based on the current components + if seasonality_type == MULTIPLICATIVE: + fitted_value = trend_level_combination * seasonality + else: # Additive seasonality, if no seasonality, then seasonality = 0 + fitted_value = trend_level_combination + seasonality + return fitted_value, damped_trend, trend_level_combination diff --git a/aeon/forecasting/_ets_fast_structtest.py b/aeon/forecasting/_ets_fast_structtest.py new file mode 100644 index 0000000000..c685ca26ba --- /dev/null +++ b/aeon/forecasting/_ets_fast_structtest.py @@ -0,0 +1,431 @@ +"""ETSForecaster class. + +An implementation of the exponential smoothing statistics forecasting algorithm. +Implements additive and multiplicative error models, +None, additive and multiplicative (including damped) trend and +None, additive and mutliplicative seasonality + +aeon enhancement proposal +https://github.com/aeon-toolkit/aeon/pull/2244/ + +""" + +__maintainer__ = [] +__all__ = ["ETSForecaster", "ModelType"] + +import numpy as np +from numba import float64, njit +from numba.experimental import jitclass + +from aeon.forecasting.base import BaseForecaster + +NONE = 0 +ADDITIVE = 1 +MULTIPLICATIVE = 2 + + +@jitclass +class ModelType: + """ + Class describing the error, trend and seasonality model of an ETS forecaster. + + Attributes + ---------- + error_type : int + The type of error model; either Additive(1) or Multiplicative(2) + trend_type : int + The type of trend model; one of None(0), additive(1) or multiplicative(2). + seasonality_type : int + The type of seasonality model; one of None(0), additive(1) or multiplicative(2). + seasonal_period : int + The period of the seasonality (m) (e.g., for quaterly data seasonal_period = 4). + """ + + error_type: int + trend_type: int + seasonality_type: int + seasonal_period: int + + def __init__( + self, + error_type=ADDITIVE, + trend_type=NONE, + seasonality_type=NONE, + seasonal_period=1, + ): + assert error_type != NONE, "Error must be either additive or multiplicative" + if seasonal_period < 1 or seasonality_type == NONE: + seasonal_period = 1 + self.error_type = error_type + self.trend_type = trend_type + self.seasonality_type = seasonality_type + self.seasonal_period = seasonal_period + + +@jitclass([("seasonality", float64[:])]) +class StateVariables: + """ + Class describing the state variables of an ETS forecaster model. + + Attributes + ---------- + level : float + The current value of the level (l) state variable + trend : float + The current value of the trend (b) state variable + seasonality : float[] + The current value of the seasonality (s) state variable + """ + + level: float + trend: float + + def __init__(self, level=0, trend=0, seasonality=None): + self.level = level + self.trend = trend + if seasonality is None: + self.seasonality = np.zeros(1, dtype=np.float64) + else: + self.seasonality = seasonality + + +@jitclass +class SmoothingParameters: + """ + Class describing the smoothing parameters of an ETS forecaster model. + + Attributes + ---------- + alpha : float, default = 0.1 + Level smoothing parameter. + beta : float, default = 0.01 + Trend smoothing parameter. + gamma : float, default = 0.01 + Seasonal smoothing parameter. + phi : float, default = 0.99 + Trend damping smoothing parameters + """ + + alpha: float + beta: float + gamma: float + phi: float + + def __init__(self, alpha=0.1, beta=0.01, gamma=0.01, phi=0.99): + self.alpha = alpha + self.beta = beta + self.gamma = gamma + self.phi = phi + + +class ETSForecaster(BaseForecaster): + """Exponential Smoothing forecaster. + + An implementation of the exponential smoothing statistics forecasting algorithm. + Implements additive and multiplicative error models, + None, additive and multiplicative (including damped) trend and + None, additive and mutliplicative seasonality[1]_. + + Parameters + ---------- + alpha : float, default = 0.1 + Level smoothing parameter. + beta : float, default = 0.01 + Trend smoothing parameter. + gamma : float, default = 0.01 + Seasonal smoothing parameter. + phi : float, default = 0.99 + Trend damping smoothing parameters + horizon : int, default = 1 + The horizon to forecast to. + model_type : ModelType, default = ModelType() + A object of type ModelType, describing the error, + trend and seasonality type of this ETS model. + + References + ---------- + .. [1] R. J. Hyndman and G. Athanasopoulos, + Forecasting: Principles and Practice. Melbourne, Australia: OTexts, 2014. + + Examples + -------- + >>> from aeon.forecasting import ETSForecaster, ModelType + >>> from aeon.datasets import load_airline + >>> y = load_airline() + >>> forecaster = ETSForecaster(alpha=0.4, beta=0.2, gamma=0.5, phi=0.8, horizon=1, + model_type=ModelType(1,2,2,4)) + >>> forecaster.fit(y) + >>> forecaster.predict() + 366.90200486015596 + """ + + default_model_type = ModelType() + default_smoothing_parameters = SmoothingParameters() + + def __init__( + self, + model_type=default_model_type, + smoothing_parameters=default_smoothing_parameters, + horizon=1, + ): + self.smoothing_parameters = smoothing_parameters + if model_type.trend_type == NONE: + self.smoothing_parameters.beta = 0 + if model_type.seasonality_type == NONE: + self.smoothing_parameters.gamma = 0 + self.forecast_val_ = 0.0 + self.states = StateVariables() + self.n_timepoints = 0 + self.avg_mean_sq_err_ = 0 + self.liklihood_ = 0 + self.residuals_ = [] + self.model_type = model_type + super().__init__(horizon=horizon, axis=1) + + def _fit(self, y, exog=None): + """Fit Exponential Smoothing forecaster to series y. + + Fit a forecaster to predict self.horizon steps ahead using y. + + Parameters + ---------- + y : np.ndarray + A time series on which to learn a forecaster to predict horizon ahead + exog : np.ndarray, default =None + Optional exogenous time series data assumed to be aligned with y + + Returns + ------- + self + Fitted BaseForecaster. + """ + data = np.array(y.squeeze(), dtype=np.float64) + self.n_timepoints = len(data) + # print(typeof(self.states.level)) + # print(typeof(data)) + # print(typeof(self.states.seasonality)) + # print(typeof(np.full(self.model_type.seasonal_period, self.states.level))) + # print(typeof(data[: self.model_type.seasonal_period])) + _initialise(self.model_type, self.states, data) + self.avg_mean_sq_err_ = 0 + self.liklihood_ = 0 + mul_liklihood_pt2 = 0 + self.residuals_ = np.zeros( + self.n_timepoints + ) # 1 Less residual than data points + for t, data_item in enumerate(data[self.model_type.seasonal_period :]): + # Calculate level, trend, and seasonal components + fitted_value, error = _update_states( + self.model_type, + self.states, + data_item, + t % self.model_type.seasonal_period, + self.smoothing_parameters, + ) + self.residuals_[t] = error + self.avg_mean_sq_err_ += (data_item - fitted_value) ** 2 + self.liklihood_ += error * error + mul_liklihood_pt2 += np.log(np.fabs(fitted_value)) + self.avg_mean_sq_err_ /= self.n_timepoints - self.model_type.seasonal_period + self.liklihood_ = ( + self.n_timepoints - self.model_type.seasonal_period + ) * np.log(self.liklihood_) + if self.model_type.error_type == MULTIPLICATIVE: + self.liklihood_ += 2 * mul_liklihood_pt2 + return self + + def _predict(self, y=None, exog=None): + """ + Predict the next horizon steps ahead. + + Parameters + ---------- + y : np.ndarray, default = None + A time series to predict the next horizon value for. If None, + predict the next horizon value after series seen in fit. + exog : np.ndarray, default =None + Optional exogenous time series data assumed to be aligned with y + + Returns + ------- + float + single prediction self.horizon steps ahead of y. + """ + y = np.array(y, dtype=np.float64) + # Generate forecasts based on the final values of level, trend, and seasonals + if self.smoothing_parameters.phi == 1: # No damping case + phi_h = float(self.horizon) + else: + # Geometric series formula for calculating phi + phi^2 + ... + phi^h + phi_h = ( + self.smoothing_parameters.phi + * (1 - self.smoothing_parameters.phi**self.horizon) + / (1 - self.smoothing_parameters.phi) + ) + seasonal_index = ( + self.n_timepoints + self.horizon + ) % self.model_type.seasonal_period + fitted_value = _predict_value( + self.model_type, self.states, seasonal_index, phi_h + )[0] + return fitted_value + + +@njit +def _initialise(model: ModelType, states: StateVariables, data): + """ + Initialize level, trend, and seasonality values for the ETS model. + + Parameters + ---------- + data : array-like + The time series data + (should contain at least two full seasons if seasonality is specified) + """ + # Initial Level: Mean of the first season + states.level = np.mean(data[: model.seasonal_period]) + # Initial Trend + if model.trend_type == ADDITIVE: + # Average difference between corresponding points in the first two seasons + states.trend = np.mean( + data[model.seasonal_period : 2 * model.seasonal_period] + - data[: model.seasonal_period] + ) + elif model.trend_type == MULTIPLICATIVE: + # Average ratio between corresponding points in the first two seasons + states.trend = np.mean( + data[model.seasonal_period : 2 * model.seasonal_period] + / data[: model.seasonal_period] + ) + else: + # No trend + states.trend = 0 + # Initial Seasonality + if model.seasonality_type == ADDITIVE: + # Seasonal component is the difference + # from the initial level for each point in the first season + states.seasonality = data[: model.seasonal_period] - states.level + elif model.seasonality_type == MULTIPLICATIVE: + # Seasonal component is the ratio of each point in the first season + # to the initial level + states.seasonality = data[: model.seasonal_period] / states.level + else: + # No seasonality + states.seasonality = np.zeros(1) + + +@njit +def _update_states( + model: ModelType, + states: StateVariables, + data_item: int, + seasonal_index: int, + parameters: SmoothingParameters, +): + """ + Update level, trend, and seasonality components. + + Using state space equations for an ETS model. + + Parameters + ---------- + data_item: float + The current value of the time series. + seasonal_index: int + The index to update the seasonal component. + """ + # Retrieve the current state values + level = states.level + seasonality = states.seasonality[seasonal_index] + fitted_value, damped_trend, trend_level_combination = _predict_value( + model, states, seasonal_index, parameters.phi + ) + # Calculate the error term (observed value - fitted value) + if model.error_type == MULTIPLICATIVE: + error = data_item / fitted_value - 1 # Multiplicative error + else: + error = data_item - fitted_value # Additive error + # Update level + if model.error_type == MULTIPLICATIVE: + states.level = trend_level_combination * (1 + parameters.alpha * error) + states.trend = damped_trend * (1 + parameters.beta * error) + states.seasonality[seasonal_index] = seasonality * ( + 1 + parameters.gamma * error + ) + if model.seasonality_type == ADDITIVE: + states.level += ( + parameters.alpha * error * seasonality + ) # Add seasonality correction + states.seasonality[seasonal_index] += ( + parameters.gamma * error * trend_level_combination + ) + if model.trend_type == ADDITIVE: + states.trend += (level + seasonality) * parameters.beta * error + else: + states.trend += seasonality / level * parameters.beta * error + elif model.trend_type == ADDITIVE: + states.trend += level * parameters.beta * error + else: + level_correction = 1 + trend_correction = 1 + seasonality_correction = 1 + if model.seasonality_type == MULTIPLICATIVE: + # Add seasonality correction + level_correction *= seasonality + trend_correction *= seasonality + seasonality_correction *= trend_level_combination + if model.trend_type == MULTIPLICATIVE: + trend_correction *= level + states.level = ( + trend_level_combination + parameters.alpha * error / level_correction + ) + states.trend = damped_trend + parameters.beta * error / trend_correction + states.seasonality[seasonal_index] = ( + seasonality + parameters.gamma * error / seasonality_correction + ) + return (fitted_value, error) + + +@njit +def _predict_value( + model: ModelType, states: StateVariables, seasonality_index: int, phi +): + """ + + Generate various useful values, including the next fitted value. + + Parameters + ---------- + trend : float + The current trend value for the model + level : float + The current level value for the model + seasonality : float + The current seasonality value for the model + phi : float + The damping parameter for the model + + Returns + ------- + fitted_value : float + single prediction based on the current state variables. + damped_trend : float + The damping parameter combined with the trend dependant on the model type + trend_level_combination : float + Combination of the trend and level based on the model type. + """ + # Apply damping parameter and + # calculate commonly used combination of trend and level components + if model.trend_type == MULTIPLICATIVE: + damped_trend = states.trend**phi + trend_level_combination = states.level * damped_trend + else: # Additive trend, if no trend, then trend = 0 + damped_trend = states.trend * phi + trend_level_combination = states.level + damped_trend + + # Calculate forecast (fitted value) based on the current components + if model.seasonality_type == MULTIPLICATIVE: + fitted_value = trend_level_combination * states.seasonality[seasonality_index] + else: # Additive seasonality, if no seasonality, then seasonality = 0 + fitted_value = trend_level_combination + states.seasonality[seasonality_index] + return fitted_value, damped_trend, trend_level_combination diff --git a/aeon/forecasting/_verify_ets.py b/aeon/forecasting/_verify_ets.py index d93b64631f..35b35546bb 100644 --- a/aeon/forecasting/_verify_ets.py +++ b/aeon/forecasting/_verify_ets.py @@ -1,10 +1,13 @@ import random import time +import timeit import numpy as np -from statsforecast.ets import etscalc +from statsforecast.ets import pegelsresid_C from statsforecast.utils import AirPassengers as ap +import aeon.forecasting._ets_fast as etsfast +import aeon.forecasting._ets_fast_structtest as ets_structtest from aeon.forecasting import ETSForecaster, ModelType NA = -99999.0 @@ -15,7 +18,6 @@ def setup(): """Generate parameters required for ETS algorithms.""" y = ap - n = len(ap) m = random.randint(1, 24) error = random.randint(1, 2) trend = random.randint(0, 2) @@ -32,26 +34,54 @@ def setup(): phi = round( random.random() * 0.18 + 0.8, 4 ) # Common constraint for phi is 0.8 < phi < 0.98 - e = np.zeros(n) - lik_fitets = np.zeros(1) - amse = np.zeros(MAX_NMSE) - nmse = 3 - return ( - y, - n, + return (y, m, error, trend, season, alpha, beta, gamma, phi) + + +def obscure_statsforecast_version( + y: np.ndarray, + m: int, + init_state: np.ndarray, + errortype: int, + trendtype: int, + seasontype: int, + alpha: float, + beta: float, + gamma: float, + phi: float, + nmse: int, +): + """Hide the differences between different statsforecast versions.""" + amse, e, x, lik = pegelsresid_C( + y[m:], m, - error, - trend, - season, + init_state, + "A" if errortype == 1 else "M", + "A" if trendtype == 1 else "M" if trendtype == 2 else "N", + "A" if seasontype == 1 else "M" if seasontype == 2 else "N", + phi != 1, alpha, beta, gamma, phi, - e, - lik_fitets, - amse, nmse, ) + # e = np.zeros(len(y)) + # amse = np.zeros(MAX_NMSE) + # lik = etscalc(y[m:], + # len(y) - m, + # init_state, + # m, + # errortype, + # trendtype, + # seasontype, + # alpha, + # beta, + # gamma, + # phi, + # e, + # amse, + # nmse) + return amse, e, lik def test_ets_comparison(setup_func, random_seed, catch_errors): @@ -59,7 +89,6 @@ def test_ets_comparison(setup_func, random_seed, catch_errors): random.seed(random_seed) ( y, - n, m, error, trend, @@ -68,53 +97,39 @@ def test_ets_comparison(setup_func, random_seed, catch_errors): beta, gamma, phi, - e, - lik_fitets, - amse, - nmse, ) = setup_func() # tsml-eval implementation - start = time.time() - f1 = ETSForecaster(alpha, beta, gamma, phi, 1, ModelType(error, trend, season, m)) + start = time.perf_counter() + f1 = ETSForecaster( + ModelType(error, trend, season, m), + alpha, + beta, + gamma, + phi, + 1, + ) f1.fit(y) - end = time.time() + end = time.perf_counter() time_fitets = end - start e_fitets = f1.residuals_ amse_fitets = f1.avg_mean_sq_err_ lik_fitets = f1.liklihood_ - # Reinitialise arrays - e.fill(0) - amse.fill(0) - f1 = ETSForecaster(alpha, beta, gamma, phi, 1, ModelType(error, trend, season, m)) + f1 = ETSForecaster(ModelType(error, trend, season, m), alpha, beta, gamma, phi, 1) f1._initialise(y) - init_states_etscalc = np.zeros(n * (1 + (trend > 0) + m * (season > 0) + 1)) + init_states_etscalc = np.zeros(len(y) * (1 + (trend > 0) + m * (season > 0) + 1)) init_states_etscalc[0] = f1.level_ init_states_etscalc[1] = f1.trend_ init_states_etscalc[1 + (trend != 0) : m + 1 + (trend != 0)] = f1.season_[::-1] if season == 0: m = 1 # Nixtla/statsforcast implementation - start = time.time() - lik_etscalc = etscalc( - y[m:], - n - m, - init_states_etscalc, - m, - error, - trend, - season, - alpha, - beta, - gamma, - phi, - e, - amse, - nmse, + start = time.perf_counter() + amse_etscalc, e_etscalc, lik_etscalc = obscure_statsforecast_version( + y, m, init_states_etscalc, error, trend, season, alpha, beta, gamma, phi, 1 ) - end = time.time() + end = time.perf_counter() time_etscalc = end - start - e_etscalc = e.copy() - amse_etscalc = amse.copy()[0] + amse_etscalc = amse_etscalc[0] if catch_errors: try: @@ -153,16 +168,135 @@ def test_ets_comparison(setup_func, random_seed, catch_errors): assert np.allclose(e_fitets, e_etscalc) assert np.allclose(amse_fitets, amse_etscalc) assert np.isclose(lik_fitets, lik_etscalc) - print(time_fitets) # noqa - print(time_etscalc) # noqa + print(f"Time for ETS: {time_fitets:0.20f}") # noqa + print(f"Time for statsforecast ETS: {time_etscalc}") # noqa return True +def time_etsfast(): + """Test function for optimised numba ets algorithm.""" + etsfast.ETSForecaster(2, 2, 2, 4).fit(ap).predict() + + +def time_ets_structtest(): + """Test function for ets algorithm using classes.""" + ets_structtest.ETSForecaster(ets_structtest.ModelType(2, 2, 2, 4)).fit(ap).predict() + + +def time_etsnoopt(): + """Test function for non-optimised ets algorithm.""" + ETSForecaster(ModelType(2, 2, 2, 4)).fit(ap).predict() + + +def time_etsfast_noclass(): + """Test function for optimised ets algorithm without the class based structure.""" + data = np.array(ap.squeeze(), dtype=np.float64) + (level, trend, seasonality, residuals_, avg_mean_sq_err_, liklihood_) = ( + etsfast._fit(data, 2, 2, 2, 4, 0.1, 0.01, 0.01, 0.99) + ) + etsfast._predict(2, 2, level, trend, seasonality, 0.99, 1, 144, 4) + + +def time_sf(): + """Test function for statsforecast ets algorithm.""" + x = np.zeros(144 * 7) + x[0:6] = [122.75, 1.123230970596215, 0.91242363, 0.96130346, 1.07535642, 1.0509165] + obscure_statsforecast_version( + ap[4:], + 4, + x, + 2, + 2, + 2, + 0.1, + 0.01, + 0.01, + 0.99, + 1, + ) + + +def time_compare(random_seed): + """Compare timings of different ets algorithms.""" + random.seed(random_seed) + (y, m, error, trend, season, alpha, beta, gamma, phi) = setup() + # etsnoopt_time = timeit.timeit(time_etsnoopt, globals={}, number=10000) + # print (f"Execution time ETS No-opt: {etsnoopt_time} seconds") + # ets_structtest_time = timeit.timeit(time_ets_structtest, globals={}, number=10000) + # print (f"Execution time ETS Structtest: {ets_structtest_time} seconds") + # Do a few iterations to remove background/overheads. Makes comparison more reliable + for _i in range(10): + time_etsfast() + time_sf() + time_etsfast_noclass() + etsfast_time = timeit.timeit(time_etsfast, globals={}, number=1000) + print(f"Execution time ETS Fast: {etsfast_time} seconds") # noqa + etsfast_noclass_time = timeit.timeit(time_etsfast_noclass, globals={}, number=1000) + print(f"Execution time ETS Fast NoClass: {etsfast_noclass_time} seconds") # noqa + statsforecast_time = timeit.timeit(time_sf, globals={}, number=1000) + print(f"Execution time StatsForecast: {statsforecast_time} seconds") # noqa + etsfast_time = timeit.timeit(time_etsfast, globals={}, number=1000) + print(f"Execution time ETS Fast: {etsfast_time} seconds") # noqa + etsfast_noclass_time = timeit.timeit(time_etsfast_noclass, globals={}, number=1000) + print(f"Execution time ETS Fast NoClass: {etsfast_noclass_time} seconds") # noqa + statsforecast_time = timeit.timeit(time_sf, globals={}, number=1000) + print(f"Execution time StatsForecast: {statsforecast_time} seconds") # noqa + # _ets_fast_nostruct implementation + start = time.perf_counter() + f3 = etsfast.ETSForecaster(error, trend, season, m, alpha, beta, gamma, phi, 1) + f3.fit(y) + end = time.perf_counter() + etsfast_time = end - start + # _ets_fast implementation + start = time.perf_counter() + f2 = ets_structtest.ETSForecaster( + ets_structtest.ModelType(error, trend, season, m), + ets_structtest.SmoothingParameters(alpha, beta, gamma, phi), + 1, + ) + f2.fit(y) + end = time.perf_counter() + ets_structtest_time = end - start + # _ets implementation + start = time.perf_counter() + f1 = ETSForecaster(ModelType(error, trend, season, m), alpha, beta, gamma, phi, 1) + f1.fit(y) + end = time.perf_counter() + etsnoopt_time = end - start + assert np.allclose(f1.residuals_, f2.residuals_) + assert np.allclose(f1.avg_mean_sq_err_, f2.avg_mean_sq_err_) + assert np.isclose(f1.liklihood_, f2.liklihood_) + assert np.allclose(f1.residuals_, f3.residuals_) + assert np.allclose(f1.avg_mean_sq_err_, f3.avg_mean_sq_err_) + assert np.isclose(f1.liklihood_, f3.liklihood_) + print( # noqa + f"ETS No-optimisation Time: {etsnoopt_time},\ + Fast Structtest time: {ets_structtest_time},\ + Fast time: {etsfast_time}" + ) + return etsnoopt_time, ets_structtest_time, etsfast_time + + if __name__ == "__main__": # np.set_printoptions(threshold=np.inf) - # test_ets_comparison(setup, 241, False) - SUCCESSES = True - for i in range(0, 30000): - SUCCESSES &= test_ets_comparison(setup, i, True) - if SUCCESSES: - print("Test Completed Successfully with no errors") # noqa + # test_ets_comparison(setup, 300, False) + # SUCCESSES = True + # for i in range(0, 30000): + # SUCCESSES &= test_ets_comparison(setup, i, True) + # if SUCCESSES: + # print("Test Completed Successfully with no errors") # noqa + time_compare(300) + # avg_ets = 0 + # avg_etsfast = 0 + # avg_etsfast_ns = 0 + # iterations = 100 + # for i in range (iterations): + # time_ets, ets_structtest_time, etsfast_time = time_compare(300) + # avg_ets += time_ets + # avg_etsfast += time_etsfast + # avg_etsfast_ns += time_etsfast_nostruct + # avg_ets/= iterations + # avg_etsfast/= iterations + # avg_etsfast_ns /= iterations + # print(f"Avg ETS Time: {avg_ets}, Avg Fast ETS time: {avg_etsfast},\ + # Avg Fast Nostruct time: {avg_etsfast_ns}") diff --git a/aeon/forecasting/tests/test_autoets_gradient_params.py b/aeon/forecasting/tests/test_autoets_gradient_params.py new file mode 100644 index 0000000000..d76673df8e --- /dev/null +++ b/aeon/forecasting/tests/test_autoets_gradient_params.py @@ -0,0 +1,21 @@ +"""Test AutoETS.""" + +# __maintainer__ = [] +# __all__ = [] +# import numpy as np + +from statsforecast.utils import AirPassengers as ap + +from aeon.forecasting._autoets_gradient_params import _fit + + +def test_autoets_forecaster(): + """TestETSForecaster.""" + parameters = _fit(ap, 1, 1, 1, 12) + print(parameters) # noqa + # assert np.allclose([parameter.item() for parameter in parameters], + # [0.1,0.05,0.05,0.98]) + + +if __name__ == "__main__": + test_autoets_forecaster()