diff --git a/docs/ADRs/023_collapsing_posterior_draws.md b/docs/ADRs/023_collapsing_posterior_draws.md index edef0bc0..ba996402 100644 --- a/docs/ADRs/023_collapsing_posterior_draws.md +++ b/docs/ADRs/023_collapsing_posterior_draws.md @@ -75,10 +75,13 @@ rented GPU — the numpy comes down first, and the conversion is re-runnable fro ### 2. The method is the one the model declares. It is never hard-coded. -`collapse_origin(..., aggregate_method=...)` implements `arithmetic_mean` and `median` — the same +`collapse_origin(..., estimator=...)` implements `arithmetic_mean` and `median` — the same two names `vhy_021` defines — and **refuses any other string** rather than falling back to a default. The converter's default is `arithmetic_mean` because all eight declare it. +*(Amended 2026-09-29: the parameter is now `estimator`, and the tool computes two more of them. +The set a model may **declare** is still exactly these two. See Amendment 1.)* + **That default and the models' declaration are two places stating one fact, so they are pinned together.** `tests/test_roster_conformance.py::test_collapse_declaration_matches_the_converter` imports `DEFAULT_AGGREGATE_METHOD` and asserts every roster member matches it, and that every @@ -207,8 +210,9 @@ python -m tools.collapse.collapse_predictions models/ --run-type calibrat python -m tools.collapse.plot_collapse_audit --draws-dir --out audit.png ``` -`--aggregate-method` exists and must be passed if a model ever stops declaring `arithmetic_mean`; -the roster test names it when that happens. +`--estimator` exists and must be passed if a model ever stops declaring `arithmetic_mean`; +the roster test names it when that happens. `--aggregate-method` remains accepted as the name it +had before Amendment 1, so runbooks and saved command lines keep working. **The plot script is part of the procedure, not a debugging aid.** The tests prove the arithmetic; they cannot tell you the field stopped looking like conflict. A person looks at the panels before @@ -252,6 +256,65 @@ a real run. --- +## Amendment 1 (2026-09-29) — an estimator this tool can compute is not a method a model may declare + +**Status:** Accepted. Extends §2; overrides nothing. + +### What changed + +`tools/collapse` gained two point estimators beyond the two `vhy_021` defines: + +| name | what it is | +|---|---| +| `q95` | the 95th percentile across draws, numpy's linear interpolation | +| `conditional_mean` | `E[y|y>0]` — the mean over the positive draws, equal to `E[y] / P(y>0)` | + +The parameter is renamed `aggregate_method` → `estimator`, because for these two the old name +was false: they are not aggregate methods and `vhy_021` does not define them. + +### Why the two sets stay separate + +`AGGREGATE_METHODS` is **contract vocabulary** — the only names a model may put in its +`aggregate_method`, checked against every roster member by +`tests/test_roster_conformance.py::test_collapse_declaration_matches_the_converter`. +`ESTIMATORS` is **what this tool can compute**, a strict superset. + +Merging them would be the whole defect: a config could declare `q95`, and the roster test — whose +job is to catch exactly that — would wave it through as a legitimate ADR-021 method. The +separation is pinned by `test_the_experiment_estimators_are_not_declarable_aggregate_methods`, +which asserts `AGGREGATE_METHODS` is a *strict* subset and names the difference. + +**The standing rule: adding an estimator means adding a key to `ESTIMATORS` and nothing else. +`AGGREGATE_METHODS` changes only when views-hydranet's ADR-021 changes.** + +### What the two extra estimators are for + +views-models#505's selector experiment. One posterior cube, collapsed three ways, submitted as +three apparent models (`_I` mean, `_II` q95, `_III` conditional mean), so the ensemble selector's +own criteria choose between a mean that under-predicts total fatalities ~5x and two estimators +that do not. Measured on the 2026-09-28 calibration run, summed over 13 origins and all 7 landed +models, `q95` and `conditional_mean` run **~4x the mean** — the right order to close that gap. + +They are **not** candidates for delivery. No model declares them, `arithmetic_mean` remains the +default and the delivered estimator, and the mean frames this amendment's code produces are +**byte-exact against the parquets delivered before it** — checked for all 7 models x 13 origins as +the acceptance condition for the change. + +### What is deliberately not recorded in the parquet + +Which estimator produced a frame. The output carries `month_id`, `priogrid_id` and three `pred_*` +columns and nothing else, by §5, and adding provenance to the frame would change a shape +`ensemble-updater` already reads. **Provenance lives in the output directory the caller picks, and +nowhere else** — so a frame moved out of its directory is a frame whose estimator is unknowable. +That is a real sharp edge and the reason the directories are named `I_mean`, `II_q95` and +`III_conditional_mean` rather than anything shorter. + +### Caveat carried forward + +At the roster's `D x K = 16`, `q95` sits between the 15th and 16th order statistics — a coarse +tail estimate that one draw moves. It is reported, not corrected for. Reading either new estimator +as a calibrated forecast rather than as a probe of the selector's criteria would be a misuse. + ## References - `tools/collapse/` — the converter and the audit plots; `tools/collapse/__init__.py` has the usage diff --git a/tests/test_collapse_predictions.py b/tests/test_collapse_predictions.py index 07fff5e4..dac74261 100644 --- a/tests/test_collapse_predictions.py +++ b/tests/test_collapse_predictions.py @@ -21,6 +21,7 @@ from tools.collapse.collapse_predictions import ( AGGREGATE_METHODS, DEFAULT_AGGREGATE_METHOD, + ESTIMATORS, TARGETS, CollapseError, collapse_origin, @@ -255,20 +256,30 @@ def test_the_default_method_is_the_one_the_roster_declares(): assert set(AGGREGATE_METHODS) == {"arithmetic_mean", "median"} +def test_the_experiment_estimators_are_not_declarable_aggregate_methods(): + """The separation is the point. `q95` and `conditional_mean` are estimators this tool can + compute; they are NOT names a model may declare. Collapsing the two sets would let a config + declare `q95` and have the roster test wave it through as an ADR-021 aggregate method.""" + assert set(AGGREGATE_METHODS) < set(ESTIMATORS), "aggregate methods must be a strict subset" + assert set(ESTIMATORS) - set(AGGREGATE_METHODS) == {"q95", "conditional_mean"} + for experiment_only in ("q95", "conditional_mean"): + assert experiment_only not in AGGREGATE_METHODS + + def test_median_is_available_and_is_not_the_mean(origin): """views-hydranet ADR-021 allows median. If a model ever declares it we must honour it.""" raw = np.load(origin / "lr_sb_best" / "y_pred.npy") - df = collapse_origin(origin, aggregate_method="median") + df = collapse_origin(origin, estimator="median") np.testing.assert_allclose( df["pred_lr_sb_best"].to_numpy(), np.median(raw.astype("float64"), axis=1), rtol=1e-12 ) assert not np.allclose(df["pred_lr_sb_best"], raw.astype("float64").mean(axis=1)) -def test_mutation_unknown_aggregate_method_is_refused(origin): +def test_mutation_unknown_estimator_is_refused(origin): """A typo must not fall back to the mean and ship an estimator nobody chose.""" - with pytest.raises(CollapseError, match="unknown aggregate_method"): - collapse_origin(origin, aggregate_method="geometric_mean") + with pytest.raises(CollapseError, match="unknown estimator"): + collapse_origin(origin, estimator="geometric_mean") def test_mutation_log_space_in_ONE_target_only_is_refused(tmp_path): @@ -319,3 +330,113 @@ def test_mutation_duplicate_identifier_rows_are_refused(origin): np.savez(d / "identifiers.npz", time=month, unit=unit) with pytest.raises(CollapseError, match="duplicate"): collapse_origin(origin) + + +# ── the views-models#505 selector experiment: q95 and E[y|y>0] ──────────────────── + + +def test_q95_is_available_and_is_not_the_mean(origin): + """Frame II. numpy's linear interpolation, pinned — a different interpolation would move + every value in the frame without changing a single row count.""" + raw = np.load(origin / "lr_sb_best" / "y_pred.npy").astype("float64") + df = collapse_origin(origin, estimator="q95") + np.testing.assert_allclose( + df["pred_lr_sb_best"].to_numpy(), np.quantile(raw, 0.95, axis=1), rtol=1e-12 + ) + assert not np.allclose(df["pred_lr_sb_best"], raw.mean(axis=1)) + + +def test_conditional_mean_is_the_mean_of_the_positive_draws(origin): + """Frame III, against the definition written the obvious way rather than the vectorised + way the implementation uses. If the two agree, the `sum / n_positive` shortcut is sound.""" + raw = np.load(origin / "lr_sb_best" / "y_pred.npy").astype("float64") + df = collapse_origin(origin, estimator="conditional_mean") + expected = np.array([row[row > 0].mean() if (row > 0).any() else 0.0 for row in raw]) + np.testing.assert_allclose(df["pred_lr_sb_best"].to_numpy(), expected, rtol=1e-12) + + +def test_conditional_mean_equals_mean_over_probability_of_positive(origin): + """The other half of the identity: `E[y|y>0] == E[y] / P(y>0)`. Stated as a test because + the docstring claims it, and a claim in a docstring is not evidence.""" + raw = np.load(origin / "lr_sb_best" / "y_pred.npy").astype("float64") + df = collapse_origin(origin, estimator="conditional_mean") + p_positive = (raw > 0).mean(axis=1) + live = p_positive > 0 + np.testing.assert_allclose( + df["pred_lr_sb_best"].to_numpy()[live], + raw.mean(axis=1)[live] / p_positive[live], + rtol=1e-12, + ) + + +def test_conditional_mean_is_never_below_the_arithmetic_mean(origin): + """A property, not an example: dividing by `n_positive <= D` cannot shrink the value. It is + what lets the scale guard's threshold stay unchanged for the experiment frames — neither + new estimator can drop a target's maximum below one the mean already cleared.""" + for target in TARGETS: + mean = collapse_origin(origin, estimator="arithmetic_mean")[f"pred_{target}"].to_numpy() + cond = collapse_origin(origin, estimator="conditional_mean")[f"pred_{target}"].to_numpy() + assert (cond >= mean - 1e-12).all() + + +def test_conditional_mean_of_an_all_silent_cell_is_zero_not_nan(tmp_path): + """The 99.6% case on real output. NaN here would be joined and scored by + `ensemble-updater` rather than refused, so the estimator must commit to 0.0.""" + o = tmp_path / "predictions_calibration_20260101_000000" / "origin_0" + draws = np.zeros((ROWS, DRAWS)) + draws[0, :] = 50.0 # one live cell, so the scale guard has something to clear + for target in TARGETS: + _write_target(o, target, draws) + df = collapse_origin(o, estimator="conditional_mean") + silent = df["pred_lr_sb_best"].to_numpy()[1:] + assert np.isfinite(silent).all(), "an all-zero cell must not produce NaN or inf" + assert (silent == 0.0).all() + + +def test_conditional_mean_equals_the_mean_when_every_draw_is_positive(tmp_path): + """The boundary of the inequality above: at `n_positive == D` the two estimators are the + same number. A frame III that differed from frame I everywhere would mean the positive mask + was wrong, and this is the case that would catch it.""" + o = tmp_path / "predictions_calibration_20260101_000000" / "origin_0" + rng = np.random.default_rng(7) + draws = rng.uniform(1.0, 100.0, size=(ROWS, DRAWS)) # no zeros at all + for target in TARGETS: + _write_target(o, target, draws) + mean = collapse_origin(o, estimator="arithmetic_mean")["pred_lr_sb_best"].to_numpy() + cond = collapse_origin(o, estimator="conditional_mean")["pred_lr_sb_best"].to_numpy() + np.testing.assert_allclose(cond, mean, rtol=1e-6) + + +def test_every_estimator_produces_the_same_shape_and_identifiers(origin): + """Whichever estimator ran, the frame the selector reads is the same shape with the same + keys. Only `pred_*` may differ — an estimator that dropped or reordered rows would make the + three frames non-comparable, which is the whole point of submitting them together.""" + frames = {name: collapse_origin(origin, estimator=name) for name in ESTIMATORS} + reference = frames[DEFAULT_AGGREGATE_METHOD] + for name, df in frames.items(): + assert list(df.columns) == list(reference.columns), name + pd.testing.assert_series_equal(df["month_id"], reference["month_id"], check_names=False) + pd.testing.assert_series_equal( + df["priogrid_id"], reference["priogrid_id"], check_names=False + ) + + +def test_convert_model_honours_the_estimator_end_to_end(tmp_path, origin): + """The CLI path, not just the function: two estimators over the same source must write + parquets that differ in value and agree in every identifier.""" + model = tmp_path / "m" + source = origin.parent + generated = model / "data" / "generated" + generated.mkdir(parents=True) + source.rename(generated / source.name) + + out_i = tmp_path / "frame_I" + out_ii = tmp_path / "frame_II" + written_i = convert_model(model, out_dir=out_i, estimator="arithmetic_mean") + written_ii = convert_model(model, out_dir=out_ii, estimator="q95") + assert len(written_i) == len(written_ii) == 1 + + a, b = pd.read_parquet(written_i[0]), pd.read_parquet(written_ii[0]) + pd.testing.assert_series_equal(a["month_id"], b["month_id"]) + pd.testing.assert_series_equal(a["priogrid_id"], b["priogrid_id"]) + assert not np.allclose(a["pred_lr_sb_best"], b["pred_lr_sb_best"]) diff --git a/tools/collapse/collapse_predictions.py b/tools/collapse/collapse_predictions.py index 5255ea83..da4f5e39 100644 --- a/tools/collapse/collapse_predictions.py +++ b/tools/collapse/collapse_predictions.py @@ -24,20 +24,32 @@ ## The collapse -Whatever the model declares in `aggregate_method` — `arithmetic_mean` or `median`, the two -views-hydranet `vhy_021` defines. Never hard-coded here; an unknown name is refused rather than -defaulted. All eight of the roster declare `arithmetic_mean`, and -`tests/test_roster_conformance.py` fails if that stops being true. +Two vocabularies, deliberately kept apart. + +**Aggregate methods** — `arithmetic_mean` and `median`, the two views-hydranet `vhy_021` +defines. These are contract vocabulary: the only names a model may *declare* in +`aggregate_method`. Never hard-coded here; an unknown name is refused rather than defaulted. +All eight of the roster declare `arithmetic_mean`, and `tests/test_roster_conformance.py` +fails if that stops being true. The mean is also the estimator the pipeline's own design points at: `feature_scaler.py:199` — *"Essential for accurate Arithmetic Mean collapse (ADR 021)"* — INVERT before COLLAPSE, so the mean is taken in count space. A better point estimate exists (`gate x mu`, ledger M70) but is unobtainable without replacing `main.py`. + +**Estimators** — the wider set this tool can compute, adding `q95` and `conditional_mean` +(`E[y|y>0]`). These are *not* aggregate methods and no model declares them. They exist for the +views-models#505 selector experiment: the same posterior cube submitted three times as three +apparent models, so the ensemble selector's own criteria decide between a mean that +under-predicts the total ~5x and two estimators that do not. Which estimator produced a frame +is not recoverable from the parquet — it is recorded by the output directory the caller picks, +and by nothing else. """ from __future__ import annotations import argparse +from collections.abc import Callable from pathlib import Path import numpy as np @@ -47,9 +59,56 @@ #: not fatalities, and are deliberately absent. TARGETS: tuple[str, ...] = ("lr_sb_best", "lr_ns_best", "lr_os_best") + +def _arithmetic_mean(draws: np.ndarray) -> np.ndarray: + return draws.mean(axis=1) + + +def _median(draws: np.ndarray) -> np.ndarray: + return np.median(draws, axis=1) + + +def _q95(draws: np.ndarray) -> np.ndarray: + """The 95th percentile across draws, linearly interpolated (numpy's default). + + At the roster's D x K = 16 this sits between the 15th and 16th order statistics, so it is + a coarse tail estimate — the upper draws are a small sample and one of them moves it. That + is a property of 16 draws, not of the quantile: it is reported, not corrected for. + """ + return np.quantile(draws, 0.95, axis=1) + + +def _conditional_mean(draws: np.ndarray) -> np.ndarray: + """`E[y|y>0]` — the mean over the positive draws only. + + Identical to `E[y] / P(y>0)`: a zero draw adds nothing to the sum, so + `sum / n_positive == (sum / D) / (n_positive / D)`. Both forms are the conditional + intensity; this one is written as the division that cannot divide by a probability of zero. + + A cell every draw calls silent has no conditional intensity to report, and gets **0.0**, + not NaN. `ensemble-updater` joins and scores these columns, so a NaN would propagate into + a metric instead of announcing itself. Where every draw is positive this equals the + arithmetic mean; it is never below it. + """ + n_positive = (draws > 0).sum(axis=1) + return np.where(n_positive > 0, draws.sum(axis=1) / np.maximum(n_positive, 1), 0.0) + + +#: Every point estimator this tool can compute, by its command-line name. +ESTIMATORS: dict[str, Callable[[np.ndarray], np.ndarray]] = { + "arithmetic_mean": _arithmetic_mean, + "median": _median, + "q95": _q95, + "conditional_mean": _conditional_mean, +} + # views-hydranet ADR-021 defines exactly these two and rejects anything else # (`volume_handler.py::collapse_to_point`). We implement the same two, under the same names, # so a model's declared `aggregate_method` can be passed straight through. +# +# `q95` and `conditional_mean` are deliberately NOT here. They are alternative estimators for +# the views-models#505 experiment, not contract vocabulary, and a model that declared one would +# be a config error `tests/test_roster_conformance.py` must keep catching. AGGREGATE_METHODS: tuple[str, ...] = ("arithmetic_mean", "median") DEFAULT_AGGREGATE_METHOD = "arithmetic_mean" @@ -101,7 +160,7 @@ def _load_target(target_dir: Path) -> tuple[np.ndarray, np.ndarray, np.ndarray]: def collapse_origin( origin_dir: Path, targets: tuple[str, ...] = TARGETS, - aggregate_method: str = DEFAULT_AGGREGATE_METHOD, + estimator: str = DEFAULT_AGGREGATE_METHOD, ) -> pd.DataFrame: """One origin directory -> one DataFrame of point predictions. @@ -109,11 +168,14 @@ def collapse_origin( identically row-for-row. A mismatch is an error, never a merge: silently joining on keys would paper over a misalignment that changes which cell a number belongs to. """ - if aggregate_method not in AGGREGATE_METHODS: + try: + collapse = ESTIMATORS[estimator] + except KeyError: raise CollapseError( - f"unknown aggregate_method {aggregate_method!r}; " - f"views-hydranet ADR-021 defines only {', '.join(AGGREGATE_METHODS)}" - ) + f"unknown estimator {estimator!r}; this tool implements " + f"{', '.join(sorted(ESTIMATORS))} — of which views-hydranet ADR-021 defines " + f"only {', '.join(AGGREGATE_METHODS)} as a declarable aggregate_method" + ) from None if not origin_dir.is_dir(): raise CollapseError(f"not a directory: {origin_dir}") @@ -152,10 +214,7 @@ def collapse_origin( ) wide = draws.astype("float64") # float32 sums make the answer depend on the draw count - if aggregate_method == "median": - frame[f"pred_{target}"] = np.median(wide, axis=1) - else: - frame[f"pred_{target}"] = wide.mean(axis=1) + frame[f"pred_{target}"] = collapse(wide) assert frame is not None # targets is non-empty by construction @@ -202,7 +261,7 @@ def convert_model( out_dir: Path | None = None, targets: tuple[str, ...] = TARGETS, min_plausible_max: float = MIN_PLAUSIBLE_MAX, - aggregate_method: str = DEFAULT_AGGREGATE_METHOD, + estimator: str = DEFAULT_AGGREGATE_METHOD, ) -> list[Path]: """Convert every origin of a model's latest prediction directory. Returns files written. @@ -230,7 +289,7 @@ def convert_model( written: list[Path] = [] for origin in origins: index = int(origin.name.split("_")[1]) - frame = collapse_origin(origin, targets, aggregate_method) + frame = collapse_origin(origin, targets, estimator) _check_scale(frame, origin, min_plausible_max) path = destination / f"predictions_{run_type}_{timestamp}_{index:02d}.parquet" frame.to_parquet(path, index=False) @@ -245,10 +304,16 @@ def main(argv: list[str] | None = None) -> int: parser.add_argument("--out-dir", type=Path, default=None) parser.add_argument("--min-plausible-max", type=float, default=MIN_PLAUSIBLE_MAX) parser.add_argument( - "--aggregate-method", + "--estimator", + "--aggregate-method", # the name before q95/conditional_mean existed; still accepted + dest="estimator", default=DEFAULT_AGGREGATE_METHOD, - choices=AGGREGATE_METHODS, - help="must match the model's own `aggregate_method` (all eight declare arithmetic_mean)", + choices=sorted(ESTIMATORS), + help=( + "arithmetic_mean or median must match the model's own `aggregate_method` (all " + "eight declare arithmetic_mean); q95 and conditional_mean are the views-models#505 " + "selector experiment and are declared by no model" + ), ) args = parser.parse_args(argv) @@ -257,7 +322,7 @@ def main(argv: list[str] | None = None) -> int: run_type=args.run_type, out_dir=args.out_dir, min_plausible_max=args.min_plausible_max, - aggregate_method=args.aggregate_method, + estimator=args.estimator, ) for path in written: print(path)