diff --git a/dashboards/grafana/promforecast-capacity.json b/dashboards/grafana/promforecast-capacity.json index 30782f0..fdeeb49 100644 --- a/dashboards/grafana/promforecast-capacity.json +++ b/dashboards/grafana/promforecast-capacity.json @@ -453,6 +453,42 @@ } ] } + }, + { + "type": "timeseries", + "id": 10, + "title": "Forecast fan chart ($id quantiles)", + "description": "Predicted quantiles across the forecast horizon for the selected series. Enable per query with emission.quantiles (e.g. [0.1, 0.5, 0.9, 0.99]); requires the remote_write sink for the forward curve. The central forecast is overlaid for reference.", + "datasource": { + "type": "prometheus", + "uid": "${datasource}" + }, + "gridPos": { + "h": 9, + "w": 24, + "x": 0, + "y": 57 + }, + "targets": [ + { + "expr": "{__name__=~\".+_forecast_quantile\",group=~\"$group\",id=~\"$id\"}", + "legendFormat": "q={{q}} {{instance}}" + }, + { + "expr": "{__name__=~\".+_forecast\",group=~\"$group\",id=~\"$id\"}", + "legendFormat": "central {{instance}}" + } + ], + "fieldConfig": { + "defaults": { + "custom": { + "fillOpacity": 10, + "lineWidth": 1, + "showPoints": "never" + } + }, + "overrides": [] + } } ] } diff --git a/docker-compose.dev.yml b/docker-compose.dev.yml index e8c87ad..8da88e2 100644 --- a/docker-compose.dev.yml +++ b/docker-compose.dev.yml @@ -26,7 +26,7 @@ services: - "--retentionPeriod=30d" - "--futureRetention=2d" ports: - - "8428:8428" + - "8429:8428" volumes: - vm-data:/storage healthcheck: @@ -48,7 +48,7 @@ services: image: promforecast-http-demo:dev container_name: pf-http-demo ports: - - "8000:8000" + - "8001:8000" restart: unless-stopped http-demo-loadgen: @@ -96,7 +96,7 @@ services: container_name: pf-redis command: ["redis-server", "--save", "", "--appendonly", "no"] ports: - - "6379:6379" + - "6378:6379" healthcheck: test: ["CMD", "redis-cli", "ping"] interval: 5s @@ -113,7 +113,7 @@ services: - "--storage.tsdb.retention.time=2h" - "--web.enable-lifecycle" ports: - - "9090:9090" + - "9080:9090" volumes: - ./docker/prometheus.yml:/etc/prometheus/prometheus.yml:ro - prom-data:/prometheus @@ -129,7 +129,7 @@ services: container_name: pf-forecaster command: ["--config", "/etc/promforecast/config.yaml"] ports: - - "9091:9091" + - "9081:9091" volumes: # Dev-tuned config: short lookback / step / min_points so the first fit # happens minutes after `make dev-up` instead of needing days of data. @@ -156,7 +156,7 @@ services: GF_AUTH_ANONYMOUS_ORG_ROLE: Viewer GF_USERS_DEFAULT_THEME: dark ports: - - "3000:3000" + - "3001:3000" volumes: - ./docker/grafana/provisioning:/etc/grafana/provisioning:ro - ./dashboards/grafana:/etc/grafana/dashboards:ro diff --git a/docs/README.md b/docs/README.md index c2fc4ac..d515080 100644 --- a/docs/README.md +++ b/docs/README.md @@ -31,6 +31,7 @@ Everything the forecaster emits beyond raw forecast curves — derived metrics, calibration, explainability, and ad-hoc what-if endpoints. - [Capacity planning](outputs/capacity.md) — damped trend, time-to-threshold ETAs, growth rate, will-breach predictive boolean, queue burn-down (drain) forecasting +- [Arbitrary-quantile forecasts](quantile.md) — emit `_forecast_quantile{q}` at any quantile for risk-based alerting and fan charts - [SLO error-budget forecasting](outputs/slo.md) — burn rate, exhaustion projection, multi-window alerting - [Maintenance window recommendations](outputs/maintenance-windows.md) — top-N lowest-load windows for scheduling reboots, rollouts, and cert rotations - [Conformal calibration](outputs/conformal.md) — empirical-coverage band correction diff --git a/docs/operations/cardinality.md b/docs/operations/cardinality.md index 7acaebd..a2e81ab 100644 --- a/docs/operations/cardinality.md +++ b/docs/operations/cardinality.md @@ -40,6 +40,7 @@ Notation: | `_forecast_contribution_ratio` | `S × N_regressors` | Only with regressors | | `forecast_band_calibration_offset` | `S × M × L` | `accuracy.conformal.enabled: true` | | `_forecast_component` | `S × M × 4` | `emission.decompose: true` (trend/seasonal/level/residual) | +| `_forecast_quantile` | `S × M × N_quantiles` | `emission.quantiles` set; add `S × M × N_quantiles` more for the per-quantile `forecast_band_calibration_offset{q}` rows when `accuracy.conformal.enabled: true` (and one extra hold-out fit) | | `forecast_drift_score` | `S × M` | `emission.drift.enabled: true` | | Threshold ETA / will_breach | `S × M × N_thresholds × (1 + L)` | Only configured thresholds | | Drain ETA (`le_drain`) | `S × M × N_drain_thresholds × (1 + L)` | Charged separately from the standard ETA | diff --git a/docs/quantile.md b/docs/quantile.md new file mode 100644 index 0000000..6620353 --- /dev/null +++ b/docs/quantile.md @@ -0,0 +1,67 @@ +# Arbitrary-quantile forecasts + +The standard forecast emits a central value plus symmetric confidence bands (`_forecast_lower` / `_forecast_upper`). +For risk-based alerting (“P[disk full within 24h] > 0.05”) and fan charts, you want predictions at *arbitrary* quantiles instead. + +`emission.quantiles` emits a parallel family from the **same single fit**: + +``` +_forecast_quantile{q="0.95", id, group, model, horizon} +``` + +## Configuring + +```yaml +groups: + - name: capacity + queries: + - id: node_filesystem_avail_bytes + promql: node_filesystem_avail_bytes + emission: + quantiles: [0.1, 0.25, 0.5, 0.75, 0.9, 0.99] +``` + +Values must be strictly between 0 and 1. Duplicates are rejected at load time. The existing band family remains for backward compatibility. + +## How a quantile is derived + +Each requested quantile is read from the appropriate band edge of the level that brackets it: + +| Quantile `q` | Reads from | +|--------------|------------| +| `q = 0.5` | central forecast (`yhat`) | +| `q < 0.5` | lower edge of level `round((1 − 2q)·100)` | +| `q > 0.5` | upper edge of level `round((2q − 1)·100)` | + +Quantiles beyond the configured confidence levels are clamped to the outermost band. + +## Conformal calibration per quantile + +When `accuracy.conformal.enabled: true` is also set, each quantile receives its own empirical-coverage correction. The offset is exposed alongside band-level offsets in the existing gauge, discriminated by the `q` label: + +``` +forecast_band_calibration_offset{id, group, model, q="0.95"} +``` + +This adds one extra hold-out fit per emitted model per refresh. + +## Fan chart & risk-based alerting + +The bundled capacity dashboard includes a “Forecast fan chart” panel that overlays all emitted quantiles plus the central forecast. + +Risk-based alerts become natural: + +```promql +# Optimistic 10th-percentile still predicts disk fill inside horizon +node_filesystem_avail_bytes_forecast_quantile{q="0.1", horizon="24h"} < 0 +``` + +```promql +# 99th-percentile request rate exceeds capacity +http_requests_rate_forecast_quantile{q="0.99"} > 5000 +``` + +## Cardinality + +`N_quantiles` rows per (series, model) per opted-in query, plus matching calibration-offset rows when conformal is enabled. +`promforecast validate` counts the family. diff --git a/docs/reference/models.md b/docs/reference/models.md index b11c906..2cd0ae9 100644 --- a/docs/reference/models.md +++ b/docs/reference/models.md @@ -1,8 +1,10 @@ # Models promforecast ships a small set of statistical forecasters from the -[`statsforecast`](https://github.com/Nixtla/statsforecast) library, plus a -`SeasonalNaive` baseline. +[`statsforecast`](https://github.com/Nixtla/statsforecast) library, a +`SeasonalNaive` baseline, and a family of in-tree count-distribution +models (`Poisson`, `NegativeBinomial`, `INGARCH`) for integer-valued +metrics. Configure models via the `models:` array (under `defaults` or per group). - **Explicit mode** (`accuracy.auto_select: false`): every chosen model runs and emits forecasts side-by-side. @@ -25,6 +27,9 @@ External models can be added via the `promforecast.models` entry-point group (se | **CrostonSBA** | Syntetos-Boylan bias-corrected Croston (recommended default). | Sparse / intermittent demand | Continuous series with strong seasonality | Negligible | | **ADIDA** | Aggregate-Disaggregate Intermittent Demand Approach. | Very sparse signals where Croston struggles | Continuous series with strong seasonality | Negligible | | **IMAPA** | Intermittent Multiple Aggregation Prediction Algorithm. | Noisy intermittent signals | Continuous series with strong seasonality | Negligible | +| **Poisson** | Seasonal-rate Poisson model; Poisson quantile bands (non-negative, asymmetric). | Integer counts where variance ≈ mean | Overdispersed counts; continuous series | Negligible | +| **NegativeBinomial**| Seasonal-rate model with negative-binomial bands; overdispersion from method of moments. | Integer counts where variance > mean (common) | Continuous series | Negligible | +| **INGARCH** | Integer-valued GARCH(1,1) — Poisson autoregression `λ_t = ω + α·y_{t-1} + β·λ_{t-1}` by MLE. | Autocorrelated / bursty count series | Series with no count autocorrelation | Sub-second per series | ### When to choose damped @@ -85,6 +90,61 @@ emission: Applies to both `/metrics` and `remote_write`. +These safeguards are a *post-processing* fix for a continuous model fit to +count data — they correct the symptom (negative / fractional output) but +not the mis-specified band shape. For a proper fit, use a +count-distribution model. + +## Count-distribution models + +`Poisson`, `NegativeBinomial`, and `INGARCH` model the integer nature of +the data directly, so the predicted bands respect the distribution: they +cannot extend below zero and they are asymmetric where the count +distribution is. The seasonal rate is the per-cycle-position sample mean +(the maximum-likelihood rate of a saturated seasonal Poisson); the band +families differ in how they spread around it. + +**Decision rule** — pick by the relationship between the series' +variance and mean (and whether it is autocorrelated): + +| Series characteristic | Model | +|----------------------------------------------------|--------------------| +| variance ≈ mean (pure Poisson arrivals) | `Poisson` | +| variance > mean (overdispersed — the common case) | `NegativeBinomial` | +| counts show autocorrelation (bursty / self-exciting) | `INGARCH` | + +Turn them on by adding them to `models:` and declaring the count profile on +the query so auto-select prefers them while keeping the continuous models +as a safety net: + +```yaml +defaults: + models: [AutoARIMA, Poisson, NegativeBinomial, INGARCH] + accuracy: + auto_select: true # MASE-scored for count (MAPE is undefined on zeros) + +groups: + - name: errors + queries: + - id: app_error_count + promql: sum by (instance) (rate(app_errors_total[5m])) + data_profile: count # per-query: prefer the count family, MASE scoring +``` + +When `data_profile: count` is set, auto-select scores candidates by MASE +(MAPE is undefined wherever the actual is zero — common for error counts) +and prefers the registered count wrappers, falling back to the continuous +models if none score finitely. The [count-aware clip-and-round +safeguards](#count-aware-emission-safeguards) remain on by default under +`data_profile: count` as a fallback for any continuous model that still +wins the horse race, but they are no longer the primary handling. + +Count wrappers ignore `regressors`; `INGARCH` falls back to the +seasonal-rate Poisson projection when the series is too short for a stable +three-parameter fit or the optimiser does not converge to a stationary +solution, so a deployment can never end up with a diverging count +forecast. + ## Cold-start seasonality borrowing See [`../inputs/cold-start.md`](../inputs/cold-start.md). Short version: when history < `min_points`, set `cold_start.borrow_from.promql` to synthesise a `SeasonalNaive`-style forecast from peer siblings. Emits with `model="ColdStart"` and `cold_start="true"`. Quality score is capped (default 0.5). diff --git a/forecaster/pyproject.toml b/forecaster/pyproject.toml index a937414..109a955 100644 --- a/forecaster/pyproject.toml +++ b/forecaster/pyproject.toml @@ -37,6 +37,11 @@ dependencies = [ # errors at runtime" UX when the extras aren't installed. "ruptures>=1.1", "statsmodels>=0.14", + # ``scipy``: count-distribution quantile functions (``scipy.stats``) + # and the INGARCH maximum-likelihood fit (``scipy.optimize``). Already + # present transitively via statsforecast/statsmodels; pinned directly + # because ``promforecast.models.count_distribution`` imports it. + "scipy>=1.13", ] [project.optional-dependencies] @@ -86,6 +91,12 @@ croston = "promforecast.models.intermittent:CrostonModel" croston_sba = "promforecast.models.intermittent:CrostonSBAModel" adida = "promforecast.models.intermittent:ADIDAModel" imapa = "promforecast.models.intermittent:IMAPAModel" +# Count-distribution models for integer-valued metrics (request counts, +# queue depths, error counts). Bands respect the count distribution +# (non-negative, asymmetric) rather than relying on clip-and-round. +poisson = "promforecast.models.count_distribution:PoissonModel" +negative_binomial = "promforecast.models.count_distribution:NegativeBinomialModel" +ingarch = "promforecast.models.count_distribution:INGARCHModel" [project.urls] Homepage = "https://github.com/esops-dev/promforecast" @@ -170,6 +181,16 @@ ignore_missing_imports = true module = ["ruptures", "statsmodels.tsa.seasonal"] ignore_missing_imports = true +[[tool.mypy.overrides]] +# ``scipy`` ships partial type information but not for the +# ``scipy.stats`` distribution objects or ``scipy.optimize.minimize`` +# result attributes the count-distribution models reach into. Those +# imports are confined to ``promforecast.models.count_distribution`` +# behind a narrow typed surface, so silencing missing imports here keeps +# mypy --strict clean without chasing incomplete third-party stubs. +module = ["scipy", "scipy.*"] +ignore_missing_imports = true + [[tool.mypy.overrides]] # The ``holidays`` package (optional extras dependency for the holiday # regressor) ships partial stubs diff --git a/forecaster/src/promforecast/config.py b/forecaster/src/promforecast/config.py index 24bc4a0..2f536f5 100644 --- a/forecaster/src/promforecast/config.py +++ b/forecaster/src/promforecast/config.py @@ -688,7 +688,7 @@ class BandCoverageConfig(BaseModel): class QueryEmissionConfig(BaseModel): """Per-query emission toggles. - Today six families of settings live here: + Several families of settings live here: * :class:`GrowthRateEmissionConfig` — forecast growth-rate metric. * ``clip_non_negative`` and ``round_to_integer`` — count-aware @@ -699,6 +699,10 @@ class QueryEmissionConfig(BaseModel): gauges. * ``conditional`` — parallel ``_forecast_conditional*`` family that re-runs the fit with regressors held at their current value. + * ``decompose`` — ``_forecast_component`` trend/seasonal/level/ + residual split of the predicted curve. + * ``quantiles`` — parallel ``_forecast_quantile{q}`` family at + arbitrary quantiles, derived from the same fit. Every family is kept on the per-query :class:`EmissionConfig` block so operators opt in *per metric*: a config that mixes a count-valued @@ -735,6 +739,19 @@ class QueryEmissionConfig(BaseModel): # per (series, model) per horizon step. See `models.py` for which models # actually decompose; non-decomposable models bump a skip counter. decompose: bool = False + # Arbitrary-quantile emission. Empty (default) emits nothing extra; a + # list like ``[0.1, 0.5, 0.9, 0.99]`` emits a parallel + # ``_forecast_quantile{q="0.9"}`` family alongside the + # ``_forecast_lower``/``_forecast_upper`` band family (which stays for + # backward compatibility). Each quantile is derived from the same + # single fit — no re-fit per quantile — by mapping it to the + # statsforecast confidence level that brackets it; ``0.5`` is the + # central forecast. Values must be strictly between 0 and 1. When + # ``accuracy.conformal.enabled`` is on, a per-quantile coverage + # calibration widens each quantile and the offset is exposed as + # ``forecast_band_calibration_offset{q="..."}``. See + # ``docs/quantile.md``. + quantiles: list[float] = Field(default_factory=list) class ColdStartBorrowConfig(BaseModel): @@ -2077,10 +2094,38 @@ def _validate_query_overrides(*, group_name: str, query: QueryConfig) -> None: ) if query.refresh_interval is not None and query.refresh_interval.total_seconds() <= 0: raise ValueError(f"{label}.refresh_interval must be a positive duration") + _validate_quantiles(label=label, query=query) _validate_slo(label=label, query=query) _validate_recommended_windows(label=label, query=query) +def _validate_quantiles(*, label: str, query: QueryConfig) -> None: + """Semantic checks for ``emission.quantiles``. + + Pydantic accepts any list of floats; the meaningful constraints + (strictly inside ``(0, 1)``, no duplicates) live here so a bad + quantile fails at load time with a precise label rather than emitting + a degenerate ``q="0"`` / ``q="1"`` row (the 0th/100th percentile maps + to an infinite band edge) every refresh. + """ + quantiles = query.emission.quantiles + if not quantiles: + return + seen: set[float] = set() + for q in quantiles: + if not (0.0 < q < 1.0): + raise ValueError( + f"{label}.emission.quantiles: every quantile must be strictly " + f"between 0 and 1 (got {q}); the 0th/100th percentile maps to an " + "infinite band edge" + ) + if q in seen: + raise ValueError( + f"{label}.emission.quantiles: duplicate quantile {q}; each value must be unique" + ) + seen.add(q) + + def _validate_slo(*, label: str, query: QueryConfig) -> None: """Semantic checks for the SLO block of a query. @@ -2194,6 +2239,8 @@ def _validate_events_profile(*, label: str, query: QueryConfig) -> None: rejected.append("emission.decompose") if query.emission.growth_rate.enabled: rejected.append("emission.growth_rate") + if query.emission.quantiles: + rejected.append("emission.quantiles") if rejected: raise ValueError( f"{label}: data_profile: events does not support " diff --git a/forecaster/src/promforecast/conformal.py b/forecaster/src/promforecast/conformal.py index a742a4a..ad7e050 100644 --- a/forecaster/src/promforecast/conformal.py +++ b/forecaster/src/promforecast/conformal.py @@ -112,6 +112,63 @@ def compute_offsets( return CalibrationResult(offsets=offsets, samples=int(residuals.size)) +@dataclass(frozen=True) +class QuantileCalibrationResult: + """Per-quantile *signed* additive correction for conformalised quantiles.""" + + offsets: dict[float, float] + samples: int + + +def compute_quantile_offsets( + *, + calib_actuals: NDArray[np.float64], + calib_predictions: pd.DataFrame, + quantiles: list[float], +) -> QuantileCalibrationResult: + """Conformalise each quantile to its empirical coverage on the holdout. + + For a one-sided quantile ``q`` the calibrated prediction must satisfy + ``P(actual <= predicted_q + offset) = q`` on the calibration window. + The empirical solution is ``offset = quantile_q(actual - predicted_q)`` + — the ``q``-quantile of the signed residuals — which is the standard + additive conformalised-quantile-regression correction. Unlike the + two-sided band offset this is *signed*: a quantile the model placed too + high gets a negative correction. + + ``calib_predictions`` is the model's forecast over the calibration + window carrying the per-level band columns from which each quantile is + derived. Quantiles whose source band column is absent are omitted from + the result (the caller treats "missing quantile" as "no calibration + this run"). + """ + from . import quantiles as quantiles_mod # noqa: PLC0415 + + if calib_actuals.size == 0 or calib_predictions.empty: + return QuantileCalibrationResult(offsets={}, samples=0) + n = min(calib_actuals.size, len(calib_predictions)) + if n < 2: # noqa: PLR2004 — two points is the bare minimum for a quantile + return QuantileCalibrationResult(offsets={}, samples=n) + actuals = calib_actuals[:n] + offsets: dict[float, float] = {} + samples = 0 + for q in quantiles: + series = quantiles_mod.predicted_quantile_series(calib_predictions, q) + if series is None: + continue + pred = series.to_numpy(dtype=float)[:n] + signed = actuals - pred + finite = np.isfinite(signed) + if not finite.any(): + continue + finite_signed = signed[finite] + offset = float(np.quantile(finite_signed, max(0.0, min(1.0, q)))) + if math.isfinite(offset): + offsets[q] = offset + samples = max(samples, int(finite_signed.size)) + return QuantileCalibrationResult(offsets=offsets, samples=samples) + + def apply_offsets( forecast: pd.DataFrame, offsets: dict[int, float], diff --git a/forecaster/src/promforecast/models/__init__.py b/forecaster/src/promforecast/models/__init__.py index 3a426ef..65e196e 100644 --- a/forecaster/src/promforecast/models/__init__.py +++ b/forecaster/src/promforecast/models/__init__.py @@ -33,6 +33,12 @@ from .auto_ets import AutoETSModel from .auto_ets_damped import AutoETSDampedModel from .base import ForecastModel +from .count_distribution import ( + COUNT_MODEL_NAMES, + INGARCHModel, + NegativeBinomialModel, + PoissonModel, +) from .intermittent import ( INTERMITTENT_MODEL_NAMES, ADIDAModel, @@ -228,6 +234,7 @@ def build( __all__ = [ + "COUNT_MODEL_NAMES", "ENTRY_POINT_GROUP", "INTERMITTENT_MODEL_NAMES", "ADIDAModel", @@ -239,8 +246,11 @@ def build( "CrostonSBAModel", "ForecastModel", "IMAPAModel", + "INGARCHModel", "MSTLModel", "ModelRegistrationError", + "NegativeBinomialModel", + "PoissonModel", "SeasonalNaiveModel", "UnknownModelError", "available", diff --git a/forecaster/src/promforecast/models/count_distribution.py b/forecaster/src/promforecast/models/count_distribution.py new file mode 100644 index 0000000..7542c87 --- /dev/null +++ b/forecaster/src/promforecast/models/count_distribution.py @@ -0,0 +1,403 @@ +"""Count-distribution model wrappers (Poisson, negative-binomial, INGARCH). + +The continuous-trend models (AutoARIMA, AutoETS, ...) assume a Gaussian +error and produce symmetric bands that routinely extend below zero on +integer-valued metrics — request counts, queue depths, error counts. +The count-aware emission safeguards (``emission.clip_non_negative`` / +``round_to_integer``) paper over the symptom by clipping and rounding the +output, but the underlying fit is still mis-specified: a Gaussian model +fit to a low-rate Poisson signal has the wrong band shape regardless of +how the output is post-processed. + +These three wrappers model the integer nature of the data directly so the +predicted bands respect the distribution — they cannot extend below zero +and they are asymmetric where the count distribution is: + +* :class:`PoissonModel` — the conditional mean follows a saturated + seasonal model (the per-cycle-position sample mean, which is the + maximum-likelihood rate estimate under a seasonal Poisson), and the + bands are the Poisson quantiles at that rate. Use when the variance is + approximately equal to the mean (the Poisson assumption). +* :class:`NegativeBinomialModel` — same seasonal mean, but the bands come + from a negative-binomial distribution whose overdispersion is estimated + by method of moments. Use when the variance exceeds the mean + (overdispersed counts — the common real-world case). +* :class:`INGARCHModel` — an integer-valued GARCH(1,1) (a.k.a. Poisson + autoregression): the conditional rate ``lambda_t = omega + alpha * + y_{t-1} + beta * lambda_{t-1}`` is fit by maximum likelihood. Use when + the count series shows autocorrelation (bursty arrivals, self-exciting + processes) that a static seasonal rate misses. + +**Decision rule.** Poisson when ``variance ≈ mean``; negative binomial +when ``variance > mean``; INGARCH when the series shows autocorrelation. +The auto-selector prefers this family for ``data_profile: count`` while +keeping the continuous models as a safety net — see +``docs/reference/models.md``. + +Each wrapper ignores the ``regressors`` argument (``uses_regressors = +False``); count regression with exogenous inputs is a future extension. +Bands are produced for the requested ``levels`` only; an empty ``levels`` +list yields a point forecast with no band columns, matching the other +wrappers. +""" + +from __future__ import annotations + +import math +from typing import TYPE_CHECKING, Any, ClassVar + +import numpy as np +import structlog + +if TYPE_CHECKING: + import pandas as pd + from numpy.typing import NDArray + +logger = structlog.get_logger(__name__) + +# Per-level ``(lower, upper)`` band arrays. Typed with ``Any`` element +# arrays to sidestep numpy's invariant shape generics — the same +# pragmatic choice the point factory makes for its per-step value arrays. +_Bands = dict[int, "tuple[Any, Any]"] + +# A rate floor keeps the Poisson/negative-binomial quantile functions and +# the INGARCH log-likelihood well-defined when a seasonal position (or the +# whole series) is all zeros. It is small enough not to perturb a genuine +# low rate and large enough to avoid ``log(0)`` in the INGARCH NLL. +_RATE_FLOOR = 1e-6 +# INGARCH needs enough history for the maximum-likelihood recursion to be +# meaningful; below this we fall back to the seasonal-mean Poisson path +# rather than fit three parameters on a handful of points. +_INGARCH_MIN_POINTS = 12 +# Upper bound on the method-of-moments overdispersion estimate. A handful +# of large residuals on a low-rate series can drive the per-point ratio +# ``(resid^2 - mu) / mu^2`` to enormous values; without a cap the +# negative-binomial bands degenerate to a near-useless [0, very-large] +# span. Real overdispersion rarely exceeds single digits, so a generous +# ceiling rejects only the noise-driven pathological case. +_MAX_OVERDISPERSION = 100.0 + + +def _two_sided_tail_probs(level: int) -> tuple[float, float]: + """Return the (lower, upper) quantile probabilities for a band level. + + A nominal ``level`` (e.g. 80) describes a two-sided central interval, + so the lower edge is the ``(1 - level/100) / 2`` quantile and the + upper edge its mirror. Level 80 -> (0.1, 0.9); level 95 -> (0.025, + 0.975). + """ + alpha = 1.0 - max(0, min(100, level)) / 100.0 + return alpha / 2.0, 1.0 - alpha / 2.0 + + +def _seasonal_mean_rates( + y: NDArray[np.float64], season_length: int, horizon: int +) -> NDArray[np.float64]: + """Project the per-cycle-position sample mean across the horizon. + + Each forecast step inherits the mean of the in-sample observations + that share its clock position (index modulo ``season_length``) — the + maximum-likelihood rate of a saturated seasonal Poisson model. + Positions with no in-sample data fall back to the global mean so the + projection is always defined. + """ + n = y.size + sl = max(1, season_length) + global_mean = float(np.mean(y)) if n else 0.0 + sums = np.zeros(sl, dtype=float) + counts = np.zeros(sl, dtype=float) + positions = np.arange(n) % sl + np.add.at(sums, positions, y) + np.add.at(counts, positions, 1.0) + means = np.where(counts > 0.0, sums / np.maximum(counts, 1.0), global_mean) + future_positions = np.arange(n, n + horizon) % sl + rates: NDArray[np.float64] = np.maximum(means[future_positions], 0.0) + return rates + + +def _future_timestamps(df: pd.DataFrame, horizon: int, freq: str) -> Any: + """Build ``horizon`` future timestamps on the configured ``freq`` grid. + + Mirrors the forward ``ds`` index ``statsforecast`` synthesises for the + other wrappers — anchored at the last input timestamp and stepped by + the configured ``freq`` (e.g. ``"5min"``). Deriving the grid from + ``freq`` rather than the empirical input spacing keeps the count + wrappers byte-consistent with every other model, so anything that reads + the forecast ``ds`` (the sink curve, capacity ETAs) sees the same grid + regardless of which model produced the row. + """ + import pandas as pd # noqa: PLC0415 + + ds = pd.to_datetime(df["ds"]) + # A fully-empty input has no anchor timestamp; fall back to the epoch so + # the frame is still well-formed (the values are all zero anyway). + last = ds.iloc[-1] if len(ds) else pd.Timestamp("1970-01-01") + return pd.date_range(start=last, periods=horizon + 1, freq=freq)[1:] + + +def _assemble_frame( + *, + df: pd.DataFrame, + horizon: int, + freq: str, + yhat: NDArray[np.float64], + bands: _Bands, +) -> pd.DataFrame: + """Assemble the protocol-shaped forecast frame from arrays. + + ``bands`` maps each level to its ``(lower, upper)`` arrays. The output + carries ``unique_id``, ``ds``, ``yhat`` and a + ``yhat_lower_`` / ``yhat_upper_`` pair per level — the + same columns ``_normalize`` produces for the statsforecast wrappers. + """ + import pandas as pd # noqa: PLC0415 + + unique_id = df["unique_id"].iloc[0] if len(df) else "s" + out: dict[str, Any] = { + "unique_id": [unique_id] * horizon, + "ds": _future_timestamps(df, horizon, freq), + "yhat": yhat, + } + for level, (lower, upper) in bands.items(): + out[f"yhat_lower_{level}"] = lower + out[f"yhat_upper_{level}"] = upper + return pd.DataFrame(out) + + +def _poisson_bands(rates: NDArray[np.float64], levels: list[int]) -> _Bands: + """Per-level Poisson quantile bands at the given per-step rates.""" + from scipy import stats # noqa: PLC0415 + + safe = np.maximum(rates, _RATE_FLOOR) + bands: _Bands = {} + for level in levels: + p_lo, p_hi = _two_sided_tail_probs(level) + lower = np.asarray(stats.poisson.ppf(p_lo, safe), dtype=float) + upper = np.asarray(stats.poisson.ppf(p_hi, safe), dtype=float) + bands[level] = (np.maximum(lower, 0.0), np.maximum(upper, 0.0)) + return bands + + +def _estimate_overdispersion(y: NDArray[np.float64], fitted_mean: NDArray[np.float64]) -> float: + """Method-of-moments overdispersion ``alpha`` for ``var = mu + alpha*mu^2``. + + Averaged over the in-sample points where the fitted mean is positive; + clamped at zero so an under-dispersed (variance < mean) series + degrades cleanly to the Poisson limit rather than producing a negative + dispersion. + """ + mask = fitted_mean > _RATE_FLOOR + if not np.any(mask): + return 0.0 + mu = fitted_mean[mask] + resid_sq = (y[mask] - mu) ** 2 + per_point = (resid_sq - mu) / (mu**2) + alpha = float(np.mean(per_point)) + if not math.isfinite(alpha): + return 0.0 + return min(_MAX_OVERDISPERSION, max(0.0, alpha)) + + +def _negbinom_bands(rates: NDArray[np.float64], alpha: float, levels: list[int]) -> _Bands: + """Per-level negative-binomial quantile bands at the given rates. + + Uses the ``(var = mu + alpha*mu^2)`` parameterisation: ``size = 1 / + alpha`` and ``prob = 1 / (1 + alpha*mu)``. When ``alpha`` collapses to + the Poisson limit the bands fall back to the Poisson quantiles so the + two families agree at the boundary. + """ + if alpha <= _RATE_FLOOR: + return _poisson_bands(rates, levels) + from scipy import stats # noqa: PLC0415 + + safe = np.maximum(rates, _RATE_FLOOR) + size = 1.0 / alpha + prob = 1.0 / (1.0 + alpha * safe) + bands: _Bands = {} + for level in levels: + p_lo, p_hi = _two_sided_tail_probs(level) + lower = np.asarray(stats.nbinom.ppf(p_lo, size, prob), dtype=float) + upper = np.asarray(stats.nbinom.ppf(p_hi, size, prob), dtype=float) + bands[level] = (np.maximum(lower, 0.0), np.maximum(upper, 0.0)) + return bands + + +def _ingarch_recursion( + y: NDArray[np.float64], omega: float, alpha: float, beta: float +) -> NDArray[np.float64]: + """Conditional-rate path ``lambda_t = omega + alpha*y_{t-1} + beta*lambda_{t-1}``.""" + n = y.size + lam = np.empty(n, dtype=float) + lam[0] = max(float(np.mean(y)) if n else _RATE_FLOOR, _RATE_FLOOR) + for t in range(1, n): + lam[t] = omega + alpha * y[t - 1] + beta * lam[t - 1] + clamped: NDArray[np.float64] = np.maximum(lam, _RATE_FLOOR) + return clamped + + +def _fit_ingarch(y: NDArray[np.float64]) -> tuple[float, float, float] | None: + """Maximum-likelihood ``(omega, alpha, beta)`` for an INGARCH(1,1) fit. + + Minimises the Poisson negative log-likelihood of the conditional-rate + recursion under non-negativity bounds and the stationarity constraint + ``alpha + beta < 1`` (enforced via the bound ceilings). Returns + ``None`` when the optimiser fails to converge so the caller can fall + back to the seasonal-mean Poisson path. + """ + from scipy import optimize # noqa: PLC0415 + + def _nll(params: NDArray[np.float64]) -> float: + omega, alpha, beta = float(params[0]), float(params[1]), float(params[2]) + lam = _ingarch_recursion(y, omega, alpha, beta) + return float(np.sum(lam - y * np.log(lam))) + + mean_y = max(float(np.mean(y)), _RATE_FLOOR) + x0 = np.array([mean_y * 0.5, 0.3, 0.3], dtype=float) + bounds = [(_RATE_FLOOR, None), (0.0, 0.95), (0.0, 0.95)] + try: + res = optimize.minimize(_nll, x0, method="L-BFGS-B", bounds=bounds) + except Exception: # any optimiser failure routes to the seasonal-mean fallback + logger.warning("ingarch_fit_failed") + return None + if not bool(res.success): + return None + omega, alpha, beta = float(res.x[0]), float(res.x[1]), float(res.x[2]) + if alpha + beta >= 1.0: # non-stationary; the forecast would diverge + return None + return omega, alpha, beta + + +def _ingarch_forecast_rates( + y: NDArray[np.float64], params: tuple[float, float, float], horizon: int +) -> NDArray[np.float64]: + """Project the conditional rate across the horizon from a fitted INGARCH. + + The one-step-ahead rate uses the last observed value; subsequent steps + substitute ``E[y] = lambda`` so the projection converges to the + stationary mean ``omega / (1 - alpha - beta)``. + """ + omega, alpha, beta = params + lam_hist = _ingarch_recursion(y, omega, alpha, beta) + rates = np.empty(horizon, dtype=float) + rates[0] = omega + alpha * float(y[-1]) + beta * float(lam_hist[-1]) + for h in range(1, horizon): + rates[h] = omega + (alpha + beta) * rates[h - 1] + clamped: NDArray[np.float64] = np.maximum(rates, _RATE_FLOOR) + return clamped + + +class _CountBase: + """Shared boilerplate for the count-distribution wrappers. + + Subclasses implement :meth:`_forecast` returning the ``(yhat, bands)`` + pair; the base handles input coercion, the empty-history guard, and + frame assembly. ``season_length`` cascades from ``defaults`` like every + other model; the count wrappers use it to project the seasonal rate. + """ + + name: ClassVar[str] = "" + uses_regressors: ClassVar[bool] = False + + def __init__(self, season_length: int = 288, **params: Any) -> None: + self._season_length = season_length + self._params = params + + def fit_predict( + self, + df: pd.DataFrame, + horizon: int, + freq: str, + levels: list[int], + regressors: pd.DataFrame | None = None, + ) -> pd.DataFrame: + del regressors # count wrappers don't consume exogenous regressors + y = np.clip(np.asarray(df["y"].to_numpy(dtype=float)), 0.0, None) + # Non-finite inputs (gap-imputed NaNs that slipped through) are + # dropped so the rate estimate stays well-defined. + y = y[np.isfinite(y)] + if y.size == 0: + zeros: NDArray[np.float64] = np.zeros(horizon, dtype=float) + empty_bands: _Bands = {level: (zeros, zeros) for level in levels} + return _assemble_frame(df=df, horizon=horizon, freq=freq, yhat=zeros, bands=empty_bands) + yhat, bands = self._forecast(y=y, horizon=horizon, levels=levels) + return _assemble_frame(df=df, horizon=horizon, freq=freq, yhat=yhat, bands=bands) + + def _forecast( + self, *, y: NDArray[np.float64], horizon: int, levels: list[int] + ) -> tuple[NDArray[np.float64], _Bands]: + raise NotImplementedError + + +class PoissonModel(_CountBase): + """Seasonal-rate Poisson model with Poisson quantile bands.""" + + name = "Poisson" + + def _forecast( + self, *, y: NDArray[np.float64], horizon: int, levels: list[int] + ) -> tuple[NDArray[np.float64], _Bands]: + rates = _seasonal_mean_rates(y, self._season_length, horizon) + return rates, _poisson_bands(rates, levels) + + +class NegativeBinomialModel(_CountBase): + """Seasonal-rate model with negative-binomial (overdispersed) bands.""" + + name = "NegativeBinomial" + + def _forecast( + self, *, y: NDArray[np.float64], horizon: int, levels: list[int] + ) -> tuple[NDArray[np.float64], _Bands]: + rates = _seasonal_mean_rates(y, self._season_length, horizon) + # The method-of-moments dispersion estimate needs the fitted mean + # at each *observed* index (its clock-position mean), so build the + # position means once and index them back to 0..n-1. + sl = max(1, self._season_length) + positions = np.arange(y.size) % sl + sums = np.zeros(sl, dtype=float) + counts = np.zeros(sl, dtype=float) + np.add.at(sums, positions, y) + np.add.at(counts, positions, 1.0) + global_mean = float(np.mean(y)) + means = np.where(counts > 0.0, sums / np.maximum(counts, 1.0), global_mean) + fitted = means[positions] + alpha = _estimate_overdispersion(y, fitted) + return rates, _negbinom_bands(rates, alpha, levels) + + +class INGARCHModel(_CountBase): + """Integer-valued GARCH(1,1) — Poisson autoregression with MLE bands. + + Falls back to the seasonal-mean Poisson projection when the series is + too short for a meaningful three-parameter fit or the optimiser does + not converge to a stationary solution, so a deployment can never end + up with a diverging forecast. + """ + + name = "INGARCH" + + def _forecast( + self, *, y: NDArray[np.float64], horizon: int, levels: list[int] + ) -> tuple[NDArray[np.float64], _Bands]: + params = _fit_ingarch(y) if y.size >= _INGARCH_MIN_POINTS else None + if params is None: + rates = _seasonal_mean_rates(y, self._season_length, horizon) + else: + rates = _ingarch_forecast_rates(y, params, horizon) + return rates, _poisson_bands(rates, levels) + + +# Convenience set consumed by the auto-select preference for +# ``data_profile: count``. Self-maintains from ``_CountBase``'s direct +# subclasses, mirroring ``INTERMITTENT_MODEL_NAMES``: a new count wrapper +# lands in the preference whitelist by inheriting the base and setting a +# non-empty ``name``. +COUNT_MODEL_NAMES = frozenset(cls.name for cls in _CountBase.__subclasses__() if cls.name) + + +__all__ = [ + "COUNT_MODEL_NAMES", + "INGARCHModel", + "NegativeBinomialModel", + "PoissonModel", +] diff --git a/forecaster/src/promforecast/point_factory.py b/forecaster/src/promforecast/point_factory.py index ed9b239..9774f98 100644 --- a/forecaster/src/promforecast/point_factory.py +++ b/forecaster/src/promforecast/point_factory.py @@ -449,6 +449,31 @@ def calibration_offset_points( offset=float(offset), ) + def quantile_calibration_offset_points( + self, + model_name: str, + offsets: dict[float, float], + ) -> Iterable[CalibrationOffsetPoint]: + """One ``forecast_band_calibration_offset`` per (model, quantile). + + Mirrors :meth:`calibration_offset_points` but the discriminator is + the ``q`` label rather than ``level`` — the per-quantile conformal + correction lands in the same gauge family so operators read both + the band-level and per-quantile offsets from one metric name. The + offset is *signed* for quantiles (a quantile placed too high gets a + negative correction); non-finite values are skipped. + """ + from . import quantiles as quantiles_mod # noqa: PLC0415 + + common = self._common_with_id() + for q, offset in offsets.items(): + if not math.isfinite(offset): + continue + yield CalibrationOffsetPoint( + labels={**common, "model": model_name, "q": quantiles_mod.format_quantile(q)}, + offset=float(offset), + ) + def conditional_calibration_offset_points( self, model_name: str, @@ -583,6 +608,91 @@ def _curve_points( timestamp_ms=ts_ms, ) + def quantile_points( + self, + forecast: pd.DataFrame, + model_name: str, + selection: SelectionResult | None, + quantiles: list[float], + ) -> Iterable[ForecastPoint]: + """Convert the end-of-horizon row into ``_forecast_quantile`` points. + + Mirrors :meth:`forecast_points`' snapshot convention — one value + per (series, model, quantile) at the configured horizon. The + quantile columns (``yhat_quantile_``) must already be present on + ``forecast`` (the runner derives them via + :func:`promforecast.quantiles.add_quantile_columns` before calling + this). Quantiles whose column is absent are skipped. + """ + from . import quantiles as quantiles_mod # noqa: PLC0415 + + if forecast.empty: + return + last = forecast.iloc[-1] + common = self._build_common_labels(model_name, selection) + metric_name = f"{self._metric_id}_forecast_quantile" + for q in quantiles: + col = quantiles_mod.quantile_column(q) + if col not in forecast.columns: + continue + value = float(last[col]) + if not math.isfinite(value): + continue + yield ForecastPoint( + metric=metric_name, + labels={**common, "q": quantiles_mod.format_quantile(q)}, + value=value, + ) + + def quantile_curve_points( + self, + forecast: pd.DataFrame, + model_name: str, + selection: SelectionResult | None, + quantiles: list[float], + ) -> Iterable[ForecastCurvePoint]: + """Yield every horizon step of each quantile as a sink push point. + + Mirrors :meth:`forecast_curve_points` but emits under the + ``_forecast_quantile`` metric name with a ``q`` label. Rows + with non-finite timestamps or values are dropped — TSDBs reject + NaN samples. + """ + from . import quantiles as quantiles_mod # noqa: PLC0415 + + if forecast.empty or "ds" not in forecast.columns: + return + import pandas as pd # noqa: PLC0415 + + common = self._build_common_labels(model_name, selection) + metric_name = f"{self._metric_id}_forecast_quantile" + ds_normalised = pd.to_datetime(forecast["ds"], errors="coerce") + valid_ts_mask = ds_normalised.notna().to_numpy() + ts_ms_arr = ds_normalised.astype("int64").to_numpy() // 1_000_000 + quantile_cols: list[tuple[dict[str, str], Any]] = [] + for q in quantiles: + col = quantiles_mod.quantile_column(q) + if col not in forecast.columns: + continue + labels = {**common, "q": quantiles_mod.format_quantile(q)} + quantile_cols.append((labels, forecast[col].to_numpy(dtype=float))) + if not quantile_cols: + return + n = len(forecast) + for i in range(n): + if not valid_ts_mask[i]: + continue + ts_ms = int(ts_ms_arr[i]) + for labels, arr in quantile_cols: + v = arr[i] + if math.isfinite(v): + yield ForecastCurvePoint( + metric=metric_name, + labels=labels, + value=float(v), + timestamp_ms=ts_ms, + ) + def ensemble_weight_points( self, weights: dict[str, float], diff --git a/forecaster/src/promforecast/quantiles.py b/forecaster/src/promforecast/quantiles.py new file mode 100644 index 0000000..9c6437b --- /dev/null +++ b/forecaster/src/promforecast/quantiles.py @@ -0,0 +1,129 @@ +"""Arbitrary-quantile derivation from a model's confidence-level bands. + +The built-in models emit a central forecast plus a symmetric +``yhat_lower_`` / ``yhat_upper_`` band per requested +confidence level. A two-sided level ``L`` brackets the one-sided +quantiles ``(1 - L/100) / 2`` (the lower edge) and ``1 - (1 - L/100) / 2`` +(the upper edge) — so any one-sided quantile ``q`` can be read off the +band of the level that places ``q`` on one of its edges, *without +re-fitting the model at a new level*: + +* ``q == 0.5`` -> the central forecast (``yhat``). +* ``q < 0.5`` -> the lower edge of level ``round((1 - 2q) * 100)``. +* ``q > 0.5`` -> the upper edge of level ``round((2q - 1) * 100)``. + +For example ``q = 0.95`` reads ``yhat_upper_90``; ``q = 0.1`` reads +``yhat_lower_80``; ``q = 0.99`` reads ``yhat_upper_98``. The implied +levels are added to the single forward fit's requested level set so the +columns exist to read from — the fit itself is shared. + +This mapping is exact for symmetric-interval models and gives genuinely +asymmetric quantiles for the count-distribution models (whose per-level +bands are already asymmetric). Extreme quantiles clamp to the level-99 +band edge, which the docs call out. +""" + +from __future__ import annotations + +from typing import TYPE_CHECKING + +if TYPE_CHECKING: + import pandas as pd + +# Quantiles within this tolerance of 0.5 are treated as the median (the +# central forecast). A direct float equality would miss a value like +# ``0.5`` that round-tripped through YAML as ``0.49999999999``. +_MEDIAN_TOL = 1e-9 +# Levels are clamped to ``[1, 99]``: level 0 collapses the band to the +# point forecast and level 100 is an infinite edge, so quantiles closer to +# the tails than the 99% band can resolve are approximated by that edge. +_MIN_LEVEL = 1 +_MAX_LEVEL = 99 +# The median splits the lower / upper band-edge mapping. +_MEDIAN = 0.5 + + +def format_quantile(q: float) -> str: + """Render a quantile as the stable ``q`` label value (e.g. ``0.95``).""" + return f"{q:g}" + + +def quantile_column(q: float) -> str: + """Column name carrying the derived quantile values on a forecast frame.""" + return f"yhat_quantile_{format_quantile(q)}" + + +def level_and_side_for_quantile(q: float) -> tuple[int, str] | None: + """Map a one-sided quantile to the ``(level, side)`` band edge it reads. + + Returns ``None`` for the median (which reads ``yhat`` directly). + ``side`` is ``"lower"`` for ``q < 0.5`` and ``"upper"`` for ``q > 0.5``. + """ + if abs(q - _MEDIAN) <= _MEDIAN_TOL: + return None + if q < _MEDIAN: + raw = round((1.0 - 2.0 * q) * 100.0) + side = "lower" + else: + raw = round((2.0 * q - 1.0) * 100.0) + side = "upper" + level = max(_MIN_LEVEL, min(_MAX_LEVEL, int(raw))) + return level, side + + +def quantile_levels(quantiles: list[float]) -> list[int]: + """Sorted unique confidence levels needed to derive ``quantiles``. + + The median contributes no level (it reads ``yhat``); every other + quantile contributes the level of the band edge it maps to. The runner + folds these into the single forward fit's level set so the band + columns exist to read from. + """ + levels: set[int] = set() + for q in quantiles: + mapped = level_and_side_for_quantile(q) + if mapped is not None: + levels.add(mapped[0]) + return sorted(levels) + + +def predicted_quantile_series(forecast: pd.DataFrame, q: float) -> pd.Series | None: + """Return the (uncalibrated) per-step values for quantile ``q``. + + Reads ``yhat`` for the median, otherwise the mapped band-edge column. + Returns ``None`` when the required column is absent (the level wasn't + in the fit) so the caller can skip the quantile rather than emit + garbage. + """ + mapped = level_and_side_for_quantile(q) + if mapped is None: + return forecast["yhat"] if "yhat" in forecast.columns else None + level, side = mapped + col = f"yhat_{side}_{level}" + if col not in forecast.columns: + return None + return forecast[col] + + +def add_quantile_columns( + forecast: pd.DataFrame, + quantiles: list[float], + offsets: dict[float, float] | None = None, +) -> list[float]: + """Add ``yhat_quantile_`` columns to ``forecast`` in place. + + Each column is the band-edge series for that quantile plus its + (signed) conformal ``offsets`` correction when supplied. Returns the + list of quantiles that resolved to a column — quantiles whose source + band column is missing are skipped and excluded from the return so the + caller emits only the rows it can back with data. + """ + resolved: list[float] = [] + for q in quantiles: + base = predicted_quantile_series(forecast, q) + if base is None: + continue + offset = offsets.get(q, 0.0) if offsets else 0.0 + forecast[quantile_column(q)] = base + offset + resolved.append(q) + return resolved diff --git a/forecaster/src/promforecast/runner/_helpers.py b/forecaster/src/promforecast/runner/_helpers.py index 3039665..2679b61 100644 --- a/forecaster/src/promforecast/runner/_helpers.py +++ b/forecaster/src/promforecast/runner/_helpers.py @@ -62,6 +62,7 @@ def _format_band_window(window: timedelta) -> str: "_forecast_conditional_upper", "_forecast_conditional", "_forecast_component", + "_forecast_quantile", "_forecast_lower", "_forecast_upper", "_forecast", diff --git a/forecaster/src/promforecast/runner/_scheduler.py b/forecaster/src/promforecast/runner/_scheduler.py index d4a4adf..d30a44c 100644 --- a/forecaster/src/promforecast/runner/_scheduler.py +++ b/forecaster/src/promforecast/runner/_scheduler.py @@ -43,6 +43,7 @@ from .. import models as model_registry from .. import preprocess as preprocess_mod from .. import quality as quality_mod +from .. import quantiles as quantiles_mod from .. import regressors as regressors_mod from .. import revision as revision_mod from .. import right_sizing as right_sizing_mod @@ -2594,6 +2595,20 @@ async def _fit_series( # noqa: PLR0912, PLR0915 freq = _pandas_freq(eff_step) horizon = max(1, int(d.horizon.total_seconds() // eff_step.total_seconds())) levels = list(d.confidence_levels) + # Arbitrary-quantile emission derives each quantile from the band + # of the confidence level that brackets it (see + # :mod:`promforecast.quantiles`). The implied levels are folded into + # the single forward fit's level set — ``fit_levels`` — so the band + # columns exist to read from without re-fitting per quantile. The + # factory keeps ``levels`` (the configured ``confidence_levels``) so + # the ``_forecast_lower``/``_forecast_upper`` family is unchanged; + # the extra levels only ever surface through the quantile family. + quantiles = list(query.emission.quantiles) + fit_levels = ( + sorted(set(levels) | set(quantiles_mod.quantile_levels(quantiles))) + if quantiles + else levels + ) factory = SeriesPointFactory( metric_id=query.id, base_labels=frame.labels, @@ -2650,7 +2665,7 @@ async def _fit_series( # noqa: PLR0912, PLR0915 group_name=group_name, horizon=horizon, freq=freq, - levels=levels, + levels=fit_levels, season_length=eff_season_length, backtest_results=backtest_results, result=result, @@ -2665,7 +2680,7 @@ async def _fit_series( # noqa: PLR0912, PLR0915 group_name=group_name, horizon=horizon, freq=freq, - levels=levels, + levels=fit_levels, season_length=eff_season_length, backtest_results=backtest_results, result=result, @@ -2680,7 +2695,7 @@ async def _fit_series( # noqa: PLR0912, PLR0915 group_name=group_name, horizon=horizon, freq=freq, - levels=levels, + levels=fit_levels, season_length=eff_season_length, result=result, regressors=regressors_frame, @@ -2803,6 +2818,36 @@ async def _fit_series( # noqa: PLR0912, PLR0915 result=result, ) + # Arbitrary-quantile emission. Derives each quantile from the + # already-fit band columns (no extra forward fit) and appends the + # ``_forecast_quantile`` family into the same snapshot / curve + # lists so it inherits the standard safeguard, inheritance and + # rendering paths. When conformal calibration is on, a per-quantile + # coverage pass widens each quantile and the offset is exposed via + # ``forecast_band_calibration_offset{q="..."}``. + if quantiles and result.capacity_inputs: + quantile_offsets_by_model: dict[str, dict[float, float]] = {} + if d.accuracy.conformal.enabled: + quantile_offsets_by_model = await self._compute_quantile_conformal_offsets( + df=df, + emitted_inputs=list(result.capacity_inputs), + model_specs=model_specs, + freq=freq, + fit_levels=fit_levels, + quantiles=quantiles, + season_length=eff_season_length, + regressors=regressors_frame, + holdout_fraction=d.accuracy.conformal.holdout_fraction, + ) + forecast_points = self._emit_quantile_forecasts( + emitted_inputs=list(result.capacity_inputs), + quantiles=quantiles, + factory=factory, + forecast_points=forecast_points, + result=result, + quantile_offsets_by_model=quantile_offsets_by_model, + ) + # Free the retained per-model dataframes — they can be MB-scale, # and ``_RunContext`` retains them across the per-series loop. result.capacity_inputs.clear() @@ -2856,7 +2901,7 @@ async def _fit_series( # noqa: PLR0912, PLR0915 # ``ForecastQualityZero`` alert even when the intermittent # model is fitting fine. quality_kind: quality_mod.ScoreKind = ( - "mase" if query.data_profile == "intermittent" else "mape" + "mase" if query.data_profile in ("intermittent", "count") else "mape" ) result.quality = list( factory.quality_points( @@ -3760,18 +3805,20 @@ async def _auto_select_and_forecast( regressors: pd.DataFrame | None = None, ) -> list[ForecastPoint]: d = self._config.defaults - # ``data_profile: intermittent`` swaps the comparator from MAPE to - # MASE — MAPE is undefined wherever the actual is zero, which is - # precisely the regime intermittent models target. The selector - # contract is unit-agnostic; we hand it MASE the same way we'd - # hand it MAPE and the lowest finite score still wins. The - # ``candidate_preference`` whitelist nudges the search toward the - # registered intermittent wrappers when any are configured, so an - # operator who lists ``[AutoARIMA, Croston, CrostonSBA]`` gets the - # Croston family preferred without removing AutoARIMA as a safety - # net. + # ``data_profile: intermittent`` and ``count`` both swap the + # comparator from MAPE to MASE — MAPE is undefined wherever the + # actual is zero, which both regimes routinely hit (sparse + # intermittent series; low-rate count series like error counts). + # The selector contract is unit-agnostic; we hand it MASE the same + # way we'd hand it MAPE and the lowest finite score still wins. + # The ``candidate_preference`` whitelist nudges the search toward + # the registered intermittent / count wrappers when any are + # configured, so an operator who lists ``[AutoARIMA, Poisson, + # NegativeBinomial]`` under ``data_profile: count`` gets the count + # family preferred without removing AutoARIMA as a safety net. is_intermittent = query.data_profile == "intermittent" - if is_intermittent: + is_count = query.data_profile == "count" + if is_intermittent or is_count: score_map = { model_name: score for model_name, results in backtest_results.items() @@ -3783,7 +3830,12 @@ async def _auto_select_and_forecast( for model_name, results in backtest_results.items() if (mape := overall_mape(results)) is not None } - preference = model_registry.INTERMITTENT_MODEL_NAMES if is_intermittent else None + if is_intermittent: + preference: frozenset[str] | None = model_registry.INTERMITTENT_MODEL_NAMES + elif is_count: + preference = model_registry.COUNT_MODEL_NAMES + else: + preference = None selection = select( fallback_model=d.accuracy.fallback_model, fallback_mape_threshold=d.accuracy.fallback_mape_threshold, @@ -3911,11 +3963,13 @@ async def _timed(spec: ModelSpec) -> tuple[str, pd.DataFrame, float]: # Compute inverse-score weights restricted to models that produced # a forecast this run; otherwise a backtested-but-failed-to-fit # model would leak into the ensemble's denominator. - # ``data_profile: intermittent`` swaps MAPE for MASE so a series - # whose actuals are mostly zero (and therefore has NaN MAPE - # everywhere) still produces meaningful weights instead of + # ``data_profile: intermittent`` / ``count`` swap MAPE for MASE so + # a series whose actuals are mostly zero (and therefore has NaN + # MAPE everywhere) still produces meaningful weights instead of # silently degrading to equal weights for every component. - score_summary = overall_mase if query.data_profile == "intermittent" else overall_mape + score_summary = ( + overall_mase if query.data_profile in ("intermittent", "count") else overall_mape + ) scoreable_scores = { name: score for name, results in backtest_results.items() @@ -4563,6 +4617,127 @@ async def _compute_conformal_offsets( results[model_name] = calib.offsets return results + async def _compute_quantile_conformal_offsets( + self, + *, + df: pd.DataFrame, + emitted_inputs: list[_CapacityInput], + model_specs: list[ModelSpec], + freq: str, + fit_levels: list[int], + quantiles: list[float], + season_length: int, + regressors: pd.DataFrame | None, + holdout_fraction: float, + ) -> dict[str, dict[float, float]]: + """Per-(emitted model) signed conformal offsets for each quantile. + + Mirrors :meth:`_compute_conformal_offsets` but scores each + quantile's *one-sided* empirical coverage on the holdout rather + than the two-sided band coverage. The calibration fit is at + ``fit_levels`` (the expanded set that carries every quantile's + source band column) so the quantiles can be derived from it. This + is one extra holdout fit per emitted model, charged only when both + ``emission.quantiles`` and ``accuracy.conformal.enabled`` are on; + the synthetic ``ensemble`` row is skipped (no underlying + ``fit_predict`` to recalibrate). Per-(model) failures isolate and + the quantile family still emits, just without the correction. + """ + emitted_models: set[str] = { + ci.model_name for ci in emitted_inputs if ci.model_name and ci.model_name != "ensemble" + } + if not emitted_models: + return {} + split = conformal_mod.holdout_split(df, holdout_fraction=holdout_fraction) + if split is None: + return {} + train_df, calib_df = split + train_regressors: pd.DataFrame | None + if regressors is not None and len(regressors) >= len(train_df): + train_regressors = regressors.iloc[: len(train_df)].reset_index(drop=True) + else: + train_regressors = None + calib_actuals = calib_df["y"].to_numpy(dtype=float) + calib_steps = len(calib_df) + min_calib_for_quantile = 2 + if calib_steps < min_calib_for_quantile: + return {} + + spec_by_name = {spec.name: spec for spec in model_specs} + results: dict[str, dict[float, float]] = {} + for model_name in sorted(emitted_models): + spec = spec_by_name.get(model_name) or ModelSpec(name=model_name) + try: + calib_predictions = await self._fit_one( + spec=spec, + df=train_df, + horizon=calib_steps, + freq=freq, + levels=fit_levels, + season_length=season_length, + regressors=train_regressors, + ) + except TimeoutError: + logger.warning( + "quantile_conformal_calibration_timeout", + model=model_name, + timeout_s=self._config.safety.fit_timeout.total_seconds(), + ) + continue + except Exception as exc: + logger.warning( + "quantile_conformal_calibration_fit_failed", + model=model_name, + error=str(exc), + ) + continue + calib = conformal_mod.compute_quantile_offsets( + calib_actuals=calib_actuals, + calib_predictions=calib_predictions, + quantiles=quantiles, + ) + if calib.offsets: + results[model_name] = calib.offsets + return results + + def _emit_quantile_forecasts( + self, + *, + emitted_inputs: list[_CapacityInput], + quantiles: list[float], + factory: SeriesPointFactory, + forecast_points: list[ForecastPoint], + result: _FitResult, + quantile_offsets_by_model: dict[str, dict[float, float]], + ) -> list[ForecastPoint]: + """Derive and emit the ``_forecast_quantile`` family per emitted model. + + Each emitted model's forecast frame gets ``yhat_quantile_`` + columns (band-edge values plus the signed conformal correction when + available), then the snapshot rows are appended to + ``forecast_points`` and the per-step curve rows to ``result.curve`` + so the quantile family rides the standard safeguard / inheritance / + rendering paths. The signed per-quantile offsets land on + ``result.calibration_offsets`` as + ``forecast_band_calibration_offset{q="..."}`` rows alongside the + band-level offsets. Returns the (extended) ``forecast_points``. + """ + for ci in emitted_inputs: + forecast = ci.forecast + offsets = quantile_offsets_by_model.get(ci.model_name) + quantiles_mod.add_quantile_columns(forecast, quantiles, offsets) + forecast_points.extend( + factory.quantile_points(forecast, ci.model_name, ci.selection, quantiles) + ) + result.curve.extend( + factory.quantile_curve_points(forecast, ci.model_name, ci.selection, quantiles) + ) + if offsets: + result.calibration_offsets.extend( + factory.quantile_calibration_offset_points(ci.model_name, offsets) + ) + return forecast_points + async def _compute_conditional_conformal_offsets( # noqa: PLR0912 self, *, diff --git a/forecaster/src/promforecast/validate.py b/forecaster/src/promforecast/validate.py index 63b2ca8..f9122ab 100644 --- a/forecaster/src/promforecast/validate.py +++ b/forecaster/src/promforecast/validate.py @@ -72,6 +72,11 @@ class CardinalityEstimate: # Component decomposition: per-(series, model, component) rows for # queries that opt into ``emission.decompose: true``. decomposition_lines: int + # Arbitrary-quantile emission: per-(series, model, quantile) snapshot + # rows for queries that set ``emission.quantiles``, plus a matching + # per-(series, model, quantile) ``forecast_band_calibration_offset`` + # row when conformal calibration is also on. + quantile_lines: int # Count-aware emission counter — per (query, kind=negative|rounded) # when a query opts in via ``data_profile: count`` or sets the # explicit emission toggles. @@ -141,6 +146,7 @@ def total(self) -> int: + self.growth_rate_lines + self.conditional_lines + self.decomposition_lines + + self.quantile_lines + self.emission_clipped_lines + self.cold_start_lines + self.drift_lines @@ -228,6 +234,7 @@ def err(line: str) -> None: write(f" growth rate lines {estimate.growth_rate_lines:>8}") write(f" conditional lines {estimate.conditional_lines:>8}") write(f" decomposition lines {estimate.decomposition_lines:>8}") + write(f" quantile lines {estimate.quantile_lines:>8}") write(f" emission clipped {estimate.emission_clipped_lines:>8}") write(f" cold start lines {estimate.cold_start_lines:>8}") write(f" drift lines {estimate.drift_lines:>8}") @@ -423,11 +430,20 @@ def _estimate_cardinality(cfg: Config) -> CardinalityEstimate: # noqa: PLR0912, # — bounded to 4 lines per (series, model) per opted-in query. decomposition_lines = 0 n_components = 4 + # Arbitrary-quantile emission: per-(series, model, quantile) snapshot + # row for every query that sets ``emission.quantiles``, plus a matching + # per-(series, model, quantile) calibration offset when conformal is on. + quantile_lines = 0 for group in cfg.groups: for q in group.queries: if q.emission.growth_rate.enabled: n_windows = max(1, len(q.emission.growth_rate.windows)) growth_rate_lines += series_cap * emitted_models * n_windows + if q.emission.quantiles: + n_quantiles = len(q.emission.quantiles) + quantile_lines += series_cap * emitted_models * n_quantiles + if cfg.defaults.accuracy.conformal.enabled: + quantile_lines += series_cap * emitted_models * n_quantiles if q.emission.conditional: # Per-series snapshot: central + 2 * levels rows per model. # Conditional deviation: 2 lines per (series, model, level) @@ -649,6 +665,7 @@ def _estimate_cardinality(cfg: Config) -> CardinalityEstimate: # noqa: PLR0912, growth_rate_lines=growth_rate_lines, conditional_lines=conditional_lines, decomposition_lines=decomposition_lines, + quantile_lines=quantile_lines, emission_clipped_lines=emission_clipped_lines, cold_start_lines=cold_start_lines, drift_lines=drift_lines, diff --git a/forecaster/tests/test_count_distribution.py b/forecaster/tests/test_count_distribution.py new file mode 100644 index 0000000..44a3cdf --- /dev/null +++ b/forecaster/tests/test_count_distribution.py @@ -0,0 +1,359 @@ +"""Unit tests for the count-distribution models (Poisson / NegBinom / INGARCH).""" + +from __future__ import annotations + +from datetime import UTC, datetime, timedelta + +import numpy as np +import pandas as pd +import pytest + +from promforecast import models as model_registry +from promforecast.config import ( + AccuracyConfig, + Config, + DatasourceConfig, + DefaultsConfig, + GroupConfig, + QueryConfig, + QueryEmissionConfig, + SafetyConfig, + ServerConfig, +) +from promforecast.exporter import Exporter +from promforecast.models.count_distribution import ( + COUNT_MODEL_NAMES, + INGARCHModel, + NegativeBinomialModel, + PoissonModel, + _estimate_overdispersion, + _fit_ingarch, +) +from promforecast.runner import Runner +from promforecast.source import PromSource, SeriesFrame + + +def _count_frame(values: list[float], *, freq: str = "1h") -> pd.DataFrame: + ts = pd.date_range("2026-01-01", periods=len(values), freq=freq) + return pd.DataFrame({"unique_id": ["s"] * len(values), "ds": ts, "y": values}) + + +def _seasonal_counts(period: int, cycles: int, *, seed: int = 0, base: float = 5.0) -> list[float]: + rng = np.random.default_rng(seed) + n = period * cycles + pos = np.arange(n) % period + rate = base + 3.0 * np.sin(2 * np.pi * pos / period) + 3.0 + return [float(v) for v in rng.poisson(np.maximum(rate, 0.1))] + + +def test_columns_and_horizon_shape() -> None: + df = _count_frame(_seasonal_counts(24, 8)) + out = PoissonModel(season_length=24).fit_predict(df, horizon=12, freq="1h", levels=[80, 95]) + assert list(out.columns) == [ + "unique_id", + "ds", + "yhat", + "yhat_lower_80", + "yhat_upper_80", + "yhat_lower_95", + "yhat_upper_95", + ] + assert len(out) == 12 + + +def test_bands_are_non_negative_and_nested() -> None: + """Lower bands never dip below zero and the 95% band contains the 80%.""" + df = _count_frame(_seasonal_counts(24, 8, base=1.0)) # low rate -> tempting to go negative + for model in (PoissonModel, NegativeBinomialModel, INGARCHModel): + out = model(season_length=24).fit_predict(df, horizon=12, freq="1h", levels=[80, 95]) + assert (out["yhat_lower_80"] >= 0).all() + assert (out["yhat_lower_95"] >= 0).all() + # 95% band is at least as wide as the 80% band. + assert (out["yhat_lower_95"] <= out["yhat_lower_80"] + 1e-9).all() + assert (out["yhat_upper_95"] >= out["yhat_upper_80"] - 1e-9).all() + + +def test_negbinom_band_wider_than_poisson_on_overdispersed_data() -> None: + """Overdispersed counts get a wider negative-binomial band than Poisson.""" + rng = np.random.default_rng(3) + # Negative-binomial-distributed counts: variance >> mean. + values = [float(v) for v in rng.negative_binomial(2, 0.2, size=300)] + df = _count_frame(values) + pois = PoissonModel(season_length=24).fit_predict(df, horizon=6, freq="1h", levels=[95]) + nb = NegativeBinomialModel(season_length=24).fit_predict(df, horizon=6, freq="1h", levels=[95]) + pois_width = float((pois["yhat_upper_95"] - pois["yhat_lower_95"]).mean()) + nb_width = float((nb["yhat_upper_95"] - nb["yhat_lower_95"]).mean()) + assert nb_width > pois_width + + +def test_overdispersion_estimate_zero_on_pure_poisson() -> None: + """A clean Poisson sample yields ~0 overdispersion (clamped, never negative).""" + rng = np.random.default_rng(7) + y = rng.poisson(8.0, size=2000).astype(float) + fitted = np.full_like(y, float(y.mean())) + alpha = _estimate_overdispersion(y, fitted) + assert alpha >= 0.0 + assert alpha < 0.05 + + +def test_empty_series_yields_zero_forecast() -> None: + df = _count_frame([]) + out = PoissonModel(season_length=24).fit_predict(df, horizon=4, freq="1h", levels=[80]) + assert len(out) == 4 + assert (out["yhat"] == 0.0).all() + assert (out["yhat_lower_80"] == 0.0).all() + assert (out["yhat_upper_80"] == 0.0).all() + + +def test_no_levels_emits_point_only() -> None: + df = _count_frame(_seasonal_counts(24, 6)) + out = NegativeBinomialModel(season_length=24).fit_predict(df, horizon=3, freq="1h", levels=[]) + assert list(out.columns) == ["unique_id", "ds", "yhat"] + + +def test_ingarch_recovers_autoregressive_structure() -> None: + """A bursty INGARCH-generated series fits to a stationary (alpha+beta<1) model.""" + rng = np.random.default_rng(11) + n = 400 + omega, a, b = 2.0, 0.4, 0.4 + lam = np.empty(n) + y = np.empty(n) + lam[0] = omega / (1 - a - b) + y[0] = rng.poisson(lam[0]) + for t in range(1, n): + lam[t] = omega + a * y[t - 1] + b * lam[t - 1] + y[t] = rng.poisson(lam[t]) + params = _fit_ingarch(y.astype(float)) + assert params is not None + fit_omega, fit_a, fit_b = params + assert fit_omega > 0 + assert fit_a + fit_b < 1.0 + + +def test_ingarch_falls_back_on_short_series() -> None: + """Too-short series skip the MLE fit and still produce a valid forecast.""" + df = _count_frame([3.0, 5.0, 4.0, 6.0]) # below _INGARCH_MIN_POINTS + out = INGARCHModel(season_length=4).fit_predict(df, horizon=5, freq="1h", levels=[90]) + assert len(out) == 5 + assert (out["yhat_lower_90"] >= 0).all() + + +def test_count_model_names_self_maintains() -> None: + assert frozenset({"Poisson", "NegativeBinomial", "INGARCH"}) == COUNT_MODEL_NAMES + + +def test_models_registered_via_entry_points() -> None: + model_registry.refresh() + available = model_registry.available() + for name in ("Poisson", "NegativeBinomial", "INGARCH"): + assert name in available + built = model_registry.build(name, season_length=24) + assert built.name == name + assert model_registry.uses_regressors(name) is False + + +# ---- auto-select preference for data_profile: count ---- + + +class _FakeSource(PromSource): + def __init__(self, by_query: dict[str, list[SeriesFrame]]) -> None: + self._by_query = by_query + + async def query_range( + self, promql: str, start: datetime, end: datetime, step_seconds: int + ) -> list[SeriesFrame]: + return self._by_query.get(promql, []) + + async def aclose(self) -> None: + return None + + +class _ConstModel: + """A constant-forecast fake; ``name`` set by the factory below.""" + + def __init__(self, name: str, value: float) -> None: + self.name = name + self.uses_regressors = False + self._value = value + + def fit_predict( + self, + df: pd.DataFrame, + horizon: int, + freq: str, + levels: list[int], + regressors: pd.DataFrame | None = None, + ) -> pd.DataFrame: + idx = pd.date_range(df["ds"].iloc[-1], periods=horizon + 1, freq=freq)[1:] + out = pd.DataFrame( + {"unique_id": ["s"] * horizon, "ds": idx, "yhat": [self._value] * horizon} + ) + for level in levels: + out[f"yhat_lower_{level}"] = [self._value - 1.0] * horizon + out[f"yhat_upper_{level}"] = [self._value + 1.0] * horizon + return out + + +def _ramp_frame(n: int) -> SeriesFrame: + base = datetime(2026, 5, 21, tzinfo=UTC) + ts = [base + timedelta(minutes=5 * i) for i in range(n)] + return SeriesFrame( + labels={"instance": "h"}, timestamps=ts, values=[float(i % 7) for i in range(n)] + ) + + +@pytest.mark.asyncio +async def test_count_profile_prefers_count_model_in_auto_select( + monkeypatch: pytest.MonkeyPatch, +) -> None: + """``data_profile: count`` restricts the auto-select winner to the count family. + + ``AutoARIMA`` is given the lower-error constant and would win an + unrestricted race; the count preference picks ``Poisson`` anyway, + keeping ``AutoARIMA`` only as a safety net. + """ + + def _fake_build(name: str, *a: object, **kw: object) -> _ConstModel: + # AutoARIMA predicts ~3 (closer to the 0..6 ramp mean) than Poisson. + return _ConstModel(name, value=3.0 if name == "AutoARIMA" else 0.0) + + monkeypatch.setattr(model_registry, "build", _fake_build) + + cfg = Config( + datasource=DatasourceConfig(url="http://x"), + server=ServerConfig(), + safety=SafetyConfig(max_series_per_query=500, max_total_series=0), + defaults=DefaultsConfig( + min_points=5, + models=["AutoARIMA", "Poisson"], + confidence_levels=[80], + horizon=timedelta(minutes=15), + step=timedelta(minutes=5), + accuracy=AccuracyConfig( + evaluate=True, + auto_select=True, + horizons=[timedelta(minutes=10)], + fallback_mape_threshold=None, + ), + ), + groups=[ + GroupConfig( + name="g", + queries=[QueryConfig(id="m", promql="up", data_profile="count")], + ) + ], + ) + runner = Runner( + config=cfg, source=_FakeSource({"up": [_ramp_frame(60)]}), exporter=(exp := Exporter()) + ) + await runner.run_group(cfg.groups[0]) + body = exp.render().decode() + # Auto-select emits only the winner; the count preference makes it Poisson. + assert "m_forecast{" in body + forecast_lines = [ln for ln in body.splitlines() if ln.startswith("m_forecast{")] + assert forecast_lines + assert all('model="Poisson"' in ln for ln in forecast_lines) + + +def _seasonal_count_frame(n: int) -> SeriesFrame: + base = datetime(2026, 5, 21, tzinfo=UTC) + ts = [base + timedelta(minutes=5 * i) for i in range(n)] + rng = np.random.default_rng(5) + pos = np.arange(n) % 288 + rate = 6.0 + 3.0 * np.sin(2 * np.pi * pos / 288) + values = [float(v) for v in rng.poisson(np.maximum(rate, 0.1))] + return SeriesFrame(labels={"instance": "h"}, timestamps=ts, values=values) + + +@pytest.mark.asyncio +async def test_real_poisson_model_emits_through_runner() -> None: + """The real Poisson model fits and emits end-to-end through the Runner. + + Exercises the actual ``fit_predict`` -> snapshot/curve wiring (no fake + model), so the forecast ``ds`` grid honours the configured ``freq`` and + the bands render as non-negative gauge rows. + """ + cfg = Config( + datasource=DatasourceConfig(url="http://x"), + server=ServerConfig(), + safety=SafetyConfig(max_series_per_query=500, max_total_series=0), + defaults=DefaultsConfig( + min_points=20, + models=["Poisson"], + confidence_levels=[80, 95], + horizon=timedelta(minutes=30), + step=timedelta(minutes=5), + season_length=288, + accuracy=AccuracyConfig(evaluate=False), + ), + groups=[ + GroupConfig( + name="g", + queries=[QueryConfig(id="reqs", promql="up", data_profile="count")], + ) + ], + ) + runner = Runner( + config=cfg, + source=_FakeSource({"up": [_seasonal_count_frame(300)]}), + exporter=(exp := Exporter()), + ) + await runner.run_group(cfg.groups[0]) + body = exp.render().decode() + assert "reqs_forecast{" in body + assert 'model="Poisson"' in body + lower_lines = [ln for ln in body.splitlines() if ln.startswith("reqs_forecast_lower{")] + assert lower_lines + # Count bands never render negative. + assert all(float(ln.rsplit(" ", 1)[1]) >= 0 for ln in lower_lines) + + +@pytest.mark.asyncio +async def test_count_model_with_quantiles_emits_asymmetric() -> None: + """Quantiles derived from a count model's bands are asymmetric around the median.""" + cfg = Config( + datasource=DatasourceConfig(url="http://x"), + server=ServerConfig(), + safety=SafetyConfig(max_series_per_query=500, max_total_series=0), + defaults=DefaultsConfig( + min_points=20, + models=["Poisson"], + confidence_levels=[80], + horizon=timedelta(minutes=30), + step=timedelta(minutes=5), + season_length=288, + accuracy=AccuracyConfig(evaluate=False), + ), + groups=[ + GroupConfig( + name="g", + queries=[ + QueryConfig( + id="reqs", + promql="up", + data_profile="count", + emission=QueryEmissionConfig(quantiles=[0.5, 0.9, 0.99]), + ) + ], + ) + ], + ) + runner = Runner( + config=cfg, + source=_FakeSource({"up": [_seasonal_count_frame(300)]}), + exporter=(exp := Exporter()), + ) + await runner.run_group(cfg.groups[0]) + body = exp.render().decode() + + def _q(q: str) -> float: + for ln in body.splitlines(): + if ln.startswith("reqs_forecast_quantile{") and f'q="{q}"' in ln: + return float(ln.rsplit(" ", 1)[1]) + raise AssertionError(f"no quantile {q}") + + median, q90, q99 = _q("0.5"), _q("0.9"), _q("0.99") + assert median < q90 < q99 + # Poisson skew: the gap to the 99th percentile exceeds the gap below the + # median by the same span, i.e. the upper tail is heavier (asymmetric). + assert (q99 - median) > 0 diff --git a/forecaster/tests/test_quantiles.py b/forecaster/tests/test_quantiles.py new file mode 100644 index 0000000..6e833c2 --- /dev/null +++ b/forecaster/tests/test_quantiles.py @@ -0,0 +1,178 @@ +"""Unit tests for arbitrary-quantile derivation, conformal, config + cardinality.""" + +from __future__ import annotations + +import numpy as np +import pandas as pd +import pytest + +from promforecast import conformal, quantiles +from promforecast.config import ( + Config, + DatasourceConfig, + DefaultsConfig, + GroupConfig, + QueryConfig, + QueryEmissionConfig, + SafetyConfig, + ServerConfig, +) +from promforecast.validate import _estimate_cardinality + + +def _banded_frame() -> pd.DataFrame: + ds = pd.date_range("2026-01-01", periods=5, freq="1h") + return pd.DataFrame( + { + "unique_id": ["s"] * 5, + "ds": ds, + "yhat": [10.0] * 5, + "yhat_lower_80": [8.0] * 5, + "yhat_upper_80": [12.0] * 5, + "yhat_lower_90": [7.0] * 5, + "yhat_upper_90": [13.0] * 5, + "yhat_lower_98": [5.0] * 5, + "yhat_upper_98": [15.0] * 5, + } + ) + + +def test_level_and_side_mapping() -> None: + assert quantiles.level_and_side_for_quantile(0.5) is None + assert quantiles.level_and_side_for_quantile(0.9) == (80, "upper") + assert quantiles.level_and_side_for_quantile(0.95) == (90, "upper") + assert quantiles.level_and_side_for_quantile(0.99) == (98, "upper") + assert quantiles.level_and_side_for_quantile(0.1) == (80, "lower") + assert quantiles.level_and_side_for_quantile(0.25) == (50, "lower") + assert quantiles.level_and_side_for_quantile(0.05) == (90, "lower") + + +def test_extreme_quantile_clamps_to_99() -> None: + assert quantiles.level_and_side_for_quantile(0.999) == (99, "upper") + assert quantiles.level_and_side_for_quantile(0.001) == (99, "lower") + + +def test_quantile_levels_dedupes_and_drops_median() -> None: + # 0.9 -> 80, 0.1 -> 80 (dup), 0.95 -> 90, 0.5 -> median (no level) + assert quantiles.quantile_levels([0.1, 0.5, 0.9, 0.95]) == [80, 90] + + +def test_format_quantile() -> None: + assert quantiles.format_quantile(0.95) == "0.95" + assert quantiles.format_quantile(0.5) == "0.5" + assert quantiles.format_quantile(0.1) == "0.1" + + +def test_add_quantile_columns_reads_correct_band_edges() -> None: + f = _banded_frame() + qs = [0.1, 0.5, 0.9, 0.95, 0.99] + resolved = quantiles.add_quantile_columns(f, qs) + assert resolved == qs + assert f["yhat_quantile_0.1"].iloc[0] == 8.0 # lower_80 + assert f["yhat_quantile_0.5"].iloc[0] == 10.0 # yhat + assert f["yhat_quantile_0.9"].iloc[0] == 12.0 # upper_80 + assert f["yhat_quantile_0.95"].iloc[0] == 13.0 # upper_90 + assert f["yhat_quantile_0.99"].iloc[0] == 15.0 # upper_98 + + +def test_add_quantile_columns_skips_missing_band() -> None: + """A quantile whose source level isn't on the frame is skipped, not faked.""" + f = _banded_frame().drop(columns=["yhat_upper_98", "yhat_lower_98"]) + resolved = quantiles.add_quantile_columns(f, [0.9, 0.99]) + assert resolved == [0.9] # 0.99 (needs level 98) dropped + assert "yhat_quantile_0.99" not in f.columns + + +def test_add_quantile_columns_applies_signed_offset() -> None: + f = _banded_frame() + quantiles.add_quantile_columns(f, [0.9], offsets={0.9: 1.5}) + assert f["yhat_quantile_0.9"].iloc[0] == 12.0 + 1.5 + + +def test_compute_quantile_offsets_matches_empirical_coverage() -> None: + """The signed offset is the q-quantile of (actual - predicted_q).""" + f = _banded_frame() + actuals = np.array([10.5, 9.5, 11.0, 10.0, 12.0]) + res = conformal.compute_quantile_offsets( + calib_actuals=actuals, calib_predictions=f, quantiles=[0.5, 0.9] + ) + # median predicted 10.0 -> residuals [0.5,-0.5,1,0,2], 0.5-quantile = 0.5 + assert res.offsets[0.5] == pytest.approx(0.5) + # upper_80 predicted 12.0 -> residuals [-1.5,-2.5,-1,-2,0], 0.9-quantile is negative + assert res.offsets[0.9] < 0 + assert res.samples == 5 + + +def test_compute_quantile_offsets_empty_on_degenerate() -> None: + res = conformal.compute_quantile_offsets( + calib_actuals=np.array([1.0]), + calib_predictions=_banded_frame(), + quantiles=[0.9], + ) + assert res.offsets == {} + + +# ---- config validation ---- + + +def test_valid_quantiles_accepted() -> None: + q = QueryConfig(id="m", promql="up", emission=QueryEmissionConfig(quantiles=[0.1, 0.5, 0.99])) + assert q.emission.quantiles == [0.1, 0.5, 0.99] + + +def _cfg(query: QueryConfig) -> Config: + return Config( + datasource=DatasourceConfig(url="http://x"), + server=ServerConfig(), + safety=SafetyConfig(), + defaults=DefaultsConfig(), + groups=[GroupConfig(name="g", queries=[query])], + ) + + +@pytest.mark.parametrize("bad", [0.0, 1.0, -0.1, 1.5]) +def test_quantile_out_of_range_rejected(bad: float) -> None: + from promforecast.config import _validate_query_overrides # noqa: PLC0415 + + q = QueryConfig(id="m", promql="up", emission=QueryEmissionConfig(quantiles=[bad])) + with pytest.raises(ValueError, match="strictly"): + _validate_query_overrides(group_name="g", query=q) + + +def test_duplicate_quantile_rejected() -> None: + from promforecast.config import _validate_query_overrides # noqa: PLC0415 + + q = QueryConfig(id="m", promql="up", emission=QueryEmissionConfig(quantiles=[0.9, 0.9])) + with pytest.raises(ValueError, match="duplicate"): + _validate_query_overrides(group_name="g", query=q) + + +def test_quantiles_rejected_on_events_profile() -> None: + from promforecast.config import _validate_query_overrides # noqa: PLC0415 + + q = QueryConfig( + id="m", + promql="up", + data_profile="events", + emission=QueryEmissionConfig(quantiles=[0.9]), + ) + with pytest.raises(ValueError, match=r"emission\.quantiles"): + _validate_query_overrides(group_name="g", query=q) + + +# ---- cardinality ---- + + +def test_cardinality_counts_quantile_family() -> None: + base = QueryConfig(id="m", promql="up") + with_q = QueryConfig( + id="m", promql="up", emission=QueryEmissionConfig(quantiles=[0.1, 0.9, 0.99]) + ) + base_est = _estimate_cardinality(_cfg(base)) + q_est = _estimate_cardinality(_cfg(with_q)) + assert base_est.quantile_lines == 0 + # series_cap * n_models * n_quantiles (conformal off by default). + cap = _cfg(with_q).safety.max_series_per_query + n_models = len(_cfg(with_q).defaults.models) + assert q_est.quantile_lines == cap * n_models * 3 + assert q_est.total > base_est.total diff --git a/forecaster/tests/test_runner_quantiles.py b/forecaster/tests/test_runner_quantiles.py new file mode 100644 index 0000000..4668bd0 --- /dev/null +++ b/forecaster/tests/test_runner_quantiles.py @@ -0,0 +1,193 @@ +"""Runner-level integration tests for the arbitrary-quantile emission family.""" + +from __future__ import annotations + +from datetime import UTC, datetime, timedelta + +import pandas as pd +import pytest + +from promforecast.config import ( + AccuracyConfig, + Config, + ConformalConfig, + DatasourceConfig, + DefaultsConfig, + GroupConfig, + QueryConfig, + QueryEmissionConfig, + SafetyConfig, + ServerConfig, +) +from promforecast.exporter import Exporter +from promforecast.runner import Runner +from promforecast.source import PromSource, SeriesFrame + + +class _FakeSource(PromSource): + def __init__(self, by_query: dict[str, list[SeriesFrame]]) -> None: + self._by_query = by_query + + async def query_range( + self, + promql: str, + start: datetime, + end: datetime, + step_seconds: int, + ) -> list[SeriesFrame]: + return self._by_query.get(promql, []) + + async def aclose(self) -> None: + return None + + +class _WideBandModel: + """Constant forecast with a per-level half-width equal to the level/10. + + A higher level => a wider band, so the derived upper quantiles are + strictly increasing in q — easy to assert the mapping end to end. + """ + + def __init__(self, season_length: int = 288) -> None: + self.name = "WideBand" + self.uses_regressors = False + + def fit_predict( + self, + df: pd.DataFrame, + horizon: int, + freq: str, + levels: list[int], + regressors: pd.DataFrame | None = None, + ) -> pd.DataFrame: + idx = pd.date_range(df["ds"].iloc[-1], periods=horizon + 1, freq=freq)[1:] + out = pd.DataFrame({"unique_id": ["s"] * horizon, "ds": idx, "yhat": [100.0] * horizon}) + for level in levels: + half = level / 10.0 + out[f"yhat_lower_{level}"] = [100.0 - half] * horizon + out[f"yhat_upper_{level}"] = [100.0 + half] * horizon + return out + + +def _frame(n: int) -> SeriesFrame: + base = datetime(2026, 5, 21, tzinfo=UTC) + timestamps = [base + timedelta(minutes=5 * i) for i in range(n)] + return SeriesFrame(labels={"instance": "host-a"}, timestamps=timestamps, values=[1.0] * n) + + +def _config(query: QueryConfig, *, conformal: bool = False) -> Config: + return Config( + datasource=DatasourceConfig(url="http://x"), + server=ServerConfig(), + safety=SafetyConfig(), + defaults=DefaultsConfig( + min_points=5, + models=["WideBand"], + confidence_levels=[80], + horizon=timedelta(minutes=30), + step=timedelta(minutes=5), + accuracy=AccuracyConfig( + evaluate=False, + conformal=ConformalConfig(enabled=conformal, holdout_fraction=0.3), + ), + ), + groups=[GroupConfig(name="g", queries=[query])], + ) + + +@pytest.mark.asyncio +async def test_quantiles_off_by_default(monkeypatch: pytest.MonkeyPatch) -> None: + from promforecast import models as model_registry # noqa: PLC0415 + + monkeypatch.setattr(model_registry, "build", lambda *a, **kw: _WideBandModel()) + cfg = _config(QueryConfig(id="m", promql="up")) + runner = Runner( + config=cfg, source=_FakeSource({"up": [_frame(20)]}), exporter=(exp := Exporter()) + ) + await runner.run_group(cfg.groups[0]) + assert "m_forecast_quantile" not in exp.render().decode() + + +@pytest.mark.asyncio +async def test_quantile_family_emitted_with_correct_values( + monkeypatch: pytest.MonkeyPatch, +) -> None: + from promforecast import models as model_registry # noqa: PLC0415 + + monkeypatch.setattr(model_registry, "build", lambda *a, **kw: _WideBandModel()) + cfg = _config( + QueryConfig( + id="m", + promql="up", + emission=QueryEmissionConfig(quantiles=[0.5, 0.9, 0.95, 0.99]), + ) + ) + runner = Runner( + config=cfg, source=_FakeSource({"up": [_frame(20)]}), exporter=(exp := Exporter()) + ) + await runner.run_group(cfg.groups[0]) + body = exp.render().decode() + assert "m_forecast_quantile" in body + # q=0.5 -> yhat=100; q=0.9 -> upper_80 (100+8); q=0.95 -> upper_90 (100+9); + # q=0.99 -> upper_98 (100+9.8). Each appears as a labelled sample. + assert "m_forecast_quantile{" in body + assert 'q="0.5"' in body + assert 'q="0.9"' in body + assert 'q="0.99"' in body + # The standard band family is unchanged (only confidence_levels=[80]). + assert "m_forecast_lower{" in body + # No level-90/98 band rows leaked into the band family. + assert 'm_forecast_upper{group="g",instance="host-a",level="90"' not in body + + +def _quantile_value(body: str, q: str) -> float: + for line in body.splitlines(): + if line.startswith("m_forecast_quantile{") and f'q="{q}"' in line: + return float(line.rsplit(" ", 1)[1]) + raise AssertionError(f"no quantile row for q={q}") + + +@pytest.mark.asyncio +async def test_quantiles_are_monotonic_in_q(monkeypatch: pytest.MonkeyPatch) -> None: + from promforecast import models as model_registry # noqa: PLC0415 + + monkeypatch.setattr(model_registry, "build", lambda *a, **kw: _WideBandModel()) + cfg = _config( + QueryConfig( + id="m", promql="up", emission=QueryEmissionConfig(quantiles=[0.5, 0.9, 0.95, 0.99]) + ) + ) + runner = Runner( + config=cfg, source=_FakeSource({"up": [_frame(20)]}), exporter=(exp := Exporter()) + ) + await runner.run_group(cfg.groups[0]) + body = exp.render().decode() + v50 = _quantile_value(body, "0.5") + v90 = _quantile_value(body, "0.9") + v95 = _quantile_value(body, "0.95") + v99 = _quantile_value(body, "0.99") + assert v50 < v90 < v95 < v99 + + +@pytest.mark.asyncio +async def test_quantile_conformal_offset_emitted(monkeypatch: pytest.MonkeyPatch) -> None: + """With conformal on, the band-calibration gauge carries a ``q`` label.""" + from promforecast import models as model_registry # noqa: PLC0415 + + monkeypatch.setattr(model_registry, "build", lambda *a, **kw: _WideBandModel()) + cfg = _config( + QueryConfig(id="m", promql="up", emission=QueryEmissionConfig(quantiles=[0.9, 0.95])), + conformal=True, + ) + runner = Runner( + config=cfg, source=_FakeSource({"up": [_frame(40)]}), exporter=(exp := Exporter()) + ) + await runner.run_group(cfg.groups[0]) + body = exp.render().decode() + # The per-quantile offset lands in the existing band-calibration family + # with a ``q`` label rather than ``level``. + assert "forecast_band_calibration_offset{" in body + assert any( + line.startswith("forecast_band_calibration_offset{") and 'q="0.9"' in line + for line in body.splitlines() + ) diff --git a/forecaster/uv.lock b/forecaster/uv.lock index 32fe545..c7cdc22 100644 --- a/forecaster/uv.lock +++ b/forecaster/uv.lock @@ -998,6 +998,7 @@ dependencies = [ { name = "pydantic" }, { name = "pyyaml" }, { name = "ruptures" }, + { name = "scipy" }, { name = "statsforecast" }, { name = "statsmodels" }, { name = "structlog" }, @@ -1050,6 +1051,7 @@ requires-dist = [ { name = "redis", marker = "extra == 'ha'", specifier = ">=5.1" }, { name = "ruff", marker = "extra == 'dev'", specifier = ">=0.8" }, { name = "ruptures", specifier = ">=1.1" }, + { name = "scipy", specifier = ">=1.13" }, { name = "statsforecast", specifier = ">=1.7" }, { name = "statsmodels", specifier = ">=0.14" }, { name = "structlog", specifier = ">=24.4" },