From 557aee9e2166a6f33b6a4c2e5cc4ce1f08901d41 Mon Sep 17 00:00:00 2001 From: "cyrille.leroux@gmail.com" Date: Wed, 30 Sep 2026 10:21:14 +0200 Subject: [PATCH 1/5] ENH gifteval: replicate the GIFT-Eval test protocol Rolling test windows follow the leaderboard (n_windows='gifteval': min(20, ceil(0.1 * min_len / pred_len)), 1 for m4_*), seasonality comes from the gluonts table the leaderboard uses, and NaN targets are imputed like the official notebooks with their original locations kept in nan_mask for the metrics to exclude them. --- benchmark_utils/forecasting_constants.py | 36 ++++++++++ datasets/gifteval.py | 87 +++++++++++++++++++++--- 2 files changed, 113 insertions(+), 10 deletions(-) diff --git a/benchmark_utils/forecasting_constants.py b/benchmark_utils/forecasting_constants.py index de64412..e019fdc 100644 --- a/benchmark_utils/forecasting_constants.py +++ b/benchmark_utils/forecasting_constants.py @@ -89,6 +89,42 @@ def from_pandas(freq_alias: str) -> tuple[str, int, int]: return freq_alias, max(1, seasonality // max(mult, 1)), default_h +# gluonts ``DEFAULT_SEASONALITIES`` (time_feature/seasonality.py): the +# table the GIFT-Eval leaderboard uses for MASE and SeasonalNaive. +# Differs from ``_BASE`` on D (1 vs 7), W (1 vs 52) and S (3600 vs 1). +_GLUONTS_SEASONALITIES = { + "S": 3600, + "T": 1440, + "H": 24, + "D": 1, + "W": 1, + "M": 12, + "Q": 4, + "Y": 1, +} + + +def gift_eval_seasonality(freq_alias: str) -> int: + """Seasonality used by the GIFT-Eval leaderboard (gluonts table). + + Mirrors ``gluonts.time_feature.get_seasonality``: look up the base + alias (anchors and modern pandas spellings normalized), divide by + the multiplier, fall back to 1 when the division is not exact or + the alias is unknown. + """ + if not freq_alias: + return 1 + m = _PANDAS_ALIAS_RE.match(freq_alias.split("-", 1)[0]) + if not m: + return 1 + mult = int(m.group(1)) if m.group(1) else 1 + base = _NORMALIZE_BASE.get(m.group(2), m.group(2)[:1].upper()) + if base not in _GLUONTS_SEASONALITIES: + return 1 + seasonality, remainder = divmod(_GLUONTS_SEASONALITIES[base], mult) + return seasonality if not remainder else 1 + + # --------------------------------------------------------------------------- # GIFT-Eval term resolution # diff --git a/datasets/gifteval.py b/datasets/gifteval.py index 0b3dc24..f75247d 100644 --- a/datasets/gifteval.py +++ b/datasets/gifteval.py @@ -71,9 +71,10 @@ Cutoffs and windows ------------------- -We don't comply with GIFT-Eval's prescribed test cutoff; we use the same -rolling-window logic as Monash via -:func:`benchmark_utils.windowing.make_forecasting_splits`. The +We replicate GIFT-Eval's prescribed test protocol: per-series windows +``W = min(20, ceil(0.1 * min_series_length / prediction_length))`` +(``m4_*`` use a single window), with non-overlapping windows ending +at the end of each series. The ``prediction_length`` for a given (freq, term) follows GIFT-Eval's canonical ``base × multiplier`` rule via :func:`benchmark_utils.forecasting_constants.gift_eval_prediction_length`. @@ -82,6 +83,7 @@ """ import csv +import math from pathlib import Path import numpy as np @@ -91,12 +93,18 @@ from benchmark_utils.download_hf import snapshot_hf_files from benchmark_utils.forecasting_constants import ( FORECASTING_METRICS, - from_pandas, + gift_eval_seasonality, gift_eval_prediction_length, ) from benchmark_utils.windowing import build_forecasting_data +# Window count of the GIFT-Eval test protocol (gift_eval.data: +# TEST_SPLIT / MAX_WINDOW). +GIFT_EVAL_TEST_SPLIT = 0.1 +GIFT_EVAL_MAX_WINDOW = 20 + + # --------------------------------------------------------------------------- # Canonical (dataset_name, term) table — derived from the GIFT-Eval # leaderboard Space (results/Seasonal_Naive/all_results.csv: 55 paths, @@ -205,6 +213,36 @@ def _skip_placeholder(reason: str) -> dict: return dict(_skip_reason=reason) +def _impute_nan(series: np.ndarray) -> np.ndarray: + """Impute NaN targets the way the official GIFT-Eval notebooks do + (gluonts LastValueImputation): forward-fill each channel, flat-fill + leading NaNs with the first valid value, 0.0 for all-NaN channels. + + 9 GIFT-Eval datasets ship NaN targets; left as is, they propagate + through solver histories (``hist[-s:]``) and metric reductions, + silently NaN-ing the scores. The original NaN locations are + returned alongside the data as ``nan_mask`` for the metrics to + exclude the imputed points. + """ + # Arrow columns can decode to zero-copy read-only views: copy + # before writing in place. + if not series.flags.writeable and np.isnan(series).any(): + series = series.copy() + for c in range(series.shape[1]): + col = series[:, c] + nan = np.isnan(col) + if not nan.any() or nan.all(): + if nan.all(): + col[:] = 0.0 + continue + idx = np.where(nan, 0, np.arange(len(col))) + np.maximum.accumulate(idx, out=idx) + col[nan] = col[idx[nan]] + first = int(np.argmax(~nan)) + col[:first] = col[first] + return series + + class Dataset(BaseDataset): """GIFT-Eval forecasting dataset (loaded from HF Salesforce/GiftEval). @@ -222,8 +260,11 @@ class Dataset(BaseDataset): prediction_length : int or None Explicit override. ``None`` → resolved from (freq, term) via :func:`benchmark_utils.forecasting_constants.gift_eval_prediction_length`. - n_windows : int - Number of rolling evaluation windows per series. + n_windows : int or "gifteval" + Number of rolling evaluation windows per series. The default + ``"gifteval"`` follows the GIFT-Eval protocol: + ``min(20, ceil(0.1 * min_series_length / pred_len))``, pinned + to 1 for ``m4_*`` datasets. max_series : int or None Optional cap on the number of series. debug : bool @@ -238,7 +279,7 @@ class Dataset(BaseDataset): "dataset_name": ["m4_weekly/W"], "term": ["short"], "prediction_length": [None], - "n_windows": [1], + "n_windows": ["gifteval"], "max_series": [None], "debug": [False], } @@ -277,6 +318,19 @@ def _snapshot(self) -> "list[str]": "data-*.arrow", ) + def _resolve_n_windows(self, series_list, pred_len): + """Number of windows for the ``"gifteval"`` sentinel; ints pass + through unchanged.""" + if self.n_windows != "gifteval": + return int(self.n_windows) + if "m4" in self.dataset_name: + return 1 + min_len = min(ts.shape[0] for ts in series_list) + return min( + GIFT_EVAL_MAX_WINDOW, + max(1, math.ceil(GIFT_EVAL_TEST_SPLIT * min_len / pred_len)), + ) + def get_data(self): # Short-circuit non-canonical combos so heavy parsing doesn't run. if self.term not in _leaderboard().get(self.dataset_name, ()): @@ -310,7 +364,7 @@ def get_data(self): # Frequency / seasonality — every series in a GIFT-Eval subset # shares the same freq, so taking it from the first entry is safe. pandas_freq = ds[0].get("freq") or "D" - freq, seasonality, _ = from_pandas(pandas_freq) + seasonality = gift_eval_seasonality(pandas_freq) pred_len = self.prediction_length if pred_len is None: @@ -335,17 +389,30 @@ def get_data(self): "had unsupported target shapes." ) + # Original NaN locations, for the metrics to exclude the + # imputed points. build_forecasting_data drops series shorter + # than pred_len + 1: filter them out here to keep the lists + # aligned. + nan_mask = [ + np.isnan(ts).any(axis=1) + for ts in series_list + if ts.shape[0] >= pred_len + 1 + ] + series_list = [_impute_nan(ts) for ts in series_list] + n_windows = self._resolve_n_windows(series_list, pred_len) + return dict( **build_forecasting_data( series_list, prediction_length=pred_len, - n_windows=self.n_windows, + n_windows=n_windows, debug=self.debug, ), covariates=Covariates(), # GIFT-Eval HF schema has no covariates task="forecasting", metrics=list(FORECASTING_METRICS), prediction_length=pred_len, - freq=freq, + freq=pandas_freq, seasonality=seasonality, + nan_mask=nan_mask, ) From 380801022fec0b1df0e83912df035a77b008bc96 Mon Sep 17 00:00:00 2001 From: "cyrille.leroux@gmail.com" Date: Wed, 30 Sep 2026 10:21:29 +0200 Subject: [PATCH 2/5] ENH metrics: gluonts MASE/sMAPE semantics, mask imputed targets MASE now scales each window by the seasonal-naive MAE of the context available at its cutoff (gluonts evaluate_forecasts aggregation) and sMAPE uses the gluonts formula 2|e|/(|y|+|yhat|) with 0/0 dropped. All forecasting metrics accept valid_mask to score the non-imputed points only; the objective builds it from the dataset nan_mask. --- benchmark_utils/metrics.py | 174 +++++++++++++++++++++++++++++++------ objective.py | 15 ++++ 2 files changed, 162 insertions(+), 27 deletions(-) diff --git a/benchmark_utils/metrics.py b/benchmark_utils/metrics.py index 4e5b8cd..3da6547 100644 --- a/benchmark_utils/metrics.py +++ b/benchmark_utils/metrics.py @@ -14,6 +14,8 @@ event_detection : metric(y_true, y_pred, **kw) -> float """ +import warnings + import numpy as np from sklearn.metrics import ( accuracy_score, @@ -58,12 +60,54 @@ def _seasonal_naive_scale(y_train, seasonality: int) -> float: return scale if scale != 0 else 1.0 -def _pinball_per_level(y_true: np.ndarray, forecast: ForecastOutput) -> np.ndarray: +def _seasonal_naive_scales_per_window( + histories, cutoff_indexes, seasonality: int, nan_masks=None +) -> np.ndarray: + """Seasonal-naive scale per window and channel, following + ``gluonts.ev.ts_stats.seasonal_error`` on the context available at + each cutoff (with its fall-back to seasonality 1 when the context + is shorter than the season). ``nan_masks``, when given, holds the + NaN locations of the original (pre-imputation) targets: diffs + touching an imputed point are excluded. Returns an (M, C) array in + the flattened series-major / cutoff-minor window order. + """ + scales = [] + if nan_masks is None: + nan_masks = [None] * len(histories) + for series, cutoffs, nan_mask in zip(histories, cutoff_indexes, nan_masks): + series = np.asarray(series) + for cutoff in cutoffs: + ctx = series[:cutoff] + s = seasonality if seasonality < ctx.shape[0] else 1 + diff = np.abs(ctx[s:] - ctx[:-s]) + if nan_mask is not None: + nm = nan_mask[:cutoff] + valid = ~(nm[s:] | nm[:-s]) + diff = np.where(valid[:, None], diff, np.nan) + # a channel with no valid diff gets a NaN scale, which + # drops its scaled errors downstream + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) + scale = np.nanmean(diff, axis=0) + else: + scale = diff.mean(axis=0) + # a constant context has a zero scale: the scaled errors are + # undefined, drop them as well + scale = np.where(scale == 0.0, np.nan, scale) + scales.append(scale) + return np.stack(scales, axis=0) + + +def _pinball_per_level( + y_true: np.ndarray, forecast: ForecastOutput, valid_mask=None +) -> np.ndarray: """Pinball loss array of shape (Q,): mean over (M, H, C) for each level.""" quants, levels = _stacked(forecast) # (M,H,C,Q), (Q,) diff = y_true[..., None] - quants # (M,H,C,Q) levels_b = levels.reshape(1, 1, 1, -1) loss = np.maximum(levels_b * diff, (levels_b - 1.0) * diff) + if valid_mask is not None: + return loss[valid_mask].mean(axis=0) # (Q,) return loss.mean(axis=(0, 1, 2)) # (Q,) @@ -72,32 +116,102 @@ def _pinball_per_level(y_true: np.ndarray, forecast: ForecastOutput) -> np.ndarr # --------------------------------------------------------------------------- -def mae(y_true, forecast: ForecastOutput, **_): +def mae(y_true, forecast: ForecastOutput, valid_mask=None, **_): """Mean Absolute Error, averaged over all windows, horizons, channels.""" - return float(np.mean(np.abs(y_true - _point_from_forecast(forecast)))) + err = np.abs(y_true - _point_from_forecast(forecast)) + return _masked_mean(err, valid_mask) -def mse(y_true, forecast: ForecastOutput, **_): +def mse(y_true, forecast: ForecastOutput, valid_mask=None, **_): """Mean Squared Error, averaged over all windows, horizons, channels.""" - return float(np.mean((y_true - _point_from_forecast(forecast)) ** 2)) + err = (y_true - _point_from_forecast(forecast)) ** 2 + return _masked_mean(err, valid_mask) + + +def rmse(y_true, forecast: ForecastOutput, valid_mask=None, **_): + return float(np.sqrt(mse(y_true, forecast, valid_mask))) + + +def mase( + y_true, + forecast: ForecastOutput, + y_train=None, + seasonality=1, + histories=None, + cutoff_indexes=None, + nan_masks=None, + valid_mask=None, + **_, +): + """Mean Absolute Scaled Error, aggregated like + ``gluonts.evaluate_forecasts(..., axis=None)``. + + The scale of each window is the seasonal-naive MAE of the context + seen at its cutoff, and the score is the mean of the per-point + scaled errors (not the ratio of two aggregated means). Requires + ``histories`` and ``cutoff_indexes``; without them, falls back to a + single scale computed on ``y_train``. + + ``nan_masks`` holds the NaN locations of the original + (pre-imputation) targets: diffs touching an imputed point are + excluded from the scale, and ``valid_mask`` (M, H, C) excludes the + imputed targets themselves. + """ + y_pred = _point_from_forecast(forecast) + err = np.abs(y_true - y_pred) + if histories is not None and cutoff_indexes is not None: + scales = _seasonal_naive_scales_per_window( + histories, cutoff_indexes, seasonality, nan_masks + ) + scaled = err / scales[:, None, :] + if valid_mask is None: + valid_mask = np.ones(err.shape, dtype=bool) + # windows with an undefined scale (NaN) are dropped + valid_mask = valid_mask & ~np.isnan(scales)[:, None, :] + with np.errstate(divide="ignore", invalid="ignore"): + return float(np.mean(scaled[valid_mask])) + scale = _seasonal_naive_scale(y_train, seasonality) + return float(np.mean(err) / scale) -def rmse(y_true, forecast: ForecastOutput, **_): - return float(np.sqrt(mse(y_true, forecast))) +def _masked_mean(values: np.ndarray, valid_mask) -> float: + """Mean over the valid points only when ``valid_mask`` is given.""" + if valid_mask is None: + return float(np.mean(values)) + return float(np.mean(values[valid_mask])) -def mase(y_true, forecast: ForecastOutput, y_train, seasonality=1, **_): - """Mean Absolute Scaled Error. Scale is naive-seasonal MAE on y_train.""" - scale = _seasonal_naive_scale(y_train, seasonality) - return float(np.mean(np.abs(y_true - _point_from_forecast(forecast))) / scale) +def windows_valid_mask(nan_masks, cutoff_indexes, target_shape) -> np.ndarray: + """Validity mask of the (M, H, C) target array: False on the points + that were imputed in the original targets. Takes the per-series NaN + masks and the jagged cutoffs, in the flattened series-major / + cutoff-minor window order; the per-step validity broadcasts over + channels.""" + valid = np.ones(target_shape[:2], dtype=bool) + i = 0 + for nan_mask, cutoffs in zip(nan_masks, cutoff_indexes): + for cutoff in cutoffs: + valid[i] = ~nan_mask[cutoff : cutoff + target_shape[1]] + i += 1 + mask = np.empty(target_shape, dtype=bool) + mask[:] = valid[:, :, None] + return mask -def smape(y_true, forecast: ForecastOutput, **_): - """Symmetric Mean Absolute Percentage Error.""" +def smape(y_true, forecast: ForecastOutput, valid_mask=None, **_): + """Symmetric Mean Absolute Percentage Error, gluonts formula: + ``2|y - yhat| / (|y| + |yhat|)``. Undefined (0/0) points are + dropped from the mean.""" y_pred = _point_from_forecast(forecast) - denom = (np.abs(y_true) + np.abs(y_pred)) / 2.0 - denom = np.where(denom == 0, 1.0, denom) - return float(np.mean(np.abs(y_true - y_pred) / denom)) + denom = np.abs(y_true) + np.abs(y_pred) + with np.errstate(divide="ignore", invalid="ignore"): + per_point = 2 * np.abs(y_true - y_pred) / denom + if valid_mask is not None: + per_point = per_point[valid_mask] + valid = per_point[~np.isnan(per_point)] + if valid.size == 0: + return float("nan") + return float(valid.mean()) def skill_score_ratio(y_true, forecast: ForecastOutput, y_train, seasonality=1, **_): @@ -117,35 +231,41 @@ def skill_score_ratio(y_true, forecast: ForecastOutput, y_train, seasonality=1, # --------------------------------------------------------------------------- -def pinball(y_true, forecast: ForecastOutput, **_): +def pinball(y_true, forecast: ForecastOutput, valid_mask=None, **_): """Mean pinball (quantile) loss, averaged over all quantile levels.""" - return float(_pinball_per_level(y_true, forecast).mean()) + return float(_pinball_per_level(y_true, forecast, valid_mask).mean()) -def crps(y_true, forecast: ForecastOutput, **_): +def crps(y_true, forecast: ForecastOutput, valid_mask=None, **_): """CRPS approximated by the quantile-score formula: 2 * mean pinball loss. Converges to the true CRPS as the quantile grid becomes dense. """ - return 2.0 * pinball(y_true, forecast) + return 2.0 * pinball(y_true, forecast, valid_mask) -def wql(y_true, forecast: ForecastOutput, **_): +def wql(y_true, forecast: ForecastOutput, valid_mask=None, **_): """Weighted Quantile Loss (Salinas et al. / Chronos). ``WQL = (1/Q) sum_q [ 2 * sum_t pinball_q(y_t) / sum_t |y_t| ]`` """ + if valid_mask is not None: + y_valid = y_true[valid_mask] + denom = float(np.sum(np.abs(y_valid))) + else: + y_valid = y_true + denom = float(np.sum(np.abs(y_true))) _, levels = _stacked(forecast) - denom = float(np.sum(np.abs(y_true))) if denom == 0: return float("nan") - # _pinball_per_level returns per-level mean; multiply back by N to get sum. - n_elem = float(y_true.size) - per_level_sum = _pinball_per_level(y_true, forecast) * n_elem + n_elem = float(y_valid.size) + # _pinball_per_level returns a mean: multiply by N to get the sum + per_level = _pinball_per_level(y_true, forecast, valid_mask) + per_level_sum = per_level * n_elem return float(2.0 * per_level_sum.sum() / (len(levels) * denom)) -def mcis(y_true, forecast: ForecastOutput, alpha=0.05, **_): +def mcis(y_true, forecast: ForecastOutput, alpha=0.05, valid_mask=None, **_): """Mean Coverage Interval Score for a ``(1 - alpha)`` prediction interval. ``IS_alpha = (U - L) + (2/alpha)(L - y) 1[y < L] + (2/alpha)(y - U) 1[y > U]`` @@ -163,7 +283,7 @@ def mcis(y_true, forecast: ForecastOutput, alpha=0.05, **_): under = np.maximum(0.0, lower - y_true) over = np.maximum(0.0, y_true - upper) score = (upper - lower) + (2.0 / alpha) * (under + over) - return float(np.mean(score)) + return _masked_mean(score, valid_mask) # --------------------------------------------------------------------------- diff --git a/objective.py b/objective.py index 5da735c..242ee07 100644 --- a/objective.py +++ b/objective.py @@ -184,10 +184,25 @@ def _eval_forecasting(self, model): # same order the flattened forecast iterates (series-major, cutoff-minor). y_true = np.concatenate([np.asarray(yt) for yt in self.y_test], axis=0) + # nan_mask carries the original NaN locations of imputed + # targets: metrics score the valid points only. + nan_masks = self.meta.get("nan_mask") + valid_mask = None + if nan_masks: + from benchmark_utils.metrics import windows_valid_mask + + valid_mask = windows_valid_mask( + nan_masks, self.cutoff_indexes, y_true.shape + ) + kwargs = dict( y_train=self.X_train, seasonality=self.meta.get("seasonality", 1), alpha=self.meta.get("mcis_alpha", 0.05), + histories=self.X_test, + cutoff_indexes=self.cutoff_indexes, + nan_masks=nan_masks, + valid_mask=valid_mask, ) result = { name: ALL_METRICS[name](y_true, forecast, **kwargs) for name in self.metrics From 78998642789ba6fb8239fa26ec27ee8ea8d47433 Mon Sep 17 00:00:00 2001 From: "cyrille.leroux@gmail.com" Date: Wed, 30 Sep 2026 10:21:52 +0200 Subject: [PATCH 3/5] ENH SeasonalNaive: auto season_length, tests for the protocol season_length='auto' resolves to meta['seasonality'] like the GIFT-Eval baseline, replacing the [1, 7, 12, 24] sweep. Tests cover the NaN imputation, the gluonts seasonality table, the per-window scales, the masked metrics and the window count resolution. --- solvers/seasonal_naive.py | 14 +- .../test_forecasting_constants.py | 30 ++++ .../test_metrics_forecasting.py | 128 ++++++++++++++++++ tests/datasets/test_gifteval.py | 74 +++++++++- 4 files changed, 234 insertions(+), 12 deletions(-) create mode 100644 tests/benchmark_utils/test_metrics_forecasting.py diff --git a/solvers/seasonal_naive.py b/solvers/seasonal_naive.py index d439975..33be74e 100644 --- a/solvers/seasonal_naive.py +++ b/solvers/seasonal_naive.py @@ -52,11 +52,10 @@ class Solver(BaseSolver): Parameters ---------- - season_length : int + season_length : int or "auto" Number of past steps to repeat. ``1`` recovers last-value - persistence; common picks are ``7`` (daily → weekly), ``12`` - (monthly → yearly), ``24`` (hourly → daily), ``52`` (weekly → - yearly). + persistence. ``"auto"`` (the default) uses the dataset's own + seasonality, ``meta["seasonality"]``. """ name = "SeasonalNaive" @@ -64,7 +63,7 @@ class Solver(BaseSolver): requirements = [] parameters = { - "season_length": [1, 7, 12, 24], + "season_length": ["auto"], } def skip(self, task, **kwargs): @@ -79,9 +78,12 @@ def set_objective(self, X_train, y_train, task, **meta): self.meta = meta def run(self, _): + season_length = self.season_length + if season_length == "auto": + season_length = self.meta.get("seasonality", 1) self._adapter = _SeasonalNaiveForecaster( prediction_length=self.meta.get("prediction_length", 1), - season_length=self.season_length, + season_length=season_length, ) def get_result(self): diff --git a/tests/benchmark_utils/test_forecasting_constants.py b/tests/benchmark_utils/test_forecasting_constants.py index b404e2d..592170d 100644 --- a/tests/benchmark_utils/test_forecasting_constants.py +++ b/tests/benchmark_utils/test_forecasting_constants.py @@ -6,9 +6,39 @@ from_aeon, from_pandas, gift_eval_prediction_length, + gift_eval_seasonality, ) +@pytest.mark.parametrize( + "alias, expected", + [ + # gluonts DEFAULT_SEASONALITIES: differs from _BASE on D, W and S + ("H", 24), + ("D", 1), + ("W", 1), + ("W-SUN", 1), + ("M", 12), + ("Q", 4), + ("Y", 1), + ("T", 1440), + ("S", 3600), + # multipliers divide the seasonality, non-exact division → 1 + ("15T", 96), + ("5T", 288), + ("6H", 4), + ("7T", 1), + ("10S", 360), + ("45T", 32), + # unknown or empty → 1 + ("", 1), + ("??", 1), + ], +) +def test_gift_eval_seasonality(alias, expected): + assert gift_eval_seasonality(alias) == expected + + @pytest.mark.parametrize( "alias, expected", [ diff --git a/tests/benchmark_utils/test_metrics_forecasting.py b/tests/benchmark_utils/test_metrics_forecasting.py new file mode 100644 index 0000000..8218164 --- /dev/null +++ b/tests/benchmark_utils/test_metrics_forecasting.py @@ -0,0 +1,128 @@ +"""Tests for the GIFT-Eval / gluonts forecasting metric semantics.""" + +import numpy as np +import pytest + +from benchmark_utils.metrics import ( + _seasonal_naive_scales_per_window, + mae, + mase, + smape, + windows_valid_mask, +) +from benchmark_utils.outputs import ForecastOutput + + +def _point_forecast(point): + """ForecastOutput with a degenerate fan around ``point`` (M, H, C).""" + point = np.asarray(point, dtype=np.float64) + quantiles = point[..., None] * np.ones((1,)) + return ForecastOutput(quantiles=[quantiles], quantile_levels=[0.5]) + + +# Two series, two channels each. A cutoff c means the context is +# series[:c]: the point at index c is the first forecast target. +_HISTORIES = [ + np.array([[1.0, 10.0], [2.0, 20.0], [3.0, 30.0], [4.0, 40.0], [5.0, 50.0]]), + np.array([[0.0, 1.0], [2.0, 3.0], [4.0, 1.0], [6.0, 3.0], [8.0, 1.0]]), +] +_CUTOFFS = [[4], [4]] + + +def test_scales_per_window_gluonts_semantics(): + # s0 diffs (s=2): |3-1|=2, |4-2|=2 -> 2.0; ch1: 20.0 + # s1 diffs (s=2): ch0: 4.0; ch1: 1,3,1,3 alternates at lag 2 -> 0.0, + # and a zero scale is *undefined* (dropped), hence NaN + scales = _seasonal_naive_scales_per_window(_HISTORIES, _CUTOFFS, 2) + np.testing.assert_allclose(scales, [[2.0, 20.0], [4.0, np.nan]]) + + +def test_scales_per_window_constant_context_is_nan(): + histories = [np.full((5, 1), 3.0)] + scales = _seasonal_naive_scales_per_window(histories, [[4]], 2) + assert np.isnan(scales).all() + + +def test_scales_per_window_short_context_falls_back_to_season_1(): + histories = [np.array([[1.0], [4.0]])] + scales = _seasonal_naive_scales_per_window(histories, [[2]], 12) + np.testing.assert_allclose(scales, [[3.0]]) + + +def test_scales_per_window_nan_mask_excludes_imputed_diffs(): + histories = [np.array([[1.0], [np.nan], [3.0], [5.0]])] + nan_masks = [np.array([False, True, False, False])] + # the NaN at position 1 kills the diff |3-1|: only |5-3| remains + scales = _seasonal_naive_scales_per_window(histories, [[4]], 1, nan_masks=nan_masks) + np.testing.assert_allclose(scales, [[2.0]]) + + +def test_mase_per_window_matches_gluonts_aggregation(): + # scales: s0 = (2, 20); s1 = (4, NaN) -> s1 ch1 points dropped + y_true = np.array([[[6.0, 60.0], [7.0, 70.0]], [[10.0, 1.0], [12.0, 3.0]]]) + forecast = _point_forecast(np.zeros_like(y_true)) + # scales: s0 = (2, 20); s1 = (4, NaN) -> s1 ch1 points dropped + value = mase( + y_true, + forecast, + histories=_HISTORIES, + cutoff_indexes=_CUTOFFS, + seasonality=2, + ) + expected = np.mean([6.0 / 2, 7.0 / 2, 60.0 / 20, 70.0 / 20, 10.0 / 4, 12.0 / 4]) + assert value == pytest.approx(expected) + + +def test_mase_legacy_fallback_without_histories(): + y_true = np.array([[[3.0], [5.0]]]) + forecast = _point_forecast(np.zeros_like(y_true)) + y_train = np.array([[1.0], [2.0]]) + # legacy scale: mean |diff season 1| of y_train = 1 -> (3 + 5) / 2 + assert mase(y_true, forecast, y_train=y_train, seasonality=1) == 4.0 + + +def test_mase_drops_undefined_scale_points(): + # s1 has a constant context: its scale is undefined, only the s0 + # points are scored + histories = [ + np.array([[1.0], [2.0], [3.0], [4.0]]), + np.array([[1.0], [1.0], [1.0], [1.0]]), + ] + cutoffs = [[4], [4]] + y_true = np.array([[[2.0]], [[5.0]]]) + forecast = _point_forecast(np.zeros_like(y_true)) + value = mase( + y_true, forecast, histories=histories, cutoff_indexes=cutoffs, seasonality=1 + ) + # s0 scale = 1 (mean |diff|), s1 scale undefined -> dropped + assert value == pytest.approx(2.0) + + +def test_windows_valid_mask(): + nan_masks = [np.array([False, False, False, True, False])] + cutoffs = [[3, 4]] + # H=1: window 0 covers index 3 (NaN -> invalid), window 1 covers 4 + mask = windows_valid_mask(nan_masks, cutoffs, (2, 1, 1)) + np.testing.assert_array_equal(mask, [[[False]], [[True]]]) + + +def test_mae_valid_mask(): + y_true = np.array([[[1.0], [2.0]]]) + forecast = _point_forecast(np.zeros_like(y_true)) + valid = np.array([[[True], [False]]]) + assert mae(y_true, forecast, valid_mask=valid) == pytest.approx(1.0) + + +def test_smape_gluonts_formula_drops_zero_zero(): + y_true = np.array([[[0.0], [2.0]]]) + forecast = _point_forecast(np.zeros_like(y_true)) + # 0/0 is undefined and dropped; the remaining point is 2*2/2 = 2 + assert smape(y_true, forecast) == pytest.approx(2.0) + + +def test_smape_valid_mask(): + y_true = np.array([[[0.0], [2.0]]]) + forecast = _point_forecast(np.zeros_like(y_true)) + valid = np.array([[[True], [False]]]) + # both points are 0/0 or masked: undefined -> NaN + assert np.isnan(smape(y_true, forecast, valid_mask=valid)) diff --git a/tests/datasets/test_gifteval.py b/tests/datasets/test_gifteval.py index e8cfef0..30fbe87 100644 --- a/tests/datasets/test_gifteval.py +++ b/tests/datasets/test_gifteval.py @@ -3,10 +3,12 @@ import inspect from pathlib import Path +import numpy as np +import pytest from benchopt.benchmark import Benchmark BENCHMARK_DIR = Path(__file__).parents[2] -Dataset, = Benchmark(BENCHMARK_DIR).check_dataset_patterns( +(Dataset,) = Benchmark(BENCHMARK_DIR).check_dataset_patterns( ["GiftEval"], class_only=True ) gifteval = inspect.getmodule(Dataset) @@ -29,9 +31,7 @@ def test_parse_leaderboard_csv(tmp_path): def test_non_canonical_combo_skips(monkeypatch): - monkeypatch.setattr( - gifteval, "_leaderboard_cache", {"m4_weekly/W": ("short",)} - ) + monkeypatch.setattr(gifteval, "_leaderboard_cache", {"m4_weekly/W": ("short",)}) ds = Dataset.get_instance(dataset_name="m4_weekly/W", term="long") data = ds.get_data() assert set(data) == {"_skip_reason"} @@ -41,6 +41,68 @@ def test_non_canonical_combo_skips(monkeypatch): def test_hf_arrow_directory(): assert gifteval._hf_arrow_directory("m4_weekly/W") == "m4_weekly" assert gifteval._hf_arrow_directory("loop_seattle/H") == "LOOP_SEATTLE/H" - assert ( - gifteval._hf_arrow_directory("car_parts/M") == "car_parts_with_missing" + assert gifteval._hf_arrow_directory("car_parts/M") == "car_parts_with_missing" + + +@pytest.mark.parametrize( + "series, expected", + [ + # forward-fill, leading flat-fill, all-NaN channel → 0.0 + ( + np.array([[np.nan], [1.0], [np.nan], [3.0]]), + np.array([[1.0], [1.0], [1.0], [3.0]]), + ), + ( + np.array([[np.nan], [np.nan], [2.0], [np.nan]]), + np.array([[2.0], [2.0], [2.0], [2.0]]), + ), + ( + np.array([[np.nan], [np.nan]]), + np.array([[0.0], [0.0]]), + ), + # multivariate: channels are filled independently + ( + np.array([[1.0, np.nan], [np.nan, 5.0], [3.0, np.nan]]), + np.array([[1.0, 5.0], [1.0, 5.0], [3.0, 5.0]]), + ), + ], +) +def test_impute_nan(series, expected): + out = gifteval._impute_nan(series.astype(np.float64)) + np.testing.assert_allclose(out, expected) + + +def test_impute_nan_read_only_view(): + # zero-copy arrow views are read-only: _impute_nan must copy + series = np.array([[np.nan], [1.0], [np.nan]]) + view = series[:] + view.flags.writeable = False + out = gifteval._impute_nan(view) + assert out is not view + np.testing.assert_allclose(out, [[1.0], [1.0], [1.0]]) + + +@pytest.mark.parametrize( + "min_len, pred_len, dataset_name, expected", + [ + (1000, 48, "electricity/15T", 3), + (1050, 48, "electricity/15T", 3), + (9600, 48, "electricity/15T", 20), + (100, 13, "m4_weekly/W", 1), + (5000, 48, "m4_hourly/H", 1), + # short series still get at least one window + (50, 48, "hospital/D", 1), + ], +) +def test_gifteval_n_windows(min_len, pred_len, dataset_name, expected): + ds = Dataset.get_instance( + dataset_name=dataset_name, term="short", n_windows="gifteval" ) + series_list = [np.zeros((min_len, 1)) for _ in range(3)] + n_windows = ds._resolve_n_windows(series_list, pred_len) + assert n_windows == expected + + +def test_gifteval_n_windows_int_passthrough(): + ds = Dataset.get_instance(dataset_name="m4_weekly/W", term="short", n_windows=4) + assert ds._resolve_n_windows([np.zeros((100, 1))], 13) == 4 From b382849191b198f443055028c0b7463f660908b2 Mon Sep 17 00:00:00 2001 From: "cyrille.leroux@gmail.com" Date: Tue, 6 Oct 2026 15:17:05 +0200 Subject: [PATCH 4/5] CHORE: comment cleanups --- benchmark_utils/metrics.py | 1 - objective.py | 3 +-- tests/benchmark_utils/test_metrics_forecasting.py | 1 - 3 files changed, 1 insertion(+), 4 deletions(-) diff --git a/benchmark_utils/metrics.py b/benchmark_utils/metrics.py index 3da6547..48a6270 100644 --- a/benchmark_utils/metrics.py +++ b/benchmark_utils/metrics.py @@ -175,7 +175,6 @@ def mase( def _masked_mean(values: np.ndarray, valid_mask) -> float: - """Mean over the valid points only when ``valid_mask`` is given.""" if valid_mask is None: return float(np.mean(values)) return float(np.mean(values[valid_mask])) diff --git a/objective.py b/objective.py index 242ee07..179af5f 100644 --- a/objective.py +++ b/objective.py @@ -55,6 +55,7 @@ ALL_METRICS, HIGHER_IS_BETTER, is_higher_better, + windows_valid_mask, ) @@ -189,8 +190,6 @@ def _eval_forecasting(self, model): nan_masks = self.meta.get("nan_mask") valid_mask = None if nan_masks: - from benchmark_utils.metrics import windows_valid_mask - valid_mask = windows_valid_mask( nan_masks, self.cutoff_indexes, y_true.shape ) diff --git a/tests/benchmark_utils/test_metrics_forecasting.py b/tests/benchmark_utils/test_metrics_forecasting.py index 8218164..0aac487 100644 --- a/tests/benchmark_utils/test_metrics_forecasting.py +++ b/tests/benchmark_utils/test_metrics_forecasting.py @@ -61,7 +61,6 @@ def test_mase_per_window_matches_gluonts_aggregation(): # scales: s0 = (2, 20); s1 = (4, NaN) -> s1 ch1 points dropped y_true = np.array([[[6.0, 60.0], [7.0, 70.0]], [[10.0, 1.0], [12.0, 3.0]]]) forecast = _point_forecast(np.zeros_like(y_true)) - # scales: s0 = (2, 20); s1 = (4, NaN) -> s1 ch1 points dropped value = mase( y_true, forecast, From 10f44ac379bc1ad3199d20f737773765a6363910 Mon Sep 17 00:00:00 2001 From: "cyrille.leroux@gmail.com" Date: Tue, 6 Oct 2026 15:18:06 +0200 Subject: [PATCH 5/5] ENH SeasonalNaive: native statsforecast quantile fan --- solvers/seasonal_naive.py | 61 ++++++++++++++++---- tests/solvers/test_seasonal_naive.py | 85 ++++++++++++++++++++++++++++ 2 files changed, 135 insertions(+), 11 deletions(-) create mode 100644 tests/solvers/test_seasonal_naive.py diff --git a/solvers/seasonal_naive.py b/solvers/seasonal_naive.py index 33be74e..084245f 100644 --- a/solvers/seasonal_naive.py +++ b/solvers/seasonal_naive.py @@ -1,17 +1,28 @@ """Seasonal-naive forecasting baseline. The forecast at horizon ``h`` is the value observed ``season_length`` steps -ago — i.e. for forecast index ``i`` (0-based within the horizon), the -prediction is ``hist[-season_length + (i mod season_length)]``. When the -available history is shorter than ``season_length``, the pattern falls -back to whatever history exists. +ago: for forecast index ``i`` (0-based within the horizon), the prediction +is ``hist[-season_length + (i mod season_length)]``. When the available +history is shorter than ``season_length``, the pattern falls back to +whatever history exists. A common, calibrated baseline for any dataset with a known seasonal period. With ``season_length=1`` it collapses to last-value persistence. + +Besides the 0.5 (median = point forecast), the adapter emits a 9-level +quantile fan (0.1 to 0.9) replicating the native prediction intervals of +``statsforecast.SeasonalNaive`` (the GIFT-Eval official baseline): per +channel, ``sigma`` is the RMS of the in-sample one-step seasonal +differences ``y[t] - y[t - season_length]``, the spread at step ``h`` is +``sigma * sqrt(k + 1)`` with ``k = floor(h / season_length)`` repeated +cycles, and level ``q`` sits at ``mean +/- norm.ppf(max(q, 1 - q)) * +sigmah``. The level 0.5 matches the native ``lo-0``/``hi-0`` columns, +i.e. the point itself. """ import numpy as np from benchopt import BaseSolver +from scipy.stats import norm from benchmark_utils.adapters.base import BaseTSFMAdapter from benchmark_utils.inputs import ForecastInput @@ -19,32 +30,60 @@ SUPPORTED_TASKS = {"forecasting"} +# GIFT-Eval leaderboard distributional metrics (gluonts +# MeanWeightedSumQuantileLoss) use a 0.1..0.9 fan. +QUANTILE_LEVELS = (0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9) + class _SeasonalNaiveForecaster(BaseTSFMAdapter): - """Repeat the last ``season_length`` observations to fill the horizon.""" + """Repeat the last ``season_length`` observations to fill the horizon + with a normal-quantile fan around each point.""" def __init__(self, prediction_length: int, season_length: int): if season_length < 1: raise ValueError(f"season_length must be >= 1, got {season_length}") self.prediction_length = prediction_length self.season_length = season_length + levels = np.asarray(QUANTILE_LEVELS) + # statsforecast native CI half-width: the interval of width + # round(200 * (max(q, 1-q) - 0.5)) percent around the mean is the + # normal quantile at max(q, 1-q); quantiles below 0.5 take the lo + # side, above 0.5 the hi side. + self._z = norm.ppf(np.maximum(levels, 1 - levels)) # (Q,) + self._sign = np.where(levels < 0.5, -1.0, 1.0)[None, None, :] # (1,1,Q) def predict(self, x: ForecastInput) -> ForecastOutput: quantiles = [] for series, cutoffs in zip(x.x, x.cutoff_indexes): series = np.asarray(series) C = series.shape[1] if series.ndim == 2 else 1 - preds = np.empty( - (len(cutoffs), self.prediction_length, C), dtype=np.float32 + fan = np.empty( + (len(cutoffs), self.prediction_length, C, len(QUANTILE_LEVELS)) ) for k, cutoff in enumerate(cutoffs): hist = series[:cutoff] - season = min(self.season_length, hist.shape[0]) + n = hist.shape[0] + season = min(self.season_length, n) pattern = hist[-season:] reps = int(np.ceil(self.prediction_length / season)) - preds[k] = np.tile(pattern, (reps, 1))[: self.prediction_length] - quantiles.append(preds[:, :, :, None]) # (n_cutoffs, H, C, 1) - return ForecastOutput(quantiles=quantiles, quantile_levels=(0.5,)) + point = np.tile(pattern, (reps, 1))[: self.prediction_length] + if n > self.season_length: + resid = hist[self.season_length :] - hist[: n - self.season_length] + sigma = np.sqrt( + np.nansum(resid * resid, axis=0) / (n - self.season_length) + ) + else: + sigma = np.zeros(C) + k_cycle = np.floor( + np.arange(self.prediction_length) / self.season_length + ) + sigmah = sigma[None, :] * np.sqrt(k_cycle + 1)[:, None] # (H, C) + fan[k] = ( + point[:, :, None] + + self._sign * self._z[None, None, :] * sigmah[:, :, None] + ) + quantiles.append(fan) + return ForecastOutput(quantiles=quantiles, quantile_levels=QUANTILE_LEVELS) class Solver(BaseSolver): diff --git a/tests/solvers/test_seasonal_naive.py b/tests/solvers/test_seasonal_naive.py new file mode 100644 index 0000000..52504e2 --- /dev/null +++ b/tests/solvers/test_seasonal_naive.py @@ -0,0 +1,85 @@ +"""Tests for the SeasonalNaive quantile fan (statsforecast-native intervals).""" + +import numpy as np +from scipy.stats import norm + +from benchmark_utils.inputs import ForecastInput +from solvers.seasonal_naive import QUANTILE_LEVELS, _SeasonalNaiveForecaster + + +def _fan(hist, prediction_length, season_length, cutoffs=None): + hist = np.asarray(hist, dtype=np.float64) + if hist.ndim == 1: + hist = hist[:, None] + cutoff_indexes = [[hist.shape[0]]] if cutoffs is None else [list(cutoffs)] + adapter = _SeasonalNaiveForecaster( + prediction_length=prediction_length, season_length=season_length + ) + out = adapter.predict(ForecastInput(x=[hist], cutoff_indexes=cutoff_indexes)) + assert out.quantile_levels == QUANTILE_LEVELS + return out.quantiles[0] # (n_cutoffs, H, C, Q) + + +Q_LO, Q_05, Q_HI = 0, 4, 8 + +# hist: [0, 1, 2, 3] with m=2 -> last season [2, 3], repeated over H=3 +_HIST = [0.0, 1.0, 2.0, 3.0] +Z10 = norm.ppf(0.9) # half-width of the 80% interval + + +def test_median_is_point_forecast(): + fan = _fan(_HIST, prediction_length=3, season_length=2) + np.testing.assert_allclose(fan[0, :, 0, Q_05], [2.0, 3.0, 2.0]) + + +def test_sigma_is_rms_seasonal_diff(): + # diffs at lag 2: (2-0)^2, (3-1)^2 -> sigma = sqrt(8/2) = 2 + # cycles over H=3: k = floor([0,1,2]/2) = [0,0,1] -> sigmah = [2, 2, 2*sqrt(2)] + fan = _fan(_HIST, prediction_length=3, season_length=2) + np.testing.assert_allclose( + fan[0, :, 0, Q_LO], + [2.0 - Z10 * 2, 3.0 - Z10 * 2, 2.0 - Z10 * 2 * np.sqrt(2)], + ) + + +def test_fan_levels_nested_around_point(): + fan = _fan(_HIST, prediction_length=3, season_length=2) + # lo(0.1) and hi(0.9) are symmetric around the point + np.testing.assert_allclose( + (fan[0, :, 0, Q_LO] + fan[0, :, 0, Q_HI]) / 2, fan[0, :, 0, Q_05] + ) + widths = fan[0, 0, 0, :] - fan[0, 0, 0, Q_05] + assert np.all(np.diff(widths[:5]) > 0) # lo side grows towards 0.5 + assert np.all(np.diff(widths[5:]) > 0) # hi side grows after 0.5 + + +def test_short_context_degenerates_to_point(): + # history shorter than the season: no seasonal diff at all -> sigma 0 + fan = _fan([1.0, 2.0], prediction_length=3, season_length=12) + np.testing.assert_allclose(fan[..., Q_LO], fan[..., Q_05]) + np.testing.assert_allclose(fan[..., Q_HI], fan[..., Q_05]) + np.testing.assert_allclose(fan[0, :, 0, Q_05], [1.0, 2.0, 1.0]) + + +def test_constant_channel_has_zero_sigma(): + hist = np.column_stack([np.arange(4.0), np.array([9.0, 9.0, 9.0, 9.0])]) + fan = _fan(hist, prediction_length=2, season_length=2) + np.testing.assert_allclose(fan[0, :, 1, Q_HI], [9.0, 9.0]) + # channel 0: sigma = 2 -> 0.9 = point + z(0.9) * 2 + np.testing.assert_allclose(fan[0, :, 0, Q_HI], [2.0 + Z10 * 2, 3.0 + Z10 * 2]) + + +def test_season_length_1_spread_grows_with_step(): + # [0, 1, 2, 3], m=1: sigma = sqrt(3/3) = 1, sigmah = 1*sqrt(1..H) + fan = _fan([0.0, 1.0, 2.0, 3.0], prediction_length=3, season_length=1) + np.testing.assert_allclose(fan[0, :, 0, Q_05], [3.0, 3.0, 3.0]) + np.testing.assert_allclose(fan[0, :, 0, Q_HI], 3.0 + Z10 * np.sqrt([1.0, 2.0, 3.0])) + + +def test_two_cutoffs_scored_separately(): + # cutoff 2: sigma 1, point 1 ; cutoff 4: sigma 1, point 3, spread + # sqrt([1, 2]) over the two steps + fan = _fan(_HIST, prediction_length=2, season_length=1, cutoffs=[2, 4]) + np.testing.assert_allclose(fan[0, 0, 0, Q_HI], 1.0 + Z10) + np.testing.assert_allclose(fan[1, 0, 0, Q_HI], 3.0 + Z10) + np.testing.assert_allclose(fan[1, 1, 0, Q_HI], 3.0 + Z10 * np.sqrt(2))