More robust detection of time series granularity. (#135)

* More robust detection of time series granularity.

Previously, we would detect the granularity of a time series as the GCD
of all timedeltas found in the time series (assuming pandas couldn't
infer the granularity on its own). However, this behavior fails for
time series with missing data that are sampled at granularities that
aren't an exact number of seconds, e.g. monthly time series would be
resampled to a daily granularity because months are of inconsistent
length.

This commit uses the most commonly observed timedelta, and it also
checks whether a k-month granularity is a better fit for the time series
than a n-day granularity.

* Simplify code.

* Introduce more flexible granularities to models.

* Retain old behavior when offsets aren't used.

* More careful handling of resampling w/ offsets.

* Fix off-by-one error.

* Add test coverage for resampling.

* Fix typo.

* Update version matrix on old docs pages.

Author: aadyotb

.github/workflows/tests.yml

@@ -38,9 +38,9 @@ jobs:
       env:
         PYTHON_VERSION: ${{ matrix.python-version }}
       with:
-        max_attempts: 3
+        max_attempts: 1
         timeout_minutes: 60
-        retry-on: error
+        retry_on: error
         command: |
           set -euxo pipefail
           # Get a comma-separated list of the directories of all python source files

benchmark_forecast.py

@@ -27,7 +27,7 @@
 from merlion.models.forecast.base import ForecasterBase
 from merlion.transform.resample import TemporalResample, granularity_str_to_seconds
 from merlion.utils import TimeSeries, UnivariateTimeSeries
-from merlion.utils.resample import get_gcd_timedelta
+from merlion.utils.resample import infer_granularity, to_pd_datetime
 
 from ts_datasets.base import BaseDataset
 from ts_datasets.forecast import *
@@ -265,12 +265,7 @@ def train_model(
             df = df.resample(dt, closed="right", label="right").mean().interpolate()
 
         vals = TimeSeries.from_pd(df)
-        # Get time-delta
-        if not is_multivariate_data:
-            dt = df.index[1] - df.index[0]
-        else:
-            dt = get_gcd_timedelta(vals.time_stamps)
-            dt = pd.to_timedelta(dt, unit="s")
+        dt = infer_granularity(vals.time_stamps)
 
         # Get the train/val split
         t = trainval.index[np.argmax(~trainval)].value // 1e9
@@ -304,7 +299,11 @@ def train_model(
         # loop over horizon conditions
         for horizon in horizons:
             horizon = granularity_str_to_seconds(horizon)
-            max_forecast_steps = math.ceil(horizon / dt.total_seconds())
+            try:
+                max_forecast_steps = int(math.ceil(horizon / dt.total_seconds()))
+            except:
+                window = TimeSeries.from_pd(test_vals.to_pd()[: to_pd_datetime(train_end_timestamp + horizon)])
+                max_forecast_steps = len(TemporalResample(granularity=dt)(window))
             logger.debug(f"horizon is {pd.Timedelta(seconds=horizon)} and max_forecast_steps is {max_forecast_steps}")
             if retrain_type == "without_retrain":
                 retrain_freq = None

docs/build_docs.sh

@@ -5,14 +5,6 @@ set -euo pipefail
 DIRNAME=$(cd "$( dirname "${BASH_SOURCE[0]}" )" &> /dev/null && pwd)
 cd "${DIRNAME}/.."
 
-# Set up virtual environment
-pip3 install --upgrade pip setuptools wheel virtualenv
-if [ ! -d venv ]; then
-  rm -f venv
-  virtualenv venv
-fi
-source venv/bin/activate
-
 # Get current git head & stash unsaved changes
 GIT_HEAD=$(git rev-parse HEAD)
 GIT_BRANCH=$(git branch --show-current)
@@ -37,19 +29,23 @@ function exit_handler {
 }
 trap exit_handler EXIT
 
-# Install Sphinx requirements. Get old Merlion docs from gh-pages branch, but only keep the version-tagged ones.
+# Set up virtual environment & install Sphinx requirements.
+pip3 install --upgrade pip setuptools wheel virtualenv
+if [ ! -d venv ]; then
+  rm -f venv
+  virtualenv venv
+fi
+source venv/bin/activate
 pip3 install -r "${DIRNAME}/requirements.txt"
+
+# Get old Merlion docs from gh-pages branch. Only keep version-tagged ones, and update the version matrix as needed.
 git checkout gh-pages && git pull && git checkout --force "${GIT_HEAD}"
 sphinx-build -M clean "${DIRNAME}/source" "${DIRNAME}/build"
 mkdir -p "${DIRNAME}/build" "${DIRNAME}/build/html"
 git --work-tree "${DIRNAME}/build/html" checkout gh-pages . && git reset --hard
-python -c \
-"import re; import os; import shutil;
-for f in [os.path.join('${DIRNAME}/build/html', f) for f in os.listdir('${DIRNAME}/build/html')]:
-    if not (os.path.isdir(f) and re.search('v([0-9].)+[0-9]$', f)):
-        shutil.rmtree(f) if os.path.isdir(f) else os.remove(f)"
+python3 "${DIRNAME}/process_old_docs.py"
 
-# Install all released versions of Merlion/ts_datasets and use them to build the appropriate API docs.
+# Install all released versions of Merlion/ts_datasets _not_ on gh-pages and use them to build the appropriate API docs.
 # Uninstall after we're done with each one.
 versions=("latest")
 for v in $(git tag --list 'v[0-9]*'); do

docs/process_old_docs.py

@@ -0,0 +1,75 @@
+#
+# Copyright (c) 2022 salesforce.com, inc.
+# All rights reserved.
+# SPDX-License-Identifier: BSD-3-Clause
+# For full license text, see the LICENSE file in the repo root or https://opensource.org/licenses/BSD-3-Clause
+#
+"""
+Script which removes redirects from the HTML API docs & updates the version matrix on old files.
+"""
+import os
+import re
+import shutil
+
+from bs4 import BeautifulSoup as bs
+from git import Repo
+
+
+def create_version_dl(soup, prefix, current_version, all_versions):
+    dl = soup.new_tag("dl")
+    dt = soup.new_tag("dt")
+    dt.string = "Versions"
+    dl.append(dt)
+    for version in all_versions:
+        # Create the href for this version & bold it if it's the current version
+        href = soup.new_tag("a", href=f"{prefix}/{version}/index.html")
+        href.string = version
+        if version == current_version:
+            strong = soup.new_tag("strong")
+            strong.append(href)
+            href = strong
+        # Create a list item & add it to the dl
+        dd = soup.new_tag("dd")
+        dd.append(href)
+        dl.append(dd)
+    return dl
+
+
+def main():
+    # Get all the versions
+    repo = Repo(search_parent_directories=True)
+    versions = sorted([tag.name for tag in repo.tags if re.match("v[0-9].*", tag.name)], reverse=True)
+    versions = ["latest", *versions]
+
+    dirname = os.path.join(os.path.dirname(os.path.abspath(__file__)), "build", "html")
+    for version in os.listdir(dirname):
+        # If this isn't a directory containing a numbered version's API docs, delete it
+        version_root = os.path.join(dirname, version)
+        if version == "latest" or version not in versions:
+            shutil.rmtree(version_root) if os.path.isdir(version_root) else os.remove(version_root)
+            continue
+
+        # Update version matrix in HTML source versioned files
+        for subdir, _, files in os.walk(version_root):
+            html_files = [os.path.join(subdir, f) for f in files if f.endswith(".html")]
+
+            # Determine how far the version root is from the files in this directory
+            prefix = ".."
+            while subdir and subdir != version_root:
+                subdir = os.path.dirname(subdir)
+                prefix += "/.."
+
+            # Create the new description list for the version & write the new file
+            for file in html_files:
+                with open(file) as f:
+                    soup = bs(f, "html.parser")
+                version_dl = [dl for dl in soup.find_all("dl") if dl.find("dt", text="Versions")]
+                if len(version_dl) == 0:
+                    continue
+                version_dl[0].replace_with(create_version_dl(soup, prefix, version, versions))
+                with open(file, "w", encoding="utf-8") as f:
+                    f.write(str(soup))
+
+
+if __name__ == "__main__":
+    main()

docs/requirements.txt

@@ -1,6 +1,7 @@
 GitPython
+beautifulsoup4
 ipykernel
-nbsphinx==0.8.7
+nbsphinx
 pandoc
 sphinx
 sphinx_autodoc_typehints

docs/source/conf.py

@@ -36,7 +36,13 @@
 # Add any Sphinx extension module names here, as strings. They can be
 # extensions coming with Sphinx (named 'sphinx.ext.*') or your custom
 # ones.
-extensions = ["nbsphinx", "sphinx.ext.autodoc", "sphinx.ext.autosummary", "sphinx_autodoc_typehints"]
+extensions = [
+    "nbsphinx",
+    "IPython.sphinxext.ipython_console_highlighting",
+    "sphinx.ext.autodoc",
+    "sphinx.ext.autosummary",
+    "sphinx_autodoc_typehints",
+]
 
 autoclass_content = "both"  # include both class docstring and __init__
 autodoc_default_options = {
@@ -91,3 +97,4 @@
     exclude_patterns = ["examples"]
 else:
     exclude_patterns = ["tutorials"]
+exclude_patterns += ["**.ipynb_checkpoints"]

merlion/models/base.py

@@ -24,7 +24,7 @@
 from merlion.transform.factory import TransformFactory
 from merlion.transform.normalize import Rescale, MeanVarNormalize
 from merlion.transform.sequence import TransformSequence
-from merlion.utils.time_series import assert_equal_timedeltas, to_pd_datetime, TimeSeries
+from merlion.utils.time_series import assert_equal_timedeltas, to_pd_datetime, infer_granularity, TimeSeries
 from merlion.utils.misc import AutodocABCMeta, ModelConfigMeta
 
 logger = logging.getLogger(__name__)
@@ -169,6 +169,7 @@ def __init__(self, config: Config):
         self.config = copy.copy(config)
         self.last_train_time = None
         self.timedelta = None
+        self.timedelta_offset = pd.to_timedelta(0)
         self.train_data = None
 
     def reset(self):
@@ -304,12 +305,10 @@ def train_pre_process(self, train_data: TimeSeries) -> TimeSeries:
 
         # Make sure timestamps are equally spaced if needed (e.g. for ARIMA)
         t = train_data.time_stamps
+        self.timedelta, self.timedelta_offset = infer_granularity(t, return_offset=True)
         if self.require_even_sampling:
-            assert_equal_timedeltas(train_data.univariates[train_data.names[0]])
+            assert_equal_timedeltas(train_data.univariates[train_data.names[0]], self.timedelta, self.timedelta_offset)
             assert train_data.is_aligned
-            self.timedelta = pd.infer_freq(to_pd_datetime(t))
-        else:
-            self.timedelta = t[1] - t[0]
         self.last_train_time = t[-1]
         return train_data.align() if self.auto_align else train_data
 

merlion/models/forecast/base.py

@@ -114,7 +114,7 @@ def resample_time_stamps(self, time_stamps: Union[int, List[int]], time_series_p
         )
 
         # Determine timedelta & initial time of forecast
-        dt = self.timedelta
+        dt, offset = self.timedelta, self.timedelta_offset
         if time_series_prev is not None and not time_series_prev.is_empty():
             t0 = to_pd_datetime(time_series_prev.tf)
         else:
@@ -124,34 +124,33 @@ def resample_time_stamps(self, time_stamps: Union[int, List[int]], time_series_p
         if isinstance(time_stamps, (int, float)):
             n = int(time_stamps)
             assert self.max_forecast_steps is None or n <= self.max_forecast_steps
-            resampled = pd.date_range(start=t0, periods=n + 1, freq=dt)[1:]
-            tf = resampled[-1]
+            resampled = pd.date_range(start=t0, periods=n + 1, freq=dt) + offset
+            resampled = resampled[1:] if resampled[0] == t0 else resampled[:-1]
             time_stamps = to_timestamp(resampled)
 
         elif not self.require_even_sampling:
             resampled = to_pd_datetime(time_stamps)
-            tf = resampled[-1]
 
         # Handle the cases where we don't have a max_forecast_steps
         elif self.max_forecast_steps is None:
             tf = to_pd_datetime(time_stamps[-1])
-            resampled = pd.date_range(start=t0, end=tf, freq=dt)[1:]
-            if resampled[-1] < tf:
-                extra = pd.date_range(start=resampled[-1], periods=2, freq=dt)[1:]
-                resampled = resampled.union(extra)
+            resampled = pd.date_range(start=t0, end=tf + 2 * dt, freq=dt) + offset
+            if resampled[0] == t0:
+                resampled = resampled[1:]
+            if len(resampled) > 1 and resampled[-2] >= tf:
+                resampled = resampled[:-1]
 
         # Handle the case where we do have a max_forecast_steps
         else:
-            resampled = pd.date_range(start=t0, periods=self.max_forecast_steps + 1, freq=dt)[1:]
-            tf = resampled[-1]
-            n = sum(t < to_pd_datetime(time_stamps[-1]) for t in resampled)
-            resampled = resampled[: n + 1]
+            resampled = pd.date_range(start=t0, periods=self.max_forecast_steps + 1, freq=dt) + offset
+            resampled = resampled[1:] if resampled[0] == t0 else resampled[:-1]
+            resampled = resampled[: 1 + sum(resampled < to_pd_datetime(time_stamps[-1]))]
 
+        tf = resampled[-1]
         assert to_pd_datetime(time_stamps[0]) >= t0 and to_pd_datetime(time_stamps[-1]) <= tf, (
             f"Expected `time_stamps` to be between {t0} and {tf}, but `time_stamps` ranges "
             f"from {to_pd_datetime(time_stamps[0])} to {to_pd_datetime(time_stamps[-1])}"
         )
-
         return to_timestamp(resampled).tolist()
 
     def train_pre_process(

merlion/models/forecast/ets.py

@@ -81,6 +81,7 @@ def __init__(self, config: ETSConfig):
         super().__init__(config)
         self.model = None
         self._last_val = None
+        self._n_train = None
 
     @property
     def require_even_sampling(self) -> bool:
@@ -136,12 +137,13 @@ def _train(self, train_data: pd.DataFrame, train_config=None):
         name = self.target_name
         train_data = train_data[name]
         times = train_data.index
-        self.model = self._instantiate_model(train_data).fit(disp=False)
+        self.model = self._instantiate_model(pd.Series(train_data.values)).fit(disp=False)
 
         # get forecast for the training data
         self._last_val = train_data[-1]
-        yhat = pd.DataFrame(self.model.fittedvalues.values.tolist(), index=times, columns=[name])
-        err = pd.DataFrame(self.model.standardized_forecasts_error.tolist(), index=times, columns=[f"{name}_err"])
+        self._n_train = len(train_data)
+        yhat = pd.DataFrame(self.model.fittedvalues.values, index=times, columns=[name])
+        err = pd.DataFrame(self.model.standardized_forecasts_error, index=times, columns=[f"{name}_err"])
         return yhat, err
 
     def _forecast(
@@ -152,10 +154,12 @@ def _forecast(
         if time_series_prev is None:
             last_val = self._last_val
             model = self.model
+            start = self._n_train
         else:
             time_series_prev = time_series_prev.iloc[:, self.target_seq_index]
-            val_prev = time_series_prev[-self._max_lookback :]
-            last_val = val_prev[-1]
+            val_prev = pd.Series(time_series_prev[-self._max_lookback :].values)
+            last_val = val_prev.iloc[-1]
+            start = len(val_prev)
 
             # the default setting of refit=False is fast and conducts exponential smoothing with given parameters,
             # while the setting of refit=True is slow and refits the model on time_series_prev.
@@ -165,16 +169,10 @@ def _forecast(
             else:
                 model = model.smooth(params=self.model.params)
 
-        # Run forecasting. Some variants of ETS model does not support prediction interval.
-        # In this case we use point forecasting and set prediction_interval as None.
-        try:
-            forecast_result = model.get_prediction(start=time_stamps[0], end=time_stamps[-1])
-            forecast = np.asarray(forecast_result.predicted_mean)
-            err = np.sqrt(np.asarray(forecast_result.var_pred_mean))
-        except (NotImplementedError, AttributeError):
-            forecast_result = model.predict(start=time_stamps[0], end=time_stamps[-1])
-            forecast = np.asarray(forecast_result)
-            err = None
+        # Run forecasting.
+        forecast_result = model.get_prediction(start=start, end=start + len(time_stamps) - 1)
+        forecast = np.asarray(forecast_result.predicted_mean)
+        err = np.sqrt(np.asarray(forecast_result.var_pred_mean))
 
         # If return_prev is True, return the forecast and error of last train window instead of time_series_prev
         if time_series_prev is not None and return_prev:

merlion/models/forecast/prophet.py

@@ -189,12 +189,6 @@ def set_seasonality(self, theta, train_data: UnivariateTimeSeries):
                 logger.debug(f"Add seasonality {str(p)} ({p * dt})")
                 self.model.add_seasonality(name=f"extra_season_{p}", period=period, fourier_order=p)
 
-    def resample_time_stamps(self, time_stamps: Union[int, List[int]], time_series_prev: TimeSeries = None):
-        if isinstance(time_stamps, (int, float)):
-            times = pd.date_range(start=self.last_train_time, freq=self.timedelta, periods=int(time_stamps + 1))[1:]
-            time_stamps = to_timestamp(times)
-        return time_stamps
-
     def _add_exog_data(self, data: pd.DataFrame, exog_data: pd.DataFrame):
         df = pd.DataFrame(data[self.target_name].rename("y"))
         if exog_data is not None:

merlion/models/forecast/smoother.py

@@ -293,7 +293,7 @@ def update(
             )
 
         new_data = TimeSeries.from_pd(new_data).univariates[name]
-        assert_equal_timedeltas(new_data, self.timedelta)
+        assert_equal_timedeltas(new_data, self.timedelta, self.timedelta_offset)
         next_train_time = self.last_train_time + self.timedelta
         if to_pd_datetime(new_data.t0) > next_train_time:
             logger.warning(