Run Inspection, Diagnostics, Scenarios & Reports

review

Monitor a running run and review a finished one — the run-inspection layer, keyed on a run_id.

Two questions, two entry points, both taking only a run_id (the run is its config, so this layer reads the run's own raw_config back to recover what it planned to do):

• monitor_run — how far along is a run in flight? Per-family job state, series done vs. the expected total, remaining, and mean fit time per family on its chosen runner. Progress is coarse: the registry has no live per-series counter — cells land when a family's write_cells runs (often at job end), so done-counts step up per job. The per-job status (from v_run_jobs) is the primary live signal; landed-cell counts refine it.

A registry row is written by the job, so a job that dies without writing leaves its row RUNNING forever and the bar simply stops moving. Two things keep that legible. Every family carries quiet_seconds — how long since its last registry signal — which is derived from rows already read, costs no runtime call, and is a fact rather than a judgement (a family that writes its cells at job end is legitimately quiet for its whole run, so a threshold here would cry wolf). And probe=True escalates the non-terminal families to their runtime via probes.reconcile.probe_run's reader, attaching a probes.reconcile.ProbeReport that says whether the job is actually still alive. Registry-first is the default deliberately: a fleet poll must never fan native calls, so escalation stays the deliberate per-run drill-down. - review_run — how did a finished run do, in data-science detail? The best model per family and overall, the full metric panel aggregated across every series (mean + p10/p50/p90), and each ensemble's lift over the best base model. - calibration_report — should you believe the numbers review_run just showed you? Whether the point-forecast correction earned its place on this run's own held-out folds, per model, and whether the prediction interval achieves the coverage it claims — broken out by horizon step, because a band that averages to nominal can still be wrong at both ends.

Same pure/I-O seam as sdk: the _assemble_* functions are pure (turn reader dicts into the result dataclasses, unit-tested offline), while monitor_run / review_run are the thin I/O callers that read the registry (registry.reads, registry.jobs) and hand off. Plotting (plot_progress, plot_leaderboard, plot_metric_distribution) is a convenience over the dataclasses, with matplotlib imported lazily so it never touches the near-instant import scale_forecasting path. For the wall-clock execution timeline, reuse sdk.build_trace_frame + sdk.plot_trace.

ArmComparison dataclass

ArmComparison(
    model_type,
    compute_engine,
    interval_calibration,
    n_series,
    n_raw_arm,
    n_auto_decided,
    n_compared,
    n_corrected_wins,
    mean_margin,
    median_margin,
)

One model's point-forecast arm choice on a finished run, and what the choice was worth.

margin is always "how much the corrected arm beat the raw one by", as a fraction of the raw arm's loss in the run's decision_metric, whichever arm the run actually shipped. Negative means the correction lost. win_rate is the share of compared series where it won — the number that matters more than the average, because a correction that helps 51% of series by a lot and hurts 49% by a lot is a different proposition from one that helps everything a little.

n_raw_arm and n_auto_decided are counts rather than a single label because under output.point_forecast="auto" the series of one model do not have to agree, and that disagreement is the feature — reporting one of them as if it spoke for all would hide it.

raw_arm_rate property

raw_arm_rate

Share of this model's series shipping the raw arm, or None when there are none.

1.0 or 0.0 under a fleetwide setting; anything between means selection actually split the model's series, which is the thing auto exists to do and the thing worth looking at.

win_rate property

win_rate

Share of compared series the corrected arm won, or None when nothing was compared.

BacktestCohort dataclass

BacktestCohort(
    n_series=0,
    n_full=0,
    n_reduced=0,
    n_unscored=0,
    n_failed=0,
    n_not_requested=0,
    fold_histogram=dict(),
    refit_modes=dict(),
    staleness_gap=None,
)

How much of the panel one model was actually scored on, split by outcome.

The context a leaderboard score is meaningless without. A ragged panel does not give every series the same backtest: backtest.make_folds drops the folds a short series cannot afford, so full (every requested fold), reduced (some), unscored (none — too short to backtest at all) and failed are all ordinary outcomes within a single run, and a model whose mean error came mostly off reduced series won an easier contest than its neighbour. n_not_requested counts series whose backtest_status is NULL, which means the run never asked for a backtest rather than that one was attempted and produced nothing.

fold_histogram maps achieved fold count → series count, so the shape of the raggedness is visible and not just its worst case. Keyed by the achieved count as an int; series with a NULL n_folds_achieved are left out of it (they are still counted in n_series).

refit_modes is the other axis of the same question — not how much of the panel was scored, but on what. It maps backtest_refit → series count: per_fold (a fresh fit at every origin), recondition (one fit carried forward on the new observations), extrapolate (one fit, never told what happened next), and unsupported (a frozen scheme was asked for and this model had no seam for it, so those series refit anyway). A model showing unsupported on a frozen run is not answering the same question as its neighbours. Ensemble cohorts inherit their members' mode, or mixed where the members disagreed — see ensemble_run.ensemble_refit_mode.

staleness_gap is what never refreshing the model cost this panel, in the run's decision metric, averaged over the series that ran a control arm and weighted by cohort size. Positive means refitting earns its keep. None on the two refit schemes, which run no control arm.

CalibrationReport dataclass

CalibrationReport(
    run_id,
    decision_metric,
    nominal_coverage,
    arms,
    coverage,
)

What the run's own data says about its point forecast and its prediction interval.

Two questions a forecaster asks of any system that corrects a model's output, answered from the run's held-out folds rather than from a claim in a doc:

• Was the correction worth applying? arms, per model — win rate and margin.
• Does the band mean what it says? coverage per horizon step against nominal, which is the width of the run's quantile set (0.8 for the shipped default of 0.1/0.5/0.9).

worst_step is the single largest deviation from nominal across every model and step: a run whose average coverage looks fine while one end of the horizon is badly wrong should not read as healthy, and an average will always say it does.

mean_coverage property

mean_coverage

Row-count-weighted achieved coverage across every model and step.

worst_step property

worst_step

The model/step furthest from nominal — the thing an average is designed to hide.

CoveragePoint dataclass

CoveragePoint(
    model_type, horizon_step, n, coverage, mean_width
)

Achieved interval coverage at one horizon step for one model.

EnsembleLift dataclass

EnsembleLift(
    model_type,
    score,
    best_base_model,
    best_base_score,
    lift,
    lift_pct,
)

How much an ensemble improved on the best base model, in the run's decision metric.

lift is best_base_score − score (positive = the ensemble is better, since lower error is better); lift_pct is that as a fraction of the base score. Compares against the single best base model overall (best_overall), the bar an ensemble has to clear to be worth keeping.

FamilyProgress dataclass

FamilyProgress(
    family,
    runtime,
    hardware,
    status,
    models,
    n_expected,
    n_done,
    fraction,
    avg_fit_seconds,
    runtime_seconds,
    last_signal_at=None,
    quiet_seconds=None,
    device_verdict=None,
)

One family's live progress on a run: its job state and how many series it has scored.

runtime / hardware name the runner this family resolved to (a Spark family's mean fit time is only comparable to another family on the same runner). n_expected is n_series × models-in-family (None when the run's series count isn't known); n_done counts landed full-fit cells; fraction is their ratio. avg_fit_seconds is the mean per-cell fit time across this family's landed cells; runtime_seconds is the job's wall-clock once it finishes.

last_signal_at is the most recent timestamp the job row carries (ended_at → started_at → created_at, first present wins) and quiet_seconds is its age at read time — both None for a family with no job row yet, or an unparseable timestamp. They are reported, never judged: how long a family may legitimately stay quiet depends on the family, and the escalation threshold that does judge it lives with the probe (probes.reconcile._DEFAULT_STALE_S), not here.

device_verdict is device_audit's answer to whether the accelerator this family paid for did any work — None for every CPU family, and for a GPU family until its job finishes, because the audit runs on the driver once the cells are written. It is the only field here that is about cost rather than progress, and it is on the progress object because this is where a reader is already looking at the per-family runtime and hardware.

ModelReview dataclass

ModelReview(
    model_type,
    family,
    ensemble_id,
    is_ensemble,
    compute_engine,
    n_series,
    score,
    metric_means=dict(),
    metric_p10=dict(),
    metric_p50=dict(),
    metric_p90=dict(),
    mean_fit_seconds=None,
    median_fit_seconds=None,
    no_artifact_rate=None,
    n_predictions=0,
    cohort=None,
    pooled_wape=None,
    n_comparable_series=None,
    base_model=None,
    scope_capability=None,
    training_mode=None,
    cov_profile=None,
    is_univariate=None,
    covariates_used=None,
    model_spec_hash=None,
    borrowed_from_run_id=None,
    n_borrowed=0,
)

One model's (or ensemble pseudo-model's) outcome on a finished run, across all its series.

score is the mean of the run's decision_metric over every series (lower = better; None when no backtest scored it). metric_means / metric_p10 / metric_p50 / metric_p90 carry the full metric panel — the cross-series mean and the 10th/50th/90th percentile of each metric in METRIC_COLUMNS — so distribution shape reads off the aggregates without pulling per-series rows. ensemble_id is None for a base model; is_ensemble is its convenience flag. n_predictions is the forecast-row count (0 flags a model that scored metadata but produced no forecasts — a fully-failed fit).

RunProgress dataclass

RunProgress(
    run_id,
    status,
    n_series,
    families,
    n_done,
    n_expected,
    fraction,
    probe=None,
)

A run's live progress snapshot: header status plus one FamilyProgress per family.

status is None when no run exists for the id yet. n_done / n_expected / fraction are the run-wide roll-up across families (n_expected and fraction are None when the series count — hence the denominator — isn't known).

probe carries the reconciled probes.reconcile.ProbeReport when monitor_run was called with probe=True, and is None for the default registry-only read — so a caller can always tell "the runtime agreed the job is alive" apart from "we never asked".

RunReview dataclass

RunReview(
    run_id,
    status,
    decision_metric,
    n_series,
    models,
    best_per_family,
    best_overall,
    ensembles,
    ensemble_lift,
)

A finished run's data-science review: the leaderboard plus derived bests and ensemble lift.

models is every model best-first (lowest score). best_per_family maps each base family to its champion; best_overall is the single best base model; ensembles are the ensemble pseudo-models; ensemble_lift scores each against best_overall.

best_overall

best_overall(models)

The single best base model (lowest score); None if no base model was scored.

Source code in src/scale_forecasting/review.py

def best_overall(models: list[ModelReview] | tuple[ModelReview, ...]) -> ModelReview | None:
    """The single best base model (lowest ``score``); ``None`` if no base model was scored."""
    scored = [m for m in models if not m.is_ensemble and m.score is not None]
    return min(scored, key=lambda m: _sort_key(m.score)) if scored else None

best_per_family

best_per_family(models)

Map each base family to its champion (lowest-score scored model in that family).

Source code in src/scale_forecasting/review.py

def best_per_family(
    models: list[ModelReview] | tuple[ModelReview, ...],
) -> dict[str, ModelReview]:
    """Map each base family to its champion (lowest-``score`` scored model in that family)."""
    best: dict[str, ModelReview] = {}
    for m in models:
        if m.is_ensemble or m.score is None:
            continue
        cur = best.get(m.family)
        if cur is None or cur.score is None or m.score < cur.score:
            best[m.family] = m
    return best

build_attributions_frame

build_attributions_frame(rows, *, level='global')

Unpack Tier 1 (level="global") or Tier 2 (level="local") attribution rows into a tidy DataFrame (pure, offline).

• level="global" unpacks fit_diagnostics["feature_attributions"] into columns ["ts_id", "model_type", "compute_engine", "feature", "importance"] sorted by ts_id, model_type, importance DESC.
• level="local" unpacks per-horizon-step explanations JSON into columns ["ts_id", "model_type", "compute_engine", "forecast_date", "yhat", "baseline_score", "feature", "attribution"] sorted by ts_id, model_type, forecast_date.

Source code in src/scale_forecasting/review.py

def build_attributions_frame(
    rows: list[dict[str, Any]],
    *,
    level: str = "global",
) -> pd.DataFrame:
    """Unpack Tier 1 (``level="global"``) or Tier 2 (``level="local"``) attribution rows into a
    tidy ``DataFrame`` (pure, offline).

    - ``level="global"`` unpacks ``fit_diagnostics["feature_attributions"]`` into columns
      ``["ts_id", "model_type", "compute_engine", "feature", "importance"]`` sorted by
      ``ts_id, model_type, importance DESC``.
    - ``level="local"`` unpacks per-horizon-step ``explanations`` JSON into columns
      ``["ts_id", "model_type", "compute_engine", "forecast_date", "yhat", "baseline_score",
      "feature", "attribution"]`` sorted by ``ts_id, model_type, forecast_date``.
    """
    import json

    import pandas as pd

    if level == "local":
        cols = [
            "ts_id",
            "model_type",
            "compute_engine",
            "forecast_date",
            "yhat",
            "baseline_score",
            "feature",
            "attribution",
        ]
        if not rows:
            return pd.DataFrame(columns=cols)
        records: list[dict[str, Any]] = []
        for r in rows:
            raw_exp = r.get("explanations")
            if isinstance(raw_exp, str) and raw_exp:
                try:
                    parsed = json.loads(raw_exp)
                except ValueError:
                    continue
            elif isinstance(raw_exp, dict):
                parsed = raw_exp
            else:
                continue
            attrs = parsed.get("attributions") if isinstance(parsed, dict) else None
            if not isinstance(attrs, dict) or not attrs:
                continue
            baseline = parsed.get("baseline_score", 0.0)
            baseline_val = (
                float(baseline)
                if isinstance(baseline, (int, float)) and math.isfinite(float(baseline))
                else 0.0
            )
            fdate_raw = _first_date(r, "forecast_date")
            if fdate_raw is None:
                continue
            fdate = fdate_raw.tz_localize(None)
            yhat_raw = r.get("yhat")
            yhat_val = (
                float(yhat_raw)
                if isinstance(yhat_raw, (int, float)) and math.isfinite(float(yhat_raw))
                else float("nan")
            )
            for feat, val in sorted(
                attrs.items(),
                key=lambda kv: -abs(float(kv[1])) if isinstance(kv[1], (int, float)) else 0.0,
            ):
                if isinstance(val, (int, float)) and math.isfinite(float(val)):
                    records.append(
                        {
                            "ts_id": str(r.get("ts_id") or ""),
                            "model_type": str(r.get("model_type") or ""),
                            "compute_engine": r.get("compute_engine"),
                            "forecast_date": fdate,
                            "yhat": yhat_val,
                            "baseline_score": baseline_val,
                            "feature": str(feat),
                            "attribution": float(val),
                        }
                    )
        if not records:
            return pd.DataFrame(columns=cols)
        return pd.DataFrame.from_records(records, columns=cols)

    cols = ["ts_id", "model_type", "compute_engine", "feature", "importance"]
    if not rows:
        return pd.DataFrame(columns=cols)
    records = []
    for r in rows:
        raw_diag = r.get("fit_diagnostics")
        if isinstance(raw_diag, str) and raw_diag:
            try:
                parsed = json.loads(raw_diag)
            except ValueError:
                continue
        elif isinstance(raw_diag, dict):
            parsed = raw_diag
        else:
            continue
        f_attrs = parsed.get("feature_attributions") if isinstance(parsed, dict) else None
        if not isinstance(f_attrs, dict) or not f_attrs:
            continue
        for feat, val in sorted(
            f_attrs.items(),
            key=lambda kv: -float(kv[1]) if isinstance(kv[1], (int, float)) else 0.0,
        ):
            if isinstance(val, (int, float)) and math.isfinite(float(val)):
                records.append(
                    {
                        "ts_id": str(r.get("ts_id") or ""),
                        "model_type": str(r.get("model_type") or ""),
                        "compute_engine": r.get("compute_engine"),
                        "feature": str(feat),
                        "importance": float(val),
                    }
                )
    if not records:
        return pd.DataFrame(columns=cols)
    return pd.DataFrame.from_records(records, columns=cols)

build_calibration_frames

build_calibration_frames(report)

Convert a CalibrationReport into (arms_df, coverage_df) pandas DataFrames (pure).

• arms_df summarizes point-forecast bias-correction arm selection (raw vs corrected win rates and relative margins per model).
• coverage_df reports empirical prediction-interval coverage and mean interval width at each horizon step ($h = 1 \dots H$) against report.nominal_coverage.

Source code in src/scale_forecasting/review.py

def build_calibration_frames(
    report: CalibrationReport,
) -> tuple[pd.DataFrame, pd.DataFrame]:
    """Convert a `CalibrationReport` into ``(arms_df, coverage_df)`` pandas ``DataFrame``s (pure).

    * ``arms_df`` summarizes point-forecast bias-correction arm selection (`raw` vs `corrected`
      win rates and relative margins per model).
    * ``coverage_df`` reports empirical prediction-interval coverage and mean interval width at
      each horizon step ($h = 1 \\dots H$) against `report.nominal_coverage`.
    """
    import pandas as pd

    arm_cols = [
        "model_type",
        "compute_engine",
        "interval_calibration",
        "n_series",
        "n_raw_arm",
        "raw_arm_rate",
        "n_auto_decided",
        "n_compared",
        "n_corrected_wins",
        "win_rate",
        "mean_margin",
        "median_margin",
    ]
    cov_cols = [
        "model_type",
        "horizon_step",
        "n",
        "coverage",
        "nominal_coverage",
        "coverage_error",
        "mean_width",
    ]
    arm_records = [
        {
            "model_type": a.model_type,
            "compute_engine": a.compute_engine,
            "interval_calibration": a.interval_calibration,
            "n_series": a.n_series,
            "n_raw_arm": a.n_raw_arm,
            "raw_arm_rate": a.raw_arm_rate,
            "n_auto_decided": a.n_auto_decided,
            "n_compared": a.n_compared,
            "n_corrected_wins": a.n_corrected_wins,
            "win_rate": a.win_rate,
            "mean_margin": a.mean_margin,
            "median_margin": a.median_margin,
        }
        for a in report.arms
    ]
    cov_records = [
        {
            "model_type": p.model_type,
            "horizon_step": p.horizon_step,
            "n": p.n,
            "coverage": p.coverage,
            "nominal_coverage": report.nominal_coverage,
            "coverage_error": (
                p.coverage - report.nominal_coverage if p.coverage is not None else None
            ),
            "mean_width": p.mean_width,
        }
        for p in report.coverage
    ]
    arms_df = (
        pd.DataFrame.from_records(arm_records, columns=arm_cols)
        if arm_records
        else pd.DataFrame(columns=arm_cols)
    )
    cov_df = (
        pd.DataFrame.from_records(cov_records, columns=cov_cols)
        if cov_records
        else pd.DataFrame(columns=cov_cols)
    )
    return arms_df, cov_df

build_cohorts_frame

build_cohorts_frame(review)

Convert a RunReview's per-model BacktestCohort and comparable holdout metrics into a tidy DataFrame (pure, offline).

Exposes v_backtest_coverage and v_model_leaderboard_comparable side-by-side so you can inspect achieved fold counts (fold_histogram), refit modes (per_fold / recondition / extrapolate / unsupported), staleness_gap (cost of never refitting), and holdout-fold pooled_wape alongside the fleet mean score.

Source code in src/scale_forecasting/review.py

def build_cohorts_frame(review: RunReview) -> pd.DataFrame:
    """Convert a `RunReview`'s per-model `BacktestCohort` and comparable holdout metrics into a
    tidy ``DataFrame`` (pure, offline).

    Exposes ``v_backtest_coverage`` and ``v_model_leaderboard_comparable`` side-by-side so you can
    inspect achieved fold counts (`fold_histogram`), refit modes (`per_fold` / `recondition` /
    `extrapolate` / `unsupported`), `staleness_gap` (cost of never refitting), and holdout-fold
    `pooled_wape` alongside the fleet mean `score`.
    """
    import pandas as pd

    cols = [
        "model_type",
        "family",
        "is_ensemble",
        "n_series",
        "score",
        "pooled_wape",
        "n_comparable_series",
        "n_full",
        "n_reduced",
        "n_unscored",
        "n_failed",
        "n_not_requested",
        "fold_histogram",
        "refit_modes",
        "staleness_gap",
    ]
    if not review.models:
        return pd.DataFrame(columns=cols)

    records: list[dict[str, Any]] = []
    for m in review.models:
        c = m.cohort or BacktestCohort()
        hist_str = ", ".join(f"{k}f:{v}" for k, v in c.fold_histogram.items()) or "none"
        modes_str = ", ".join(f"{k}:{v}" for k, v in c.refit_modes.items()) or "none"
        records.append(
            {
                "model_type": m.model_type,
                "family": m.family,
                "is_ensemble": m.is_ensemble,
                "n_series": m.n_series if m.n_series is not None else c.n_series,
                "score": m.score,
                "pooled_wape": m.pooled_wape,
                "n_comparable_series": m.n_comparable_series,
                "n_full": c.n_full,
                "n_reduced": c.n_reduced,
                "n_unscored": c.n_unscored,
                "n_failed": c.n_failed,
                "n_not_requested": c.n_not_requested,
                "fold_histogram": hist_str,
                "refit_modes": modes_str,
                "staleness_gap": c.staleness_gap,
            }
        )
    return pd.DataFrame.from_records(records, columns=cols)

build_ensemble_weights_frame

build_ensemble_weights_frame(best_params_rows)

Unpack learned ensemble stacking weights from read_best_params rows or a build_best_params_frame DataFrame into a long-form DataFrame (pure, offline).

Columns: ts_id, ensemble_model, strategy, base_model, weight, wape.

Source code in src/scale_forecasting/review.py

def build_ensemble_weights_frame(
    best_params_rows: list[dict[str, Any]] | pd.DataFrame,
) -> pd.DataFrame:
    """Unpack learned ensemble stacking weights from `read_best_params` rows or a
    `build_best_params_frame` ``DataFrame`` into a long-form ``DataFrame`` (pure, offline).

    Columns: ``ts_id``, ``ensemble_model``, ``strategy``, ``base_model``, ``weight``, ``wape``.
    """
    import json

    import pandas as pd

    cols = ["ts_id", "ensemble_model", "strategy", "base_model", "weight", "wape"]
    raw_list: Sequence[Mapping[Any, Any]]
    if isinstance(best_params_rows, pd.DataFrame):
        if best_params_rows.empty:
            return pd.DataFrame(columns=cols)
        raw_list = best_params_rows.to_dict(orient="records")
    else:
        raw_list = list(best_params_rows)
    if not raw_list:
        return pd.DataFrame(columns=cols)

    records: list[dict[str, Any]] = []
    for r in raw_list:
        m_type = str(r.get("model_type") or "")
        is_ens = (
            bool(r.get("is_ensemble"))
            or (r.get("ensemble_id") is not None)
            or m_type.startswith("ensemble_")
        )
        if not is_ens:
            continue
        raw_bp = r.get("best_params")
        parsed: Any = None
        if isinstance(raw_bp, dict):
            parsed = raw_bp
        elif isinstance(raw_bp, str) and raw_bp:
            try:
                parsed = json.loads(raw_bp)
            except ValueError:
                parsed = None
        if not isinstance(parsed, dict) or not parsed:
            continue
        strategy = m_type.removeprefix("ensemble_")
        for base_model, weight in sorted(parsed.items()):
            w_val = _num(weight)
            if w_val is None:
                continue
            records.append(
                {
                    "ts_id": str(r.get("ts_id") or ""),
                    "ensemble_model": m_type,
                    "strategy": strategy,
                    "base_model": str(base_model),
                    "weight": w_val,
                    "wape": _num(r.get("wape")),
                }
            )
    if not records:
        return pd.DataFrame(columns=cols)
    return pd.DataFrame.from_records(records, columns=cols)

build_hierarchy_frame

build_hierarchy_frame(pred_rows)

Summarize hierarchical predictions by model and level, verifying additive coherence ($\max_t |\hat{y}{\text{total},t} - \sum|$) (pure).}} \hat{y}_{b,t

Source code in src/scale_forecasting/review.py

def build_hierarchy_frame(pred_rows: list[dict[str, Any]]) -> pd.DataFrame:
    """Summarize hierarchical predictions by model and level, verifying additive coherence
    ($\\max_t |\\hat{y}_{\\text{total},t} - \\sum_{b \\in \\text{bottom}} \\hat{y}_{b,t}|$) (pure).
    """
    import pandas as pd

    cols = [
        "model_type",
        "level",
        "n_series",
        "n_points",
        "mean_yhat",
        "sum_yhat",
        "max_coherence_residual",
    ]
    if not pred_rows:
        return pd.DataFrame(columns=cols)

    df = pd.DataFrame.from_records(pred_rows)
    df["ts_id"] = df["ts_id"].astype(str)
    df["level"] = df["ts_id"].map(_hierarchy_level_label)
    date_col = "forecast_date" if "forecast_date" in df.columns else "ds"

    records: list[dict[str, Any]] = []
    for model_type, m_grp in df.groupby("model_type", sort=True):
        tot = m_grp[m_grp["level"] == "total"].groupby(date_col)["yhat"].sum()
        bot = m_grp[m_grp["level"] == "bottom"].groupby(date_col)["yhat"].sum()
        residual: float | None = None
        if not tot.empty and not bot.empty:
            aligned = pd.concat([tot.rename("total"), bot.rename("bottom")], axis=1).dropna()
            if not aligned.empty:
                residual = float((aligned["total"] - aligned["bottom"]).abs().max())
        for lvl in ("total", "aggregate", "bottom"):
            l_grp = m_grp[m_grp["level"] == lvl]
            if l_grp.empty:
                continue
            records.append(
                {
                    "model_type": str(model_type),
                    "level": lvl,
                    "n_series": int(l_grp["ts_id"].nunique()),
                    "n_points": len(l_grp),
                    "mean_yhat": float(l_grp["yhat"].astype(float).mean()),
                    "sum_yhat": float(l_grp["yhat"].astype(float).sum()),
                    "max_coherence_residual": residual,
                }
            )
    return pd.DataFrame.from_records(records, columns=cols)

build_leaderboard_frame

build_leaderboard_frame(review, *, all_metrics=False)

Convert a RunReview into a ranked pandas DataFrame (pure, offline).

Includes core ranking columns (rank, model_type, family, is_ensemble, compute_engine, n_series, score, pooled_wape, n_comparable_series), key accuracy and interval metrics (mean_wape, mean_smape, mean_mase, mean_rmsse, mean_mae, mean_rmse, mean_coverage, mean_interval_score), runtime/reliability (mean_fit_seconds, median_fit_seconds, no_artifact_rate, n_predictions), and ensemble lift (lift_vs_best_base, lift_pct). When all_metrics=True, appends mean_<metric> and p50_<metric> for every metric in METRIC_COLUMNS.

Source code in src/scale_forecasting/review.py

def build_leaderboard_frame(review: RunReview, *, all_metrics: bool = False) -> pd.DataFrame:
    """Convert a `RunReview` into a ranked pandas ``DataFrame`` (pure, offline).

    Includes core ranking columns (`rank`, `model_type`, `family`, `is_ensemble`, `compute_engine`,
    `n_series`, `score`, `pooled_wape`, `n_comparable_series`), key accuracy and interval metrics
    (`mean_wape`, `mean_smape`, `mean_mase`, `mean_rmsse`, `mean_mae`, `mean_rmse`,
    `mean_coverage`, `mean_interval_score`), runtime/reliability (`mean_fit_seconds`,
    `median_fit_seconds`, `no_artifact_rate`, `n_predictions`), and ensemble lift
    (`lift_vs_best_base`, `lift_pct`). When ``all_metrics=True``, appends ``mean_<metric>`` and
    ``p50_<metric>`` for every metric in `METRIC_COLUMNS`.
    """
    import pandas as pd

    from .models import base_model_name

    lift_map = {(e.model_type): e for e in review.ensemble_lift}
    records: list[dict[str, Any]] = []
    for idx, m in enumerate(review.models, start=1):
        lift = lift_map.get(m.model_type)
        row: dict[str, Any] = {
            "rank": idx,
            "model_type": m.model_type,
            "base_model": m.base_model
            or (m.model_type if m.is_ensemble else base_model_name(m.model_type)),
            "family": m.family,
            "scope_capability": m.scope_capability,
            "training_mode": m.training_mode,
            "cov_profile": m.cov_profile,
            "is_univariate": m.is_univariate,
            "covariates_used": m.covariates_used,
            "model_spec_hash": m.model_spec_hash,
            "borrowed_from_run_id": m.borrowed_from_run_id,
            "n_borrowed": m.n_borrowed,
            "is_ensemble": m.is_ensemble,
            "ensemble_id": m.ensemble_id,
            "compute_engine": m.compute_engine,
            "n_series": m.n_series,
            "decision_metric": review.decision_metric,
            "score": m.score,
            "pooled_wape": m.pooled_wape,
            "n_comparable_series": m.n_comparable_series,
            "mean_wape": m.metric_means.get("wape"),
            "mean_smape": m.metric_means.get("smape"),
            "mean_mase": m.metric_means.get("mase"),
            "mean_rmsse": m.metric_means.get("rmsse"),
            "mean_mae": m.metric_means.get("mae"),
            "mean_rmse": m.metric_means.get("rmse"),
            "mean_coverage": m.metric_means.get("coverage"),
            "mean_interval_score": m.metric_means.get("interval_score"),
            "mean_fit_seconds": m.mean_fit_seconds,
            "median_fit_seconds": m.median_fit_seconds,
            "no_artifact_rate": m.no_artifact_rate,
            "n_predictions": m.n_predictions,
            "lift_vs_best_base": lift.lift if lift is not None else None,
            "lift_pct": lift.lift_pct if lift is not None else None,
        }
        if all_metrics:
            for col in METRIC_COLUMNS:
                row[f"mean_{col}"] = m.metric_means.get(col)
                row[f"p50_{col}"] = m.metric_p50.get(col)
        records.append(row)
    if not records:
        return pd.DataFrame(
            columns=[
                "rank",
                "model_type",
                "base_model",
                "family",
                "scope_capability",
                "training_mode",
                "cov_profile",
                "is_univariate",
                "covariates_used",
                "model_spec_hash",
                "borrowed_from_run_id",
                "n_borrowed",
                "is_ensemble",
                "ensemble_id",
                "compute_engine",
                "n_series",
                "decision_metric",
                "score",
                "pooled_wape",
                "n_comparable_series",
                "mean_wape",
                "mean_smape",
                "mean_mase",
                "mean_rmsse",
                "mean_mae",
                "mean_rmse",
                "mean_coverage",
                "mean_interval_score",
                "mean_fit_seconds",
                "median_fit_seconds",
                "no_artifact_rate",
                "n_predictions",
                "lift_vs_best_base",
                "lift_pct",
            ]
        )
    return pd.DataFrame.from_records(records)

build_predictions_frame

build_predictions_frame(
    pred_rows, *, oof_rows=None, history_rows=None
)

Stack historical observations, out-of-fold backtests, and forward predictions into one tidy DataFrame (pure, offline).

Columns: ts_id, segment ("history" / "oof" / "forecast"), model_type, ds, y_true, yhat, yhat_lower, yhat_upper, fold_id.

Source code in src/scale_forecasting/review.py

def build_predictions_frame(
    pred_rows: list[dict[str, Any]],
    *,
    oof_rows: list[dict[str, Any]] | None = None,
    history_rows: list[dict[str, Any]] | None = None,
) -> pd.DataFrame:
    """Stack historical observations, out-of-fold backtests, and forward predictions into one tidy
    ``DataFrame`` (pure, offline).

    Columns: ``ts_id``, ``segment`` (``"history"`` / ``"oof"`` / ``"forecast"``), ``model_type``,
    ``ds``, ``y_true``, ``yhat``, ``yhat_lower``, ``yhat_upper``, ``fold_id``.
    """
    import pandas as pd

    cols = [
        "ts_id",
        "segment",
        "model_type",
        "ds",
        "y_true",
        "yhat",
        "yhat_lower",
        "yhat_upper",
        "fold_id",
    ]
    records: list[dict[str, Any]] = []
    for r in history_rows or []:
        records.append(
            {
                "ts_id": str(r["ts_id"]),
                "segment": "history",
                "model_type": "actual",
                "ds": _first_date(r, "ds", "forecast_date"),
                "y_true": _num(r.get("y") if "y" in r else r.get("y_true")),
                "yhat": None,
                "yhat_lower": None,
                "yhat_upper": None,
                "fold_id": None,
            }
        )
    for r in oof_rows or []:
        records.append(
            {
                "ts_id": str(r["ts_id"]),
                "segment": "oof",
                "model_type": str(r["model_type"]),
                "ds": _first_date(r, "forecast_date", "ds"),
                "y_true": _num(r.get("y_true")),
                "yhat": _num(r.get("yhat")),
                "yhat_lower": _num(r.get("yhat_lower")),
                "yhat_upper": _num(r.get("yhat_upper")),
                "fold_id": r.get("fold_id"),
            }
        )
    for r in pred_rows:
        records.append(
            {
                "ts_id": str(r["ts_id"]),
                "segment": "forecast",
                "model_type": str(r["model_type"]),
                "ds": _first_date(r, "forecast_date", "ds"),
                "y_true": None,
                "yhat": _num(r.get("yhat")),
                "yhat_lower": _num(r.get("yhat_lower")),
                "yhat_upper": _num(r.get("yhat_upper")),
                "fold_id": None,
            }
        )
    if not records:
        return pd.DataFrame(columns=cols)
    df = pd.DataFrame.from_records(records, columns=cols)
    return df.sort_values(["ts_id", "segment", "model_type", "ds"]).reset_index(drop=True)

calibration_report

calibration_report(run_id, *, settings=None)

Read a finished run's point-forecast and interval diagnostic.

Reads the config (for the decision metric and the quantile set that defines nominal coverage), the per-model arm rollup (registry.reads.read_arm_comparison) and the per-step coverage panel (registry.reads.read_coverage_by_step), then composes via _assemble_calibration.

Source code in src/scale_forecasting/review.py

def calibration_report(run_id: str, *, settings: Settings | None = None) -> CalibrationReport:
    """Read a finished run's point-forecast and interval diagnostic.

    Reads the config (for the decision metric and the quantile set that defines nominal coverage),
    the per-model arm rollup (`registry.reads.read_arm_comparison`) and the per-step coverage panel
    (`registry.reads.read_coverage_by_step`), then composes via `_assemble_calibration`.
    """
    from .registry.reads import read_arm_comparison, read_coverage_by_step, read_run_config

    raw = read_run_config(run_id, settings=settings)
    cfg = RunConfig.model_validate(raw) if raw else None
    decision_metric = cfg.backtest.decision_metric if cfg else "wape"
    # Nominal is the span of the shipped quantile set, not a constant: a run that widened its
    # quantiles is not under-covering just because it left the default behind.
    from .models.base_model import DEFAULT_QUANTILES

    nominal = max(DEFAULT_QUANTILES) - min(DEFAULT_QUANTILES)
    return _assemble_calibration(
        run_id,
        decision_metric,
        nominal,
        read_arm_comparison(run_id, settings=settings),
        read_coverage_by_step(run_id, settings=settings),
    )

ensemble_lift

ensemble_lift(models)

Score each ensemble's improvement over the best base model, best lift first.

Empty when there is no scored base model to compare against, or no scored ensemble.

Source code in src/scale_forecasting/review.py

def ensemble_lift(
    models: list[ModelReview] | tuple[ModelReview, ...],
) -> list[EnsembleLift]:
    """Score each ensemble's improvement over the best base model, best lift first.

    Empty when there is no scored base model to compare against, or no scored ensemble.
    """
    champ = best_overall(models)
    if champ is None or champ.score is None:
        return []
    champ_score = champ.score
    lifts = [
        EnsembleLift(
            model_type=m.model_type,
            score=m.score,
            best_base_model=champ.model_type,
            best_base_score=champ_score,
            lift=champ_score - m.score,
            lift_pct=((champ_score - m.score) / champ_score if champ_score else None),
        )
        for m in models
        if m.is_ensemble and m.score is not None
    ]
    return sorted(lifts, key=lambda x: x.lift, reverse=True)

explain_forecast_frame

explain_forecast_frame(
    frame,
    *,
    ts_id=None,
    model_type=None,
    seasonal_period=7,
    covariate_df=None,
)

Decompose a single series' historical actuals, backtest OOF trajectory, and future forecast into interpretable structural components (pure, offline).

Given a frame from build_predictions_frame (or Forecaster.predictions_df), extracts: 1. trend_baseline & level_shift: Smooth rolling/piecewise trend plus abrupt regime-jump detection via features.level_shift_step. 2. seasonal_effect: Repeating periodic cycle (period seasonal_period, default 7) estimated from detrended history and projected across the forecast horizon. 3. covariate_effect (plus cov_<name> columns when covariate_df is provided): Ridge-stabilized linear attribution of numeric exogenous covariates across history and the future forecast horizon. 4. oof_residual & interval_width: Out-of-fold error ($y - \hat{y}{\text{OOF}}$) on backtest dates and prediction-interval width ($\hat{y} - \hat{y}_{\text{lower}}$) across the future horizon.}

Source code in src/scale_forecasting/review.py

def explain_forecast_frame(
    frame: pd.DataFrame,
    *,
    ts_id: str | None = None,
    model_type: str | None = None,
    seasonal_period: int = 7,
    covariate_df: pd.DataFrame | None = None,
) -> pd.DataFrame:
    """Decompose a single series' historical actuals, backtest OOF trajectory, and future forecast
    into interpretable structural components (pure, offline).

    Given a ``frame`` from `build_predictions_frame` (or `Forecaster.predictions_df`), extracts:
    1. ``trend_baseline`` & ``level_shift``: Smooth rolling/piecewise trend plus abrupt regime-jump
       detection via `features.level_shift_step`.
    2. ``seasonal_effect``: Repeating periodic cycle (period ``seasonal_period``, default ``7``)
       estimated from detrended history and projected across the forecast horizon.
    3. ``covariate_effect`` (plus ``cov_<name>`` columns when ``covariate_df`` is provided):
       Ridge-stabilized linear attribution of numeric exogenous covariates across history and the
       future forecast horizon.
    4. ``oof_residual`` & ``interval_width``: Out-of-fold error ($y - \\hat{y}_{\\text{OOF}}$) on
       backtest dates and prediction-interval width ($\\hat{y}_{\\text{upper}} -
       \\hat{y}_{\\text{lower}}$) across the future horizon.
    """
    import numpy as np
    import pandas as pd

    from .features import level_shift_step

    base_cols = [
        "ts_id",
        "model_type",
        "segment",
        "ds",
        "y_true",
        "yhat",
        "yhat_lower",
        "yhat_upper",
        "trend_baseline",
        "level_shift",
        "seasonal_effect",
        "covariate_effect",
        "oof_residual",
        "interval_width",
    ]
    if frame.empty:
        return pd.DataFrame(columns=base_cols)

    df = frame.copy()
    if "ds" not in df.columns and "forecast_date" in df.columns:
        df["ds"] = pd.to_datetime(df["forecast_date"])
    else:
        df["ds"] = pd.to_datetime(df["ds"])
    if "segment" not in df.columns:
        df["segment"] = "forecast"

    all_ts = list(dict.fromkeys(df["ts_id"].astype(str)))
    chosen_ts = ts_id if (ts_id is not None and ts_id in set(all_ts)) else all_ts[0]
    sub = df[df["ts_id"].astype(str) == chosen_ts]

    fc_models = list(dict.fromkeys(sub[sub["segment"] != "history"]["model_type"].astype(str)))
    chosen_model = (
        model_type
        if (model_type is not None and model_type in set(fc_models))
        else (fc_models[0] if fc_models else "actual")
    )

    hist = sub[sub["segment"] == "history"].sort_values("ds").drop_duplicates(subset=["ds"])
    fc = (
        sub[(sub["segment"] == "forecast") & (sub["model_type"].astype(str) == chosen_model)]
        .sort_values("ds")
        .drop_duplicates(subset=["ds"])
    )
    oof = (
        sub[(sub["segment"] == "oof") & (sub["model_type"].astype(str) == chosen_model)]
        .sort_values("ds")
        .drop_duplicates(subset=["ds"])
    )

    timeline_parts: list[pd.DataFrame] = []
    if not hist.empty:
        timeline_parts.append(
            pd.DataFrame(
                {
                    "ds": hist["ds"].to_numpy(),
                    "segment": "history",
                    "y_true": hist["y_true"].astype(float).to_numpy(),
                    "yhat": np.nan,
                    "yhat_lower": np.nan,
                    "yhat_upper": np.nan,
                }
            )
        )
    if not fc.empty:
        timeline_parts.append(
            pd.DataFrame(
                {
                    "ds": fc["ds"].to_numpy(),
                    "segment": "forecast",
                    "y_true": np.nan,
                    "yhat": fc["yhat"].astype(float).to_numpy(),
                    "yhat_lower": (
                        fc["yhat_lower"].astype(float).to_numpy()
                        if "yhat_lower" in fc.columns
                        else np.nan
                    ),
                    "yhat_upper": (
                        fc["yhat_upper"].astype(float).to_numpy()
                        if "yhat_upper" in fc.columns
                        else np.nan
                    ),
                }
            )
        )
    if not timeline_parts:
        return pd.DataFrame(columns=base_cols)

    tl = pd.concat(timeline_parts, ignore_index=True).sort_values("ds").reset_index(drop=True)
    tl["ts_id"] = chosen_ts
    tl["model_type"] = chosen_model

    # Overlay OOF yhat & residuals on matching history dates (or append if history was omitted)
    oof_map: dict[pd.Timestamp, tuple[float | None, float | None, float | None, float | None]] = {}
    for rec in oof.to_dict(orient="records"):
        oof_map[pd.Timestamp(rec["ds"])] = (
            _num(rec.get("y_true")),
            _num(rec.get("yhat")),
            _num(rec.get("yhat_lower")),
            _num(rec.get("yhat_upper")),
        )

    oof_residuals = np.full(len(tl), np.nan)
    for idx, row_ds in enumerate(tl["ds"]):
        ts_key = pd.Timestamp(row_ds)
        if ts_key in oof_map:
            y_t, yh_oof, lo_oof, hi_oof = oof_map[ts_key]
            if tl.loc[idx, "segment"] == "history" and yh_oof is not None:
                tl.loc[idx, "yhat"] = yh_oof
                if lo_oof is not None:
                    tl.loc[idx, "yhat_lower"] = lo_oof
                if hi_oof is not None:
                    tl.loc[idx, "yhat_upper"] = hi_oof
                actual_val = _num(tl.loc[idx, "y_true"])
                if actual_val is None:
                    actual_val = y_t
                if actual_val is not None:
                    oof_residuals[idx] = actual_val - yh_oof

    # Combined continuous signal (y_true on history, yhat on forecast)
    signal = np.where(
        tl["segment"] == "history",
        tl["y_true"].to_numpy(dtype=float),
        tl["yhat"].to_numpy(dtype=float),
    )
    n_hist = int((tl["segment"] == "history").sum())

    # 1. Level-shift detection on history + rolling trend baseline
    level_shift_arr = np.zeros(len(tl), dtype=float)
    if n_hist >= 16:
        hist_series = pd.Series(signal[:n_hist])
        step_dummy = level_shift_step(hist_series)
        if step_dummy.max() > 0:
            pre_mean = float(hist_series[step_dummy == 0].mean())
            post_mean = float(hist_series[step_dummy == 1].mean())
            jump = post_mean - pre_mean
            level_shift_arr[:n_hist] = step_dummy * jump
            level_shift_arr[n_hist:] = jump

    adj_signal = signal - level_shift_arr
    win = max(3, int(seasonal_period))
    trend_adj = (
        pd.Series(adj_signal)
        .rolling(window=win, center=True, min_periods=1)
        .mean()
        .to_numpy(dtype=float)
    )
    trend_baseline = trend_adj + level_shift_arr
    detrended = signal - trend_baseline

    # 2. Periodic seasonal profile (period = seasonal_period)
    seasonal_effect = np.zeros(len(tl), dtype=float)
    p = max(1, int(seasonal_period))
    if p > 1 and len(tl) >= p:
        ref_len = n_hist if n_hist >= p else len(tl)
        phase_means = np.zeros(p, dtype=float)
        for rem in range(p):
            vals = detrended[:ref_len][np.arange(ref_len) % p == rem]
            vals = vals[np.isfinite(vals)]
            phase_means[rem] = float(vals.mean()) if vals.size else 0.0
        phase_means -= float(phase_means.mean())
        for idx in range(len(tl)):
            seasonal_effect[idx] = phase_means[idx % p]

    # 3. Exogenous covariate attribution (when covariate_df is supplied)
    covariate_effect = np.zeros(len(tl), dtype=float)
    cov_added_cols: list[str] = []
    if covariate_df is not None and not covariate_df.empty:
        cov_work = covariate_df.copy()
        if "ts_id" in cov_work.columns:
            cov_work = cov_work[cov_work["ts_id"].astype(str) == chosen_ts]
        date_c = (
            "ds"
            if "ds" in cov_work.columns
            else ("forecast_date" if "forecast_date" in cov_work.columns else None)
        )
        if date_c is not None and not cov_work.empty:
            cov_work["ds"] = pd.to_datetime(cov_work[date_c])
            cov_work = cov_work.sort_values("ds").drop_duplicates(subset=["ds"])
            exclude = {
                "ts_id",
                "ds",
                "forecast_date",
                "y",
                "y_true",
                "yhat",
                "yhat_lower",
                "yhat_upper",
                "segment",
                "fold_id",
            }
            num_cols = [
                c
                for c in cov_work.columns
                if c not in exclude and pd.api.types.is_numeric_dtype(cov_work[c])
            ]
            if num_cols:
                merged_cov = pd.merge(tl[["ds"]], cov_work[["ds", *num_cols]], on="ds", how="left")
                X_raw = merged_cov[num_cols].ffill().bfill().fillna(0.0).to_numpy(dtype=float)
                fit_n = n_hist if n_hist >= len(num_cols) + 2 else len(tl)
                col_means = X_raw[:fit_n].mean(axis=0)
                X_centered = X_raw - col_means
                target_res = (detrended - seasonal_effect)[:fit_n]
                valid_mask = np.isfinite(target_res)
                if int(valid_mask.sum()) >= 2:
                    X_fit = X_centered[:fit_n][valid_mask]
                    y_fit = target_res[valid_mask]
                    # Small ridge penalty so collinear covariates decompose cleanly
                    gram = X_fit.T @ X_fit + 1e-3 * np.eye(X_fit.shape[1])
                    betas = np.linalg.solve(gram, X_fit.T @ y_fit)
                    for c_idx, c_name in enumerate(num_cols):
                        col_contrib = X_centered[:, c_idx] * float(betas[c_idx])
                        out_col = f"cov_{c_name}"
                        tl[out_col] = col_contrib
                        cov_added_cols.append(out_col)
                        covariate_effect += col_contrib

    tl["trend_baseline"] = trend_baseline
    tl["level_shift"] = level_shift_arr
    tl["seasonal_effect"] = seasonal_effect
    tl["covariate_effect"] = covariate_effect
    tl["oof_residual"] = oof_residuals
    tl["interval_width"] = np.where(
        tl["yhat_upper"].notna() & tl["yhat_lower"].notna(),
        tl["yhat_upper"].astype(float) - tl["yhat_lower"].astype(float),
        np.nan,
    )
    return tl[[*base_cols, *cov_added_cols]]

family_of

family_of(model_type, ensemble_id=None)

The family a result row belongs to: "ensemble" if ensemble_id is set, else the model's registered family ("unknown" for a name this build no longer registers).

Source code in src/scale_forecasting/review.py

def family_of(model_type: str, ensemble_id: str | None = None) -> str:
    """The family a result row belongs to: ``"ensemble"`` if ``ensemble_id`` is set, else the
    model's registered ``family`` (``"unknown"`` for a name this build no longer registers)."""
    if ensemble_id is not None:
        return "ensemble"
    from .errors import ModelError
    from .models import get_model

    try:
        return get_model(model_type).family
    except ModelError:
        return "unknown"

monitor_run

monitor_run(
    run_id,
    *,
    probe=False,
    stale_after_s=None,
    settings=None,
)

Read a run's live progress: header status + per-family job state + series done vs. expected.

Reads the run's header (registry.reads.read_run_summary), its config (registry.reads.read_run_config, for the expected-work denominator), its jobs (registry.jobs.read_run_jobs) and its landed-cell counts (registry.reads.read_progress), then composes them via _assemble_progress. Poll it while a run is in flight; returns a status-only snapshot when the run id has never run. Every family carries quiet_seconds either way — the "is this bar frozen or just coarse" signal, free because it comes off rows already read.

probe=True additionally escalates the run's non-terminal jobs to their runtime and attaches the reconciled probes.reconcile.ProbeReport as RunProgress.probe — the answer to is this job still alive, which the registry alone cannot give. It shares one pass of reads with the registry side (probes.reconcile.read_and_probe), so probing costs the native calls and not a second set of queries; an already-terminal run short-circuits and touches no runtime at all. stale_after_s overrides the probe's startup grace (see probes.reconcile.probe_run) and is ignored when probe is False.

Source code in src/scale_forecasting/review.py

def monitor_run(
    run_id: str,
    *,
    probe: bool = False,
    stale_after_s: float | None = None,
    settings: Settings | None = None,
) -> RunProgress:  # pragma: no cover - GCP I/O
    """Read a run's live progress: header status + per-family job state + series done vs. expected.

    Reads the run's header (`registry.reads.read_run_summary`), its config
    (`registry.reads.read_run_config`, for the expected-work denominator), its jobs
    (`registry.jobs.read_run_jobs`) and its landed-cell counts (`registry.reads.read_progress`),
    then composes them via `_assemble_progress`. Poll it while a run is in flight; returns a
    status-only snapshot when the run id has never run. Every family carries ``quiet_seconds``
    either way — the "is this bar frozen or just coarse" signal, free because it comes off rows
    already read.

    ``probe=True`` additionally escalates the run's non-terminal jobs to their runtime and attaches
    the reconciled `probes.reconcile.ProbeReport` as ``RunProgress.probe`` — the answer to *is this
    job still alive*, which the registry alone cannot give. It shares one pass of reads with the
    registry side (`probes.reconcile.read_and_probe`), so probing costs the native calls and not a
    second set of queries; an already-terminal run short-circuits and touches no runtime at all.
    ``stale_after_s`` overrides the probe's startup grace (see `probes.reconcile.probe_run`) and
    is ignored when ``probe`` is ``False``.
    """
    from .registry.jobs import read_run_jobs
    from .registry.reads import read_progress, read_run_config, read_run_summary

    if probe:
        from .probes.reconcile import read_and_probe
        from .settings import Settings as _Settings

        s = settings if settings is not None else _Settings.resolve()
        progress, report, _rows = read_and_probe(
            run_id, job=None, settings=s, stale_after_s=stale_after_s
        )
        return replace(progress, probe=report)

    summary = read_run_summary(run_id, settings=settings)
    raw = read_run_config(run_id, settings=settings)
    cfg = RunConfig.model_validate(raw) if raw else None
    if cfg is None:
        return _assemble_progress(run_id, summary, None, [], [])
    job_rows = read_run_jobs(run_id, settings=settings)
    progress_rows = read_progress(run_id, settings=settings)
    return _assemble_progress(run_id, summary, cfg, job_rows, progress_rows)

plot_attributions

plot_attributions(
    attributions_df,
    *,
    ts_id=None,
    model_type=None,
    top_k=12,
    ax=None,
    title=None,
)

Plot Tier 1 global driver importance or Tier 2 per-horizon-step local feature attributions (pure, offline).

Source code in src/scale_forecasting/review.py

def plot_attributions(
    attributions_df: pd.DataFrame,
    *,
    ts_id: str | None = None,
    model_type: str | None = None,
    top_k: int = 12,
    ax: Any = None,
    title: str | None = None,
) -> Any:
    """Plot Tier 1 global driver importance or Tier 2 per-horizon-step local feature attributions
    (pure, offline).
    """
    plt = _pyplot()
    import numpy as np

    if ax is None:
        _, ax = plt.subplots(figsize=(9.5, 4.5))
    if attributions_df.empty:
        ax.set_title(f"{title or 'Feature Attributions'} (no attribution rows)")
        return ax

    df = attributions_df.copy()
    if ts_id is not None and "ts_id" in df.columns:
        df = df[df["ts_id"].astype(str) == str(ts_id)]
    if model_type is not None and "model_type" in df.columns:
        df = df[df["model_type"].astype(str) == str(model_type)]
    if df.empty:
        ax.set_title(f"{title or 'Feature Attributions'} (no matching rows)")
        return ax

    is_local = "attribution" in df.columns and "forecast_date" in df.columns
    if is_local:
        # Aggregate mean |attribution| to pick top_k features, then plot signed mean attribution
        # or horizon trajectory.
        top_feats = (
            df["attribution"]
            .abs()
            .groupby(df["feature"])
            .mean()
            .sort_values(ascending=False)
            .head(max(1, top_k))
            .index.tolist()
        )
        sub = df[df["feature"].isin(top_feats)]
        agg = sub.groupby("feature")["attribution"].mean().sort_values(ascending=True)
        ys = np.arange(len(agg))
        vals = agg.to_numpy(dtype=float)
        colors = ["#0072B2" if v >= 0 else "#D55E00" for v in vals]
        ax.barh(ys, vals, color=colors, height=0.6)
        ax.axvline(0.0, color="#666666", linestyle=":", linewidth=1.0)
        ax.set_yticks(list(ys))
        ax.set_yticklabels([str(f) for f in agg.index])
        ax.set_xlabel("mean signed local attribution across horizon")
        ax.set_title(title or "Tier 2 Local Feature Attributions (Horizon Mean)")
        return ax

    agg = (
        df.groupby("feature")["importance"]
        .mean()
        .sort_values(ascending=False)
        .head(max(1, top_k))
        .iloc[::-1]
    )
    ys = np.arange(len(agg))
    vals = agg.to_numpy(dtype=float)
    ax.barh(ys, vals, color="#0072B2", height=0.6)
    ax.set_yticks(list(ys))
    ax.set_yticklabels([str(f) for f in agg.index])
    ax.set_xlabel("normalized driver importance")
    ax.set_title(title or "Tier 1 Global Feature Attributions")
    return ax

plot_calibration

plot_calibration(report, *, ax=None, title=None)

Plot empirical prediction-interval coverage and mean interval width by horizon step from a CalibrationReport (pure, offline).

When ax=None, creates a two-panel side-by-side figure: * Left panel: Empirical coverage by horizon step ($h = 1 \dots H$) vs. nominal target (dashed horizontal reference line at report.nominal_coverage). * Right panel: Mean prediction-interval width ($\hat{y}{\text{upper}} - \hat{y}$) by horizon step.}

Source code in src/scale_forecasting/review.py

def plot_calibration(
    report: CalibrationReport,
    *,
    ax: Any = None,
    title: str | None = None,
) -> Any:
    """Plot empirical prediction-interval coverage and mean interval width by horizon step from a
    `CalibrationReport` (pure, offline).

    When ``ax=None``, creates a two-panel side-by-side figure:
    * Left panel: Empirical coverage by horizon step ($h = 1 \\dots H$) vs. nominal target
      (dashed horizontal reference line at ``report.nominal_coverage``).
    * Right panel: Mean prediction-interval width ($\\hat{y}_{\\text{upper}} -
      \\hat{y}_{\\text{lower}}$) by horizon step.
    """
    plt = _pyplot()

    _, cov_df = build_calibration_frames(report)
    palette = ["#0072B2", "#E69F00", "#009E73", "#D55E00", "#CC79A7", "#56B4E9"]
    heading = (
        title or f"{report.run_id} — Interval Calibration (nominal {report.nominal_coverage:.0%})"
    )

    if ax is not None:
        fig = ax.figure
        ax_cov = ax
        ax_width = None
    elif not cov_df.empty and cov_df["mean_width"].notna().any():
        fig, (ax_cov, ax_width) = plt.subplots(1, 2, figsize=(12, 4.2))
    else:
        fig, ax_cov = plt.subplots(figsize=(8, 4.2))
        ax_width = None

    if cov_df.empty:
        ax_cov.set_title(f"{heading} (no OOF coverage rows)")
        return ax_cov

    for idx, (m_name, grp) in enumerate(cov_df.groupby("model_type", sort=True)):
        ordered = grp.sort_values("horizon_step")
        color = palette[idx % len(palette)]
        ax_cov.plot(
            ordered["horizon_step"],
            ordered["coverage"],
            marker="o",
            markersize=4,
            linewidth=1.8,
            color=color,
            label=str(m_name),
        )
        if ax_width is not None and ordered["mean_width"].notna().any():
            ax_width.plot(
                ordered["horizon_step"],
                ordered["mean_width"],
                marker="s",
                markersize=3.5,
                linewidth=1.8,
                color=color,
                label=str(m_name),
            )

    ax_cov.axhline(
        report.nominal_coverage,
        color="#222222",
        linestyle="--",
        linewidth=1.4,
        label=f"nominal ({report.nominal_coverage:.0%})",
    )
    ax_cov.set_ylim(-0.02, 1.05)
    ax_cov.set_xlabel("horizon step (h)")
    ax_cov.set_ylabel("empirical coverage")
    ax_cov.legend(loc="best", fontsize=8)

    if ax_width is not None:
        ax_cov.set_title(f"Empirical Coverage by Step (nominal {report.nominal_coverage:.0%})")
        ax_width.set_xlabel("horizon step (h)")
        ax_width.set_ylabel("mean interval width (yhat_upper − yhat_lower)")
        ax_width.set_title("Mean Interval Width by Step")
        ax_width.legend(loc="best", fontsize=8)
        fig.suptitle(heading, fontsize=11, y=1.02)
        fig.tight_layout()
        return (ax_cov, ax_width)

    ax_cov.set_title(heading)
    return ax_cov

plot_ensemble_weights

plot_ensemble_weights(weights_df, *, ax=None, title=None)

Plot base-model stacking weights across learned ensemble strategies (nnls, ridge, xgb) as a horizontal stacked bar chart (pure, offline).

Source code in src/scale_forecasting/review.py

def plot_ensemble_weights(
    weights_df: list[dict[str, Any]] | pd.DataFrame,
    *,
    ax: Any = None,
    title: str | None = None,
) -> Any:
    """Plot base-model stacking weights across learned ensemble strategies (`nnls`, `ridge`, `xgb`)
    as a horizontal stacked bar chart (pure, offline).
    """
    plt = _pyplot()
    import numpy as np
    import pandas as pd

    df = (
        weights_df
        if isinstance(weights_df, pd.DataFrame) and "base_model" in weights_df.columns
        else build_ensemble_weights_frame(weights_df)
    )
    if ax is None:
        _, ax = plt.subplots(figsize=(10, 3.8))
    if df.empty:
        ax.set_title(f"{title or 'Ensemble Stacking Weights'} (no learned weight rows)")
        return ax

    pivot = (
        df.groupby(["ensemble_model", "base_model"], as_index=False)["weight"]
        .mean()
        .pivot(index="ensemble_model", columns="base_model", values="weight")
        .fillna(0.0)
    )
    # Normalize each strategy row by sum(|w|) so unconstrained (ridge) negative weights
    # never shift stacked bar segments backward or exceed 1.0
    abs_pivot = pivot.abs()
    row_sums = abs_pivot.sum(axis=1).replace(0.0, 1.0)
    norm_pivot = abs_pivot.div(row_sums, axis=0)

    palette = ["#0072B2", "#E69F00", "#009E73", "#D55E00", "#CC79A7", "#56B4E9", "#F0E442"]
    strategies = list(norm_pivot.index)
    base_models = list(norm_pivot.columns)
    ys = np.arange(len(strategies))
    left = np.zeros(len(strategies))

    for idx, bm in enumerate(base_models):
        vals = norm_pivot[bm].to_numpy(dtype=float)
        color = palette[idx % len(palette)]
        ax.barh(ys, vals, left=left, height=0.55, color=color, label=bm)
        for y_idx, v in enumerate(vals):
            if v >= 0.08:
                ax.text(
                    left[y_idx] + v / 2.0,
                    y_idx,
                    f"{v:.0%}",
                    ha="center",
                    va="center",
                    color="white" if idx in (0, 2, 3) else "#111111",
                    fontsize=8.5,
                    fontweight="bold",
                )
        left += vals

    ax.set_yticks(list(ys))
    ax.set_yticklabels(strategies)
    ax.set_xlim(0, 1.05)
    ax.set_xlabel("mean relative weight share (|w| / sum |w|) across series")
    ax.set_title(title or "Learned Ensemble Stacking Weights by Strategy")
    ax.legend(
        loc="upper center", bbox_to_anchor=(0.5, -0.18), ncol=min(4, len(base_models)), fontsize=8.5
    )
    return ax

plot_forecast_explanation

plot_forecast_explanation(
    explanation_df, *, history_tail=None, title=None
)

Render a 4-panel forecast explainability & decomposition chart from an explain_forecast_frame DataFrame (pure, offline).

Panels: 1. Trajectory, Trend Baseline & Regime Shift: History, OOF backtest, future forecast + 80% prediction interval, and structural trend baseline. 2. Periodic Seasonal Cycle: Extracted seasonal wave across history and forecast horizon. 3. Exogenous Covariate Attribution (or Detrended Residual Wave when univariate): Individual cov_<name> contributions and total covariate effect across history and the future horizon. 4. Backtest OOF Residuals & Forecast Horizon Uncertainty: Out-of-fold errors ($y - \hat{y}{\text{OOF}}$) transitioning into future 80% interval width ($\hat{y}$).}} - \hat{y}_{\text{lower}

Source code in src/scale_forecasting/review.py

def plot_forecast_explanation(
    explanation_df: pd.DataFrame,
    *,
    history_tail: int | None = None,
    title: str | None = None,
) -> Any:
    """Render a 4-panel forecast explainability & decomposition chart from an
    `explain_forecast_frame` ``DataFrame`` (pure, offline).

    Panels:
    1. **Trajectory, Trend Baseline & Regime Shift**: History, OOF backtest, future forecast + 80%
       prediction interval, and structural trend baseline.
    2. **Periodic Seasonal Cycle**: Extracted seasonal wave across history and forecast horizon.
    3. **Exogenous Covariate Attribution** (or Detrended Residual Wave when univariate): Individual
       ``cov_<name>`` contributions and total covariate effect across history and the future
       horizon.
    4. **Backtest OOF Residuals & Forecast Horizon Uncertainty**: Out-of-fold errors ($y -
       \\hat{y}_{\\text{OOF}}$) transitioning into future 80% interval width
       ($\\hat{y}_{\\text{upper}} - \\hat{y}_{\\text{lower}}$).
    """
    plt = _pyplot()
    import pandas as pd

    if explanation_df.empty:
        _, ax = plt.subplots(figsize=(10, 3))
        ax.set_title(f"{title or 'Forecast Explanation'} (no rows)")
        return ax

    df = explanation_df.copy()
    df["ds"] = pd.to_datetime(df["ds"])
    if history_tail is not None and history_tail > 0:
        hist_part = df[df["segment"] == "history"].tail(history_tail)
        fc_part = df[df["segment"] != "history"]
        df = pd.concat([hist_part, fc_part], ignore_index=True)
    tid = str(df["ts_id"].iloc[0])
    m_name = str(df["model_type"].iloc[0])
    heading = title or f"Forecast Decomposition & Attribution — series={tid} · model={m_name}"

    _, axes = plt.subplots(4, 1, figsize=(11, 9.2), sharex=True)
    ax_main, ax_seas, ax_cov, ax_err = axes

    hist = df[df["segment"] == "history"]
    fc = df[df["segment"] == "forecast"]

    # Panel 1: Main trajectory + trend + level shift
    if not hist.empty:
        ax_main.plot(
            hist["ds"],
            hist["y_true"],
            color="#222222",
            linewidth=1.5,
            label="actual (history)",
        )
        oof_valid = hist[hist["yhat"].notna()]
        if not oof_valid.empty:
            ax_main.plot(
                oof_valid["ds"],
                oof_valid["yhat"],
                color="#E69F00",
                linestyle="--",
                linewidth=1.4,
                label=f"{m_name} (OOF)",
            )
    if not fc.empty:
        ax_main.plot(
            fc["ds"],
            fc["yhat"],
            color="#0072B2",
            linewidth=2.1,
            label=f"{m_name} (forecast)",
        )
        if fc["yhat_lower"].notna().any() and fc["yhat_upper"].notna().any():
            ax_main.fill_between(
                fc["ds"],
                fc["yhat_lower"].astype(float),
                fc["yhat_upper"].astype(float),
                color="#0072B2",
                alpha=0.18,
                label="80% interval",
            )
        ax_main.axvline(fc["ds"].iloc[0], color="#666666", linestyle=":", linewidth=1.2)

    ax_main.plot(
        df["ds"],
        df["trend_baseline"],
        color="#009E73",
        linewidth=1.8,
        alpha=0.9,
        label="trend + regime baseline",
    )
    if df["level_shift"].abs().max() > 0:
        shift_idx = df.index[df["level_shift"].abs() > 0][0]
        ax_main.axvline(
            df.loc[shift_idx, "ds"],
            color="#D55E00",
            linestyle="-.",
            linewidth=1.3,
            label="detected level shift",
        )
    ax_main.set_title(heading)
    ax_main.set_ylabel("value")
    ax_main.legend(loc="best", fontsize=8, ncol=3)

    # Panel 2: Seasonal effect
    ax_seas.plot(
        df["ds"], df["seasonal_effect"], color="#0072B2", linewidth=1.4, label="seasonal cycle"
    )
    ax_seas.axhline(0.0, color="#999999", linestyle=":", linewidth=0.9)
    if not fc.empty:
        ax_seas.axvline(fc["ds"].iloc[0], color="#666666", linestyle=":", linewidth=1.2)
    ax_seas.set_ylabel("seasonal")
    ax_seas.legend(loc="upper right", fontsize=8)

    # Panel 3: Exogenous covariate attribution
    cov_cols = [c for c in df.columns if c.startswith("cov_")]
    palette = ["#D55E00", "#009E73", "#CC79A7", "#E69F00", "#56B4E9"]
    if cov_cols:
        for idx, c_col in enumerate(cov_cols):
            ax_cov.plot(
                df["ds"],
                df[c_col],
                linewidth=1.3,
                color=palette[idx % len(palette)],
                label=c_col.removeprefix("cov_"),
            )
        ax_cov.plot(
            df["ds"],
            df["covariate_effect"],
            color="#222222",
            linestyle="--",
            linewidth=1.4,
            label="total covariate effect",
        )
    else:
        ax_cov.plot(
            df["ds"],
            df["covariate_effect"],
            color="#666666",
            linewidth=1.2,
            label="covariate effect (none configured)",
        )
    ax_cov.axhline(0.0, color="#999999", linestyle=":", linewidth=0.9)
    if not fc.empty:
        ax_cov.axvline(fc["ds"].iloc[0], color="#666666", linestyle=":", linewidth=1.2)
    ax_cov.set_ylabel("covariates")
    ax_cov.legend(loc="upper right", fontsize=8, ncol=min(4, max(1, len(cov_cols) + 1)))

    # Panel 4: OOF residuals & future interval width
    oof_res = df[df["oof_residual"].notna()]
    if not oof_res.empty:
        ax_err.bar(
            oof_res["ds"],
            oof_res["oof_residual"],
            width=0.8,
            color="#E69F00",
            alpha=0.75,
            label="OOF residual (y − yhat)",
        )
    if not fc.empty and fc["interval_width"].notna().any():
        ax_err.plot(
            fc["ds"],
            fc["interval_width"],
            color="#0072B2",
            marker="o",
            markersize=3.5,
            linewidth=1.6,
            label="forecast 80% interval width",
        )
        ax_err.axvline(fc["ds"].iloc[0], color="#666666", linestyle=":", linewidth=1.2)
    ax_err.axhline(0.0, color="#999999", linestyle=":", linewidth=0.9)
    ax_err.set_ylabel("error / width")
    ax_err.set_xlabel("date")
    ax_err.legend(loc="best", fontsize=8, ncol=2)
    return axes

plot_forecasts_frame

plot_forecasts_frame(
    frame,
    *,
    ts_id=None,
    ts_ids=None,
    models=None,
    max_series=3,
    history_tail=None,
    ax=None,
    title=None,
)

Plot historical actuals, out-of-fold backtests, and forward predictions with 80% intervals from a build_predictions_frame DataFrame (pure, offline).

Renders up to max_series subplots (one per ts_id) and returns the matplotlib Axes (or array of Axes when multiple series are plotted).

Source code in src/scale_forecasting/review.py

def plot_forecasts_frame(
    frame: pd.DataFrame,
    *,
    ts_id: str | None = None,
    ts_ids: Sequence[str] | None = None,
    models: Sequence[str] | None = None,
    max_series: int = 3,
    history_tail: int | None = None,
    ax: Any = None,
    title: str | None = None,
) -> Any:
    """Plot historical actuals, out-of-fold backtests, and forward predictions with 80% intervals
    from a `build_predictions_frame` ``DataFrame`` (pure, offline).

    Renders up to ``max_series`` subplots (one per ``ts_id``) and returns the matplotlib ``Axes``
    (or array of ``Axes`` when multiple series are plotted).
    """
    plt = _pyplot()
    import pandas as pd

    if frame.empty:
        if ax is None:
            _, ax = plt.subplots(figsize=(10, 3))
        ax.set_title(f"{title or 'Forecasts'} (no rows)")
        return ax

    df = frame.copy()
    if "ds" not in df.columns and "forecast_date" in df.columns:
        df["ds"] = pd.to_datetime(df["forecast_date"])
    if "segment" not in df.columns:
        df["segment"] = "forecast"

    effective_ts_ids = [ts_id] if ts_id is not None else (list(ts_ids) if ts_ids else None)
    all_ts = list(dict.fromkeys(df["ts_id"].astype(str)))
    chosen_ts = [t for t in all_ts if not effective_ts_ids or t in set(effective_ts_ids)][
        : max(1, max_series)
    ]
    if not chosen_ts:
        if ax is None:
            _, ax = plt.subplots(figsize=(10, 3))
        ax.set_title(f"{title or 'Forecasts'} (no matching series)")
        return ax

    palette = ["#0072B2", "#E69F00", "#009E73", "#D55E00", "#CC79A7", "#56B4E9"]
    if ax is None:
        _, axes_raw = plt.subplots(
            len(chosen_ts),
            1,
            figsize=(11, 3.4 * len(chosen_ts)),
            squeeze=False,
        )
        axes = [axes_raw[i, 0] for i in range(len(chosen_ts))]
    else:
        axes = [ax]
        chosen_ts = chosen_ts[:1]

    for sub_ax, tid in zip(axes, chosen_ts, strict=True):
        sub = df[df["ts_id"].astype(str) == tid]
        hist = sub[sub["segment"] == "history"].sort_values("ds")
        if history_tail is not None and history_tail > 0:
            hist = hist.tail(history_tail)
        if not hist.empty:
            sub_ax.plot(
                hist["ds"],
                hist["y_true"],
                color="#222222",
                linewidth=1.6,
                label="actual (history)",
            )
        oof = sub[sub["segment"] == "oof"]
        if not oof.empty and hist.empty and "y_true" in oof.columns:
            actual_oof = (
                oof.dropna(subset=["y_true"]).drop_duplicates(subset=["ds"]).sort_values("ds")
            )
            if not actual_oof.empty:
                sub_ax.plot(
                    actual_oof["ds"],
                    actual_oof["y_true"],
                    color="#222222",
                    linewidth=1.5,
                    label="actual (OOF)",
                )
        fc_models = [
            m
            for m in dict.fromkeys(sub[sub["segment"] != "history"]["model_type"].astype(str))
            if not models or m in set(models)
        ]
        for m_idx, m_name in enumerate(fc_models):
            color = palette[m_idx % len(palette)]
            m_oof = oof[oof["model_type"].astype(str) == m_name].sort_values("ds")
            if not m_oof.empty:
                sub_ax.plot(
                    m_oof["ds"],
                    m_oof["yhat"],
                    linestyle="--",
                    linewidth=1.3,
                    color=color,
                    alpha=0.85,
                    label=f"{m_name} (OOF)",
                )
            m_fc = sub[
                (sub["segment"] == "forecast") & (sub["model_type"].astype(str) == m_name)
            ].sort_values("ds")
            if not m_fc.empty:
                sub_ax.plot(
                    m_fc["ds"],
                    m_fc["yhat"],
                    linestyle="-",
                    linewidth=2.0,
                    color=color,
                    label=f"{m_name} (forecast)",
                )
                if (
                    "yhat_lower" in m_fc.columns
                    and "yhat_upper" in m_fc.columns
                    and m_fc["yhat_lower"].notna().any()
                    and m_fc["yhat_upper"].notna().any()
                ):
                    sub_ax.fill_between(
                        m_fc["ds"],
                        m_fc["yhat_lower"].astype(float),
                        m_fc["yhat_upper"].astype(float),
                        color=color,
                        alpha=0.16,
                    )
        sub_ax.set_title(f"{title + ' — ' if title else ''}series: {tid}")
        sub_ax.set_ylabel("value")
        sub_ax.legend(loc="best", fontsize=8, ncol=2)
    axes[-1].set_xlabel("date")
    return axes[0] if len(axes) == 1 else axes

plot_hierarchy_frame

plot_hierarchy_frame(
    pred_rows, *, model_type=None, ax=None, title=None
)

Plot top-level (__total__) forecast against the sum of bottom-level leaf forecasts (or a level summary from build_hierarchy_frame) to visually verify hierarchical coherence.

Source code in src/scale_forecasting/review.py

def plot_hierarchy_frame(
    pred_rows: list[dict[str, Any]] | pd.DataFrame,
    *,
    model_type: str | None = None,
    ax: Any = None,
    title: str | None = None,
) -> Any:
    """Plot top-level (``__total__``) forecast against the sum of bottom-level leaf forecasts (or
    a level summary from `build_hierarchy_frame`) to visually verify hierarchical coherence.
    """
    plt = _pyplot()
    import pandas as pd

    df = (
        pred_rows.copy()
        if isinstance(pred_rows, pd.DataFrame)
        else pd.DataFrame.from_records(pred_rows)
    )
    if ax is None:
        _, ax = plt.subplots(figsize=(10, 4.5))
    if df.empty:
        ax.set_title(f"{title or 'Hierarchical Coherence'} (no rows)")
        return ax

    # Support passing the summary DataFrame from build_hierarchy_frame directly
    if {"level", "sum_yhat", "model_type"} <= set(df.columns) and "ts_id" not in df.columns:
        chosen_model = model_type or _default_hierarchy_model(df["model_type"])
        sub = df[df["model_type"].astype(str) == chosen_model]
        ax.bar(sub["level"].astype(str), sub["sum_yhat"].astype(float), color="#0072B2", alpha=0.85)
        ax.set_title(f"{title or 'Hierarchical Rollup Totals'}: {chosen_model}")
        ax.set_xlabel("hierarchy level")
        ax.set_ylabel("sum of yhat")
        return ax

    df["ts_id"] = df["ts_id"].astype(str)
    df["level"] = df["ts_id"].map(_hierarchy_level_label)
    date_col = "forecast_date" if "forecast_date" in df.columns else "ds"
    df[date_col] = pd.to_datetime(df[date_col])

    chosen_model = model_type or _default_hierarchy_model(df["model_type"])
    sub = df[df["model_type"].astype(str) == chosen_model]
    if sub.empty:
        ax.set_title(f"{title or 'Hierarchical Coherence'} (no rows for {chosen_model})")
        return ax

    tot = sub[sub["level"] == "total"].groupby(date_col)["yhat"].sum().sort_index()
    bot = sub[sub["level"] == "bottom"].groupby(date_col)["yhat"].sum().sort_index()
    if not tot.empty:
        ax.plot(
            tot.index,
            tot.to_numpy(),
            color="#0072B2",
            linewidth=2.4,
            label=f"Top-level total ({chosen_model})",
        )
    if not bot.empty:
        ax.plot(
            bot.index,
            bot.to_numpy(),
            color="#E69F00",
            linestyle="--",
            linewidth=2.0,
            label=f"Σ bottom-level leaves ({chosen_model})",
        )
    residual_str = ""
    if not tot.empty and not bot.empty:
        aligned = pd.concat([tot.rename("t"), bot.rename("b")], axis=1).dropna()
        if not aligned.empty:
            max_err = float((aligned["t"] - aligned["b"]).abs().max())
            residual_str = f" — max |total − Σ bottom| = {max_err:.2e}"
    ax.set_title(f"{title or 'Hierarchical Coherence'}: {chosen_model}{residual_str}")
    ax.set_xlabel("forecast date")
    ax.set_ylabel("forecast value")
    ax.tick_params(axis="x", labelrotation=15)
    ax.legend(loc="best", fontsize=9)
    return ax

plot_leaderboard

plot_leaderboard(review, *, ax=None, top=None, title=None)

Render a RunReview's scored models as a ranked bar chart; return the matplotlib Axes.

One horizontal bar per model, best (lowest decision-metric error) on top, coloured base vs. ensemble (_BASE_COLOR/_ENSEMBLE_COLOR) with a legend when both are present and the score labelled at each bar end. Unscored models (no backtest) are dropped. top caps the bar count. matplotlib imports lazily; a review with no scored model renders an empty titled axes.

Source code in src/scale_forecasting/review.py

def plot_leaderboard(
    review: RunReview, *, ax: Any = None, top: int | None = None, title: str | None = None
) -> Any:
    """Render a `RunReview`'s scored models as a ranked bar chart; return the matplotlib ``Axes``.

    One horizontal bar per model, best (lowest decision-metric error) on top, coloured base vs.
    ensemble (`_BASE_COLOR`/`_ENSEMBLE_COLOR`) with a legend when both are present and the score
    labelled at each bar end. Unscored models (no backtest) are dropped. ``top`` caps the bar count.
    matplotlib imports lazily; a review with no scored model renders an empty titled axes.
    """
    plt = _pyplot()

    scored = [m for m in review.models if m.score is not None]
    if top is not None:
        scored = scored[:top]
    heading = title or f"{review.run_id} — model leaderboard ({review.decision_metric})"
    if ax is None:
        _, ax = plt.subplots(figsize=(10, max(2.0, 0.5 * len(scored) + 1)))
    ax.set_title(heading if scored else f"{heading} (no scored models)")
    if not scored:
        return ax

    ranked = list(reversed(scored))  # best on top
    ys = range(len(ranked))
    ax.barh(
        list(ys),
        [m.score for m in ranked],
        height=0.6,
        color=[_ENSEMBLE_COLOR if m.is_ensemble else _BASE_COLOR for m in ranked],
    )
    for y, m in zip(ys, ranked, strict=True):
        ax.text(m.score, y, f" {m.score:.4g}", va="center", fontsize=9)
    max_score = max((float(m.score) for m in ranked if m.score is not None), default=0.0)
    if max_score > 0:
        ax.set_xlim(0, max_score * 1.32)
    ax.set_yticks(list(ys))
    ax.set_yticklabels([m.model_type for m in ranked])
    ax.set_xlabel(f"mean {review.decision_metric} (lower is better)")
    if any(m.is_ensemble for m in ranked) and any(not m.is_ensemble for m in ranked):
        from matplotlib.patches import Patch

        ax.legend(
            handles=[
                Patch(color=_BASE_COLOR, label="base model"),
                Patch(color=_ENSEMBLE_COLOR, label="ensemble"),
            ],
            loc="lower right",
            fontsize=9,
        )
    return ax

plot_metric_distribution

plot_metric_distribution(
    review, *, metric=None, ax=None, title=None
)

Render each model's cross-series spread for one metric (p10–p90 range, p50 dot) as Axes.

One row per scored model: a thin line from the 10th to the 90th cross-series percentile with a marker at the median, coloured base vs. ensemble — the distribution shape (not just the mean) of a metric across every series, read straight off the server-side aggregates so it holds at scale. metric defaults to the run's decision metric. matplotlib imports lazily; a review with no aggregated percentiles renders an empty titled axes.

Source code in src/scale_forecasting/review.py

def plot_metric_distribution(
    review: RunReview, *, metric: str | None = None, ax: Any = None, title: str | None = None
) -> Any:
    """Render each model's cross-series spread for one metric (p10–p90 range, p50 dot) as ``Axes``.

    One row per scored model: a thin line from the 10th to the 90th cross-series percentile with a
    marker at the median, coloured base vs. ensemble — the distribution shape (not just the mean) of
    a metric across every series, read straight off the server-side aggregates so it holds at scale.
    ``metric`` defaults to the run's decision metric. matplotlib imports lazily; a review with no
    aggregated percentiles renders an empty titled axes.
    """
    plt = _pyplot()

    chosen = metric or review.decision_metric
    rows = [m for m in review.models if m.metric_p50.get(chosen) is not None]
    heading = title or f"{review.run_id} — {chosen} across series (p10–p50–p90)"
    if ax is None:
        _, ax = plt.subplots(figsize=(10, max(2.0, 0.5 * len(rows) + 1)))
    ax.set_title(heading if rows else f"{heading} (no aggregated percentiles)")
    if not rows:
        return ax

    ordered = sorted(rows, key=lambda m: _sort_key(m.metric_p50.get(chosen)), reverse=True)
    for y, m in enumerate(ordered):
        color = _ENSEMBLE_COLOR if m.is_ensemble else _BASE_COLOR
        lo = m.metric_p10.get(chosen)
        hi = m.metric_p90.get(chosen)
        mid = m.metric_p50[chosen]
        if lo is not None and hi is not None:
            ax.hlines(y, lo, hi, color=color, linewidth=2)
        ax.plot(mid, y, "o", color=color, markersize=8)
    ax.set_yticks(range(len(ordered)))
    ax.set_yticklabels([m.model_type for m in ordered])
    ax.set_xlabel(f"{chosen} (p10–p90 range, dot = median)")
    if any(m.is_ensemble for m in ordered) and any(not m.is_ensemble for m in ordered):
        from matplotlib.lines import Line2D

        ax.legend(
            handles=[
                Line2D([0], [0], color=_BASE_COLOR, marker="o", label="base model"),
                Line2D([0], [0], color=_ENSEMBLE_COLOR, marker="o", label="ensemble"),
            ],
            loc="lower right",
            fontsize=9,
        )
    return ax

plot_progress

plot_progress(progress, *, ax=None, title=None)

Render a RunProgress as a per-family progress bar chart and return the matplotlib Axes.

One horizontal bar per family — length is the fraction of expected cells that have landed, colour is the family's job status (_STATUS_COLORS) — with the done/expected count and status labelled at the bar end (the label doubles as the status's secondary encoding). Families keep DAG order, ensemble last. matplotlib imports lazily so it never touches the package import; an empty run renders an empty titled axes rather than raising.

A non-terminal family also gets its quiet_seconds in the label (quiet 22m): the bar of a job that died mid-run stops moving and its status stays RUNNING, so without this a dead run and a slow one are pixel-identical. It is reported as an age, not flagged against a threshold — a family that writes its cells at job end is legitimately quiet the whole time, and a cry-wolf marker teaches the reader to ignore it. When a probes.reconcile.ProbeReport is attached (monitor_run(probe=True)), its verdict replaces the age for the families it covers, since a live reading beats an inference from silence.

A family that finished with a device verdict gets two more words at the end of its label (gpu used / gpu idle / no gpu). A GPU family's bar is otherwise indistinguishable from a CPU family's, which is the whole reason the accelerator went twenty-one jobs without anyone noticing it was doing nothing.

Source code in src/scale_forecasting/review.py

def plot_progress(progress: RunProgress, *, ax: Any = None, title: str | None = None) -> Any:
    """Render a `RunProgress` as a per-family progress bar chart and return the matplotlib ``Axes``.

    One horizontal bar per family — length is the fraction of expected cells that have landed,
    colour is the family's job status (`_STATUS_COLORS`) — with the ``done/expected`` count and
    status labelled at the bar end (the label doubles as the status's secondary encoding). Families
    keep DAG order, ensemble last. matplotlib imports lazily so it never touches the package import;
    an empty run renders an empty titled axes rather than raising.

    A non-terminal family also gets its ``quiet_seconds`` in the label (``quiet 22m``): the bar of a
    job that died mid-run stops moving and its status stays ``RUNNING``, so without this a dead run
    and a slow one are pixel-identical. It is reported as an age, not flagged against a threshold —
    a family that writes its cells at job end is legitimately quiet the whole time, and a
    cry-wolf marker teaches the reader to ignore it. When a `probes.reconcile.ProbeReport` is
    attached (``monitor_run(probe=True)``), its verdict replaces the age for the families it
    covers, since a live reading beats an inference from silence.

    A family that finished with a device verdict gets two more words at the end of its label
    (``gpu used`` / ``gpu idle`` / ``no gpu``). A GPU family's bar is otherwise indistinguishable
    from a CPU family's, which is the whole reason the accelerator went twenty-one jobs without
    anyone noticing it was doing nothing.
    """
    plt = _pyplot()

    heading = title or f"{progress.run_id} — {progress.status or 'unknown'}"
    if progress.fraction is not None:
        heading += f" — {progress.fraction:.0%} of cells landed"
    if ax is None:
        _, ax = plt.subplots(figsize=(10, max(2.0, 0.6 * len(progress.families) + 1)))
    ax.set_title(heading if progress.families else f"{heading} (no families)")
    if not progress.families:
        return ax

    fams = list(reversed(progress.families))  # first family on top
    ys = range(len(fams))
    ax.barh(
        list(ys),
        [f.fraction if f.fraction is not None else 0.0 for f in fams],
        height=0.6,
        color=[_STATUS_COLORS.get(f.status, _STATUS_DEFAULT) for f in fams],
    )
    verdicts = _reportable_verdicts(progress)
    for y, f in zip(ys, fams, strict=True):
        expected = f.n_expected if f.n_expected is not None else "?"
        label = f"{f.n_done}/{expected} · {f.status or 'pending'}"
        if (verdict := verdicts.get(f.family)) is not None:
            label += f" · {verdict.lower().replace('_', ' ')}"
        elif (f.status or "").upper() == "RUNNING" and f.quiet_seconds is not None:
            label += f" · quiet {_human_age(f.quiet_seconds)}"
        # Appended rather than chained into the ladder above: liveness and cost are different
        # questions, and a finished GPU family has an answer to both.
        if (device := verdict_label(f.device_verdict)) is not None:
            label += f" · {device}"
        ax.text(
            (f.fraction if f.fraction is not None else 0.0) + 0.01,
            y,
            label,
            va="center",
            fontsize=9,
        )
    ax.set_yticks(list(ys))
    ax.set_yticklabels([f.family for f in fams])
    ax.set_xlim(0, 1.15)
    ax.set_xlabel("fraction of expected cells landed")
    return ax

review_run

review_run(run_id, *, settings=None)

Read a finished run's data-science review: bests per family/overall + ensemble lift + panel.

Reads the header (registry.reads.read_run_summary), the config (for the decision metric and series count), the leaderboard (registry.reads.read_leaderboard), the cross-series aggregates (registry.reads.read_metric_aggregates), per-model prediction counts (registry.reads.read_prediction_counts), and the two cohort reads that say whether the ranking is comparable at all — registry.reads.read_backtest_coverage and registry.reads.read_comparable_leaderboard — then composes via _assemble_review.

Source code in src/scale_forecasting/review.py

def review_run(
    run_id: str, *, settings: Settings | None = None
) -> RunReview:  # pragma: no cover - GCP I/O
    """Read a finished run's data-science review: bests per family/overall + ensemble lift + panel.

    Reads the header (`registry.reads.read_run_summary`), the config (for the decision metric and
    series count), the leaderboard (`registry.reads.read_leaderboard`), the cross-series aggregates
    (`registry.reads.read_metric_aggregates`), per-model prediction counts
    (`registry.reads.read_prediction_counts`), and the two cohort reads that say whether the
    ranking is comparable at all — `registry.reads.read_backtest_coverage` and
    `registry.reads.read_comparable_leaderboard` — then composes via `_assemble_review`.
    """
    from .registry.reads import (
        read_backtest_coverage,
        read_comparable_leaderboard,
        read_leaderboard,
        read_metric_aggregates,
        read_prediction_counts,
        read_run_config,
        read_run_summary,
    )

    summary = read_run_summary(run_id, settings=settings)
    raw = read_run_config(run_id, settings=settings)
    cfg = RunConfig.model_validate(raw) if raw else None
    decision_metric = cfg.backtest.decision_metric if cfg else "wape"
    n_series = (summary or {}).get("n_series") or (cfg.data.series_limit if cfg else None)
    leaderboard_rows = read_leaderboard(run_id, settings=settings)
    aggregate_rows = read_metric_aggregates(run_id, settings=settings)
    prediction_counts = read_prediction_counts(run_id, settings=settings)
    coverage_rows = read_backtest_coverage(run_id, settings=settings)
    comparable_rows = read_comparable_leaderboard(run_id, settings=settings)
    return _assemble_review(
        run_id,
        summary,
        decision_metric,
        n_series,
        leaderboard_rows,
        aggregate_rows,
        prediction_counts,
        coverage_rows,
        comparable_rows,
    )

diagnostics

Econometric & operational forecast diagnostics — pure DataFrame analytics (zero GCP imports).

Provides the post-backtest and post-run diagnostic suite for evaluating model value, lead-time accuracy decay, directional bias, and temporal/lead-time aggregation:

• Forecast Value Added (FVA) & Statistical Significance (compute_fva, diebold_mariano_test) — step-by-step accuracy waterfall across modeling tiers (Baseline -> Statistical -> ML -> Deep Learning -> AutoML -> Ensemble) with Harvey, Leybourne & Newbold (1997) modified Diebold–Mariano and paired Wilcoxon signed-rank $p$-values on out-of-fold loss differentials.
• Lead-Time Horizon Decay (compute_horizon_decay) — accuracy and empirical interval coverage at each forecast step $h = 1 \dots H$, plus the crossover horizon step $h^*$ where errors exceed the baseline.
• Directional Bias & Tracking Signal Audit (compute_bias_report) — per-(ts_id, model_type) audit of signed percentage bias, Trigg / APICS Tracking Signal ($\text{TS} = \sum e_t / \text{MAD}$), and consecutive over/under-forecasting fold streaks.
• Autocorrelation-Aware Temporal & Lead-Time Rollup (rollup_time) — aggregates point forecasts and actuals to Weekly ("W"), Monthly ("MS"/"M"), Quarterly ("QS"/"Q"), or fixed $L$-step supplier lead-time windows (freq="7D" or window=7) while computing period-level prediction intervals from empirical multi-step error covariance rather than naively summing daily interval bounds.
• Operational Inventory & Replenishment Policy Translators (compute_inventory_policy, plot_inventory_policy) — translates calibrated multi-step forecasts, AR(1) error covariance, and supplier lead-time variability ($\mu_L, \sigma_L$) into continuous-review $(s, S)$ and periodic-review $(R, S)$ Reorder Points (ROP), Safety Stock (SS), Order-Up-To Levels ($S$), Recommended Order Quantities ($Q$), and Newsvendor Expected Holding vs. Stockout Cost tables under Cycle Service Level ($\alpha$), Unit Fill Rate ($\beta$), or critical-fractile cost optimization.

compute_anomalies

compute_anomalies(
    oof_df=None,
    *,
    series_df=None,
    predictions_df=None,
    actuals_df=None,
    model_name=None,
    anomaly_prob_threshold=0.95,
    confidence_level=0.8,
    anomalies_only=False,
    freq="D",
)

Universal time-series anomaly detection across any of the 37 models, 9 ensembles, or STL.

Computes the exact two-sided Gaussian / conformal tail probability used by BigQuery ML (ML.DETECT_ANOMALIES):

$$p_{\text{anom}}(t) = 2\,\Phi\left(|y_t - \hat{y}_t| / \hat{\sigma}_t\right) - 1$$

where $\hat{\sigma}t = (U_t - L_t) / (2\,z)$ is recovered from each model's prediction interval $[L_t, U_t]$ at confidence_level (falling back to the per-series robust MAD scale $1.4826 \cdot \text{median}(|y - \hat{y}|)$ when intervals are absent or degenerate), and rescales lower_bound and upper_bound to anomaly_prob_threshold ($\tau \in (0, 1)$).

Supports three input modes: 1. Out-of-fold backtest predictions (oof_df): Evaluates any Python, AutoML, BigQuery, or Ensemble model (or model_name="all") over its out-of-fold predictions (backtest_oof). 2. Published forecasts + new actuals (predictions_df + actuals_df): Flags live anomalies as new actuals arrive against published forecast_predictions. 3. Model-free historical series (series_df or model_name="stl"): Decomposes raw historical series via seasonal-trend baseline + robust MAD residuals when no model OOF frame is provided.

Returns a tidy DataFrame with columns ["ts_id", "model_type", "ds", "y", "yhat", "is_anomaly", "anomaly_probability", "lower_bound", "upper_bound"].

Source code in src/scale_forecasting/diagnostics.py

def compute_anomalies(
    oof_df: pd.DataFrame | None = None,
    *,
    series_df: pd.DataFrame | None = None,
    predictions_df: pd.DataFrame | None = None,
    actuals_df: pd.DataFrame | None = None,
    model_name: str | None = None,
    anomaly_prob_threshold: float = 0.95,
    confidence_level: float = 0.80,
    anomalies_only: bool = False,
    freq: str = "D",
) -> pd.DataFrame:
    """Universal time-series anomaly detection across any of the 37 models, 9 ensembles, or STL.

    Computes the exact two-sided Gaussian / conformal tail probability used by BigQuery ML
    (`ML.DETECT_ANOMALIES`):

    $$p_{\\text{anom}}(t) = 2\\,\\Phi\\left(|y_t - \\hat{y}_t| / \\hat{\\sigma}_t\\right) - 1$$

    where $\\hat{\\sigma}_t = (U_t - L_t) / (2\\,z_{1-\\alpha/2})$ is recovered from each model's
    prediction interval $[L_t, U_t]$ at ``confidence_level`` (falling back to the per-series
    robust MAD scale $1.4826 \\cdot \\text{median}(|y - \\hat{y}|)$ when intervals are absent or
    degenerate), and rescales ``lower_bound`` and ``upper_bound`` to
    ``anomaly_prob_threshold`` ($\\tau \\in (0, 1)$).

    Supports three input modes:
    1. **Out-of-fold backtest predictions (``oof_df``):** Evaluates any Python, AutoML, BigQuery, or
       Ensemble model (or ``model_name="all"``) over its out-of-fold predictions (`backtest_oof`).
    2. **Published forecasts + new actuals (``predictions_df`` + ``actuals_df``):** Flags live
       anomalies as new actuals arrive against published `forecast_predictions`.
    3. **Model-free historical series (``series_df`` or ``model_name="stl"``):** Decomposes raw
       historical series via seasonal-trend baseline + robust MAD residuals when no model OOF frame
       is provided.

    Returns a tidy ``DataFrame`` with columns ``["ts_id", "model_type", "ds", "y", "yhat",
    "is_anomaly", "anomaly_probability", "lower_bound", "upper_bound"]``.
    """
    from .errors import ConfigError
    from .features import local_seasonal_baseline
    from .seasonality import seasonal_period

    if not (0.0 < float(anomaly_prob_threshold) < 1.0):
        raise ConfigError(
            f"anomaly_prob_threshold={anomaly_prob_threshold!r} must lie strictly in (0.0, 1.0)."
        )
    eff_conf = float(confidence_level) if 0.0 < float(confidence_level) < 1.0 else 0.80
    z_conf = float(stats.norm.ppf(0.5 + 0.5 * eff_conf))
    z_thresh = float(stats.norm.ppf(0.5 + 0.5 * float(anomaly_prob_threshold)))

    # 1. Join predictions_df + actuals_df if provided.
    work_df: pd.DataFrame | None = None
    if (
        predictions_df is not None
        and not predictions_df.empty
        and actuals_df is not None
        and not actuals_df.empty
    ):
        p = predictions_df.copy()
        a = actuals_df.copy()
        p_date = (
            "forecast_date"
            if "forecast_date" in p.columns
            else ("target_date" if "target_date" in p.columns else "ds")
        )
        a_date = (
            "ds"
            if "ds" in a.columns
            else ("target_date" if "target_date" in a.columns else "forecast_date")
        )
        a_val = "y" if "y" in a.columns else "y_true"
        p["ds"] = pd.to_datetime(p[p_date]).dt.date
        a["ds"] = pd.to_datetime(a[a_date]).dt.date
        if "ts_id" not in p.columns:
            p["ts_id"] = "series_0"
        if "ts_id" not in a.columns:
            a["ts_id"] = "series_0"
        if "model_type" not in p.columns:
            p["model_type"] = model_name or "model"
        merged = p.merge(
            a[["ts_id", "ds", a_val]].rename(columns={a_val: "y"}),
            on=["ts_id", "ds"],
            how="inner",
        )
        if not merged.empty:
            work_df = merged

    # 2. Otherwise use oof_df when available (unless model_name == "stl" explicitly).
    if (
        work_df is None
        and oof_df is not None
        and not oof_df.empty
        and (model_name is None or str(model_name).lower() != "stl")
    ):
        work_df = oof_df.copy()

    if work_df is not None and not work_df.empty:
        df = work_df.copy()
        if "ts_id" not in df.columns:
            df["ts_id"] = "series_0"
        if "model_type" not in df.columns:
            df["model_type"] = model_name or "model"
        date_col = (
            "ds"
            if "ds" in df.columns
            else ("target_date" if "target_date" in df.columns else "forecast_date")
        )
        y_col = "y" if "y" in df.columns else "y_true"
        yhat_col = "yhat" if "yhat" in df.columns else "forecast"
        df["ds"] = pd.to_datetime(df[date_col]).dt.date
        df["y"] = pd.to_numeric(df[y_col], errors="coerce")
        df["yhat"] = pd.to_numeric(df[yhat_col], errors="coerce")
        df["yhat_lower"] = (
            pd.to_numeric(df["yhat_lower"], errors="coerce")
            if "yhat_lower" in df.columns
            else np.nan
        )
        df["yhat_upper"] = (
            pd.to_numeric(df["yhat_upper"], errors="coerce")
            if "yhat_upper" in df.columns
            else np.nan
        )

        if model_name is not None and str(model_name).lower() not in ("all", "*"):
            req = str(model_name).strip()
            exact = df[df["model_type"] == req]
            if exact.empty and req == "arima_plus_xreg":
                exact = df[df["model_type"] == "arima_plus@xreg"]
            df = exact

        df = df.dropna(subset=["ds", "y", "yhat"])
        if df.empty:
            return pd.DataFrame(columns=_ANOMALY_COLUMNS)

        # Deduplicate overlapping backtest fold dates by keeping the latest fold.
        if "fold" in df.columns:
            df = (
                df.sort_values(["ts_id", "model_type", "ds", "fold"])
                .drop_duplicates(subset=["ts_id", "model_type", "ds"], keep="last")
                .reset_index(drop=True)
            )
        else:
            df = (
                df.sort_values(["ts_id", "model_type", "ds"])
                .drop_duplicates(subset=["ts_id", "model_type", "ds"], keep="last")
                .reset_index(drop=True)
            )

        out_frames: list[pd.DataFrame] = []
        for (tid, mtype), grp in df.groupby(["ts_id", "model_type"], sort=False):
            g = grp.sort_values("ds").copy()
            y_arr = g["y"].to_numpy(dtype=float)
            yhat_arr = g["yhat"].to_numpy(dtype=float)
            lo_arr = g["yhat_lower"].to_numpy(dtype=float)
            hi_arr = g["yhat_upper"].to_numpy(dtype=float)

            abs_err = np.abs(y_arr - yhat_arr)
            med_err = float(np.nanmedian(abs_err))
            mad_sigma = max(1.4826 * med_err, float(np.nanstd(y_arr - yhat_arr)), 1e-6)

            width = hi_arr - lo_arr
            sigma = np.where(
                np.isfinite(width) & (width > 1e-9),
                width / (2.0 * max(z_conf, 1e-6)),
                mad_sigma,
            )
            sigma = np.maximum(sigma, 1e-6)

            z_val = abs_err / sigma
            prob = np.clip(2.0 * stats.norm.cdf(z_val) - 1.0, 0.0, 1.0)
            lower_b = yhat_arr - z_thresh * sigma
            upper_b = yhat_arr + z_thresh * sigma
            is_anom = (
                (prob >= float(anomaly_prob_threshold)) | (y_arr < lower_b) | (y_arr > upper_b)
            )

            out_frames.append(
                pd.DataFrame(
                    {
                        "ts_id": str(tid),
                        "model_type": str(mtype),
                        "ds": g["ds"].to_numpy(),
                        "y": y_arr,
                        "yhat": yhat_arr,
                        "is_anomaly": is_anom.astype(bool),
                        "anomaly_probability": np.round(prob, 6),
                        "lower_bound": lower_b,
                        "upper_bound": upper_b,
                    }
                )
            )
        if not out_frames:
            return pd.DataFrame(columns=_ANOMALY_COLUMNS)
        res = pd.concat(out_frames, ignore_index=True)
        if anomalies_only:
            res = res.loc[res["is_anomaly"] == True].reset_index(drop=True)  # noqa: E712
        return res

    # 3. Model-free seasonal-trend + MAD decomposition on series_df.
    if series_df is None or series_df.empty:
        return pd.DataFrame(columns=_ANOMALY_COLUMNS)

    s = series_df.copy()
    if "ts_id" not in s.columns:
        s["ts_id"] = "series_0"
    date_c = "ds" if "ds" in s.columns else ("date" if "date" in s.columns else s.columns[1])
    y_c = "y" if "y" in s.columns else ("y_true" if "y_true" in s.columns else s.columns[-1])
    s["ds"] = pd.to_datetime(s[date_c]).dt.date
    s["y"] = pd.to_numeric(s[y_c], errors="coerce")
    s = s.dropna(subset=["ds", "y"]).sort_values(["ts_id", "ds"])
    if s.empty:
        return pd.DataFrame(columns=_ANOMALY_COLUMNS)

    period = max(1, int(seasonal_period(freq)))
    label_model = str(model_name) if model_name and str(model_name).lower() != "all" else "stl"
    stl_frames: list[pd.DataFrame] = []
    for tid, grp in s.groupby("ts_id", sort=False):
        g = grp.sort_values("ds")
        y_arr = g["y"].to_numpy(dtype=float)
        baseline = local_seasonal_baseline(y_arr, period)
        resid = y_arr - baseline
        med = float(np.nanmedian(resid))
        yhat_arr = baseline + med
        abs_dev = np.abs(y_arr - yhat_arr)
        mad_sigma = max(1.4826 * float(np.nanmedian(abs_dev)), float(np.nanstd(resid)), 1e-6)
        z_val = abs_dev / mad_sigma
        prob = np.clip(2.0 * stats.norm.cdf(z_val) - 1.0, 0.0, 1.0)
        lower_b = yhat_arr - z_thresh * mad_sigma
        upper_b = yhat_arr + z_thresh * mad_sigma
        is_anom = (prob >= float(anomaly_prob_threshold)) | (y_arr < lower_b) | (y_arr > upper_b)
        stl_frames.append(
            pd.DataFrame(
                {
                    "ts_id": str(tid),
                    "model_type": label_model,
                    "ds": g["ds"].to_numpy(),
                    "y": y_arr,
                    "yhat": yhat_arr,
                    "is_anomaly": is_anom.astype(bool),
                    "anomaly_probability": np.round(prob, 6),
                    "lower_bound": lower_b,
                    "upper_bound": upper_b,
                }
            )
        )
    res_stl = pd.concat(stl_frames, ignore_index=True)
    if anomalies_only:
        res_stl = res_stl.loc[res_stl["is_anomaly"] == True].reset_index(drop=True)  # noqa: E712
    return res_stl

compute_bias_report

compute_bias_report(
    oof_df,
    *,
    threshold_pct=0.05,
    tracking_signal_limit=4.0,
    min_consecutive_folds=2,
)

Audit every (ts_id, model_type) cell for chronic over- or under-forecasting bias.

Computes signed error $\hat{y}_t - y_t$, percentage bias $\sum(\hat{y}_t - y_t)/\sum|y_t|$, the Trigg / APICS Tracking Signal $\text{TS} = \sum(\hat{y}_t - y_t) / \text{MAD}$, and the longest streak of consecutive folds with positive (over-forecasting) or negative (under-forecasting) fold-mean error.

Source code in src/scale_forecasting/diagnostics.py

def compute_bias_report(
    oof_df: pd.DataFrame,
    *,
    threshold_pct: float = 0.05,
    tracking_signal_limit: float = 4.0,
    min_consecutive_folds: int = 2,
) -> pd.DataFrame:
    """Audit every ``(ts_id, model_type)`` cell for chronic over- or under-forecasting bias.

    Computes signed error $\\hat{y}_t - y_t$, percentage bias $\\sum(\\hat{y}_t - y_t)/\\sum|y_t|$,
    the **Trigg / APICS Tracking Signal** $\\text{TS} = \\sum(\\hat{y}_t - y_t) / \\text{MAD}$, and
    the longest streak of consecutive folds with positive (over-forecasting) or negative
    (under-forecasting) fold-mean error.
    """
    cols = [
        "ts_id",
        "model_type",
        "n_obs",
        "n_folds",
        "mean_y_true",
        "mean_yhat",
        "mean_signed_error",
        "bias_pct",
        "mad",
        "tracking_signal",
        "max_consecutive_over_folds",
        "max_consecutive_under_folds",
        "bias_status",
    ]
    if oof_df is None or oof_df.empty:
        return pd.DataFrame(columns=cols)

    df = _ensure_oof_columns(oof_df)
    records: list[dict[str, Any]] = []
    for (tid, mtype), grp in df.groupby(["ts_id", "model_type"], sort=True):
        y_true = pd.to_numeric(grp["y_true"], errors="coerce").to_numpy(dtype=float)
        yhat = pd.to_numeric(grp["yhat"], errors="coerce").to_numpy(dtype=float)
        keep = np.isfinite(y_true) & np.isfinite(yhat)
        if not keep.any():
            continue
        yt, yh = y_true[keep], yhat[keep]
        signed_err = yh - yt
        n_obs = int(keep.sum())
        mean_yt = float(np.mean(yt))
        mean_yh = float(np.mean(yh))
        mean_se = float(np.mean(signed_err))
        abs_sum = float(np.sum(np.abs(yt)))
        bias_pct = float(np.sum(signed_err) / abs_sum) if abs_sum > 0 else 0.0
        mad = float(np.mean(np.abs(signed_err)))
        tracking_sig = float(np.sum(signed_err) / mad) if mad > 1e-12 else 0.0

        # Consecutive fold streaks.
        sub_valid = grp.loc[keep].copy()
        sub_valid["_se"] = signed_err
        fold_means = sub_valid.groupby("fold_id", sort=True)["_se"].mean().to_numpy(dtype=float)
        max_over, max_under, cur_over, cur_under = 0, 0, 0, 0
        for fm in fold_means:
            if fm > 0:
                cur_over += 1
                cur_under = 0
                max_over = max(max_over, cur_over)
            elif fm < 0:
                cur_under += 1
                cur_over = 0
                max_under = max(max_under, cur_under)
            else:
                cur_over, cur_under = 0, 0

        exceeds_mag = abs(bias_pct) >= threshold_pct
        exceeds_ts = abs(tracking_sig) >= tracking_signal_limit
        if exceeds_mag and (exceeds_ts or max_over >= min_consecutive_folds) and bias_pct > 0:
            status = "over_forecasting"
        elif exceeds_mag and (exceeds_ts or max_under >= min_consecutive_folds) and bias_pct < 0:
            status = "under_forecasting"
        else:
            status = "unbiased"

        records.append(
            {
                "ts_id": str(tid),
                "model_type": str(mtype),
                "n_obs": n_obs,
                "n_folds": len(fold_means),
                "mean_y_true": mean_yt,
                "mean_yhat": mean_yh,
                "mean_signed_error": mean_se,
                "bias_pct": bias_pct,
                "mad": mad,
                "tracking_signal": tracking_sig,
                "max_consecutive_over_folds": max_over,
                "max_consecutive_under_folds": max_under,
                "bias_status": status,
            }
        )

    out = pd.DataFrame.from_records(records, columns=cols)
    return out.sort_values(
        "tracking_signal", key=lambda s: s.abs(), ascending=False, ignore_index=True
    )

compute_change_points

compute_change_points(
    series_df=None,
    *,
    oof_df=None,
    model_name=None,
    exclude_gap_edges=True,
    min_gap_steps=1,
    freq="D",
    sigma_threshold=3.0,
    min_segment=7,
)

Detect sustained structural change-point windows on raw series or model residuals (pure).

• When model_name is None (or "raw") and series_df is provided, segments the raw target series $y_t$ to identify sustained level shifts (model_type = "raw").
• When oof_df is provided with model_name (e.g. "xgboost", "tide", "arima_plus@xreg", "ensemble_champion", or "all"), segments each model's residuals ($e_t = y_t - \hat{y}_t$) to identify structural regime shifts that the model failed to absorb (window_avg, window_min, window_max, and window_stddev are reported in the series' original $y_t$ units over the detected break window so they remain directly comparable to BigQuery ML.DETECT_CHANGE_POINTS).

Annotates each detected window with gap_edge_missing_steps and is_gap_edge (filtering out ingestion outage boundaries when exclude_gap_edges=True).

Returns a tidy DataFrame with columns ["ts_id", "model_type", "begin_date", "end_date", "window_count", "window_avg", "window_min", "window_max", "window_stddev", "gap_edge_missing_steps", "is_gap_edge", "status"].

Source code in src/scale_forecasting/diagnostics.py

def compute_change_points(
    series_df: pd.DataFrame | None = None,
    *,
    oof_df: pd.DataFrame | None = None,
    model_name: str | None = None,
    exclude_gap_edges: bool = True,
    min_gap_steps: int = 1,
    freq: str = "D",
    sigma_threshold: float = 3.0,
    min_segment: int = 7,
) -> pd.DataFrame:
    """Detect sustained structural change-point windows on raw series or model residuals (pure).

    - When ``model_name`` is ``None`` (or ``"raw"``) and ``series_df`` is provided, segments the
      raw target series $y_t$ to identify sustained level shifts (`model_type = "raw"`).
    - When ``oof_df`` is provided with ``model_name`` (e.g. ``"xgboost"``, ``"tide"``,
      ``"arima_plus@xreg"``, ``"ensemble_champion"``, or ``"all"``), segments each model's
      **residuals ($e_t = y_t - \\hat{y}_t$)** to identify structural regime shifts that the model
      failed to absorb (`window_avg`, `window_min`, `window_max`, and `window_stddev` are reported
      in the series' original $y_t$ units over the detected break window so they remain directly
      comparable to BigQuery `ML.DETECT_CHANGE_POINTS`).

    Annotates each detected window with ``gap_edge_missing_steps`` and ``is_gap_edge`` (filtering
    out ingestion outage boundaries when ``exclude_gap_edges=True``).

    Returns a tidy ``DataFrame`` with columns ``["ts_id", "model_type", "begin_date", "end_date",
    "window_count", "window_avg", "window_min", "window_max", "window_stddev",
    "gap_edge_missing_steps", "is_gap_edge", "status"]``.
    """
    from .errors import ConfigError

    if isinstance(min_gap_steps, bool) or not isinstance(min_gap_steps, int) or min_gap_steps < 1:
        raise ConfigError(f"min_gap_steps={min_gap_steps!r} must be an integer >= 1.")

    expected_step = _freq_timedelta(freq)
    half_win = max(2, min(int(min_segment) // 2, 5))
    records: list[dict[str, Any]] = []

    def _scan_group(
        tid: str,
        mtype: str,
        dates_ts: pd.Series,
        signal_arr: np.ndarray,
        raw_y_arr: np.ndarray,
    ) -> None:
        n = int(signal_arr.size)
        if n < 2 * min_segment:
            return
        breaks = _segment_change_indices(
            signal_arr,
            min_segment=max(3, int(min_segment)),
            sigma_threshold=float(sigma_threshold),
        )
        if not breaks:
            return
        dt_vals = pd.to_datetime(dates_ts).reset_index(drop=True)
        diffs = dt_vals.diff()
        step_days = max(1.0, float(expected_step.total_seconds()) / 86400.0)
        missing_steps_arr = np.maximum(
            0,
            np.floor(
                diffs.dt.total_seconds().fillna(0.0).to_numpy() / (step_days * 86400.0) + 1e-6
            ).astype(int)
            - 1,
        )

        for b_idx in breaks:
            w_start = max(0, b_idx - half_win)
            w_end = min(n - 1, b_idx + half_win)
            win_y = raw_y_arr[w_start : w_end + 1]
            # Check whether any gap edge falls inside [w_start, w_end]
            win_gaps = missing_steps_arr[w_start : w_end + 1]
            qualifying_gaps = win_gaps[win_gaps >= int(min_gap_steps)]
            max_gap = int(qualifying_gaps.max()) if qualifying_gaps.size > 0 else 0
            is_gap = bool(max_gap > 0)
            if exclude_gap_edges and is_gap:
                continue
            records.append(
                {
                    "ts_id": str(tid),
                    "model_type": str(mtype),
                    "begin_date": dt_vals.iloc[w_start].date(),
                    "end_date": dt_vals.iloc[w_end].date(),
                    "window_count": int(win_y.size),
                    "window_avg": float(np.mean(win_y)),
                    "window_min": float(np.min(win_y)),
                    "window_max": float(np.max(win_y)),
                    "window_stddev": float(np.std(win_y, ddof=1)) if win_y.size > 1 else 0.0,
                    "gap_edge_missing_steps": max_gap,
                    "is_gap_edge": is_gap,
                    "status": "",
                }
            )

    # 1. Model-residual mode when oof_df is provided and (model_name is set or series_df is None)
    if oof_df is not None and not oof_df.empty and (model_name is not None or series_df is None):
        df = oof_df.copy()
        if "ts_id" not in df.columns:
            df["ts_id"] = "series_0"
        if "model_type" not in df.columns:
            df["model_type"] = model_name or "model"
        date_col = (
            "ds"
            if "ds" in df.columns
            else ("target_date" if "target_date" in df.columns else "forecast_date")
        )
        y_col = "y" if "y" in df.columns else "y_true"
        yhat_col = "yhat" if "yhat" in df.columns else "forecast"
        df["ds"] = pd.to_datetime(df[date_col])
        df["y"] = pd.to_numeric(df[y_col], errors="coerce")
        df["yhat"] = pd.to_numeric(df[yhat_col], errors="coerce")
        if model_name is not None and str(model_name).lower() not in ("all", "*", "raw"):
            req = str(model_name).strip()
            exact = df[df["model_type"] == req]
            if exact.empty and req == "arima_plus_xreg":
                exact = df[df["model_type"] == "arima_plus@xreg"]
            df = exact
        df = df.dropna(subset=["ds", "y", "yhat"])
        if "fold" in df.columns:
            df = (
                df.sort_values(["ts_id", "model_type", "ds", "fold"])
                .drop_duplicates(subset=["ts_id", "model_type", "ds"], keep="last")
                .reset_index(drop=True)
            )
        else:
            df = (
                df.sort_values(["ts_id", "model_type", "ds"])
                .drop_duplicates(subset=["ts_id", "model_type", "ds"], keep="last")
                .reset_index(drop=True)
            )
        for (tid, mtype), grp in df.groupby(["ts_id", "model_type"], sort=False):
            g = grp.sort_values("ds").reset_index(drop=True)
            y_arr = g["y"].to_numpy(dtype=float)
            yhat_arr = g["yhat"].to_numpy(dtype=float)
            resid_arr = y_arr - yhat_arr
            _scan_group(str(tid), str(mtype), g["ds"], resid_arr, y_arr)
        if not records:
            return pd.DataFrame(columns=_CHANGE_POINT_COLUMNS)
        return pd.DataFrame.from_records(records, columns=_CHANGE_POINT_COLUMNS).sort_values(
            ["ts_id", "model_type", "begin_date"], ignore_index=True
        )

    # 2. Raw series mode
    if series_df is None or series_df.empty:
        return pd.DataFrame(columns=_CHANGE_POINT_COLUMNS)
    s = series_df.copy()
    if "ts_id" not in s.columns:
        s["ts_id"] = "series_0"
    date_c = "ds" if "ds" in s.columns else ("date" if "date" in s.columns else s.columns[1])
    y_c = "y" if "y" in s.columns else ("y_true" if "y_true" in s.columns else s.columns[-1])
    s["ds"] = pd.to_datetime(s[date_c])
    s["y"] = pd.to_numeric(s[y_c], errors="coerce")
    s = s.dropna(subset=["ds", "y"]).sort_values(["ts_id", "ds"])
    label_m = str(model_name) if model_name else "raw"
    for tid, grp in s.groupby("ts_id", sort=False):
        g = grp.sort_values("ds").reset_index(drop=True)
        y_arr = g["y"].to_numpy(dtype=float)
        _scan_group(str(tid), label_m, g["ds"], y_arr, y_arr)
    if not records:
        return pd.DataFrame(columns=_CHANGE_POINT_COLUMNS)
    return pd.DataFrame.from_records(records, columns=_CHANGE_POINT_COLUMNS).sort_values(
        ["ts_id", "model_type", "begin_date"], ignore_index=True
    )

compute_exceptions_queue

compute_exceptions_queue(
    *,
    oof_df=None,
    predictions_df=None,
    realized_df=None,
    stability_df=None,
    floor=None,
    cap=None,
    metric="wape",
    baseline=None,
    top_n=50,
)

Build the Top-$N$ Exception Triage Queue across 5 operational forecasting rules.

Evaluates each (ts_id, model_type) cell and assigns a composite severity_score and plain-English exception_reasons across:

1. worse_than_naive: Model's held-out backtest loss exceeds the naive benchmark on this series.
2. chronic_bias: Trigg Tracking Signal $|\text{TS}| \ge 4.0$ or $\ge 2$ consecutive biased folds (chronic over- or under-forecasting).
3. realized_drift: Realized live error exceeds $1.5\times$ expected backtest error.
4. whipsaw_jump: Cross-run forecast revision jump $\ge 30\%$ on an overlapping date.
5. capacity_ceiling_hit: Forward forecast hits or clips against physical floor / cap bounds on $\ge 20\%$ of horizon steps.

Source code in src/scale_forecasting/diagnostics.py

def compute_exceptions_queue(
    *,
    oof_df: pd.DataFrame | None = None,
    predictions_df: pd.DataFrame | None = None,
    realized_df: pd.DataFrame | None = None,
    stability_df: pd.DataFrame | None = None,
    floor: float | None = None,
    cap: float | None = None,
    metric: str = "wape",
    baseline: str | None = None,
    top_n: int = 50,
) -> pd.DataFrame:
    """Build the Top-$N$ Exception Triage Queue across 5 operational forecasting rules.

    Evaluates each ``(ts_id, model_type)`` cell and assigns a composite ``severity_score`` and
    plain-English ``exception_reasons`` across:

    1. ``worse_than_naive``: Model's held-out backtest loss exceeds the naive benchmark on this
       series.
    2. ``chronic_bias``: Trigg Tracking Signal $|\\text{TS}| \\ge 4.0$ or $\\ge 2$ consecutive
       biased folds (chronic over- or under-forecasting).
    3. ``realized_drift``: Realized live error exceeds $1.5\\times$ expected backtest error.
    4. ``whipsaw_jump``: Cross-run forecast revision jump $\\ge 30\\%$ on an overlapping date.
    5. ``capacity_ceiling_hit``: Forward forecast hits or clips against physical ``floor`` / ``cap``
       bounds on $\\ge 20\\%$ of horizon steps.
    """
    cols = [
        "ts_id",
        "model_type",
        "severity_score",
        "n_flags",
        "exception_flags",
        "exception_reasons",
        "metric_loss",
        "baseline_loss",
        "bias_pct",
        "tracking_signal",
        "drift_ratio",
        "max_jump_pct",
    ]
    cell_flags: dict[tuple[str, str], list[str]] = {}
    cell_reasons: dict[tuple[str, str], list[str]] = {}
    cell_sev: dict[tuple[str, str], float] = {}
    cell_meta: dict[tuple[str, str], dict[str, float]] = {}

    def _touch(key: tuple[str, str]) -> None:
        if key not in cell_flags:
            cell_flags[key] = []
            cell_reasons[key] = []
            cell_sev[key] = 0.0
            cell_meta[key] = {
                "metric_loss": float("nan"),
                "baseline_loss": float("nan"),
                "bias_pct": float("nan"),
                "tracking_signal": float("nan"),
                "drift_ratio": float("nan"),
                "max_jump_pct": float("nan"),
            }

    # 1 & 2. Backtest OOF checks: worse_than_naive + chronic_bias
    if oof_df is not None and not oof_df.empty:
        oof_norm = _ensure_oof_columns(oof_df)
        all_models = list(dict.fromkeys(oof_norm["model_type"].astype(str)))
        base_m = _resolve_baseline_model(all_models, baseline)

        losses: dict[tuple[str, str], float] = {}
        for (tid, mtype), grp in oof_norm.groupby(["ts_id", "model_type"], sort=False):
            losses[(str(tid), str(mtype))] = _cell_metric_loss(grp, metric)

        for (tid_s, mtype_s), m_loss in losses.items():
            b_loss = losses.get((tid_s, base_m), float("nan"))
            key = (tid_s, mtype_s)
            _touch(key)
            cell_meta[key]["metric_loss"] = m_loss
            cell_meta[key]["baseline_loss"] = b_loss
            if (
                mtype_s != base_m
                and math.isfinite(m_loss)
                and math.isfinite(b_loss)
                and m_loss > b_loss * 1.02
            ):
                ratio = m_loss / max(b_loss, 1e-6)
                cell_flags[key].append("worse_than_naive")
                cell_reasons[key].append(
                    f"Backtest {metric.upper()} ({m_loss:.3f}) is worse than {base_m} "
                    f"({b_loss:.3f}, {ratio:.2f}x)."
                )
                cell_sev[key] += min(40.0, 15.0 * ratio)

        bias_df = compute_bias_report(oof_norm)
        for b_row in bias_df.to_dict(orient="records"):
            key = (str(b_row["ts_id"]), str(b_row["model_type"]))
            _touch(key)
            b_pct = float(b_row["bias_pct"])
            t_sig = float(b_row["tracking_signal"])
            cell_meta[key]["bias_pct"] = b_pct
            cell_meta[key]["tracking_signal"] = t_sig
            if str(b_row["bias_status"]) != "unbiased":
                direction = (
                    "over-forecasting"
                    if b_row["bias_status"] == "over_forecasting"
                    else "under-forecasting"
                )
                cell_flags[key].append("chronic_bias")
                cell_reasons[key].append(
                    f"Chronic {direction}: bias={b_pct * 100:+.1f}%, Trigg TS={t_sig:+.1f}."
                )
                cell_sev[key] += min(35.0, 5.0 * abs(t_sig))

    # 3. Realized drift check
    if realized_df is not None and not realized_df.empty:
        for row in realized_df.itertuples(index=False):
            key = (str(row.ts_id), str(row.model_type))
            _touch(key)
            d_ratio = float(getattr(row, "drift_ratio", float("nan")))
            cell_meta[key]["drift_ratio"] = d_ratio
            if bool(getattr(row, "drift_flagged", False)):
                cell_flags[key].append("realized_drift")
                cell_reasons[key].append(
                    f"Live accuracy drift: realized loss is {d_ratio:.2f}x backtest expectation."
                )
                cell_sev[key] += min(45.0, 20.0 * d_ratio)

    # 4. Whipsaw revision stability check
    if stability_df is not None and not stability_df.empty:
        for row in stability_df.itertuples(index=False):
            key = (str(row.ts_id), str(row.model_type))
            _touch(key)
            jump = float(getattr(row, "max_single_step_jump_pct", float("nan")))
            cell_meta[key]["max_jump_pct"] = jump
            if bool(getattr(row, "whipsaw_flagged", False)):
                cell_flags[key].append("whipsaw_jump")
                cell_reasons[key].append(
                    f"Forecast whipsaw: {jump * 100:.1f}% single-step jump across consecutive runs."
                )
                cell_sev[key] += min(30.0, 40.0 * jump)

    # 5. Physical capacity ceiling / floor hit check
    if (
        predictions_df is not None
        and not predictions_df.empty
        and (floor is not None or cap is not None)
    ):
        pred = predictions_df.copy()
        if "segment" in pred.columns and (pred["segment"] == "forecast").any():
            pred = pred[pred["segment"] == "forecast"].copy()
        yhat_col = "yhat" if "yhat" in pred.columns else "forecast"
        if "ts_id" not in pred.columns:
            pred["ts_id"] = "series_0"
        if "model_type" not in pred.columns:
            pred["model_type"] = "model"
        for (tid, mtype), grp in pred.groupby(["ts_id", "model_type"], sort=False):
            yh = pd.to_numeric(grp[yhat_col], errors="coerce").dropna().to_numpy(dtype=float)
            if yh.size == 0:
                continue
            hit_mask = np.zeros(yh.size, dtype=bool)
            if cap is not None:
                hit_mask |= yh >= (float(cap) - 1e-6)
            if floor is not None and float(floor) > 0.0:
                hit_mask |= yh <= (float(floor) + 1e-6)
            hit_share = float(np.mean(hit_mask))
            if hit_share >= 0.20:
                key = (str(tid), str(mtype))
                _touch(key)
                cell_flags[key].append("capacity_ceiling_hit")
                cell_reasons[key].append(
                    f"Capacity bound reached on {hit_share * 100:.0f}% of horizon steps."
                )
                cell_sev[key] += 25.0 * hit_share

    records: list[dict[str, Any]] = []
    for key, flags in cell_flags.items():
        if not flags:
            continue
        tid, mtype = key
        meta = cell_meta[key]
        records.append(
            {
                "ts_id": tid,
                "model_type": mtype,
                "severity_score": round(cell_sev[key], 2),
                "n_flags": len(flags),
                "exception_flags": ",".join(flags),
                "exception_reasons": " ".join(cell_reasons[key]),
                "metric_loss": meta["metric_loss"],
                "baseline_loss": meta["baseline_loss"],
                "bias_pct": meta["bias_pct"],
                "tracking_signal": meta["tracking_signal"],
                "drift_ratio": meta["drift_ratio"],
                "max_jump_pct": meta["max_jump_pct"],
            }
        )

    if not records:
        return pd.DataFrame(columns=cols)
    out = pd.DataFrame.from_records(records, columns=cols)
    out = out.sort_values(
        ["severity_score", "n_flags"], ascending=[False, False], ignore_index=True
    )
    return out.head(max(1, int(top_n)))

compute_fva

compute_fva(
    oof_df, *, baseline=None, metric="wape", series_df=None
)

Compute the Forecast Value Added (FVA) Waterfall & Statistical Significance table.

Evaluates each model_type in oof_df against baseline (defaulting to "naive_seasonal" or the simplest benchmark present) and reports:

• fva_stage_rank / fva_stage: position on the FVA ladder (Baseline, Statistical, Machine Learning, Deep Learning, AutoML, Ensemble).
• metric: the decision metric name (e.g. "wape").
• mean_metric / pooled_wape / mean_mae: accuracy across all scored series.
• baseline_model / baseline_metric: benchmark reference loss.
• fva_abs ($L_{\text{baseline}} - L_{\text{model}}$) and fva_pct ($(L_{\text{baseline}} - L_{\text{model}}) / L_{\text{baseline}}$): positive values mean the model added value over the baseline; negative values mean it destroyed value.
• win_rate_vs_baseline: fraction of series where model_type achieved strictly lower loss than baseline.
• dm_stat, p_value_dm, p_value_wilcoxon, significant_at_05: Harvey–Leybourne– Newbold (1997) modified Diebold–Mariano statistic/p-value and paired Wilcoxon signed-rank p-value testing whether the accuracy gain over baseline is statistically significant.

Source code in src/scale_forecasting/diagnostics.py

def compute_fva(
    oof_df: pd.DataFrame,
    *,
    baseline: str | None = None,
    metric: str = "wape",
    series_df: pd.DataFrame | None = None,
) -> pd.DataFrame:
    """Compute the Forecast Value Added (FVA) Waterfall & Statistical Significance table.

    Evaluates each ``model_type`` in ``oof_df`` against ``baseline`` (defaulting to
    ``"naive_seasonal"`` or the simplest benchmark present) and reports:

    - ``fva_stage_rank`` / ``fva_stage``: position on the FVA ladder (`Baseline`, `Statistical`,
      `Machine Learning`, `Deep Learning`, `AutoML`, `Ensemble`).
    - ``metric``: the decision metric name (e.g. ``"wape"``).
    - ``mean_metric`` / ``pooled_wape`` / ``mean_mae``: accuracy across all scored series.
    - ``baseline_model`` / ``baseline_metric``: benchmark reference loss.
    - ``fva_abs`` ($L_{\\text{baseline}} - L_{\\text{model}}$) and ``fva_pct``
      ($(L_{\\text{baseline}} - L_{\\text{model}}) / L_{\\text{baseline}}$): positive values mean
      the model added value over the baseline; negative values mean it destroyed value.
    - ``win_rate_vs_baseline``: fraction of series where ``model_type`` achieved strictly lower
      loss than ``baseline``.
    - ``dm_stat``, ``p_value_dm``, ``p_value_wilcoxon``, ``significant_at_05``: Harvey–Leybourne–
      Newbold (1997) modified Diebold–Mariano statistic/p-value and paired Wilcoxon signed-rank
      p-value testing whether the accuracy gain over ``baseline`` is statistically significant.
    """
    from .models import base_model_name

    cols = [
        "fva_stage_rank",
        "fva_stage",
        "model_type",
        "base_model",
        "demand_class",
        "n_series",
        "n_points",
        "metric",
        "mean_metric",
        "pooled_wape",
        "mean_mae",
        "baseline_model",
        "baseline_metric",
        "fva_abs",
        "fva_pct",
        "win_rate_vs_baseline",
        "dm_stat",
        "p_value_dm",
        "p_value_wilcoxon",
        "significant_at_05",
        "univariate_baseline_model",
        "covariate_fva_abs",
        "covariate_fva_pct",
        "covariate_p_value_dm",
        "covariate_significant_at_05",
        "local_baseline_model",
        "pooling_fva_abs",
        "pooling_fva_pct",
        "pooling_p_value_dm",
        "pooling_significant_at_05",
    ]
    if oof_df is None or oof_df.empty:
        return pd.DataFrame(columns=cols)

    df = _ensure_oof_columns(oof_df)
    if series_df is not None and not series_df.empty and "demand_class" not in df.columns:
        from .features import compute_demand_profile

        ts_col = "ts_id" if "ts_id" in series_df.columns else None
        y_col = "y" if "y" in series_df.columns else None
        if ts_col and y_col:
            regime_map = {
                str(tid): str(compute_demand_profile(grp[y_col])["demand_class"])
                for tid, grp in series_df.groupby(ts_col, sort=False)
            }
            df["demand_class"] = df["ts_id"].astype(str).map(regime_map).fillna("all")
    if "demand_class" not in df.columns:
        df["demand_class"] = "all"

    all_models = list(dict.fromkeys(df["model_type"].astype(str)))
    base_model = _resolve_baseline_model(all_models, baseline)

    # Compute per-(ts_id, model_type) loss and per-(ts_id, fold_id, forecast_date, model_type)
    # pointwise absolute errors so DM tests have power even on single-series panels.
    series_losses: dict[tuple[str, str], float] = {}
    fold_losses: dict[tuple[str, int, str], float] = {}
    for (tid, mtype), grp in df.groupby(["ts_id", "model_type"], sort=False):
        series_losses[(str(tid), str(mtype))] = _cell_metric_loss(grp, metric)
    for (tid, fid, mtype), grp in df.groupby(["ts_id", "fold_id", "model_type"], sort=False):
        fold_losses[(str(tid), int(str(fid)), str(mtype))] = _cell_metric_loss(grp, metric)

    df["_abs_err"] = np.abs(
        pd.to_numeric(df["y_true"], errors="coerce") - pd.to_numeric(df["yhat"], errors="coerce")
    )
    point_pivot = df.pivot_table(
        index=["ts_id", "fold_id", "forecast_date"],
        columns="model_type",
        values="_abs_err",
        aggfunc="mean",
    )

    all_ts = sorted({str(t) for t in df["ts_id"].unique()})
    all_folds = sorted({(str(t), int(f)) for t, f in zip(df["ts_id"], df["fold_id"], strict=False)})

    def _compare_pair(
        cand_model: str,
        ref_model: str | None,
        cand_series_vec: np.ndarray,
        cand_mean: float,
        max_h: int,
    ) -> tuple[float, float, float, float, float, float, float, bool]:
        if ref_model is None:
            return (
                float("nan"),
                float("nan"),
                float("nan"),
                float("nan"),
                float("nan"),
                float("nan"),
                float("nan"),
                False,
            )
        ref_series_vec = np.array(
            [series_losses.get((tid, ref_model), float("nan")) for tid in all_ts], dtype=float
        )
        ref_mean = (
            float(np.nanmean(ref_series_vec)) if np.isfinite(ref_series_vec).any() else float("nan")
        )
        both_series = np.isfinite(cand_series_vec) & np.isfinite(ref_series_vec)
        if cand_model == ref_model:
            win_rate = 0.5
            dm_stat, p_dm, p_wilcox = 0.0, 1.0, 1.0
        else:
            win_rate = (
                float(np.mean(cand_series_vec[both_series] < ref_series_vec[both_series]))
                if both_series.any()
                else float("nan")
            )
            if int(both_series.sum()) >= 5:
                dm_stat, p_dm = diebold_mariano_test(cand_series_vec, ref_series_vec, h=1)
                p_wilcox = wilcoxon_pvalue(cand_series_vec, ref_series_vec)
            else:
                m_fold_vec = np.array(
                    [
                        fold_losses.get((tid, fid, cand_model), float("nan"))
                        for tid, fid in all_folds
                    ],
                    dtype=float,
                )
                b_fold_vec = np.array(
                    [
                        fold_losses.get((tid, fid, ref_model), float("nan"))
                        for tid, fid in all_folds
                    ],
                    dtype=float,
                )
                both_folds = np.isfinite(m_fold_vec) & np.isfinite(b_fold_vec)
                if int(both_folds.sum()) >= 5:
                    dm_stat, p_dm = diebold_mariano_test(m_fold_vec, b_fold_vec, h=1)
                    p_wilcox = wilcoxon_pvalue(m_fold_vec, b_fold_vec)
                elif cand_model in point_pivot.columns and ref_model in point_pivot.columns:
                    mp = point_pivot[cand_model].to_numpy(dtype=float)
                    bp = point_pivot[ref_model].to_numpy(dtype=float)
                    dm_stat, p_dm = diebold_mariano_test(mp, bp, h=max(1, min(max_h, 7)))
                    p_wilcox = wilcoxon_pvalue(mp, bp)
                else:
                    dm_stat, p_dm, p_wilcox = float("nan"), float("nan"), float("nan")
        fva_abs = (
            float(ref_mean - cand_mean)
            if math.isfinite(ref_mean) and math.isfinite(cand_mean)
            else float("nan")
        )
        fva_pct = (
            float(fva_abs / ref_mean) if math.isfinite(fva_abs) and ref_mean > 0 else float("nan")
        )
        sig_05 = bool(
            cand_model != ref_model
            and math.isfinite(fva_abs)
            and fva_abs > 0.0
            and math.isfinite(p_dm)
            and p_dm < 0.05
        )
        return ref_mean, fva_abs, fva_pct, win_rate, dm_stat, p_dm, p_wilcox, sig_05

    records: list[dict[str, Any]] = []
    for mtype in all_models:
        stage_rank, stage_label = _classify_model_stage(mtype, base_model)
        m_df = df[df["model_type"] == mtype]
        y_true = pd.to_numeric(m_df["y_true"], errors="coerce").to_numpy(dtype=float)
        yhat = pd.to_numeric(m_df["yhat"], errors="coerce").to_numpy(dtype=float)
        keep = np.isfinite(y_true) & np.isfinite(yhat)
        n_points = int(keep.sum())
        n_series = int(m_df["ts_id"].nunique())

        denom = float(np.sum(np.abs(y_true[keep]))) if n_points else 0.0
        pooled_wape = (
            float(np.sum(np.abs(y_true[keep] - yhat[keep])) / denom) if denom > 0 else float("nan")
        )
        mean_mae = float(np.mean(np.abs(y_true[keep] - yhat[keep]))) if n_points else float("nan")

        m_series_vec = np.array(
            [series_losses.get((tid, mtype), float("nan")) for tid in all_ts], dtype=float
        )
        mean_metric = (
            float(np.nanmean(m_series_vec)) if np.isfinite(m_series_vec).any() else float("nan")
        )
        max_h = int(pd.to_numeric(m_df["horizon_step"], errors="coerce").max() or 1)

        (
            base_mean_metric,
            fva_abs,
            fva_pct,
            win_rate,
            dm_stat,
            p_dm,
            p_wilcox,
            sig_05,
        ) = _compare_pair(mtype, base_model, m_series_vec, mean_metric, max_h)

        uni_base, local_base = _match_sibling_baselines(mtype, all_models)
        (
            _,
            cov_fva_abs,
            cov_fva_pct,
            _,
            _,
            cov_p_dm,
            _,
            cov_sig_05,
        ) = _compare_pair(mtype, uni_base, m_series_vec, mean_metric, max_h)
        (
            _,
            pool_fva_abs,
            pool_fva_pct,
            _,
            _,
            pool_p_dm,
            _,
            pool_sig_05,
        ) = _compare_pair(mtype, local_base, m_series_vec, mean_metric, max_h)

        dom_regime = (
            str(m_df["demand_class"].mode().iloc[0])
            if not m_df["demand_class"].dropna().empty
            else "all"
        )
        records.append(
            {
                "fva_stage_rank": stage_rank,
                "fva_stage": stage_label,
                "model_type": mtype,
                "base_model": base_model_name(mtype),
                "demand_class": dom_regime,
                "n_series": n_series,
                "n_points": n_points,
                "metric": metric,
                "mean_metric": mean_metric,
                "pooled_wape": pooled_wape,
                "mean_mae": mean_mae,
                "baseline_model": base_model,
                "baseline_metric": base_mean_metric,
                "fva_abs": fva_abs,
                "fva_pct": fva_pct,
                "win_rate_vs_baseline": win_rate,
                "dm_stat": dm_stat,
                "p_value_dm": p_dm,
                "p_value_wilcoxon": p_wilcox,
                "significant_at_05": sig_05,
                "univariate_baseline_model": uni_base,
                "covariate_fva_abs": cov_fva_abs,
                "covariate_fva_pct": cov_fva_pct,
                "covariate_p_value_dm": cov_p_dm,
                "covariate_significant_at_05": cov_sig_05,
                "local_baseline_model": local_base,
                "pooling_fva_abs": pool_fva_abs,
                "pooling_fva_pct": pool_fva_pct,
                "pooling_p_value_dm": pool_p_dm,
                "pooling_significant_at_05": pool_sig_05,
            }
        )

    out = pd.DataFrame.from_records(records, columns=cols)
    return out.sort_values(
        ["fva_stage_rank", "mean_metric"], ascending=[True, True], ignore_index=True
    )

compute_horizon_decay

compute_horizon_decay(
    oof_df, *, metric="wape", baseline=None
)

Evaluate forecast accuracy, bias, and interval coverage by lead-time step ($h = 1 \dots H$).

Also identifies crossover_step ($h^*$) for each model: the first horizon step where the model's loss exceeds the baseline model's loss (or None if it beats the baseline across the entire horizon).

Source code in src/scale_forecasting/diagnostics.py

def compute_horizon_decay(
    oof_df: pd.DataFrame,
    *,
    metric: str = "wape",
    baseline: str | None = None,
) -> pd.DataFrame:
    """Evaluate forecast accuracy, bias, and interval coverage by lead-time step ($h = 1 \\dots H$).

    Also identifies ``crossover_step`` ($h^*$) for each model: the first horizon step where the
    model's loss exceeds the ``baseline`` model's loss (or ``None`` if it beats the baseline across
    the entire horizon).
    """
    cols = [
        "model_type",
        "horizon_step",
        "n_obs",
        "metric",
        "metric_value",
        "mae",
        "rmse",
        "bias",
        "coverage_80",
        "baseline_model",
        "baseline_metric_value",
        "fva_vs_baseline_pct",
        "crossover_step",
    ]
    if oof_df is None or oof_df.empty:
        return pd.DataFrame(columns=cols)

    df = _ensure_oof_columns(oof_df)
    all_models = list(dict.fromkeys(df["model_type"].astype(str)))
    base_model = _resolve_baseline_model(all_models, baseline)

    step_loss: dict[tuple[str, int], float] = {}
    partial_rows: list[dict[str, Any]] = []
    for (mtype, step_raw), grp in df.groupby(["model_type", "horizon_step"], sort=True):
        step = int(str(step_raw))
        y_true = pd.to_numeric(grp["y_true"], errors="coerce").to_numpy(dtype=float)
        yhat = pd.to_numeric(grp["yhat"], errors="coerce").to_numpy(dtype=float)
        keep = np.isfinite(y_true) & np.isfinite(yhat)
        n_obs = int(keep.sum())
        if n_obs == 0:
            continue
        yt, yh = y_true[keep], yhat[keep]
        m_val = _cell_metric_loss(grp, metric)
        step_loss[(str(mtype), step)] = m_val
        err = yh - yt
        mae = float(np.mean(np.abs(err)))
        rmse = float(np.sqrt(np.mean(err**2)))
        bias = float(np.mean(err))

        cov80 = float("nan")
        if "yhat_lower" in grp.columns and "yhat_upper" in grp.columns:
            lo = pd.to_numeric(grp["yhat_lower"], errors="coerce").to_numpy(dtype=float)[keep]
            hi = pd.to_numeric(grp["yhat_upper"], errors="coerce").to_numpy(dtype=float)[keep]
            b_keep = np.isfinite(lo) & np.isfinite(hi)
            if b_keep.any():
                cov80 = float(np.mean((yt[b_keep] >= lo[b_keep]) & (yt[b_keep] <= hi[b_keep])))

        partial_rows.append(
            {
                "model_type": str(mtype),
                "horizon_step": step,
                "n_obs": n_obs,
                "metric": metric,
                "metric_value": m_val,
                "mae": mae,
                "rmse": rmse,
                "bias": bias,
                "coverage_80": cov80,
                "baseline_model": base_model,
            }
        )

    # Determine crossover step h* per model relative to base_model.
    crossover_by_model: dict[str, int | None] = {}
    for mtype in all_models:
        if mtype == base_model:
            crossover_by_model[mtype] = None
            continue
        steps_for_m = sorted(s for (m, s) in step_loss if m == mtype)
        cross: int | None = None
        for s in steps_for_m:
            ml = step_loss.get((mtype, s), float("nan"))
            bl = step_loss.get((base_model, s), float("nan"))
            if math.isfinite(ml) and math.isfinite(bl) and ml > bl:
                cross = s
                break
        crossover_by_model[mtype] = cross

    for row in partial_rows:
        mtype = row["model_type"]
        step = row["horizon_step"]
        ml = row["metric_value"]
        bl = step_loss.get((base_model, step), float("nan"))
        row["baseline_metric_value"] = bl
        row["fva_vs_baseline_pct"] = (
            float((bl - ml) / bl)
            if math.isfinite(ml) and math.isfinite(bl) and bl > 0
            else float("nan")
        )
        row["crossover_step"] = crossover_by_model.get(mtype)

    out = pd.DataFrame.from_records(partial_rows, columns=cols)
    return out.sort_values(["model_type", "horizon_step"], ignore_index=True)

compute_inventory_policy

compute_inventory_policy(
    predictions_df,
    *,
    oof_df=None,
    inventory_df=None,
    model_name=None,
    model_type=None,
    lead_time=7.0,
    lead_time_std=0.0,
    review_period=0.0,
    service_level=0.95,
    fill_rate=None,
    on_hand=0.0,
    on_order=0.0,
    backorders=0.0,
    holding_cost=1.0,
    stockout_cost=9.0,
    min_order_qty=0.0,
    order_multiple=0.0,
    confidence_level=0.8,
)

Translate calibrated multi-step forecasts into inventory & replenishment policies.

Combines non-stationary forward forecasts ($\hat{y}_1 \dots \hat{y}_H$), empirical AR(1) multi-step forecast error covariance ($\hat{\rho}_1$ estimated from oof_df), and supplier lead-time variability ($\mu_L, \sigma_L$) to compute continuous-review $(s, S)$ and periodic-review $(R, S)$ policies:

$$\sigma_P = \sqrt{\sum_{i=1}^{\lceil P \rceil} \sigma_i^2 + 2\sum_{1 \le i < j \le \lceil P \rceil} \hat{\rho}1^{j-i}\,\sigma_i\,\sigma_j + \sigma_L^2\,\bar{\mu}$$}}^2

Supports three service-level targeting modes (globally or per-SKU via inventory_df): 1. Cycle Service Level (service_level, Type-1 CSL $\alpha \in (0, 1)$): Sets $z = \Phi^{-1}(\alpha)$ so the probability of zero stockout during the protection window $P = \mu_L + R$ equals $\alpha$. 2. Unit Fill Rate (fill_rate, Type-2 Fill Rate $\beta \in (0, 1)$): Solves the standard normal loss equation $G(z) = (1 - \beta)\,\mu_P / \sigma_P$ for $z$ so the expected fraction of demand satisfied immediately from stock equals $\beta$. 3. Newsvendor Critical Fractile (service_level=None and fill_rate=None): Computes $\alpha^* = c_u / (c_u + c_o)$ from stockout_cost ($c_u$) and holding_cost ($c_o$) to minimize total expected holding + stockout cost per cycle.

Also reports naive_iid_safety_stock ($\rho = 0$) and naive_linear_safety_stock ($\rho = 1$) alongside safety_stock so supply-chain planners can audit the buffer impact of autocorrelation-aware interval aggregation.

Source code in src/scale_forecasting/diagnostics.py

def compute_inventory_policy(
    predictions_df: pd.DataFrame,
    *,
    oof_df: pd.DataFrame | None = None,
    inventory_df: pd.DataFrame | Mapping[str, Mapping[str, Any]] | None = None,
    model_name: str | None = None,
    model_type: str | None = None,
    lead_time: float = 7.0,
    lead_time_std: float = 0.0,
    review_period: float = 0.0,
    service_level: float | None = 0.95,
    fill_rate: float | None = None,
    on_hand: float = 0.0,
    on_order: float = 0.0,
    backorders: float = 0.0,
    holding_cost: float = 1.0,
    stockout_cost: float = 9.0,
    min_order_qty: float = 0.0,
    order_multiple: float = 0.0,
    confidence_level: float = 0.80,
) -> pd.DataFrame:
    """Translate calibrated multi-step forecasts into inventory & replenishment policies.

    Combines non-stationary forward forecasts ($\\hat{y}_1 \\dots \\hat{y}_H$), empirical AR(1)
    multi-step forecast error covariance ($\\hat{\\rho}_1$ estimated from ``oof_df``), and supplier
    lead-time variability ($\\mu_L, \\sigma_L$) to compute continuous-review $(s, S)$ and
    periodic-review $(R, S)$ policies:

    $$\\sigma_P = \\sqrt{\\sum_{i=1}^{\\lceil P \\rceil} \\sigma_i^2 +
      2\\sum_{1 \\le i < j \\le \\lceil P \\rceil} \\hat{\\rho}_1^{j-i}\\,\\sigma_i\\,\\sigma_j +
      \\sigma_L^2\\,\\bar{\\mu}_{\\text{step}}^2}$$

    Supports three service-level targeting modes (globally or per-SKU via ``inventory_df``):
    1. **Cycle Service Level (``service_level``, Type-1 CSL $\\alpha \\in (0, 1)$):**
       Sets $z = \\Phi^{-1}(\\alpha)$ so the probability of zero stockout during the protection
       window $P = \\mu_L + R$ equals $\\alpha$.
    2. **Unit Fill Rate (``fill_rate``, Type-2 Fill Rate $\\beta \\in (0, 1)$):**
       Solves the standard normal loss equation $G(z) = (1 - \\beta)\\,\\mu_P / \\sigma_P$ for $z$
       so the expected fraction of demand satisfied immediately from stock equals $\\beta$.
    3. **Newsvendor Critical Fractile (``service_level=None`` and ``fill_rate=None``):**
       Computes $\\alpha^* = c_u / (c_u + c_o)$ from ``stockout_cost`` ($c_u$) and ``holding_cost``
       ($c_o$) to minimize total expected holding + stockout cost per cycle.

    Also reports ``naive_iid_safety_stock`` ($\\rho = 0$) and ``naive_linear_safety_stock``
    ($\\rho = 1$) alongside ``safety_stock`` so supply-chain planners can audit the buffer impact
    of autocorrelation-aware interval aggregation.
    """
    from .errors import ConfigError

    if float(lead_time) <= 0.0:
        raise ConfigError(f"lead_time={lead_time!r} must be > 0.")
    if float(lead_time_std) < 0.0:
        raise ConfigError(f"lead_time_std={lead_time_std!r} must be >= 0.")
    if float(review_period) < 0.0:
        raise ConfigError(f"review_period={review_period!r} must be >= 0.")
    if service_level is not None and not (0.0 < float(service_level) < 1.0):
        raise ConfigError(f"service_level={service_level!r} must lie strictly in (0.0, 1.0).")
    if fill_rate is not None and not (0.0 < float(fill_rate) < 1.0):
        raise ConfigError(f"fill_rate={fill_rate!r} must lie strictly in (0.0, 1.0).")
    if float(holding_cost) < 0.0 or float(stockout_cost) < 0.0:
        raise ConfigError("holding_cost and stockout_cost must be >= 0.")
    if (
        service_level is None
        and fill_rate is None
        and (float(holding_cost) <= 0.0 or float(stockout_cost) <= 0.0)
    ):
        raise ConfigError(
            "When both service_level and fill_rate are None (Newsvendor mode), "
            "both holding_cost and stockout_cost must be > 0."
        )
    if float(min_order_qty) < 0.0 or float(order_multiple) < 0.0:
        raise ConfigError("min_order_qty and order_multiple must be >= 0.")

    if predictions_df is None or predictions_df.empty:
        return pd.DataFrame(columns=_INVENTORY_POLICY_COLUMNS)

    df = predictions_df.copy()
    if "segment" in df.columns and (df["segment"] == "forecast").any():
        df = df[df["segment"] == "forecast"].copy()
    if df.empty:
        return pd.DataFrame(columns=_INVENTORY_POLICY_COLUMNS)

    date_col = (
        "forecast_date"
        if "forecast_date" in df.columns
        else (
            "target_date" if "target_date" in df.columns else ("ds" if "ds" in df.columns else None)
        )
    )
    if date_col is None:
        raise ValueError(
            "predictions_df must contain a 'ds', 'forecast_date', or 'target_date' column"
        )
    yhat_col = "yhat" if "yhat" in df.columns else "forecast"
    lo_col = (
        "yhat_lower"
        if "yhat_lower" in df.columns
        else ("lower_bound" if "lower_bound" in df.columns else None)
    )
    hi_col = (
        "yhat_upper"
        if "yhat_upper" in df.columns
        else ("upper_bound" if "upper_bound" in df.columns else None)
    )

    eff_model = model_name or model_type
    if "ts_id" not in df.columns:
        df["ts_id"] = "series_0"
    if "model_type" not in df.columns:
        df["model_type"] = eff_model or "model"

    if eff_model is not None and str(eff_model).lower() not in ("all", "*"):
        req = str(eff_model).strip()
        exact = df[df["model_type"] == req]
        if exact.empty and req == "arima_plus_xreg":
            exact = df[df["model_type"] == "arima_plus@xreg"]
        df = exact
        if df.empty:
            return pd.DataFrame(columns=_INVENTORY_POLICY_COLUMNS)

    df["_dt"] = pd.to_datetime(df[date_col])
    df = df.sort_values(["ts_id", "model_type", "_dt"])

    rho_map = _estimate_oof_error_rho(oof_df)
    oof_sigma_map: dict[tuple[str, str], float] = {}
    if oof_df is not None and not oof_df.empty:
        oof_norm = _ensure_oof_columns(oof_df)
        for (tid_o, mtype_o), grp_o in oof_norm.groupby(["ts_id", "model_type"], sort=False):
            err_o = (
                (
                    pd.to_numeric(grp_o["y_true"], errors="coerce")
                    - pd.to_numeric(grp_o["yhat"], errors="coerce")
                )
                .dropna()
                .to_numpy(dtype=float)
            )
            err_o = err_o[np.isfinite(err_o)]
            if err_o.size > 0:
                mad_o = 1.4826 * float(np.median(np.abs(err_o)))
                rmse_o = float(np.sqrt(np.mean(err_o**2)))
                oof_sigma_map[(str(tid_o), str(mtype_o))] = max(mad_o, rmse_o, 1e-6)

    eff_conf = float(confidence_level) if 0.0 < float(confidence_level) < 1.0 else 0.80
    z_conf = max(float(stats.norm.ppf(0.5 + 0.5 * eff_conf)), 1e-6)
    overrides = _normalize_inventory_overrides(inventory_df)

    records: list[dict[str, Any]] = []
    for (tid, mtype), grp in df.groupby(["ts_id", "model_type"], sort=True):
        tid_s, mtype_s = str(tid), str(mtype)
        sku = overrides.get(tid_s, {})

        lt = float(sku.get("lead_time", lead_time))
        lt_std = float(sku.get("lead_time_std", lead_time_std))
        rev_p = float(sku.get("review_period", review_period))
        oh = float(sku.get("on_hand", on_hand))
        oo = float(sku.get("on_order", on_order))
        bo = float(sku.get("backorders", backorders))
        h_cost = float(sku.get("holding_cost", holding_cost))
        s_cost = float(sku.get("stockout_cost", stockout_cost))
        moq = float(sku.get("min_order_qty", min_order_qty))
        mult = float(sku.get("order_multiple", order_multiple))
        sku_fr = sku.get("fill_rate", fill_rate)
        sku_sl = sku.get("service_level", service_level)

        if lt <= 0.0 or lt_std < 0.0 or rev_p < 0.0:
            raise ConfigError(f"Invalid lead_time/review_period override for ts_id={tid_s!r}.")

        yh = np.maximum(
            0.0,
            np.nan_to_num(
                pd.to_numeric(grp[yhat_col], errors="coerce").to_numpy(dtype=float),
                nan=0.0,
            ),
        )
        n_avail = int(yh.size)
        if n_avail == 0:
            continue

        # Recover step-level upper-tail standard deviations sigma_h.
        sigma_hi = np.zeros(n_avail, dtype=float)
        if hi_col is not None:
            yu = pd.to_numeric(grp[hi_col], errors="coerce").to_numpy(dtype=float)
            if lo_col is not None:
                yl = pd.to_numeric(grp[lo_col], errors="coerce").to_numpy(dtype=float)
                half_w = np.maximum(yu - yh, 0.5 * (yu - yl))
            else:
                half_w = yu - yh
            valid_w = np.isfinite(half_w) & (half_w > 1e-9)
            sigma_hi = np.where(valid_w, np.maximum(half_w, 0.0) / z_conf, 0.0)

        # If upper bounds were absent, check `quantiles` column if present.
        if not (sigma_hi > 1e-9).any() and "quantiles" in grp.columns:
            q_sigmas: list[float] = []
            for y_val, q_raw in zip(yh, grp["quantiles"], strict=False):
                q_dict: dict[str, Any] = {}
                if isinstance(q_raw, str) and q_raw.startswith("{"):
                    try:
                        q_dict = json.loads(q_raw)
                    except Exception:  # noqa: BLE001 - malformed JSON falls back to empty dict
                        q_dict = {}
                elif isinstance(q_raw, Mapping):
                    q_dict = dict(q_raw)
                hi_q_candidates: list[tuple[float, float]] = []
                for k_q, v_q in q_dict.items():
                    k_clean = str(k_q).lstrip("pqPQ")
                    try:
                        qf = float(k_clean)
                        if qf > 1.0:
                            qf /= 100.0
                        vf = float(v_q)
                        if 0.55 <= qf < 1.0 and math.isfinite(vf):
                            hi_q_candidates.append((qf, vf))
                    except (TypeError, ValueError):
                        continue
                if hi_q_candidates:
                    qf_best, vf_best = max(hi_q_candidates, key=lambda item: item[0])
                    zq = max(float(stats.norm.ppf(qf_best)), 1e-6)
                    q_sigmas.append(max(0.0, vf_best - float(y_val)) / zq)
                else:
                    q_sigmas.append(0.0)
            sigma_hi = np.asarray(q_sigmas, dtype=float)

        fallback_sigma = oof_sigma_map.get((tid_s, mtype_s), float(np.nanmean(yh)) * 0.15)
        sigma_hi = np.where(sigma_hi > 1e-9, sigma_hi, max(fallback_sigma, 1e-6))

        prot_steps = lt + rev_p
        lt_demand = _window_sum(yh, lt)
        prot_demand = _window_sum(yh, prot_steps)
        mean_step = prot_demand / max(prot_steps, 1e-9)

        rho = rho_map.get((tid_s, mtype_s), 0.35)
        prot_sigmas = _protection_sigma_vector(sigma_hi, prot_steps)
        fcst_std = _ar1_sum_scale(prot_sigmas, rho)
        iid_fcst_std = _ar1_sum_scale(prot_sigmas, 0.0)
        lin_fcst_std = float(np.sum(prot_sigmas))

        lt_var_std = lt_std * max(mean_step, 0.0)
        total_std = math.sqrt(fcst_std**2 + lt_var_std**2)
        total_iid_std = math.sqrt(iid_fcst_std**2 + lt_var_std**2)
        total_lin_std = math.sqrt(lin_fcst_std**2 + lt_var_std**2)

        # Resolve safety factor z, target_mode, and effective cycle service level alpha.
        if sku_fr is not None:
            fr_val = float(sku_fr)
            if not (0.0 < fr_val < 1.0):
                raise ConfigError(f"fill_rate={sku_fr!r} must lie strictly in (0.0, 1.0).")
            z_factor = _solve_z_for_fill_rate(fr_val, prot_demand, total_std)
            target_csl = float(stats.norm.cdf(z_factor))
            target_mode = "fill_rate"
        elif sku_sl is not None and str(sku_sl).lower() != "newsvendor":
            sl_val = float(sku_sl)
            if not (0.0 < sl_val < 1.0):
                raise ConfigError(f"service_level={sku_sl!r} must lie strictly in (0.0, 1.0).")
            target_csl = sl_val
            z_factor = max(0.0, float(stats.norm.ppf(target_csl)))
            target_mode = "cycle_service_level"
        else:
            if h_cost <= 0.0 or s_cost <= 0.0:
                raise ConfigError(
                    f"Newsvendor service level for ts_id={tid_s!r} requires "
                    "holding_cost > 0 and stockout_cost > 0."
                )
            target_csl = float(s_cost / (s_cost + h_cost))
            z_factor = max(0.0, float(stats.norm.ppf(target_csl)))
            target_mode = "newsvendor_critical_fractile"

        ss = max(0.0, z_factor * total_std)
        ss_iid = max(0.0, z_factor * total_iid_std)
        ss_lin = max(0.0, z_factor * total_lin_std)

        rop = lt_demand + ss
        out_level = prot_demand + ss
        inv_pos = oh + oo - bo
        pol_type = "periodic_R_S" if rev_p > 0.0 else "continuous_s_S"
        triggered = bool(inv_pos < out_level) if rev_p > 0.0 else bool(inv_pos <= rop)

        raw_q = max(0.0, out_level - inv_pos) if triggered else 0.0
        rec_q = raw_q
        if rec_q > 0.0:
            if moq > 0.0:
                rec_q = max(rec_q, moq)
            if mult > 0.0:
                rec_q = math.ceil(rec_q / mult - 1e-12) * mult

        exp_short = total_std * _normal_loss(z_factor)
        imp_fill = float(np.clip(1.0 - exp_short / max(prot_demand, 1e-9), 0.0, 1.0))
        exp_overage = ss + exp_short
        exp_hold_cost = h_cost * exp_overage
        exp_stock_cost = s_cost * exp_short
        exp_tot_cost = exp_hold_cost + exp_stock_cost

        records.append(
            {
                "ts_id": tid_s,
                "model_type": mtype_s,
                "policy_type": pol_type,
                "target_mode": target_mode,
                "lead_time": lt,
                "lead_time_std": lt_std,
                "review_period": rev_p,
                "protection_steps": prot_steps,
                "mean_step_demand": mean_step,
                "lead_time_demand": lt_demand,
                "protection_demand": prot_demand,
                "error_autocorr_rho": rho,
                "forecast_error_std": fcst_std,
                "lead_time_var_std": lt_var_std,
                "total_protection_std": total_std,
                "target_service_level": target_csl,
                "implied_fill_rate": imp_fill,
                "safety_factor_z": z_factor,
                "safety_stock": ss,
                "naive_iid_safety_stock": ss_iid,
                "naive_linear_safety_stock": ss_lin,
                "reorder_point": rop,
                "order_up_to_level": out_level,
                "inventory_position": inv_pos,
                "reorder_triggered": triggered,
                "recommended_order_qty": rec_q,
                "expected_units_short": exp_short,
                "expected_holding_cost": exp_hold_cost,
                "expected_stockout_cost": exp_stock_cost,
                "expected_total_cost": exp_tot_cost,
            }
        )

    if not records:
        return pd.DataFrame(columns=_INVENTORY_POLICY_COLUMNS)
    return pd.DataFrame.from_records(records, columns=_INVENTORY_POLICY_COLUMNS)

compute_realized_accuracy

compute_realized_accuracy(
    predictions_df,
    actuals_df,
    *,
    oof_df=None,
    metric="wape",
    drift_threshold=1.5,
)

Join published forward forecasts with newly arrived actuals on (ts_id, date) and score realized accuracy, interval coverage, and accuracy drift relative to backtest expectations.

Parameters:

Name	Type	Description	Default
predictions_df	DataFrame	Forward predictions frame containing ts_id, model_type, forecast_date (or ds), yhat (or forecast), and optional yhat_lower / yhat_upper (or lower_bound / upper_bound).	required
actuals_df	DataFrame	Newly arrived actuals containing ts_id, ds (or forecast_date), and y (or y_true / actual).	required
oof_df	DataFrame | None	Optional backtest OOF frame from the same run to compare realized error against expected backtest error (drift_ratio = realized_metric / backtest_metric).	None
metric	str	Primary accuracy metric (default "wape").	'wape'
drift_threshold	float	Ratio threshold above which drift_flagged is set True (default 1.5).	1.5

Source code in src/scale_forecasting/diagnostics.py

def compute_realized_accuracy(
    predictions_df: pd.DataFrame,
    actuals_df: pd.DataFrame,
    *,
    oof_df: pd.DataFrame | None = None,
    metric: str = "wape",
    drift_threshold: float = 1.5,
) -> pd.DataFrame:
    """Join published forward forecasts with newly arrived actuals on ``(ts_id, date)`` and score
    realized accuracy, interval coverage, and accuracy drift relative to backtest expectations.

    Args:
        predictions_df: Forward predictions frame containing ``ts_id``, ``model_type``,
            ``forecast_date`` (or ``ds``), ``yhat`` (or ``forecast``), and optional ``yhat_lower``
            / ``yhat_upper`` (or ``lower_bound`` / ``upper_bound``).
        actuals_df: Newly arrived actuals containing ``ts_id``, ``ds`` (or ``forecast_date``), and
            ``y`` (or ``y_true`` / ``actual``).
        oof_df: Optional backtest OOF frame from the same run to compare realized error against
            expected backtest error (`drift_ratio = realized_metric / backtest_metric`).
        metric: Primary accuracy metric (default ``"wape"``).
        drift_threshold: Ratio threshold above which ``drift_flagged`` is set True (default 1.5).
    """
    cols = [
        "ts_id",
        "model_type",
        "n_realized",
        "metric",
        "realized_metric",
        "realized_mae",
        "realized_rmse",
        "realized_bias",
        "realized_coverage_80",
        "backtest_metric",
        "drift_ratio",
        "drift_flagged",
    ]
    if predictions_df is None or predictions_df.empty or actuals_df is None or actuals_df.empty:
        return pd.DataFrame(columns=cols)

    pred = predictions_df.copy()
    if "segment" in pred.columns and (pred["segment"] == "forecast").any():
        pred = pred[pred["segment"] == "forecast"].copy()
    p_date = "forecast_date" if "forecast_date" in pred.columns else "ds"
    p_yhat = "yhat" if "yhat" in pred.columns else "forecast"
    p_lo = (
        "yhat_lower"
        if "yhat_lower" in pred.columns
        else ("lower_bound" if "lower_bound" in pred.columns else None)
    )
    p_hi = (
        "yhat_upper"
        if "yhat_upper" in pred.columns
        else ("upper_bound" if "upper_bound" in pred.columns else None)
    )
    if "ts_id" not in pred.columns:
        pred["ts_id"] = "series_0"
    if "model_type" not in pred.columns:
        pred["model_type"] = "model"

    act = actuals_df.copy()
    a_date = "ds" if "ds" in act.columns else "forecast_date"
    a_y = "y" if "y" in act.columns else ("y_true" if "y_true" in act.columns else "actual")
    if "ts_id" not in act.columns:
        act["ts_id"] = "series_0"

    pred["_dt"] = pd.to_datetime(pred[p_date])
    act["_dt"] = pd.to_datetime(act[a_date])
    act_dedup = (
        act[["ts_id", "_dt", a_y]]
        .drop_duplicates(subset=["ts_id", "_dt"], keep="last")
        .rename(columns={a_y: "_y_actual"})
    )
    merged = pred.merge(act_dedup, on=["ts_id", "_dt"], how="inner")
    if merged.empty:
        return pd.DataFrame(columns=cols)

    bt_loss_map: dict[tuple[str, str], float] = {}
    if oof_df is not None and not oof_df.empty:
        oof_norm = _ensure_oof_columns(oof_df)
        for (tid, mtype), grp in oof_norm.groupby(["ts_id", "model_type"], sort=False):
            bt_loss_map[(str(tid), str(mtype))] = _cell_metric_loss(grp, metric)

    records: list[dict[str, Any]] = []
    for (tid, mtype), grp in merged.groupby(["ts_id", "model_type"], sort=True):
        yt = pd.to_numeric(grp["_y_actual"], errors="coerce").to_numpy(dtype=float)
        yh = pd.to_numeric(grp[p_yhat], errors="coerce").to_numpy(dtype=float)
        keep = np.isfinite(yt) & np.isfinite(yh)
        n_real = int(keep.sum())
        if n_real == 0:
            continue
        yt_k, yh_k = yt[keep], yh[keep]
        lo_k = (
            pd.to_numeric(grp[p_lo], errors="coerce").to_numpy(dtype=float)[keep]
            if p_lo is not None
            else None
        )
        hi_k = (
            pd.to_numeric(grp[p_hi], errors="coerce").to_numpy(dtype=float)[keep]
            if p_hi is not None
            else None
        )
        m_panel = compute_metrics(
            yt_k,
            yh_k,
            lower=lo_k if lo_k is not None and np.isfinite(lo_k).any() else None,
            upper=hi_k if hi_k is not None and np.isfinite(hi_k).any() else None,
        )
        real_loss = float(loss_of(metric, m_panel.get(metric, float("nan"))))
        bt_loss = bt_loss_map.get((str(tid), str(mtype)), float("nan"))
        drift = (
            float(real_loss / bt_loss)
            if math.isfinite(real_loss) and math.isfinite(bt_loss) and bt_loss > 1e-12
            else float("nan")
        )
        flagged = bool(math.isfinite(drift) and drift >= float(drift_threshold))
        records.append(
            {
                "ts_id": str(tid),
                "model_type": str(mtype),
                "n_realized": n_real,
                "metric": metric,
                "realized_metric": real_loss,
                "realized_mae": float(m_panel.get("mae", float("nan"))),
                "realized_rmse": float(m_panel.get("rmse", float("nan"))),
                "realized_bias": float(m_panel.get("bias", float("nan"))),
                "realized_coverage_80": float(m_panel.get("coverage", float("nan"))),
                "backtest_metric": bt_loss,
                "drift_ratio": drift,
                "drift_flagged": flagged,
            }
        )

    out = pd.DataFrame.from_records(records, columns=cols)
    return out.sort_values(
        ["drift_flagged", "realized_metric"], ascending=[False, False], ignore_index=True
    )

compute_revision_stability

compute_revision_stability(
    predictions_by_run, *, jump_threshold_pct=0.3
)

Measure cross-run forecast revision stability ("whipsaw") across overlapping target dates.

Accepts either a mapping {run_id: predictions_df} (ordered oldest-to-newest) or a single concatenated DataFrame containing a run_id column. Computes per (ts_id, model_type):

• n_runs / n_overlapping_dates
• mean_forecast: average forecast magnitude across overlapping dates
• mar: Mean Absolute Revision $|\hat{y}^{(r)}_t - \hat{y}^{(r-1)}_t|$
• revision_cov: Coefficient of Variation of revisions ($\text{std}(\Delta \hat{y}) / |\overline{\hat{y}}|$)
• max_single_step_jump_pct: largest single-date percentage jump between consecutive runs
• whipsaw_flagged: True when max_single_step_jump_pct >= jump_threshold_pct

Source code in src/scale_forecasting/diagnostics.py

def compute_revision_stability(
    predictions_by_run: Mapping[str, pd.DataFrame] | pd.DataFrame,
    *,
    jump_threshold_pct: float = 0.30,
) -> pd.DataFrame:
    """Measure cross-run forecast revision stability ("whipsaw") across overlapping target dates.

    Accepts either a mapping ``{run_id: predictions_df}`` (ordered oldest-to-newest) or a single
    concatenated DataFrame containing a ``run_id`` column. Computes per ``(ts_id, model_type)``:

    - ``n_runs`` / ``n_overlapping_dates``
    - ``mean_forecast``: average forecast magnitude across overlapping dates
    - ``mar``: Mean Absolute Revision $|\\hat{y}^{(r)}_t - \\hat{y}^{(r-1)}_t|$
    - ``revision_cov``: Coefficient of Variation of revisions ($\\text{std}(\\Delta \\hat{y}) /
      |\\overline{\\hat{y}}|$)
    - ``max_single_step_jump_pct``: largest single-date percentage jump between consecutive runs
    - ``whipsaw_flagged``: True when ``max_single_step_jump_pct >= jump_threshold_pct``
    """
    cols = [
        "ts_id",
        "model_type",
        "n_runs",
        "n_overlapping_dates",
        "mean_forecast",
        "mar",
        "revision_cov",
        "max_single_step_jump_pct",
        "whipsaw_flagged",
    ]
    frames: list[pd.DataFrame] = []
    run_order: list[str] = []
    if isinstance(predictions_by_run, Mapping):
        for idx, (rid, f_df) in enumerate(predictions_by_run.items()):
            if f_df is None or f_df.empty:
                continue
            sub = f_df.copy()
            sub["_run_order"] = idx
            sub["run_id"] = str(rid)
            frames.append(sub)
            run_order.append(str(rid))
    elif isinstance(predictions_by_run, pd.DataFrame) and not predictions_by_run.empty:
        raw_df = predictions_by_run.copy()
        if "run_id" not in raw_df.columns:
            return pd.DataFrame(columns=cols)
        for idx, rid in enumerate(dict.fromkeys(raw_df["run_id"].astype(str))):
            sub = raw_df[raw_df["run_id"].astype(str) == rid].copy()
            sub["_run_order"] = idx
            frames.append(sub)
            run_order.append(rid)

    if len(frames) < 2:
        return pd.DataFrame(columns=cols)

    combined = pd.concat(frames, ignore_index=True)
    if "segment" in combined.columns and (combined["segment"] == "forecast").any():
        combined = combined[combined["segment"] == "forecast"].copy()
    date_col = "forecast_date" if "forecast_date" in combined.columns else "ds"
    yhat_col = "yhat" if "yhat" in combined.columns else "forecast"
    if "ts_id" not in combined.columns:
        combined["ts_id"] = "series_0"
    if "model_type" not in combined.columns:
        combined["model_type"] = "model"

    combined["_dt"] = pd.to_datetime(combined[date_col])
    combined["_yh"] = pd.to_numeric(combined[yhat_col], errors="coerce")

    records: list[dict[str, Any]] = []
    for (tid, mtype), grp in combined.groupby(["ts_id", "model_type"], sort=True):
        piv = grp.pivot_table(index="_dt", columns="_run_order", values="_yh", aggfunc="mean")
        if piv.shape[1] < 2:
            continue
        diffs: list[float] = []
        pct_jumps: list[float] = []
        overlap_dates: set[Any] = set()
        run_cols = sorted(piv.columns)
        for c_prev, c_curr in itertools.pairwise(run_cols):
            pair = piv[[c_prev, c_curr]].dropna()
            if pair.empty:
                continue
            overlap_dates.update(pair.index.tolist())
            prev_v = pair[c_prev].to_numpy(dtype=float)
            curr_v = pair[c_curr].to_numpy(dtype=float)
            delta = curr_v - prev_v
            diffs.extend(np.abs(delta).tolist())
            denom = np.maximum(np.abs(prev_v), 1e-6)
            pct_jumps.extend((np.abs(delta) / denom).tolist())

        if not diffs:
            continue
        mean_fc = float(np.nanmean(piv.to_numpy(dtype=float)))
        mar = float(np.mean(diffs))
        rev_std = float(np.std(diffs))
        rev_cov = rev_std / max(abs(mean_fc), 1e-6)
        max_jump = float(np.max(pct_jumps))
        flagged = bool(max_jump >= float(jump_threshold_pct))
        records.append(
            {
                "ts_id": str(tid),
                "model_type": str(mtype),
                "n_runs": int(piv.shape[1]),
                "n_overlapping_dates": len(overlap_dates),
                "mean_forecast": mean_fc,
                "mar": mar,
                "revision_cov": rev_cov,
                "max_single_step_jump_pct": max_jump,
                "whipsaw_flagged": flagged,
            }
        )

    out = pd.DataFrame.from_records(records, columns=cols)
    return out.sort_values(
        ["whipsaw_flagged", "max_single_step_jump_pct"],
        ascending=[False, False],
        ignore_index=True,
    )

diebold_mariano_test

diebold_mariano_test(loss_model, loss_baseline, *, h=1)

Harvey, Leybourne & Newbold (1997) modified Diebold–Mariano test on paired losses.

Tests $H_0: \mathbb{E}[L(\text{model}) - L(\text{baseline})] = 0$ against the two-sided alternative using Newey–West / Bartlett-style autocovariance up to lag $h - 1$ and the HLN small-sample Student-$t_{n-1}$ adjustment:

$$\kappa_{\text{HLN}} = \sqrt{\frac{n + 1 - 2h + h(h - 1)/n}{n}}$$

Returns (dm_stat, p_value). Negative dm_stat means loss_model < loss_baseline (the candidate model improved on the baseline).

Source code in src/scale_forecasting/diagnostics.py

def diebold_mariano_test(
    loss_model: Sequence[float] | np.ndarray,
    loss_baseline: Sequence[float] | np.ndarray,
    *,
    h: int = 1,
) -> tuple[float, float]:
    """Harvey, Leybourne & Newbold (1997) modified Diebold–Mariano test on paired losses.

    Tests $H_0: \\mathbb{E}[L(\\text{model}) - L(\\text{baseline})] = 0$ against the two-sided
    alternative using Newey–West / Bartlett-style autocovariance up to lag $h - 1$ and the HLN
    small-sample Student-$t_{n-1}$ adjustment:

    $$\\kappa_{\\text{HLN}} = \\sqrt{\\frac{n + 1 - 2h + h(h - 1)/n}{n}}$$

    Returns ``(dm_stat, p_value)``. Negative ``dm_stat`` means ``loss_model < loss_baseline``
    (the candidate model improved on the baseline).
    """
    lm = np.asarray(loss_model, dtype=float)
    lb = np.asarray(loss_baseline, dtype=float)
    mask = np.isfinite(lm) & np.isfinite(lb)
    if int(mask.sum()) < 3:
        return float("nan"), float("nan")

    d = lm[mask] - lb[mask]
    n = int(d.size)
    d_bar = float(np.mean(d))
    if np.allclose(d, 0.0):
        return 0.0, 1.0

    centered = d - d_bar
    gamma_0 = float(np.dot(centered, centered) / n)
    if gamma_0 <= 1e-18:
        return (float("-inf") if d_bar < 0 else float("inf")), 0.0

    eff_h = max(1, min(int(h), n // 2))
    var_d = gamma_0
    for lag in range(1, eff_h):
        weight = 1.0 - lag / eff_h
        gamma_k = float(np.dot(centered[lag:], centered[:-lag]) / n)
        var_d += 2.0 * weight * gamma_k

    # Guard against negative finite-sample spectral density estimates.
    var_d_bar = max(var_d, gamma_0 * 0.1) / n
    raw_dm = d_bar / math.sqrt(var_d_bar)

    hln_num = n + 1.0 - 2.0 * eff_h + (eff_h * (eff_h - 1.0)) / n
    hln_factor = math.sqrt(max(hln_num / n, 1.0 / n))
    dm_stat = float(raw_dm * hln_factor)
    p_val = float(2.0 * stats.t.sf(abs(dm_stat), df=max(1, n - 1)))
    return dm_stat, p_val

plot_fva

plot_fva(fva_df, *, ax=None, title=None)

Plot the Forecast Value Added (FVA) accuracy waterfall by model and modeling tier.

Source code in src/scale_forecasting/diagnostics.py

def plot_fva(
    fva_df: pd.DataFrame,
    *,
    ax: Any = None,
    title: str | None = None,
) -> Any:
    """Plot the Forecast Value Added (FVA) accuracy waterfall by model and modeling tier."""
    from .errors import require_extra

    require_extra("notebook", purpose="Plotting FVA waterfall")
    import matplotlib.pyplot as plt

    if ax is None:
        _, ax = plt.subplots(figsize=(10, 5))
    if fva_df is None or fva_df.empty:
        ax.set_title(title or "Forecast Value Added (no backtest OOF rows)")
        return ax

    ordered = fva_df.sort_values(["fva_stage_rank", "mean_metric"], ascending=[True, False])
    labels = [
        f"{r.model_type} ({r.fva_stage})" + (" *" if r.significant_at_05 else "")
        for r in ordered.itertuples(index=False)
    ]
    vals = ordered["fva_pct"].fillna(0.0).to_numpy(dtype=float) * 100.0
    colors = ["#2E7D32" if v > 0 else ("#757575" if v == 0 else "#C62828") for v in vals]

    y_pos = np.arange(len(labels))
    ax.barh(y_pos, vals, color=colors, alpha=0.88)
    ax.axvline(0.0, color="#333333", linewidth=1.0, linestyle="--")
    ax.set_yticks(y_pos)
    ax.set_yticklabels(labels)
    ax.set_xlabel(f"FVA % Improvement vs. {ordered['baseline_model'].iloc[0]} (* p < 0.05)")
    ax.set_title(title or "Forecast Value Added (FVA) Ladder")
    ax.grid(axis="x", alpha=0.25)
    return ax

plot_horizon_decay

plot_horizon_decay(decay_df, *, ax=None, title=None)

Plot forecast error decay across lead-time steps $h = 1 \dots H$ per model.

Source code in src/scale_forecasting/diagnostics.py

def plot_horizon_decay(
    decay_df: pd.DataFrame,
    *,
    ax: Any = None,
    title: str | None = None,
) -> Any:
    """Plot forecast error decay across lead-time steps $h = 1 \\dots H$ per model."""
    from .errors import require_extra

    require_extra("notebook", purpose="Plotting horizon decay")
    import matplotlib.pyplot as plt

    if ax is None:
        _, ax = plt.subplots(figsize=(10, 5))
    if decay_df is None or decay_df.empty:
        ax.set_title(title or "Horizon Error Decay (no backtest OOF rows)")
        return ax

    metric_name = str(decay_df["metric"].iloc[0]) if "metric" in decay_df.columns else "loss"
    for mtype, grp in decay_df.groupby("model_type", sort=False):
        ordered = grp.sort_values("horizon_step")
        ax.plot(
            ordered["horizon_step"],
            ordered["metric_value"],
            marker="o",
            markersize=4,
            label=str(mtype),
        )
    ax.set_xlabel("Forecast Horizon Step (h)")
    ax.set_ylabel(metric_name.upper())
    ax.set_title(title or f"Lead-Time Horizon Accuracy Decay ({metric_name.upper()})")
    ax.legend(loc="best", frameon=True)
    ax.grid(alpha=0.25)
    return ax

plot_inventory_policy

plot_inventory_policy(
    policy_df, *, max_rows=15, ax=None, title=None
)

Plot Reorder Points (ROP = Lead-Time Demand + Safety Stock), Order-Up-To Levels ($S$), and Current Inventory Position per (ts_id, model_type).

Source code in src/scale_forecasting/diagnostics.py

def plot_inventory_policy(
    policy_df: pd.DataFrame,
    *,
    max_rows: int = 15,
    ax: Any = None,
    title: str | None = None,
) -> Any:
    """Plot Reorder Points (`ROP` = Lead-Time Demand + Safety Stock), Order-Up-To Levels ($S$),
    and Current Inventory Position per ``(ts_id, model_type)``.
    """
    from .errors import require_extra

    require_extra("notebook", purpose="Plotting inventory policy")
    import matplotlib.pyplot as plt

    if ax is None:
        _, ax = plt.subplots(figsize=(10, 5))
    if policy_df is None or policy_df.empty:
        ax.set_title(title or "Inventory & Replenishment Policy (no rows)")
        return ax

    sub = policy_df.head(max(1, int(max_rows))).copy()
    labels = [f"{r.ts_id} ({r.model_type})" for r in sub.itertuples(index=False)]
    lt_dem = sub["lead_time_demand"].to_numpy(dtype=float)
    ss = sub["safety_stock"].to_numpy(dtype=float)
    rev_extra = np.maximum(0.0, sub["protection_demand"].to_numpy(dtype=float) - lt_dem)
    inv_pos = sub["inventory_position"].to_numpy(dtype=float)
    rec_q = sub["recommended_order_qty"].to_numpy(dtype=float)

    y_pos = np.arange(len(labels))
    ax.barh(y_pos, lt_dem, color="#2563eb", alpha=0.85, label="Lead-Time Demand (μ_L)")
    ax.barh(
        y_pos,
        ss,
        left=lt_dem,
        color="#f59e0b",
        alpha=0.88,
        label="Calibrated Safety Stock (SS)",
    )
    if (rev_extra > 1e-6).any():
        ax.barh(
            y_pos,
            rev_extra,
            left=lt_dem + ss,
            color="#9333ea",
            alpha=0.75,
            label="Review-Period Demand (R)",
        )
    ax.scatter(
        inv_pos,
        y_pos,
        color="#16a34a",
        marker="D",
        s=48,
        zorder=5,
        label="Inventory Position (IP)",
    )
    for idx, q_val in enumerate(rec_q):
        if q_val > 0.0:
            total_bar = lt_dem[idx] + ss[idx] + rev_extra[idx]
            ax.text(
                total_bar * 1.01,
                y_pos[idx],
                f" Order Q={q_val:,.0f}",
                va="center",
                fontsize=8.5,
                color="#b91c1c",
                fontweight="bold",
            )

    ax.set_yticks(y_pos)
    ax.set_yticklabels(labels)
    ax.invert_yaxis()
    ax.set_xlabel("Units")
    ax.set_title(title or "Operational Inventory Policy: ROP, Safety Stock & Order Quantities")
    ax.legend(loc="best", frameon=True)
    ax.grid(axis="x", alpha=0.25)
    return ax

rollup_time

rollup_time(
    predictions_df, *, freq="W", oof_df=None, level=0.8
)

Aggregate forecasts (and optional actuals) to coarser calendar buckets or lead-time windows with autocorrelation-aware prediction intervals.

Why this exists: daily point forecasts (yhat) sum linearly into Weekly ("W") or Monthly ("MS" / "M") totals, but daily yhat_lower and yhat_upper cannot be summed linearly (summing daily bounds assumes perfect $\rho = 1.0$ error correlation across days, inflating a 7-day interval by up to $2.6\times$). rollup_time estimates the lag-1 residual autocorrelation $\hat{\rho}1$ from oof_df (defaulting to a conservative $\rho = 0.35$ when oof_df is omitted) and computes period-level uncertainty via the AR(1) sum variance $\text{Var}(\sum \sigma_i\sigma_j$, reporting both the calibrated bounds (}^k e_i) = \sum_i \sigma_i^2 + 2\sum_{i<j} \hat{\rho}_1^{j-iyhat_lower, yhat_upper) and the naive linear sum (naive_sum_lower, naive_sum_upper) for side-by-side comparison.

Source code in src/scale_forecasting/diagnostics.py

def rollup_time(
    predictions_df: pd.DataFrame,
    *,
    freq: str = "W",
    oof_df: pd.DataFrame | None = None,
    level: float = 0.80,
) -> pd.DataFrame:
    """Aggregate forecasts (and optional actuals) to coarser calendar buckets or lead-time windows
    with autocorrelation-aware prediction intervals.

    Why this exists: daily point forecasts (`yhat`) sum linearly into Weekly (`"W"`) or Monthly
    (`"MS"` / `"M"`) totals, **but daily `yhat_lower` and `yhat_upper` cannot be summed linearly**
    (summing daily bounds assumes perfect $\\rho = 1.0$ error correlation across days, inflating a
    7-day interval by up to $2.6\\times$). `rollup_time` estimates the lag-1 residual
    autocorrelation $\\hat{\\rho}_1$ from ``oof_df`` (defaulting to a conservative $\\rho = 0.35$
    when ``oof_df`` is omitted) and computes period-level uncertainty via the AR(1) sum variance
    $\\text{Var}(\\sum_{i=1}^k e_i) = \\sum_i \\sigma_i^2 + 2\\sum_{i<j} \\hat{\\rho}_1^{j-i}
    \\sigma_i\\sigma_j$, reporting both the calibrated bounds (``yhat_lower``, ``yhat_upper``) and
    the naive linear sum (``naive_sum_lower``, ``naive_sum_upper``) for side-by-side comparison.
    """
    cols = [
        "ts_id",
        "model_type",
        "period_start",
        "period_end",
        "n_steps",
        "y_true",
        "yhat",
        "yhat_lower",
        "yhat_upper",
        "naive_sum_lower",
        "naive_sum_upper",
        "error_autocorr_rho",
    ]
    if predictions_df is None or predictions_df.empty:
        return pd.DataFrame(columns=cols)

    df = predictions_df.copy()
    if "segment" in df.columns and (df["segment"] == "forecast").any():
        df = df[df["segment"] == "forecast"].copy()
    date_col = (
        "ds" if "ds" in df.columns else ("forecast_date" if "forecast_date" in df.columns else None)
    )
    if date_col is None:
        raise ValueError("predictions_df must contain a 'ds' or 'forecast_date' column")
    if "ts_id" not in df.columns:
        df["ts_id"] = "series_0"
    if "model_type" not in df.columns:
        df["model_type"] = "model"

    df["_dt"] = pd.to_datetime(df[date_col])
    df = df.sort_values(["ts_id", "model_type", "_dt"])

    norm_freq = {"M": "ME", "Q": "QE"}.get(freq, freq)
    df["_period"] = df["_dt"].dt.to_period(
        norm_freq if norm_freq not in ("MS", "QS", "ME", "QE") else norm_freq[0]
    )

    rho_map = _estimate_oof_error_rho(oof_df)
    alpha_tail = 0.5 + 0.5 * float(level)
    z_val = float(stats.norm.ppf(alpha_tail))
    z_80 = float(stats.norm.ppf(0.90))

    records: list[dict[str, Any]] = []
    for (tid, mtype, _per), grp in df.groupby(["ts_id", "model_type", "_period"], sort=True):
        yh = pd.to_numeric(grp["yhat"], errors="coerce").to_numpy(dtype=float)
        yt = (
            pd.to_numeric(grp["y_true"], errors="coerce").to_numpy(dtype=float)
            if "y_true" in grp.columns
            else np.full(len(grp), np.nan)
        )
        yl = (
            pd.to_numeric(grp["yhat_lower"], errors="coerce").to_numpy(dtype=float)
            if "yhat_lower" in grp.columns
            else yh.copy()
        )
        yu = (
            pd.to_numeric(grp["yhat_upper"], errors="coerce").to_numpy(dtype=float)
            if "yhat_upper" in grp.columns
            else yh.copy()
        )

        sum_yh = float(np.nansum(yh))
        sum_yt = float(np.nansum(yt)) if np.isfinite(yt).any() else float("nan")
        naive_lo = float(np.nansum(yl))
        naive_hi = float(np.nansum(yu))

        # Step-level sigma_i from (yhat_upper - yhat_lower) / (2 * z_80), preserving asymmetry.
        half_lo = np.maximum(yh - yl, 0.0)
        half_hi = np.maximum(yu - yh, 0.0)
        sigma_lo = np.where(np.isfinite(half_lo), half_lo / z_80, 0.0)
        sigma_hi = np.where(np.isfinite(half_hi), half_hi / z_80, 0.0)

        rho = rho_map.get((str(tid), str(mtype)), 0.35)
        scale_lo = _ar1_sum_scale(sigma_lo, rho)
        scale_hi = _ar1_sum_scale(sigma_hi, rho)

        records.append(
            {
                "ts_id": str(tid),
                "model_type": str(mtype),
                "period_start": grp["_dt"].min(),
                "period_end": grp["_dt"].max(),
                "n_steps": len(grp),
                "y_true": sum_yt,
                "yhat": sum_yh,
                "yhat_lower": float(sum_yh - z_val * scale_lo),
                "yhat_upper": float(sum_yh + z_val * scale_hi),
                "naive_sum_lower": naive_lo,
                "naive_sum_upper": naive_hi,
                "error_autocorr_rho": rho,
            }
        )

    return pd.DataFrame.from_records(records, columns=cols)

wilcoxon_pvalue

wilcoxon_pvalue(loss_model, loss_baseline)

Two-sided paired Wilcoxon signed-rank $p$-value on finite non-tied loss differentials.

Source code in src/scale_forecasting/diagnostics.py

def wilcoxon_pvalue(
    loss_model: Sequence[float] | np.ndarray,
    loss_baseline: Sequence[float] | np.ndarray,
) -> float:
    """Two-sided paired Wilcoxon signed-rank $p$-value on finite non-tied loss differentials."""
    lm = np.asarray(loss_model, dtype=float)
    lb = np.asarray(loss_baseline, dtype=float)
    mask = np.isfinite(lm) & np.isfinite(lb)
    if int(mask.sum()) < 3:
        return float("nan")
    diff = lm[mask] - lb[mask]
    non_zero = diff[np.abs(diff) > 1e-12]
    if non_zero.size < 3:
        return 1.0 if diff.size >= 3 else float("nan")
    try:
        res = stats.wilcoxon(non_zero, zero_method="wilcox", alternative="two-sided")
        return float(res.pvalue)
    except Exception:  # noqa: BLE001 - degenerate rank distributions return NaN gracefully
        return float("nan")

scenarios

Counterfactual "What-If" scenario & price-elasticity simulator — pure, zero GCP imports.

Allows planners and AI agents to perturb future exogenous covariates (promotions, pricing, weather, macro indices, capacity caps) across the entire forecast horizon or within specific date windows (e.g., a 2-week promotion followed by post-promo dip) and compare baseline vs. counterfactual unit volume, prediction intervals, and dollar revenue (price * yhat).

apply_scenario_shocks

apply_scenario_shocks(df, shocks, *, date_col='ds')

Apply counterfactual covariate shocks to a copy of df (pure).

Supports full-horizon or date-windowed perturbations:

• Multiplicative: {"price_index": "*0.85"}, {"price_index": "-15%"}, or {"price_index": {"op": "mul", "value": 0.85, "start_date": "2026-03-01", "end_date": "2026-03-07"}}
• Additive: {"temperature": "+5.0"} or {"temperature": {"op": "add", "value": 5.0}}
• Override: {"promo": "=1"}, {"promo": 1}, or {"promo": {"op": "set", "value": 1.0}}

Source code in src/scale_forecasting/scenarios.py

def apply_scenario_shocks(
    df: pd.DataFrame,
    shocks: Mapping[str, Any],
    *,
    date_col: str = "ds",
) -> pd.DataFrame:
    """Apply counterfactual covariate shocks to a copy of ``df`` (pure).

    Supports full-horizon or date-windowed perturbations:

    - Multiplicative: ``{"price_index": "*0.85"}``, ``{"price_index": "-15%"}``, or
      ``{"price_index": {"op": "mul", "value": 0.85, "start_date": "2026-03-01",
      "end_date": "2026-03-07"}}``
    - Additive: ``{"temperature": "+5.0"}`` or ``{"temperature": {"op": "add", "value": 5.0}}``
    - Override: ``{"promo": "=1"}``, ``{"promo": 1}``, or ``{"promo": {"op": "set", "value": 1.0}}``
    """
    if df is None or df.empty or not shocks:
        return df.copy() if df is not None else pd.DataFrame()

    out = df.copy()
    eff_date_col = (
        date_col
        if date_col in out.columns
        else (
            "ds"
            if "ds" in out.columns
            else ("forecast_date" if "forecast_date" in out.columns else None)
        )
    )
    dt_series = pd.to_datetime(out[eff_date_col]) if eff_date_col is not None else None

    for col, raw_spec in shocks.items():
        if col not in out.columns:
            out[col] = 0.0
        col_vals = pd.to_numeric(out[col], errors="coerce").fillna(0.0).to_numpy(dtype=float)
        for rule in _parse_single_shock(raw_spec):
            mask = np.ones(len(out), dtype=bool)
            if dt_series is not None:
                if rule["start_date"] is not None:
                    mask &= (dt_series >= pd.Timestamp(rule["start_date"])).to_numpy(dtype=bool)
                if rule["end_date"] is not None:
                    mask &= (dt_series <= pd.Timestamp(rule["end_date"])).to_numpy(dtype=bool)
            op = rule["op"]
            val = float(rule["value"])
            if op == "mul":
                col_vals[mask] = col_vals[mask] * val
            elif op == "add":
                col_vals[mask] = col_vals[mask] + val
            elif op == "set":
                col_vals[mask] = val
        out[col] = col_vals
    return out

plot_scenarios

plot_scenarios(
    sim_df,
    *,
    ts_id=None,
    model_type=None,
    show_revenue=False,
    ax=None,
    title=None,
)

Plot baseline vs. counterfactual scenario forecast trajectories (or revenue).

Source code in src/scale_forecasting/scenarios.py

def plot_scenarios(
    sim_df: pd.DataFrame,
    *,
    ts_id: str | None = None,
    model_type: str | None = None,
    show_revenue: bool = False,
    ax: Any = None,
    title: str | None = None,
) -> Any:
    """Plot baseline vs. counterfactual scenario forecast trajectories (or revenue)."""
    from .errors import require_extra

    require_extra("notebook", purpose="Plotting scenario simulations")
    import matplotlib.pyplot as plt

    if ax is None:
        _, ax = plt.subplots(figsize=(10, 5))
    if sim_df is None or sim_df.empty:
        ax.set_title(title or "Scenario Simulation (empty)")
        return ax

    sub = sim_df.copy()
    chosen_ts = ts_id or str(sub["ts_id"].iloc[0])
    sub = sub[sub["ts_id"].astype(str) == str(chosen_ts)]
    if sub.empty:
        ax.set_title(title or f"Scenario Simulation ({chosen_ts}: no rows)")
        return ax
    chosen_model = model_type or str(sub["model_type"].iloc[0])
    sub = sub[sub["model_type"].astype(str) == str(chosen_model)]

    val_col = "scenario_revenue" if show_revenue and "scenario_revenue" in sub.columns else "yhat"
    for scen, grp in sub.groupby("scenario", sort=False):
        ordered = grp.sort_values("ds")
        style = "--" if scen == "baseline" else "-"
        lw = 2.2 if scen == "baseline" else 1.8
        ax.plot(ordered["ds"], ordered[val_col], linestyle=style, linewidth=lw, label=str(scen))
        if val_col == "yhat" and scen == "baseline":
            ax.fill_between(
                ordered["ds"],
                ordered["yhat_lower"],
                ordered["yhat_upper"],
                alpha=0.15,
                label="baseline 80% PI",
            )

    ax.set_xlabel("Forecast Date")
    ax.set_ylabel("Revenue" if val_col == "scenario_revenue" else "Forecast Units (yhat)")
    ax.set_title(title or f"What-If Scenarios — {chosen_ts} ({chosen_model})")
    ax.legend(loc="best", frameon=True)
    ax.grid(alpha=0.25)
    return ax

simulate_scenarios

simulate_scenarios(
    data,
    cfg,
    scenarios,
    *,
    models=None,
    future_covariates_df=None,
    price_col=None,
)

Run counterfactual covariate scenarios across data and models (pure, offline).

Fits each candidate model once on each series' historical trajectory, then evaluates both the unshocked "baseline" future covariate trajectory and each named scenario in scenarios. When price_col is provided (or automatically detected from "price" / "price_index"), computes baseline_revenue, scenario_revenue, revenue_delta, revenue_delta_pct, and implied_elasticity ($\frac{\%\Delta \hat{y}}{\%\Delta p}$).

Source code in src/scale_forecasting/scenarios.py

def simulate_scenarios(
    data: pd.DataFrame,
    cfg: RunConfig,
    scenarios: Mapping[str, Mapping[str, Any]],
    *,
    models: Sequence[str] | None = None,
    future_covariates_df: pd.DataFrame | None = None,
    price_col: str | None = None,
) -> pd.DataFrame:
    """Run counterfactual covariate scenarios across ``data`` and ``models`` (pure, offline).

    Fits each candidate model once on each series' historical trajectory, then evaluates both the
    unshocked ``"baseline"`` future covariate trajectory and each named scenario in ``scenarios``.
    When ``price_col`` is provided (or automatically detected from ``"price"`` / ``"price_index"``),
    computes ``baseline_revenue``, ``scenario_revenue``, ``revenue_delta``, ``revenue_delta_pct``,
    and ``implied_elasticity`` ($\\frac{\\%\\Delta \\hat{y}}{\\%\\Delta p}$).
    """
    shock_cols = sorted({str(k) for shocks in scenarios.values() for k in shocks})
    known_covs = list(dict.fromkeys([*cfg.features.known_future_covariates, *shock_cols]))
    if price_col and price_col not in known_covs:
        known_covs.append(price_col)

    work_cfg = cfg
    if set(known_covs) != set(cfg.features.known_future_covariates):
        work_cfg = cfg.model_copy(
            update={"features": cfg.features.model_copy(update={"future_covariates": known_covs})}
        )

    eff_price_col = _infer_price_col(price_col, known_covs, scenarios)

    candidate_models = list(models) if models is not None else list(work_cfg.models)
    py_models = [
        m
        for m in candidate_models
        if get_model(m).runtime == "python" and get_model(m).is_available()
    ]
    if not py_models:
        py_models = ["naive_seasonal"]

    d = work_cfg.data
    panel = data.copy()
    if d.ts_id_col not in panel.columns:
        panel[d.ts_id_col] = "series_0"
    all_covs = list(dict.fromkeys([*work_cfg.features.all_covariates, *known_covs]))
    for c in all_covs:
        if c not in panel.columns:
            panel[c] = 1.0 if c == eff_price_col else 0.0

    all_scenarios: dict[str, Mapping[str, Any]] = {"baseline": {}, **dict(scenarios)}
    target_qs = work_cfg.output.effective_quantiles
    fallback_arm = corrected_arm_for(work_cfg.backtest.decision_metric)
    horizon = d.horizon

    rows: list[dict[str, Any]] = []
    for raw_tid, s_group in panel.groupby(d.ts_id_col, sort=True):
        tid = str(raw_tid)
        s_clean = regularize_series_grid(s_group, work_cfg)
        raw_target = s_clean[d.target_col].astype(float)
        lam = fit_transform_lambda(raw_target, work_cfg.features.transform)
        base_fut = _build_baseline_future_frame(
            s_clean, work_cfg, known_covs, future_covariates_df, tid
        )
        base_prices = (
            pd.to_numeric(base_fut[eff_price_col], errors="coerce").to_numpy(dtype=float)
            if eff_price_col and eff_price_col in base_fut.columns
            else None
        )

        for mname in py_models:
            m_cls = get_model(mname)
            spec = resolve_model_variant(work_cfg, mname, model_cls=m_cls)
            model_cfg = spec.effective_cfg
            m_params: dict[str, Any] = dict(spec.clean_params)
            ctx = ModelContext(
                freq=d.freq,
                horizon=model_cfg.max_horizon,
                transform=model_cfg.features.transform,
                transform_lambda=lam,
            )
            y_hist, X_hist = build_features(
                s_clean, model_cfg, lam, m_cls.lags_covariates_internally
            )
            est = m_cls(m_params, ctx)
            est.fit(y_hist, X_hist)

            cal: Calibration | None = None
            oof_df: pd.DataFrame | None = None
            if model_cfg.backtest.enabled:
                oof_raw, _, _ = backtest_cell(
                    s_clean,
                    partial(m_cls, m_params, ctx),
                    model_cfg,
                    lam,
                    None,
                    m_cls.lags_covariates_internally,
                )
                oof_df = oof_raw if not oof_raw.empty else None
                if oof_df is not None:
                    cal = calibrate_from_oof(oof_df, target_qs)

            arm, _ = resolve_arm(
                model_cfg.output.point_forecast,
                oof_df,
                model_cfg.backtest.decision_metric,
                fallback_arm,
            )

            baseline_yhat_arr: np.ndarray | None = None
            scen_frames: list[tuple[str, pd.DataFrame, np.ndarray | None]] = []

            for scen_name, shocks in all_scenarios.items():
                fut_cov_df = apply_scenario_shocks(base_fut, shocks, date_col=d.date_col)
                X_fut = build_future_features(
                    y_hist,
                    X_hist,
                    model_cfg,
                    horizon=horizon,
                    future_covariates_df=fut_cov_df,
                )
                raw_pred = est.predict(horizon, X_fut, quantiles=target_qs)
                cal_pred, _ = apply_calibration(raw_pred, cal, arm)
                cal_pred, _ = apply_quantile_and_cost_rules(cal_pred, model_cfg.output)
                if has_output_constraints(model_cfg.output):
                    cal_pred = apply_output_constraints(
                        cal_pred, model_cfg.output, bounds_df=fut_cov_df
                    )
                scen_prices = (
                    pd.to_numeric(fut_cov_df[eff_price_col], errors="coerce").to_numpy(dtype=float)
                    if eff_price_col and eff_price_col in fut_cov_df.columns
                    else None
                )
                if scen_name == "baseline":
                    baseline_yhat_arr = pd.to_numeric(cal_pred["yhat"], errors="coerce").to_numpy(
                        dtype=float
                    )
                scen_frames.append((scen_name, cal_pred, scen_prices))

            assert baseline_yhat_arr is not None
            for scen_name, cal_pred, scen_prices in scen_frames:
                yh = pd.to_numeric(cal_pred["yhat"], errors="coerce").to_numpy(dtype=float)
                yl = pd.to_numeric(cal_pred["yhat_lower"], errors="coerce").to_numpy(dtype=float)
                yu = pd.to_numeric(cal_pred["yhat_upper"], errors="coerce").to_numpy(dtype=float)
                ds_vals = pd.to_datetime(cal_pred["ds"]).tolist()
                for idx in range(len(cal_pred)):
                    b_yh = float(baseline_yhat_arr[idx])
                    s_yh = float(yh[idx])
                    d_yh = s_yh - b_yh
                    d_pct = (d_yh / abs(b_yh)) if abs(b_yh) > 1e-12 else 0.0
                    rec: dict[str, Any] = {
                        "scenario": scen_name,
                        "ts_id": tid,
                        "model_type": mname,
                        "ds": ds_vals[idx],
                        "horizon_step": idx + 1,
                        "yhat": s_yh,
                        "yhat_lower": float(yl[idx]),
                        "yhat_upper": float(yu[idx]),
                        "baseline_yhat": b_yh,
                        "delta_yhat": d_yh,
                        "delta_pct": d_pct,
                    }
                    if base_prices is not None and scen_prices is not None:
                        bp = float(base_prices[idx])
                        sp = float(scen_prices[idx])
                        brev = bp * b_yh
                        srev = sp * s_yh
                        rev_delta = srev - brev
                        rev_pct = (rev_delta / abs(brev)) if abs(brev) > 1e-12 else 0.0
                        p_pct = ((sp - bp) / abs(bp)) if abs(bp) > 1e-12 else 0.0
                        elasticity = (d_pct / p_pct) if abs(p_pct) > 1e-9 else float("nan")
                        rec.update(
                            {
                                "baseline_price": bp,
                                "price": sp,
                                "baseline_revenue": brev,
                                "scenario_revenue": srev,
                                "revenue_delta": rev_delta,
                                "revenue_delta_pct": rev_pct,
                                "implied_elasticity": elasticity,
                            }
                        )
                    rows.append(rec)

    return pd.DataFrame.from_records(rows)

summarize_scenarios

summarize_scenarios(sim_df)

Roll up step-level simulate_scenarios output into horizon totals per (scenario, ts_id, model_type).

Source code in src/scale_forecasting/scenarios.py

def summarize_scenarios(sim_df: pd.DataFrame) -> pd.DataFrame:
    """Roll up step-level `simulate_scenarios` output into horizon totals per
    ``(scenario, ts_id, model_type)``.
    """
    base_cols = [
        "scenario",
        "ts_id",
        "model_type",
        "n_steps",
        "baseline_units",
        "scenario_units",
        "delta_units",
        "delta_pct",
    ]
    rev_cols = [
        "baseline_revenue",
        "scenario_revenue",
        "revenue_delta",
        "revenue_delta_pct",
        "mean_implied_elasticity",
    ]
    if sim_df is None or sim_df.empty:
        return pd.DataFrame(columns=base_cols)

    has_rev = "scenario_revenue" in sim_df.columns
    records: list[dict[str, Any]] = []
    for (scen, tid, mtype), grp in sim_df.groupby(["scenario", "ts_id", "model_type"], sort=False):
        b_u = float(pd.to_numeric(grp["baseline_yhat"], errors="coerce").sum())
        s_u = float(pd.to_numeric(grp["yhat"], errors="coerce").sum())
        d_u = s_u - b_u
        d_pct = (d_u / abs(b_u)) if abs(b_u) > 1e-12 else 0.0
        rec: dict[str, Any] = {
            "scenario": str(scen),
            "ts_id": str(tid),
            "model_type": str(mtype),
            "n_steps": len(grp),
            "baseline_units": b_u,
            "scenario_units": s_u,
            "delta_units": d_u,
            "delta_pct": d_pct,
        }
        if has_rev:
            b_r = float(pd.to_numeric(grp["baseline_revenue"], errors="coerce").sum())
            s_r = float(pd.to_numeric(grp["scenario_revenue"], errors="coerce").sum())
            d_r = s_r - b_r
            r_pct = (d_r / abs(b_r)) if abs(b_r) > 1e-12 else 0.0
            el_vals = pd.to_numeric(grp["implied_elasticity"], errors="coerce").dropna()
            finite_el = el_vals[np.isfinite(el_vals)]
            rec.update(
                {
                    "baseline_revenue": b_r,
                    "scenario_revenue": s_r,
                    "revenue_delta": d_r,
                    "revenue_delta_pct": r_pct,
                    "mean_implied_elasticity": (
                        float(finite_el.mean()) if not finite_el.empty else float("nan")
                    ),
                }
            )
        records.append(rec)

    cols = base_cols + (rev_cols if has_rev else [])
    return pd.DataFrame.from_records(records, columns=cols)

report

Self-contained interactive HTML & JSON executive forecast report — pure, zero GCP imports.

Generates a single-file, zero-server, offline-capable .html report (plus optional .json summary payload) combining:

1. Executive Summary & KPI Cards — run identity, series count, model count, winning model, headline Forecast Value Added (FVA %), and active exception count.
2. Model Leaderboard & FVA Waterfall — full metric table alongside an inline SVG Forecast Value Added bar chart with Harvey–Leybourne–Newbold (1997) modified Diebold–Mariano significance indicators (p < 0.05).
3. Top-N Exception Triage Queue — interactive filterable table of flagged (ts_id, model_type) cells across worse_than_naive, chronic_bias, realized_drift, whipsaw_jump, and capacity_ceiling_hit with plain-English explanations.
4. Lead-Time Horizon Decay & Bias Audit — inline SVG step-by-step error curve ($h = 1 \dots H$) and Trigg Tracking Signal table.
5. Interactive Series & Multi-Quantile Fan Explorer — client-side ts_id and model_type selector rendering an interactive SVG chart of historical actuals, point forecast, $80\%$ calibrated interval, and multi-percentile fan (p10–p95).
6. Inventory & Replenishment Policy (ROP / Safety Stock) — dynamic safety stock, reorder point ($s$), order-up-to target ($S$), and recommended order quantity combining forward forecasts, AR(1) multi-step error covariance, and supplier lead-time variability.

build_report_payload

build_report_payload(
    *,
    title="Scale-Forecasting Executive Workbench Report",
    run_id="local-workbench",
    leaderboard_df=None,
    oof_df=None,
    predictions_df=None,
    history_df=None,
    fva_df=None,
    horizon_decay_df=None,
    bias_df=None,
    exceptions_df=None,
    realized_df=None,
    stability_df=None,
    inventory_policy_df=None,
    metric="wape",
)

Assemble the structured dictionary backing both HTML and JSON report exports.

Source code in src/scale_forecasting/report.py

def build_report_payload(
    *,
    title: str = "Scale-Forecasting Executive Workbench Report",
    run_id: str = "local-workbench",
    leaderboard_df: pd.DataFrame | None = None,
    oof_df: pd.DataFrame | None = None,
    predictions_df: pd.DataFrame | None = None,
    history_df: pd.DataFrame | None = None,
    fva_df: pd.DataFrame | None = None,
    horizon_decay_df: pd.DataFrame | None = None,
    bias_df: pd.DataFrame | None = None,
    exceptions_df: pd.DataFrame | None = None,
    realized_df: pd.DataFrame | None = None,
    stability_df: pd.DataFrame | None = None,
    inventory_policy_df: pd.DataFrame | None = None,
    metric: str = "wape",
) -> dict[str, Any]:
    """Assemble the structured dictionary backing both HTML and JSON report exports."""
    if fva_df is None and oof_df is not None and not oof_df.empty:
        fva_df = compute_fva(oof_df, metric=metric, series_df=history_df)
    if horizon_decay_df is None and oof_df is not None and not oof_df.empty:
        horizon_decay_df = compute_horizon_decay(oof_df, metric=metric)
    if bias_df is None and oof_df is not None and not oof_df.empty:
        bias_df = compute_bias_report(oof_df)
    if exceptions_df is None:
        exceptions_df = compute_exceptions_queue(
            oof_df=oof_df,
            predictions_df=predictions_df,
            realized_df=realized_df,
            stability_df=stability_df,
            metric=metric,
        )
    if inventory_policy_df is None and predictions_df is not None and not predictions_df.empty:
        inventory_policy_df = compute_inventory_policy(predictions_df, oof_df=oof_df)

    fva_recs = _df_to_records(fva_df)
    lb_recs = _df_to_records(leaderboard_df)
    exc_recs = _df_to_records(exceptions_df)
    decay_recs = _df_to_records(horizon_decay_df)
    bias_recs = _df_to_records(bias_df)
    real_recs = _df_to_records(realized_df)
    stab_recs = _df_to_records(stability_df)
    inv_recs = _df_to_records(inventory_policy_df)
    explorer = _build_explorer_series_map(predictions_df, history_df)

    best_model = "n/a"
    best_fva_pct: float | None = None
    if fva_recs:
        best_row = min(
            fva_recs,
            key=lambda r: float(r["mean_metric"]) if r.get("mean_metric") is not None else 1e9,
        )
        best_model = str(best_row.get("model_type", "n/a"))
        if best_row.get("fva_pct") is not None:
            best_fva_pct = round(float(best_row["fva_pct"]) * 100.0, 2)
    elif lb_recs:
        first_row = lb_recs[0]
        best_model = str(first_row.get("model_type") or first_row.get("model") or "n/a")

    return {
        "title": title,
        "run_id": run_id,
        "metric": metric,
        "kpis": {
            "n_series": len(explorer["ts_ids"]),
            "n_models": len(explorer["models"]) or len(lb_recs) or len(fva_recs),
            "champion_model": best_model,
            "best_fva_pct": best_fva_pct,
            "n_exceptions": len(exc_recs),
        },
        "leaderboard": lb_recs,
        "fva": fva_recs,
        "exceptions": exc_recs,
        "horizon_decay": decay_recs,
        "bias_report": bias_recs,
        "realized_accuracy": real_recs,
        "revision_stability": stab_recs,
        "inventory_policy": inv_recs,
        "explorer": explorer,
    }

export_report

export_report(
    target,
    out_path,
    *,
    title=None,
    metric=None,
    write_json=False,
)

Export a self-contained .html report (and optional .json sidecar) from a BakeOff, PlaygroundRun, Forecaster, or dictionary of DataFrames.

Source code in src/scale_forecasting/report.py

def export_report(
    target: Any,
    out_path: str | Path,
    *,
    title: str | None = None,
    metric: str | None = None,
    write_json: bool = False,
) -> Path:
    """Export a self-contained `.html` report (and optional `.json` sidecar) from a `BakeOff`,
    `PlaygroundRun`, `Forecaster`, or dictionary of DataFrames.
    """
    out = Path(out_path)
    out.parent.mkdir(parents=True, exist_ok=True)

    run_id = "local-workbench"
    eff_metric = metric or "wape"
    leaderboard_df: pd.DataFrame | None = None
    oof_df: pd.DataFrame | None = None
    predictions_df: pd.DataFrame | None = None
    history_df: pd.DataFrame | None = None
    fva_df: pd.DataFrame | None = None
    horizon_decay_df: pd.DataFrame | None = None
    bias_df: pd.DataFrame | None = None
    exceptions_df: pd.DataFrame | None = None
    inventory_policy_df: pd.DataFrame | None = None

    # 1. BakeOff duck-typing
    if hasattr(target, "leaderboard") and hasattr(target, "predictions") and hasattr(target, "oof"):
        leaderboard_df = target.leaderboard
        predictions_df = target.predictions
        oof_df = target.oof
        history_df = getattr(target, "series", None)
        if hasattr(target, "inventory_policy_df"):
            inventory_policy_df = target.inventory_policy_df()
        if hasattr(target, "config"):
            run_id = getattr(target.config, "run_name", run_id)
            if metric is None:
                eff_metric = target.config.backtest.decision_metric
    # 2. PlaygroundRun duck-typing
    elif hasattr(target, "result") and hasattr(target, "series"):
        res = target.result
        run_id = getattr(res, "run_id", run_id)
        predictions_df = res.predictions.assign(ts_id=res.ts_id, model_type=res.model_type)
        if res.oof is not None and not res.oof.empty:
            oof_df = res.oof.assign(ts_id=res.ts_id, model_type=res.model_type)
        history_df = target.series
        if hasattr(target, "inventory_policy_df"):
            inventory_policy_df = target.inventory_policy_df()
    # 3. Forecaster duck-typing
    elif hasattr(target, "leaderboard_df") and hasattr(target, "fva_df"):
        run_id = getattr(target, "run_id", run_id)
        leaderboard_df = target.leaderboard_df(run_id)
        fva_df = target.fva_df(run_id, metric=metric)
        horizon_decay_df = target.horizon_decay_df(run_id, metric=metric)
        bias_df = target.bias_report_df(run_id)
        predictions_df = target.predictions_df(run_id)
        if hasattr(target, "exceptions_df"):
            exceptions_df = target.exceptions_df(run_id)
        if hasattr(target, "inventory_policy_df"):
            inventory_policy_df = target.inventory_policy_df(run_id)
    # 4. Mapping of DataFrames
    elif isinstance(target, dict):
        run_id = str(target.get("run_id", run_id))
        leaderboard_df = target.get("leaderboard_df")
        oof_df = target.get("oof_df")
        predictions_df = target.get("predictions_df")
        history_df = target.get("history_df")
        fva_df = target.get("fva_df")
        horizon_decay_df = target.get("horizon_decay_df")
        bias_df = target.get("bias_df")
        exceptions_df = target.get("exceptions_df")
        inventory_policy_df = target.get("inventory_policy_df")

    eff_title = title or f"Scale-Forecasting Executive Report ({run_id})"
    html_str = render_html_report(
        title=eff_title,
        run_id=run_id,
        leaderboard_df=leaderboard_df,
        oof_df=oof_df,
        predictions_df=predictions_df,
        history_df=history_df,
        fva_df=fva_df,
        horizon_decay_df=horizon_decay_df,
        bias_df=bias_df,
        exceptions_df=exceptions_df,
        inventory_policy_df=inventory_policy_df,
        metric=eff_metric,
    )
    out.write_text(html_str, encoding="utf-8")

    if write_json:
        payload = build_report_payload(
            title=eff_title,
            run_id=run_id,
            leaderboard_df=leaderboard_df,
            oof_df=oof_df,
            predictions_df=predictions_df,
            history_df=history_df,
            fva_df=fva_df,
            horizon_decay_df=horizon_decay_df,
            bias_df=bias_df,
            exceptions_df=exceptions_df,
            inventory_policy_df=inventory_policy_df,
            metric=eff_metric,
        )
        json_path = out.with_suffix(".json")
        json_path.write_text(json.dumps(payload, indent=2), encoding="utf-8")

    return out

main

main(argv=None)

CLI entry point: python -m scale_forecasting.report --out report.html.

Source code in src/scale_forecasting/report.py

def main(argv: list[str] | None = None) -> int:
    """CLI entry point: ``python -m scale_forecasting.report --out report.html``."""
    parser = argparse.ArgumentParser(
        prog="python -m scale_forecasting.report",
        description="Generate a self-contained HTML Forecaster Workbench report.",
    )
    parser.add_argument(
        "--out",
        required=True,
        help="Output .html file path.",
    )
    parser.add_argument(
        "--models",
        nargs="+",
        default=["naive_seasonal", "theta", "holtwinters"],
        help="Models to include when running an offline bake-off report.",
    )
    parser.add_argument(
        "--horizon",
        type=int,
        default=14,
        help="Forecast horizon for offline bake-off report.",
    )
    parser.add_argument(
        "--n-folds",
        type=int,
        default=2,
        help="Backtest folds for offline bake-off report.",
    )
    parser.add_argument(
        "--json",
        action="store_true",
        help="Also write a companion .json summary payload next to --out.",
    )
    args = parser.parse_args(argv)

    from .playground import bakeoff

    bo = bakeoff(models=list(args.models), horizon=int(args.horizon), n_folds=int(args.n_folds))
    written = export_report(bo, args.out, write_json=bool(args.json))
    print(f"Wrote interactive HTML report to {written}")
    return 0

render_html_report

render_html_report(
    *,
    title="Scale-Forecasting Executive Workbench Report",
    run_id="local-workbench",
    leaderboard_df=None,
    oof_df=None,
    predictions_df=None,
    history_df=None,
    fva_df=None,
    horizon_decay_df=None,
    bias_df=None,
    exceptions_df=None,
    realized_df=None,
    stability_df=None,
    inventory_policy_df=None,
    metric="wape",
)

Render a complete, self-contained single-file interactive HTML report.

Source code in src/scale_forecasting/report.py

def render_html_report(
    *,
    title: str = "Scale-Forecasting Executive Workbench Report",
    run_id: str = "local-workbench",
    leaderboard_df: pd.DataFrame | None = None,
    oof_df: pd.DataFrame | None = None,
    predictions_df: pd.DataFrame | None = None,
    history_df: pd.DataFrame | None = None,
    fva_df: pd.DataFrame | None = None,
    horizon_decay_df: pd.DataFrame | None = None,
    bias_df: pd.DataFrame | None = None,
    exceptions_df: pd.DataFrame | None = None,
    realized_df: pd.DataFrame | None = None,
    stability_df: pd.DataFrame | None = None,
    inventory_policy_df: pd.DataFrame | None = None,
    metric: str = "wape",
) -> str:
    """Render a complete, self-contained single-file interactive HTML report."""
    payload = build_report_payload(
        title=title,
        run_id=run_id,
        leaderboard_df=leaderboard_df,
        oof_df=oof_df,
        predictions_df=predictions_df,
        history_df=history_df,
        fva_df=fva_df,
        horizon_decay_df=horizon_decay_df,
        bias_df=bias_df,
        exceptions_df=exceptions_df,
        realized_df=realized_df,
        stability_df=stability_df,
        inventory_policy_df=inventory_policy_df,
        metric=metric,
    )
    kpis = payload["kpis"]
    fva_svg = _render_fva_svg(payload["fva"])
    decay_svg = _render_decay_svg(payload["horizon_decay"])
    lb_html = _html_table(
        payload["leaderboard"] or payload["fva"],
        empty_msg="No leaderboard rows available.",
    )
    fva_html = _html_table(payload["fva"], empty_msg="No backtest OOF rows available for FVA.")
    exc_html = _html_table(
        payload["exceptions"],
        table_id="exceptions-table",
        empty_msg="Zero exceptions flagged — all series passed operational checks.",
    )
    decay_html = _html_table(payload["horizon_decay"], empty_msg="No horizon decay rows.")
    bias_html = _html_table(payload["bias_report"], empty_msg="No bias report rows.")
    inv_html = _html_table(
        payload["inventory_policy"],
        empty_msg="No forward predictions available for inventory policy translation.",
    )

    fva_kpi_str = f"{kpis['best_fva_pct']:+.1f}%" if kpis["best_fva_pct"] is not None else "n/a"
    explorer_json = _script_safe_json(payload["explorer"])

    return f"""<!DOCTYPE html>
<html lang="en">
<head>
<meta charset="utf-8" />
<meta name="viewport" content="width=device-width, initial-scale=1" />
<title>{html.escape(title)} — {html.escape(run_id)}</title>
<style>
  :root {{
    --bg: #f8fafc;
    --card: #ffffff;
    --text: #0f172a;
    --muted: #64748b;
    --border: #e2e8f0;
    --primary: #2563eb;
    --success: #16a34a;
    --warn: #d97706;
    --danger: #dc2626;
  }}
  * {{ box-sizing: border-box; }}
  body {{
    margin: 0; padding: 24px 32px;
    font-family: -apple-system, BlinkMacSystemFont, "Segoe UI", Roboto, sans-serif;
    background: var(--bg); color: var(--text); line-height: 1.5;
  }}
  header {{
    display: flex; justify-content: space-between; align-items: baseline;
    border-bottom: 1px solid var(--border); padding-bottom: 16px; margin-bottom: 20px;
  }}
  h1 {{ margin: 0; font-size: 1.45rem; }}
  .badge-run {{
    font-family: ui-monospace, SFMono-Regular, Menlo, monospace;
    background: #e0e7ff; color: #1e40af; padding: 4px 10px; border-radius: 6px; font-size: 0.85rem;
  }}
  .kpi-grid {{
    display: grid; grid-template-columns: repeat(auto-fit, minmax(170px, 1fr));
    gap: 14px; margin-bottom: 22px;
  }}
  .kpi-card {{
    background: var(--card); border: 1px solid var(--border); border-radius: 10px;
    padding: 14px 16px; box-shadow: 0 1px 2px rgba(0,0,0,0.03);
  }}
  .kpi-label {{ font-size: 0.78rem; text-transform: uppercase; color: var(--muted); }}
  .kpi-val {{ font-size: 1.35rem; font-weight: 700; margin-top: 4px; }}
  .tabs {{
    display: flex; gap: 8px; border-bottom: 1px solid var(--border); margin-bottom: 20px;
    flex-wrap: wrap;
  }}
  .tab-btn {{
    background: none; border: none; border-bottom: 3px solid transparent;
    padding: 10px 16px; cursor: pointer; font-weight: 600; color: var(--muted);
    font-size: 0.92rem;
  }}
  .tab-btn.active {{ color: var(--primary); border-bottom-color: var(--primary); }}
  .tab-panel {{ display: none; }}
  .tab-panel.active {{ display: block; }}
  .card {{
    background: var(--card); border: 1px solid var(--border); border-radius: 10px;
    padding: 18px 20px; margin-bottom: 20px;
  }}
  .card h2 {{ margin-top: 0; font-size: 1.1rem; }}
  .table-wrap {{ overflow-x: auto; }}
  table.data-table {{
    width: 100%; border-collapse: collapse; font-size: 0.85rem;
  }}
  table.data-table th, table.data-table td {{
    padding: 8px 10px; border-bottom: 1px solid var(--border); text-align: left;
  }}
  table.data-table th {{ background: #f1f5f9; font-weight: 600; color: #334155; }}
  .chart-svg {{ width: 100%; max-width: 820px; height: auto; display: block; margin: 8px 0; }}
  .svg-label {{ font-size: 12px; fill: #1e293b; }}
  .svg-val {{ font-size: 11px; font-weight: 600; fill: #0f172a; }}
  .svg-axis {{ font-size: 12px; fill: #475569; font-weight: 600; }}
  .filter-bar {{ display: flex; gap: 8px; margin-bottom: 12px; flex-wrap: wrap; }}
  .pill {{
    border: 1px solid var(--border); background: #f8fafc; border-radius: 999px;
    padding: 4px 12px; font-size: 0.8rem; cursor: pointer;
  }}
  .pill.active {{ background: var(--primary); color: #fff; border-color: var(--primary); }}
  .controls {{
    display: flex; gap: 14px; align-items: center; margin-bottom: 14px; flex-wrap: wrap;
  }}
  select {{
    padding: 6px 10px; border-radius: 6px; border: 1px solid var(--border); font-size: 0.9rem;
  }}
  .muted {{ color: var(--muted); font-style: italic; }}
</style>
</head>
<body>
<header>
  <h1>{html.escape(title)}</h1>
  <span class="badge-run">run_id: {html.escape(run_id)}</span>
</header>

<section class="kpi-grid">
  <div class="kpi-card">
    <div class="kpi-label">Series Evaluated</div>
    <div class="kpi-val">{kpis["n_series"]}</div>
  </div>
  <div class="kpi-card">
    <div class="kpi-label">Models &amp; Ensembles</div>
    <div class="kpi-val">{kpis["n_models"]}</div>
  </div>
  <div class="kpi-card">
    <div class="kpi-label">Champion Model ({html.escape(metric.upper())})</div>
    <div class="kpi-val">{html.escape(str(kpis["champion_model"]))}</div>
  </div>
  <div class="kpi-card">
    <div class="kpi-label">Headline FVA vs. Baseline</div>
    <div class="kpi-val">{fva_kpi_str}</div>
  </div>
  <div class="kpi-card">
    <div class="kpi-label">Flagged Exceptions</div>
    <div class="kpi-val">{kpis["n_exceptions"]}</div>
  </div>
</section>

<nav class="tabs">
  <button class="tab-btn active" onclick="switchTab('tab-fva', this)">
    Leaderboard &amp; FVA
  </button>
  <button class="tab-btn" onclick="switchTab('tab-exceptions', this)">
    Exception Triage ({kpis["n_exceptions"]})
  </button>
  <button class="tab-btn" onclick="switchTab('tab-decay', this)">
    Horizon Decay &amp; Bias
  </button>
  <button class="tab-btn" onclick="switchTab('tab-inventory', this)">
    Inventory Policy (ROP / SS)
  </button>
  <button class="tab-btn" onclick="switchTab('tab-explorer', this)">
    Interactive Series &amp; Fan Explorer
  </button>
</nav>

<section id="tab-fva" class="tab-panel active">
  <div class="card">
    <h2>Forecast Value Added (FVA) Waterfall &amp; Significance</h2>
    {fva_svg}
    {fva_html}
  </div>
  <div class="card">
    <h2>Model Leaderboard</h2>
    {lb_html}
  </div>
</section>

<section id="tab-exceptions" class="tab-panel">
  <div class="card">
    <h2>Top-N Exception Triage Queue</h2>
    <div class="filter-bar">
      <button class="pill active" onclick="filterExceptions('all', this)">All</button>
      <button class="pill" onclick="filterExceptions('worse_than_naive', this)">
        Worse Than Naive
      </button>
      <button class="pill" onclick="filterExceptions('chronic_bias', this)">
        Chronic Bias
      </button>
      <button class="pill" onclick="filterExceptions('realized_drift', this)">
        Realized Drift
      </button>
      <button class="pill" onclick="filterExceptions('whipsaw_jump', this)">
        Whipsaw Jump
      </button>
      <button class="pill" onclick="filterExceptions('capacity_ceiling_hit', this)">
        Capacity Hit
      </button>
    </div>
    {exc_html}
  </div>
</section>

<section id="tab-decay" class="tab-panel">
  <div class="card">
    <h2>Lead-Time Horizon Accuracy Decay</h2>
    {decay_svg}
    {decay_html}
  </div>
  <div class="card">
    <h2>Directional Bias &amp; Trigg Tracking Signal Audit</h2>
    {bias_html}
  </div>
</section>

<section id="tab-inventory" class="tab-panel">
  <div class="card">
    <h2>Inventory &amp; Replenishment Policy (ROP / Safety Stock / Order-Up-To)</h2>
    {inv_html}
  </div>
</section>

<section id="tab-explorer" class="tab-panel">
  <div class="card">
    <h2>Interactive Series &amp; Multi-Quantile Fan Explorer</h2>
    <div class="controls">
      <label>Series (<code>ts_id</code>):
        <select id="sel-ts" onchange="renderFanChart()"></select>
      </label>
      <label>Model (<code>model_type</code>):
        <select id="sel-model" onchange="renderFanChart()"></select>
      </label>
    </div>
    <div id="fan-chart-container"></div>
  </div>
</section>

<script>
const EXPLORER_DATA = {explorer_json};

function escHtml(s) {{
  return String(s)
    .replace(/&/g, '&amp;').replace(/</g, '&lt;').replace(/>/g, '&gt;')
    .replace(/"/g, '&quot;').replace(/'/g, '&#39;');
}}

function switchTab(panelId, btn) {{
  document.querySelectorAll('.tab-panel').forEach(p => p.classList.remove('active'));
  document.querySelectorAll('.tab-btn').forEach(b => b.classList.remove('active'));
  document.getElementById(panelId).classList.add('active');
  btn.classList.add('active');
}}

function filterExceptions(flag, btn) {{
  document.querySelectorAll('.filter-bar .pill').forEach(p => p.classList.remove('active'));
  btn.classList.add('active');
  const tbl = document.getElementById('exceptions-table');
  if (!tbl) return;
  tbl.querySelectorAll('tbody tr').forEach(tr => {{
    const flags = tr.getAttribute('data-flags') || '';
    tr.style.display = (flag === 'all' || flags.indexOf(flag) !== -1) ? '' : 'none';
  }});
}}

function initExplorer() {{
  const tsSel = document.getElementById('sel-ts');
  const modSel = document.getElementById('sel-model');
  if (!tsSel || !modSel) return;
  (EXPLORER_DATA.ts_ids || []).forEach(t => {{
    const o = document.createElement('option');
    o.value = t; o.textContent = t; tsSel.appendChild(o);
  }});
  (EXPLORER_DATA.models || []).forEach(m => {{
    const o = document.createElement('option');
    o.value = m; o.textContent = m; modSel.appendChild(o);
  }});
  renderFanChart();
}}

function renderFanChart() {{
  const box = document.getElementById('fan-chart-container');
  const tsSel = document.getElementById('sel-ts');
  const modSel = document.getElementById('sel-model');
  if (!box || !tsSel || !modSel) return;
  const tid = tsSel.value;
  const mtype = modSel.value;
  const sObj = (EXPLORER_DATA.series || {{}})[tid];
  if (!sObj || !sObj.models || !sObj.models[mtype]) {{
    box.innerHTML = '<p class="muted">No forecast trajectory available for this selection.</p>';
    return;
  }}
  const hist = sObj.history || [];
  const fc = sObj.models[mtype] || [];
  const allY = [];
  hist.forEach(p => {{ if (p.y !== null && p.y !== undefined) allY.push(p.y); }});
  fc.forEach(p => {{
    if (p.yhat !== null) allY.push(p.yhat);
    if (p.lower !== null) allY.push(p.lower);
    if (p.upper !== null) allY.push(p.upper);
    const q = p.quantiles || {{}};
    if (q.p10 !== undefined && q.p10 !== null) allY.push(q.p10);
    if (q.p90 !== undefined && q.p90 !== null) allY.push(q.p90);
    if (q.p95 !== undefined && q.p95 !== null) allY.push(q.p95);
  }});
  if (!allY.length) {{
    box.innerHTML = '<p class="muted">Empty series coordinates.</p>';
    return;
  }}
  const W = 780, H = 310, LP = 58, RP = 30, TP = 22, BP = 38;
  const PW = W - LP - RP, PH = H - TP - BP;
  const totalN = hist.length + fc.length;
  const minY = Math.min(...allY), maxY = Math.max(...allY);
  const spanY = (maxY - minY) || 1.0;
  const sx = i => LP + (totalN <= 1 ? PW * 0.5 : (i / (totalN - 1)) * PW);
  const sy = v => TP + PH - ((v - minY) / spanY) * PH;

  let polyBand = '';
  const bandTop = [], bandBot = [];
  fc.forEach((p, idx) => {{
    const x = sx(hist.length + idx);
    const lo = (p.lower !== null && p.lower !== undefined) ? p.lower : p.yhat;
    const hi = (p.upper !== null && p.upper !== undefined) ? p.upper : p.yhat;
    bandTop.push(`${{x.toFixed(1)}},${{sy(hi).toFixed(1)}}`);
    bandBot.unshift(`${{x.toFixed(1)}},${{sy(lo).toFixed(1)}}`);
  }});
  if (bandTop.length) {{
    const pts = bandTop.concat(bandBot).join(' ');
    polyBand = `<polygon points="${{pts}}" fill="#3b82f6" opacity="0.20" />`;
  }}

  const histPts = hist.map((p, idx) => `${{sx(idx).toFixed(1)}},${{sy(p.y).toFixed(1)}}`).join(' ');
  const fcPts = fc.map(
    (p, idx) => `${{sx(hist.length + idx).toFixed(1)}},${{sy(p.yhat).toFixed(1)}}`
  ).join(' ');
  const splitX = sx(Math.max(0, hist.length - 1));
  const sxStr = splitX.toFixed(1);
  const histLine = histPts
    ? `<polyline fill="none" stroke="#0f172a" stroke-width="2" points="${{histPts}}" />`
    : '';
  const fcLine = fcPts
    ? `<polyline fill="none" stroke="#2563eb" stroke-width="2.4" points="${{fcPts}}" />`
    : '';
  const lblRange = `min=${{minY.toFixed(1)}}, max=${{maxY.toFixed(1)}}`;
  const lblTxt = escHtml(`${{tid}} · ${{mtype}} (${{lblRange}})`);

  box.innerHTML = `
    <svg viewBox="0 0 ${{W}} ${{H}}" class="chart-svg">
      <rect x="${{LP}}" y="${{TP}}" width="${{PW}}" height="${{PH}}"
        fill="#f8fafc" stroke="#cbd5e1" />
      <line x1="${{sxStr}}" y1="${{TP}}" x2="${{sxStr}}" y2="${{TP + PH}}"
        stroke="#94a3b8" stroke-dasharray="4,3" />
      ${{polyBand}}
      ${{histLine}}
      ${{fcLine}}
      <text x="${{LP + 8}}" y="${{TP + 16}}" class="svg-label">${{lblTxt}}</text>
    </svg>`;
}}
window.addEventListener('DOMContentLoaded', initExplorer);
</script>
</body>
</html>
"""