feat: quantile interpolation for arbitrary levels and quantiles

docs/examples/index.md

@@ -21,4 +21,5 @@ For model API details, see the [Model Hub](../model-hub.md).
 | [TiRex Family](tirex-family.ipynb) | Forecast with TiRex 1.0 and 2.0 | Python 3.11+ |
 | [TabPFN Family](tabpfn-family.ipynb) | Forecast with TabPFN-2 and TabPFN-3 (LOCAL) | Python 3.10–3.12; `TABPFN_TOKEN` |
 | [Toto Family](toto-family.ipynb) | Forecast with Toto 1.0 and 2.0 | Python 3.10+ |
+| [PatchTST-FM Family](patchtst-fm-family.ipynb) | Forecast with IBM Research and Granite PatchTST-FM r1/r2 | Python 3.11–3.13; GPU recommended |
 | [Finetuning](finetuning.ipynb) | Adapt Chronos 2 and TimeGPT to your data | Python 3.10+; GPU recommended |

docs/forecasting-parameters.md

@@ -57,3 +57,60 @@ anomalies_df = ff.detect_anomalies(df=df, freq="MS", level=99)
 * **`freq`** must match your data's timestamp spacing. Irregular or gapped series may fail inference.
 * **`h`** should cover the horizon you care about for evaluation or deployment.
 * For **anomaly detection**, use a horizon aligned with the seasonal cycle when possible.
+
+
+## Probabilistic forecasting (`level` and `quantiles`)
+
+Foundation models support two mutually exclusive ways to request probabilistic
+forecasts in `forecast()` and `cross_validation()`:
+
+### `level` — prediction intervals
+
+Pass confidence levels as percentages, e.g. `level=[80, 95]`. Each level `L`
+maps to symmetric quantiles `(α/2, 1 − α/2)` with `α = 1 − L/100`. The output
+includes `{model}-lo-{L}` and `{model}-hi-{L}` columns.
+
+```python
+fcst_df = ff.forecast(df=df, h=12, freq="MS", level=[80, 95])
+```
+
+### `quantiles` — direct quantile forecasts
+
+Pass quantile levels in `(0, 1)`, e.g. `quantiles=[0.1, 0.5, 0.9]`. The output
+includes `{model}-q-{pct}` columns where `pct = int(100 × quantile)`.
+
+```python
+fcst_df = ff.forecast(df=df, h=12, freq="MS", quantiles=[0.1, 0.5, 0.9])
+```
+
+### Point forecasts
+
+The point forecast is **always** returned in the `{model}` column, regardless of
+whether you use `level` or `quantiles`. You do not need a special level to
+request the median.
+
+### Interpolation on fixed-knot models
+
+Some models (TiRex, TimesFM, TabPFN, FlowState, Tafsut, Toto 2.0, …) predict a
+fixed set of native quantile **knots** internally. When you request levels or
+quantiles that do not exactly match those knots, the library linearly
+interpolates between adjacent knots.
+
+### Edge clamping
+
+When a requested quantile falls **outside** the range the model supports, the
+forecast is computed at the nearest in-range quantile (same semantics as
+`numpy.interp` on fixed knots). Values are not extrapolated beyond that range.
+Output column names still reflect **your** request (`level` / `quantiles`).
+
+Fixed-knot models interpolate between native knots and clamp to the knot edges.
+Native quantile backends (for example PatchTST-FM, T0, Chronos) clamp to their
+documented quantile range before calling the model, then map results back to
+your requested levels.
+
+For example, on a model with native knots `0.1` through `0.9`:
+
+- `level=[95]` maps to quantiles `0.025` and `0.975`, which are clamped to
+  `0.1` and `0.9` respectively, with columns `{model}-lo-95` and `{model}-hi-95`.
+- `quantiles=[0.01]` returns the same forecast values as `quantiles=[0.1]`, with
+  column `{model}-q-1`.

foundationforecast/core/__init__.py

@@ -8,6 +8,12 @@
 )
 from .gluonts_forecaster import GluonTSForecaster
 from .multi_model import MultiModelForecasterMixin
+from .quantiles import (
+    FIXED_KNOT_QUANTILES_NOTE,
+    interpolate_quantiles,
+    resolve_quantile_values,
+    validate_levels,
+)
 from .utils import (
     PanelData,
     TimeSeriesDataset,
@@ -17,15 +23,19 @@
 
 __all__ = [
     "Forecaster",
+    "FIXED_KNOT_QUANTILES_NOTE",
     "GluonTSForecaster",
     "MultiModelForecasterMixin",
     "PanelData",
     "QuantileConverter",
+    "interpolate_quantiles",
+    "resolve_quantile_values",
     "TimeSeriesDataset",
     "grouped_std_by_id",
     "process_panel_from_df",
     "_DataProcessor",
     "get_seasonality",
     "maybe_convert_col_to_datetime",
     "maybe_infer_freq",
+    "validate_levels",
 ]

foundationforecast/core/forecaster.py

@@ -26,6 +26,7 @@
     vertical_concat,
 )
 
+from .quantiles import _LEVEL_ZERO_ERROR
 from .utils import PanelData, TimeSeriesDataset, grouped_std_by_id
 
 T = TypeVar("T")
@@ -237,6 +238,11 @@ def _assign_quantile_forecasts(
         fcsts_quantiles_np: np.ndarray,
     ) -> pd.DataFrame:
         q_cols = [f"{alias}-q-{int(q * 100)}" for q in quantiles]
+        if len(q_cols) != len(set(q_cols)):
+            raise ValueError(
+                "Requested quantiles map to duplicate output column names "
+                f"(using int(100 × quantile) suffixes): {quantiles}"
+            )
         q_vals = [fcsts_quantiles_np[..., i].reshape(-1) for i in range(len(quantiles))]
         for q_col, q_val in zip(q_cols, q_vals, strict=True):
             fcst_df = ufp.assign_columns(fcst_df, q_col, q_val)
@@ -381,7 +387,23 @@ def detect_anomalies(
 
 
 class QuantileConverter:
-    """Handles inputs and outputs for probabilistic forecasts."""
+    """Handles inputs and outputs for probabilistic forecasts.
+
+    Pass either ``level`` (confidence intervals as percentages) or ``quantiles``
+    (values in ``(0, 1)``), but not both. The point forecast is always returned
+    in the model alias column regardless of which probabilistic API is used.
+
+    When ``level`` is provided, each level ``L`` maps to symmetric quantiles
+    ``(α/2, 1 − α/2)`` with ``α = 1 − L/100``. The output uses
+    ``{alias}-lo-{L}`` and ``{alias}-hi-{L}`` columns.
+
+    When ``quantiles`` is provided, the output uses ``{alias}-q-{pct}`` columns
+    where ``pct = int(100 × quantile)``.
+
+    Fixed-knot models interpolate linearly from their native quantile knots to
+    the requested levels or quantiles, clamping to edge knots when needed. See
+    ``foundationforecast.core.quantiles`` for edge-clamping semantics.
+    """
 
     def __init__(
         self,
@@ -405,6 +427,8 @@ def _prepare_level_and_quantiles(
                 "You must not provide both `level` and `quantiles` simultaneously."
             )
         if quantiles is None and level is not None:
+            if 0 in level:
+                raise ValueError(_LEVEL_ZERO_ERROR)
             _quantiles = []
             for lv in level:
                 q_lo, q_hi = QuantileConverter._level_to_quantiles(lv)
@@ -463,12 +487,6 @@ def maybe_convert_quantiles_to_level(
         out_cols = [c for c in df.columns if "-q-" not in c]
         df = ufp.copy_if_pandas(df, deep=False)
         for model in models:
-            if 0 in self.level:
-                mid_col = f"{model}-q-50"
-                if mid_col in df:
-                    df = ufp.assign_columns(df, model, df[mid_col])
-                    if model not in out_cols:
-                        out_cols.append(model)
             for lv in self.level:
                 q_lo, q_hi = self._level_to_quantiles(lv)
                 lo_src = f"{model}-q-{int(q_lo * 100)}"

foundationforecast/core/quantiles.py

@@ -0,0 +1,177 @@
+"""Quantile interpolation utilities for fixed-knot forecast models.
+
+Fixed-knot models predict a discrete set of native quantile levels (knots).
+User-requested ``level`` or ``quantiles`` values that do not match those knots
+are obtained by piecewise linear interpolation along the quantile axis.
+
+Edge clamping
+-------------
+When a requested quantile falls below the lowest native knot (or above the
+highest), the forecast at the nearest edge knot is returned. This matches
+``numpy.interp`` semantics: values are not extrapolated beyond the model's
+native quantile range. For example, requesting ``q=0.01`` on a model with
+knots ``[0.1, ..., 0.9]`` returns the same values as ``q=0.1``.
+"""
+
+from __future__ import annotations
+
+from collections.abc import Sequence
+
+import numpy as np
+from scipy.interpolate import interp1d
+
+_LEVEL_ZERO_ERROR = (
+    "`level=0` is not supported. Point forecasts are always returned in the "
+    "model column. Pass standard confidence levels (e.g. `[80, 95]`) or use "
+    "`quantiles=` directly."
+)
+
+FIXED_KNOT_QUANTILES_NOTE = (
+    "This model predicts a fixed set of native quantile knots internally, then "
+    "linearly interpolates (with edge clamping) to any requested ``level`` or "
+    "``quantiles``. See ``foundationforecast.core.quantiles`` for details."
+)
+
+# Documented native quantile ranges for backends that evaluate quantiles directly.
+PATCHTST_FM_QUANTILE_RANGE = (0.01, 0.99)
+T0_ALPHA_QUANTILE_RANGE = (0.1, 0.9)
+T0_BETA_QUANTILE_RANGE = (0.01, 0.99)
+DEFAULT_NATIVE_QUANTILE_RANGE = (0.01, 0.99)
+
+
+def validate_levels(level: Sequence[int | float] | None) -> list[int | float] | None:
+    """Validate levels, rejecting the legacy ``level=0`` sentinel."""
+    if level is None:
+        return None
+    levels = list(level)
+    if any(lv == 0 for lv in levels):
+        raise ValueError(_LEVEL_ZERO_ERROR)
+    return levels
+
+
+def _match_quantile_level(levels: Sequence[float], q: float) -> int:
+    arr = np.asarray(levels, dtype=np.float64)
+    matches = np.where(np.isclose(arr, q, rtol=0.0, atol=1e-9))[0]
+    if matches.size == 0:
+        raise ValueError(
+            f"Quantile level {q} missing from backend output levels {list(levels)}"
+        )
+    return int(matches[0])
+
+
+def clip_quantiles_to_range(
+    quantiles: Sequence[float],
+    q_min: float,
+    q_max: float,
+) -> np.ndarray:
+    """Clip quantile levels to ``[q_min, q_max]`` for backend evaluation."""
+    return np.clip(np.asarray(quantiles, dtype=np.float64), q_min, q_max)
+
+
+def backend_quantile_levels(
+    requested: Sequence[float],
+    *,
+    q_min: float,
+    q_max: float,
+    include_median: bool = True,
+) -> list[float]:
+    """Unique sorted quantile levels to request from a native quantile backend."""
+    levels = {float(q) for q in clip_quantiles_to_range(requested, q_min, q_max)}
+    if include_median:
+        levels.add(float(np.clip(0.5, q_min, q_max)))
+    return sorted(levels)
+
+
+def select_clipped_quantile_values(
+    backend_levels: Sequence[float],
+    values: np.ndarray,
+    requested: Sequence[float],
+    *,
+    q_min: float,
+    q_max: float,
+    axis: int = -1,
+) -> np.ndarray:
+    """Map backend forecasts to user-requested quantiles with edge clamping."""
+    clipped = clip_quantiles_to_range(requested, q_min, q_max)
+    indices = [_match_quantile_level(backend_levels, float(q)) for q in clipped]
+    return np.take(values, indices, axis=axis)
+
+
+def _knot_indices(
+    knot_quantiles: np.ndarray,
+    requested: Sequence[float],
+) -> list[int] | None:
+    indices: list[int] = []
+    for q in requested:
+        matches = np.where(np.isclose(knot_quantiles, q))[0]
+        if len(matches) == 0:
+            return None
+        indices.append(int(matches[0]))
+    return indices
+
+
+def interpolate_quantiles(
+    knot_quantiles: Sequence[float],
+    knot_values: np.ndarray,
+    requested_quantiles: Sequence[float],
+    *,
+    axis: int = -1,
+) -> np.ndarray:
+    """Linearly interpolate knot forecasts onto ``requested_quantiles``.
+
+    Args:
+        knot_quantiles: Native quantile levels predicted by the model, in
+            ascending order.
+        knot_values: Forecast array with native quantiles along ``axis``.
+        requested_quantiles: Quantile levels to return.
+        axis: Axis index of ``knot_quantiles`` in ``knot_values``.
+
+    Returns:
+        Array with the same shape as ``knot_values`` except ``axis`` is replaced
+        by ``len(requested_quantiles)``.
+    """
+    knot_qs = np.asarray(knot_quantiles, dtype=np.float64)
+    if knot_qs.ndim != 1 or len(knot_qs) < 1:
+        raise ValueError("`knot_quantiles` must be a non-empty 1-D sequence.")
+    if len(knot_qs) == 1:
+        v = np.moveaxis(knot_values, axis, -1)
+        out = np.repeat(v[..., :1], len(requested_quantiles), axis=-1)
+        return np.moveaxis(out, -1, axis) if axis != -1 else out
+
+    requested = np.clip(
+        np.asarray(requested_quantiles, dtype=np.float64),
+        knot_qs[0],
+        knot_qs[-1],
+    )
+    values = np.moveaxis(knot_values, axis, -1)
+    orig_shape = values.shape[:-1]
+    flat = values.reshape(-1, len(knot_qs))
+    interp = interp1d(
+        knot_qs,
+        flat,
+        axis=1,
+        kind="linear",
+        assume_sorted=True,
+    )
+    out = interp(requested).reshape(*orig_shape, len(requested_quantiles))
+    return np.moveaxis(out, -1, axis) if axis != -1 else out
+
+
+def resolve_quantile_values(
+    knot_quantiles: Sequence[float],
+    knot_values: np.ndarray,
+    requested_quantiles: Sequence[float],
+    *,
+    axis: int = -1,
+) -> np.ndarray:
+    """Select or interpolate ``knot_values`` onto ``requested_quantiles``."""
+    knot_qs = np.asarray(knot_quantiles, dtype=np.float64)
+    indices = _knot_indices(knot_qs, requested_quantiles)
+    if indices is not None:
+        return np.take(knot_values, indices, axis=axis)
+    return interpolate_quantiles(
+        knot_quantiles,
+        knot_values,
+        requested_quantiles,
+        axis=axis,
+    )