diff --git a/docs/EVAL.md b/docs/EVAL.md index eeba100..eccad52 100644 --- a/docs/EVAL.md +++ b/docs/EVAL.md @@ -1120,3 +1120,79 @@ LightGBM result — needs narrowing to a statement about decisions rather than a No claim was contradicted outright. The Q13 negative result is independently predicted by the forecastability literature, which is the strongest single piece of evidence that the measurement apparatus here is working correctly. + +--- + +# Racing the predictability measures (2026-08-08) + +Q13 showed the shipped diagnostic does not generalise, and the forecastability literature +offered a specific reason and a specific fix: a single lagged correlation reads one +frequency, spectral entropy reads all of them. **The fix was tried and it did not work.** + +## What was raced + +Three training-free measures, scored against the same target — did forecasting pay, by total +economic cost, at a 6 h and 12 h commitment — on the same 43 serverless workloads from Q13 +plus the 5 non-GPU fleet workloads: + +1. **Daily autocorrelation** — what shipped. +2. **Spectral predictability**, `1 − normalised spectral entropy` over a Welch-averaged + periodogram — the field's standard forecastability feature. +3. **Low-frequency power fraction** — the share of spectral power at periods *longer than + the commitment window*. Not from the literature; motivated by this project's own Q1 + result, that a decision fixed for N hours can only exploit structure slower than N hours. + Alone among the three it is a function of the decision horizon, not of the series only. + +Comparison is by **AUC**, which needs no cutoff. Reporting each measure at its own best +threshold would let a measure win by fitting a cutoff to 43 points, so best-threshold +accuracy is shown only as a clearly-labelled optimistic bound. + +## Result + +Serverless workloads, n = 43, with a 20,000-draw permutation test on each AUC: + +| measure | AUC @ 6 h | p | AUC @ 12 h | p | +|---|---:|---:|---:|---:| +| daily autocorrelation | 0.539 | 0.67 | 0.438 | 0.50 | +| spectral predictability | 0.396 | 0.26 | 0.456 | 0.63 | +| low-frequency power | 0.595 | 0.31 | **0.700** | **0.027** | + +**Spectral entropy failed.** At 0.396 and 0.456 it is at or below a coin flip, and on the +fleet workloads it is actively anti-predictive (AUC 0.167 at 6 h, though n = 5 there makes +that number nearly meaningless). The literature's recommended instrument did not rescue the +diagnostic on this population, and predicting that it would was wrong. + +**The horizon-relative measure is the only one showing signal**, and it is not from the +literature — it comes from this project's own Q1 finding. At a twelve-hour commitment it +reaches AUC 0.700 with a nominal p = 0.027. + +## Why that 0.700 is not being shipped + +Six tests were run (3 measures x 2 windows). A Bonferroni-corrected threshold is p < 0.0083, +and 0.027 does not clear it. Picking the largest of six AUCs and quoting its uncorrected +p-value is precisely how a result gets manufactured, and this project's own doctrine says an +implausibly clean number is an artefact until proven otherwise. + +So the honest status is: **a promising lead on one cohort, at one horizon, that does not +survive correction for the number of things tried.** It is not a finding, and the shipped +diagnostic is unchanged. + +## What this means for the diagnostic + +The product keeps daily autocorrelation and keeps the scope caveat Q13 forced onto it. That +is not because autocorrelation is good — on serverless it is a coin flip — but because +nothing tested beat it well enough to justify swapping the claim a visitor reads. + +The measures are implemented, unit-tested against signals with known answers, and available +in `delphi.forecast.predictability` for the follow-up, which is now well defined: test +low-frequency power on a much larger cohort with the horizon fixed in advance, as a single +pre-registered hypothesis rather than one of six. On a fresh cohort with one test, p < 0.05 +would mean something. + +## Added to the negatives ledger + +- **Spectral entropy did not improve on daily autocorrelation** for predicting whether + forecasting pays, on either population — contradicting the expectation drawn from the + forecastability literature. +- **The one measure that showed signal does not survive multiple-comparison correction**, and + is reported as a lead rather than promoted into the product. diff --git a/scripts/compare_predictability.py b/scripts/compare_predictability.py new file mode 100644 index 0000000..1511b89 --- /dev/null +++ b/scripts/compare_predictability.py @@ -0,0 +1,164 @@ +"""Which training-free measure actually predicts whether forecasting pays? + +Q13 established that the shipped diagnostic — daily autocorrelation — does not generalise +past the fleet-aggregate traces it was derived on. The forecastability literature says why: +a single lagged correlation reads one frequency, while spectral entropy reads all of them, +and spectral measures are the established indicator for exactly this question. + +This races three candidates against the same target on the same workloads: + +* ``daily_autocorrelation`` — what shipped; +* ``spectral_predictability`` (1 - spectral entropy) — the field's standard; +* ``low_frequency_power_fraction`` — the share of power at periods longer than the + commitment window, motivated by this project's own Q1 result that what a fixed-for-N-hours + decision can exploit is structure slower than N hours. + +**Comparison is threshold-free.** Reporting each measure at its own best cutoff would let a +measure win by overfitting a cutoff to 43 points. AUC — the probability that a randomly +chosen paying workload scores above a randomly chosen non-paying one — needs no cutoff and +is the honest way to ask which measure carries more signal. Best-threshold accuracy is +printed alongside as an optimistic upper bound, clearly labelled as one. + +Run: `python scripts/compare_predictability.py --cohort 400 --sample 60` +""" + +from __future__ import annotations + +import argparse +from dataclasses import dataclass + +import numpy as np +from evaluate_diagnostic import build_workloads # type: ignore[import-not-found] +from evaluate_threshold import ( # type: ignore[import-not-found] + BIN_SECONDS, + build_candidates, + stratify, + strict_dominance, +) + +from delphi.forecast.predictability import ( + daily_autocorrelation, + low_frequency_power_fraction, + spectral_predictability, +) + + +@dataclass(frozen=True) +class Scored: + name: str + population: str + measures: dict[str, float] + paid: bool + + +def auc(scores: list[float], labels: list[bool]) -> float: + """Mann-Whitney U / ROC AUC. 0.5 is a coin flip; below 0.5 is anti-predictive.""" + positive = [s for s, y in zip(scores, labels, strict=True) if y] + negative = [s for s, y in zip(scores, labels, strict=True) if not y] + if not positive or not negative: + return float("nan") + order = np.argsort(np.asarray(scores, dtype=np.float64), kind="mergesort") + ranks = np.empty(len(scores), dtype=np.float64) + ranks[order] = np.arange(1, len(scores) + 1, dtype=np.float64) + # Average ranks within ties so a constant measure scores exactly 0.5, not 1.0. + values = np.asarray(scores, dtype=np.float64) + for value in np.unique(values): + mask = values == value + if mask.sum() > 1: + ranks[mask] = ranks[mask].mean() + positive_rank_sum = float(ranks[np.asarray(labels, dtype=bool)].sum()) + n_pos, n_neg = len(positive), len(negative) + return (positive_rank_sum - n_pos * (n_pos + 1) / 2.0) / (n_pos * n_neg) + + +def best_threshold_accuracy(scores: list[float], labels: list[bool]) -> tuple[float, float]: + """Optimistic upper bound: the best cutoff chosen *on this same data*.""" + best = (float("nan"), 0.0) + for cut in sorted(set(scores)): + predicted = [s >= cut for s in scores] + acc = sum(p == y for p, y in zip(predicted, labels, strict=True)) / len(labels) + if acc > best[1]: + best = (cut, acc) + return best + + +def measures_for(values: np.ndarray, day_steps: int, window_steps: int) -> dict[str, float]: + return { + "daily autocorr": daily_autocorrelation(values, day_steps), + "spectral predictability": spectral_predictability(values), + "low-freq power": low_frequency_power_fraction(values, window_steps), + } + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("--cohort", type=int, default=400) + parser.add_argument("--sample", type=int, default=60) + parser.add_argument("--windows", nargs="+", type=int, default=[6, 12]) + parser.add_argument("--seed", type=int, default=20260808) + args = parser.parse_args() + + print("# Which predictability measure actually works?\n") + print( + "Threshold-free comparison by AUC — the probability a paying workload scores above " + "a non-paying one. 0.50 is a coin flip. Best-threshold accuracy is an optimistic " + "upper bound, since the cutoff is fitted on the same data it is scored on.\n" + ) + + serverless = stratify(build_candidates(args.cohort, args.seed), args.sample, args.seed) + fleet = build_workloads(include_gpu=False) + + for window_hours in args.windows: + rows: list[Scored] = [] + + window_steps = window_hours * 3600 // BIN_SECONDS + for candidate in serverless: + workload = candidate.workload + wins, _ = strict_dominance(workload, window_steps) + rows.append( + Scored( + workload.name, + "serverless", + measures_for(workload.series.values, workload.day_steps, window_steps), + wins > 0, + ) + ) + + for workload in fleet: + steps = window_hours * 3600 // workload.series.step_seconds + wins, _ = strict_dominance(workload, steps) + rows.append( + Scored( + workload.name, + "fleet", + measures_for(workload.series.values, workload.day_steps, steps), + wins > 0, + ) + ) + + print(f"\n## {window_hours}-hour commitment\n") + for population in ("serverless", "fleet", "both"): + subset = [r for r in rows if population in (r.population, "both")] + labels = [r.paid for r in subset] + if len(set(labels)) < 2: + print( + f"### {population}: all outcomes identical, AUC undefined (n={len(subset)})\n" + ) + continue + print( + f"### {population} (n={len(subset)}, " + f"{sum(labels)} paid / {len(labels) - sum(labels)} did not)\n" + ) + print("| measure | AUC | best-threshold accuracy (optimistic) |") + print("|---|---:|---:|") + for key in ("daily autocorr", "spectral predictability", "low-freq power"): + scores = [r.measures[key] for r in subset] + area = auc(scores, labels) + cut, acc = best_threshold_accuracy(scores, labels) + print(f"| {key} | {area:.3f} | {acc:.0%} at >= {cut:.3f} |") + baseline = max(sum(labels), len(labels) - sum(labels)) / len(labels) + print(f"\nAlways-guess-the-majority baseline: {baseline:.0%}.\n") + + +if __name__ == "__main__": + main() diff --git a/src/delphi/api/app.py b/src/delphi/api/app.py index 350e473..4ac169d 100644 --- a/src/delphi/api/app.py +++ b/src/delphi/api/app.py @@ -32,6 +32,7 @@ ) from delphi.config import get_settings from delphi.control.newsvendor import CostRatio +from delphi.forecast.predictability import optional_lag_autocorrelation SNAPSHOT_ENV = "DELPHI_SNAPSHOT_PATH" DEFAULT_SNAPSHOT = Path("data/snapshot.json") @@ -190,12 +191,8 @@ def score(request: Annotated[DiagnosticRequest, Body()]) -> DiagnosticResponse: status_code=422, detail="a constant series has no structure to diagnose" ) - centred = values - values.mean() - def autocorrelation(lag: int) -> float | None: - if lag < 1 or len(values) <= lag + 8: - return None - return float(np.corrcoef(centred[:-lag], centred[lag:])[0, 1]) + return optional_lag_autocorrelation(values, lag) per_day = max(round(86400 / request.step_seconds), 1) daily = autocorrelation(per_day) diff --git a/src/delphi/forecast/predictability.py b/src/delphi/forecast/predictability.py new file mode 100644 index 0000000..53cd6c2 --- /dev/null +++ b/src/delphi/forecast/predictability.py @@ -0,0 +1,165 @@ +"""Training-free measures of whether a demand series is worth forecasting. + +The project shipped a diagnostic built on **daily autocorrelation**, and Q13 established +that it does not generalise: on individual Azure Functions workloads the 0.50 cutoff is a +coin flip and the relationship is not even monotone. The forecastability literature predicts +exactly that failure. A single lagged correlation interrogates one frequency; a workload can +be highly structured at a period the chosen lag does not look at, or carry a strong day-lag +correlation whose power sits at frequencies too fast to help a six-hour commitment. + +The established alternative is **spectral entropy** — the Shannon entropy of the normalised +power spectral density, which reads all frequencies at once. Near 0 means the power is +concentrated in a few components (a near-sinusoid, highly predictable); near 1 means it is +spread evenly (white noise, unforecastable). It is the standard forecastability feature and +recent work proposes it specifically as a fast, training-free indicator of whether +forecasting will beat a simple baseline — the same job DELPHI's diagnostic does. + +This module also carries a third measure that the project's own Q1 result motivates. +Forecasting pays over a commitment window only when demand carries structure *slower* than +that window; power at frequencies faster than the window averages out inside it and cannot +be exploited by a decision that is fixed for its duration. ``low_frequency_power_fraction`` +measures that share directly, and unlike the other two it is a function of the horizon being +committed to rather than of the series alone. + +Numpy only, deliberately: a diagnostic that must run inside a web request should not pull a +transitive scipy import, and a periodogram is fifteen lines. +""" + +import numpy as np +import numpy.typing as npt + +FloatArray = npt.NDArray[np.float64] + + +def _validate(values: FloatArray, minimum: int) -> FloatArray: + series = np.asarray(values, dtype=np.float64) + if series.ndim != 1: + raise ValueError("demand must be one-dimensional") + if len(series) < minimum: + raise ValueError(f"need at least {minimum} points, got {len(series)}") + if not np.isfinite(series).all(): + raise ValueError("demand must be finite") + if float(series.std()) <= 0: + raise ValueError("a constant series has no spectrum to measure") + return series + + +def lag_autocorrelation(values: FloatArray, lag: int) -> float: + """Pearson correlation between the series and itself at ``lag`` steps. + + The single definition of the project's headline diagnostic. It previously existed twice + — once in the API and once in the evaluation script — which is the same drift risk that + let the shipped verdict contradict `docs/EVAL.md` for a day. + """ + if lag < 1: + raise ValueError("lag must be positive") + series = _validate(values, lag + 9) + lagged, leading = series[:-lag], series[lag:] + if lagged.std() <= 0 or leading.std() <= 0: + raise ValueError("one half of the lagged pair is constant") + return float(np.corrcoef(lagged, leading)[0, 1]) + + +def optional_lag_autocorrelation(values: FloatArray, lag: int) -> float | None: + """``lag_autocorrelation``, returning ``None`` where it is not measurable. + + The API reports "could not be computed" rather than "no structure" for a series too + short to carry the lag, and those two must never be confused: one is missing evidence, + the other is evidence of absence. + """ + try: + return lag_autocorrelation(values, lag) + except ValueError: + return None + + +def daily_autocorrelation(values: FloatArray, day_steps: int) -> float: + """Correlation at a one-day lag — the measure the shipped diagnostic uses.""" + return lag_autocorrelation(values, day_steps) + + +def power_spectrum(values: FloatArray, *, segments: int = 8) -> tuple[FloatArray, FloatArray]: + """Welch-style averaged periodogram: returns (frequencies in cycles/step, power). + + A raw periodogram is an inconsistent estimator — its variance does not fall as the + series grows — so entropy computed from one is dominated by noise. Averaging Hann- + windowed half-overlapping segments trades frequency resolution for the variance + reduction that makes the entropy stable. The DC term is dropped: a series' mean level + says nothing about whether its *shape* is predictable. + """ + series = _validate(values, 32) + if segments < 1: + raise ValueError("segments must be positive") + + # Half-overlapping segments, each at least 16 points; fall back to one segment if short. + length = max(len(series) // max(segments, 1) * 2, 16) + length = min(length, len(series)) + step = max(length // 2, 1) + starts = range(0, len(series) - length + 1, step) + + window = np.hanning(length) + correction = float(np.sum(window**2)) + accumulated = np.zeros(length // 2 + 1, dtype=np.float64) + count = 0 + for start in starts: + chunk = series[start : start + length] + chunk = chunk - chunk.mean() + spectrum = np.abs(np.fft.rfft(chunk * window)) ** 2 / correction + accumulated += spectrum + count += 1 + if count == 0: # pragma: no cover - length is clamped to len(series) above + raise ValueError("series too short to segment") + + power = accumulated[1:] / count + freqs = np.fft.rfftfreq(length, d=1.0)[1:] + return freqs, power + + +def spectral_entropy(values: FloatArray, *, segments: int = 8) -> float: + """Normalised Shannon entropy of the power spectrum, in [0, 1]. + + 0 is a pure sinusoid; 1 is white noise. This is the field's standard forecastability + feature, and it is reported here in its raw orientation — **lower means more + predictable** — so it cannot be silently confused with an autocorrelation. + """ + _, power = power_spectrum(values, segments=segments) + total = float(power.sum()) + if total <= 0: + raise ValueError("spectrum carries no power") + density = power / total + density = density[density > 0] + if len(density) < 2: + return 0.0 + entropy = float(-(density * np.log(density)).sum()) + return entropy / float(np.log(len(power))) + + +def spectral_predictability(values: FloatArray, *, segments: int = 8) -> float: + """``1 - spectral_entropy``, so that higher means more predictable. + + Provided because every other number in the diagnostic points the same way, and a metric + whose direction is inverted relative to its neighbours is a reporting bug waiting to + happen. + """ + return 1.0 - spectral_entropy(values, segments=segments) + + +def low_frequency_power_fraction( + values: FloatArray, window_steps: int, *, segments: int = 8 +) -> float: + """Share of spectral power at periods **longer than** ``window_steps``. + + Motivated by this project's own Q1 result: what a commitment controller can exploit is + structure that persists across the window it is committing to. Variation faster than the + window averages out inside it, and a capacity level fixed for the whole window cannot + track it however well it was predicted. Unlike entropy and autocorrelation this is a + function of the decision horizon, not of the series alone — which is the point. + """ + if window_steps < 2: + raise ValueError("window_steps must be at least two") + freqs, power = power_spectrum(values, segments=segments) + total = float(power.sum()) + if total <= 0: + raise ValueError("spectrum carries no power") + slow = freqs < (1.0 / float(window_steps)) + return float(power[slow].sum() / total) diff --git a/tests/test_predictability.py b/tests/test_predictability.py new file mode 100644 index 0000000..2021c83 --- /dev/null +++ b/tests/test_predictability.py @@ -0,0 +1,102 @@ +"""Training-free predictability measures, checked against signals with known answers.""" + +import numpy as np +import pytest + +from delphi.forecast.predictability import ( + daily_autocorrelation, + low_frequency_power_fraction, + power_spectrum, + spectral_entropy, + spectral_predictability, +) + +DAY = 288 # 5-minute bins + + +def _sine(periods: float = 10.0, day_steps: int = DAY, noise: float = 0.0) -> np.ndarray: + rng = np.random.default_rng(20260808) + t = np.arange(int(periods * day_steps), dtype=np.float64) + values = 50 + 20 * np.sin(2 * np.pi * t / day_steps) + if noise: + values = values + rng.normal(0.0, noise, size=len(t)) + return values + + +def _white_noise(length: int = DAY * 10) -> np.ndarray: + rng = np.random.default_rng(7) + return 50 + rng.normal(0.0, 10.0, size=length) + + +def test_spectral_entropy_separates_a_sinusoid_from_white_noise() -> None: + """The defining property: near 0 for one frequency, near 1 for all of them.""" + assert spectral_entropy(_sine()) < 0.25 + assert spectral_entropy(_white_noise()) > 0.90 + + +def test_spectral_entropy_is_bounded_and_ordered_by_noise_level() -> None: + previous = spectral_entropy(_sine(noise=0.0)) + for noise in (2.0, 5.0, 10.0, 20.0): + current = spectral_entropy(_sine(noise=noise)) + assert 0.0 <= current <= 1.0 + assert current > previous, f"entropy must rise with noise (at sigma={noise})" + previous = current + + +def test_spectral_entropy_is_invariant_to_shift_and_scale() -> None: + """A forecastability measure must not depend on the units demand is reported in.""" + values = _sine(noise=4.0) + base = spectral_entropy(values) + assert spectral_entropy(values * 7.5) == pytest.approx(base, abs=1e-9) + assert spectral_entropy(values + 1000.0) == pytest.approx(base, abs=1e-9) + + +def test_spectral_predictability_is_the_complement() -> None: + values = _sine(noise=4.0) + assert spectral_predictability(values) == pytest.approx(1.0 - spectral_entropy(values)) + + +def test_low_frequency_fraction_tracks_the_commitment_window() -> None: + """A daily cycle is slow relative to six hours and fast relative to three days.""" + daily = _sine(noise=3.0) + assert low_frequency_power_fraction(daily, window_steps=DAY // 4) > 0.8 + assert low_frequency_power_fraction(daily, window_steps=DAY * 3) < 0.2 + + +def test_low_frequency_fraction_separates_fast_from_slow_structure() -> None: + """Two series, equally predictable, differing only in whether a window can use it.""" + t = np.arange(DAY * 10, dtype=np.float64) + slow = 50 + 20 * np.sin(2 * np.pi * t / DAY) + fast = 50 + 20 * np.sin(2 * np.pi * t / (DAY // 24)) + window = DAY // 4 + + # Both are near-pure tones, so spectral entropy cannot tell them apart... + assert abs(spectral_entropy(slow) - spectral_entropy(fast)) < 0.15 + # ...but only one carries structure a six-hour commitment can exploit. + assert low_frequency_power_fraction(slow, window) > 0.8 + assert low_frequency_power_fraction(fast, window) < 0.05 + + +def test_daily_autocorrelation_matches_the_shipped_definition() -> None: + values = _sine(noise=4.0) + expected = float(np.corrcoef(values[:-DAY], values[DAY:])[0, 1]) + assert daily_autocorrelation(values, DAY) == pytest.approx(expected) + + +def test_power_spectrum_drops_dc_and_returns_positive_frequencies() -> None: + freqs, power = power_spectrum(_sine(noise=2.0)) + assert len(freqs) == len(power) + assert (freqs > 0).all(), "the mean level is not a statement about predictability" + assert (power >= 0).all() + + +def test_measures_reject_degenerate_input() -> None: + for bad in (np.zeros(DAY * 2), np.full(DAY * 2, 5.0)): + with pytest.raises(ValueError): + spectral_entropy(bad) + with pytest.raises(ValueError): + spectral_entropy(np.array([1.0, 2.0, 3.0])) + with pytest.raises(ValueError): + spectral_entropy(np.array([1.0, np.nan] * 100)) + with pytest.raises(ValueError): + low_frequency_power_fraction(_sine(), window_steps=1)