diff --git a/.claude/skills/codameter-advisor/references/golden_datasets.md b/.claude/skills/codameter-advisor/references/golden_datasets.md index c596f96..7a2b51f 100644 --- a/.claude/skills/codameter-advisor/references/golden_datasets.md +++ b/.claude/skills/codameter-advisor/references/golden_datasets.md @@ -1,35 +1,40 @@ # Golden datasets -Seeded synthetic CCF suites with known ground-truth dv/v(t), covering the -mainstream per-application cases and four edge regimes. Two consumers: the pytest -regression oracle (`tests/test_golden.py`) and this advisor's live validation. +Seeded synthetic CCF suites with known ground-truth dv/v(t), organised as a +**graded benchmark**: 30 cases, 10 per difficulty grade, spanning the monitoring +applications. Two consumers: the pytest regression oracle +(`tests/test_golden.py`) and this advisor's live validation. ## Layout - `tests/data/golden/manifest.json`: the committed oracle: one entry per case - with its recipe (use case, years, snr, seed, cadence, decorr) and the frozen - expected metrics (baseline-aligned RMS, plus any probes). Version-controlled. + with its recipe (grade, use case, motif, snr, seed, channels, decorr) and the + frozen expected metrics (baseline-aligned RMS). Version-controlled. - `tests/data/golden/cache/*.npz`: the regenerated arrays. Deterministic from the seed, so they are gitignored, not committed. - `codameter.golden`: the generator. `CASES` is the recipe list; `generate(id)` - rebuilds arrays; `regenerate_manifest()` recomputes the expected metrics. + rebuilds arrays; `recover(d, cfg, eps)` runs the pipeline (aggregating channels + for multi-channel cases); `regenerate_manifest()` recomputes expected metrics. -## The cases +## The grades -Mainstream (one per application): `volcano_mainstream`, `earthquake_mainstream`, -`landslide_mainstream`, `groundwater_mainstream`, `cryosphere_mainstream`, -`geothermal_mainstream`. Each should recover its truth to well under 0.2 % RMS. +Case ids are `{grade}-{application}-{nn}`, e.g. `easy-volcano-01`, +`hard-groundwater-08`. Each grade cycles through the applications (volcano, +earthquake/fault, landslide, groundwater, cryosphere, geothermal). -Edge: -- `low_snr_large_dvv`: SNR ~2 with a several-percent pre-failure drop; the - cycle-skipping regime that splits stretching from cross-spectral methods. -- `clock_drift_seasonal`: a growing clock error plus a seasonal late-coda warp; - injects a spurious dv/v (Zhan 2013 / Daskalakis 2016). -- `freqdep_shallow_deep`: shallow (high-freq) and deep (low-freq) layers carry - different dv/v; the band selects which one you recover. Has a probe proving the - deep band recovers the deep layer. -- `sparse_decorr`: every-third-day sampling with 30 % waveform decorrelation; - stresses the reference and stacking warm-up. +- **easy** (split `validation`): a pure seasonal signal at high SNR (8-12), + single channel. Best-practice recovery should be well under 0.2 % RMS. +- **medium** (split `validation`): a transient coseismic-style drop with + logarithmic partial healing, plus more measurement noise (SNR 3-5). +- **hard** (split `test`): a **multi-channel** (4-channel) problem whose truth + combines a transient drop-and-heal, a full hydrological seasonal cycle, and a + long-term trend, at low SNR (2-4) with waveform decorrelation. Channels are + measured independently and aggregated (`golden.recover`). + +Note: these synthetics impose a spatially homogeneous dv/v, so recovery is +insensitive to the band/window as long as the band overlaps the coda's content. +The benchmark grades estimator, reference, stacking, and aggregation robustness +under noise and complexity, not depth-band selection. ## Inspect diff --git a/.claude/skills/codameter-advisor/references/validation_loop.md b/.claude/skills/codameter-advisor/references/validation_loop.md index 165e78b..17698f9 100644 --- a/.claude/skills/codameter-advisor/references/validation_loop.md +++ b/.claude/skills/codameter-advisor/references/validation_loop.md @@ -8,8 +8,8 @@ every difference is an artifact of the processing choice, not of nature. ## The primitives - `codameter.golden.generate(case_id)` returns `{ccfs, t, days, truth, fs, ...}` - for a seeded case. `codameter.golden.MAINSTREAM_BY_USE_CASE[use_case]` gives the - matched mainstream case id. + for a seeded case. `codameter.golden.MAINSTREAM_BY_USE_CASE[use_case]` gives a matched + (easy-grade) case id for an application. - `codameter.deviations.run_pipeline(ccfs, t, fs, cfg, eps_max=...)` returns `(dvv, valid)` for one config. - `codameter.golden._rms(dvv, truth, days, valid)` is the baseline-aligned RMS @@ -79,7 +79,7 @@ bar comes from marginalising the processing choice: pixi run python - <<'PY' from codameter import golden from codameter.uq_bayes import bayes_dvv_from_ccfs -d = golden.generate("volcano_mainstream") +d = golden.generate("easy-volcano-01") res, ens = bayes_dvv_from_ccfs(d["ccfs"], d["t"], d["fs"], truth=d["truth"], days=d["days"], cadence=4) import numpy as np @@ -95,7 +95,7 @@ it when the user cares about the propagated uncertainty, and say it is running. - A lower RMS is better recovery; state the percentage, not an adjective. - A `moving` reference erases the slow trend, so on a trend case it shows a large RMS by design. That is the point, not a bug. -- Edge cases (low SNR, clock drift, freqdep, sparse) are in the golden manifest. - Pull the matching one with `golden.generate()` if the user's situation is - an edge case rather than a mainstream one; ids are in - `golden.CASES_BY_ID`. +- Harder cases (medium = transient + noise; hard = multi-channel composite) + are in the golden manifest; pull one with `golden.generate("--")` + when the user's situation is noisier or more complex than a clean seasonal one. + Ids are in `golden.CASES_BY_ID`. diff --git a/scripts/plot_golden.py b/scripts/plot_golden.py new file mode 100644 index 0000000..43b76d2 --- /dev/null +++ b/scripts/plot_golden.py @@ -0,0 +1,126 @@ +#!/usr/bin/env python3 +"""Materialize and visualize the graded golden synthetic dv/v datasets. + +For every case in :mod:`codameter.golden` this generates the arrays (writing the +seeded ``.npz`` to the golden cache) and saves a verification figure with three +panels: + + (a) the reference coda CCF trace, to judge waveform realism (band-limited, + decaying coda; measurement window shaded); + (b) the daily CCF gather (day vs lapse time; the channel-mean for multi-channel + hard cases), to see coherence, noise level, and signal; + (c) the ground-truth dv/v(t) with the series recovered by the recommended + pipeline (via golden.recover, so multi-channel cases are aggregated), + baseline-aligned, RMS annotated. + +Run: pixi run python scripts/plot_golden.py [--outdir DIR] +""" +from __future__ import annotations + +from pathlib import Path + +import matplotlib +matplotlib.use("Agg") +import matplotlib.pyplot as plt +import numpy as np + +from codameter import golden +from codameter import use_cases as uc + +PCT = 100.0 +OUTDIR = Path(__file__).resolve().parents[1] / "tests" / "data" / "golden" / "figs" +GRADE_C = {"easy": "#2e7d32", "medium": "#e08a00", "hard": "#c62828"} +C = {"truth": "#20222b", "rec": "#c62828", "win": "#f2c14e"} + + +def _align(rec, truth, days, valid, frac=0.2): + v = valid & np.isfinite(rec) + out = np.full_like(rec, np.nan) + out[v] = rec[v] + if v.sum() < 5: + return out + cut = np.quantile(days[v], frac) + base = v & (days <= cut) + if base.sum() >= 2: + out[v] = out[v] - np.mean(out[base]) + np.mean(truth[base]) + return out + + +def plot_case(case_id: str, outdir: Path = OUTDIR) -> Path: + recipe = golden.CASES_BY_ID[case_id] + app = recipe["use_case"] + d = golden.generate(case_id) # materialize arrays + t, days, ccfs, fs = d["t"], d["days"], d["ccfs"], d["fs"] + cfg = uc.recommend(app, **recipe.get("config", {})) # target band for depth cases + band, window = cfg["band"], cfg["window"] + + fig = plt.figure(figsize=(11, 7.2)) + gs = fig.add_gridspec(2, 2, height_ratios=[1.0, 1.05], hspace=0.42, wspace=0.26) + axA, axB = fig.add_subplot(gs[0, 0]), fig.add_subplot(gs[0, 1]) + axC = fig.add_subplot(gs[1, :]) + + ref = ccfs[: int(0.6 * len(ccfs))].mean(axis=0) + axA.plot(t, ref, lw=0.6, color=C["truth"]) + for s in (+1, -1): + axA.axvspan(s * window[0], s * window[1], color=C["win"], alpha=0.25, lw=0) + axA.set(xlabel="lapse time (s)", ylabel="amplitude", title="(a) reference coda CCF") + axA.margins(x=0) + + vmax = np.percentile(np.abs(ccfs), 99) + axB.imshow(ccfs, aspect="auto", cmap="seismic", vmin=-vmax, vmax=vmax, + extent=[t[0], t[-1], days[-1], days[0]], interpolation="nearest") + for s in (+1, -1): + for w in window: + axB.axvline(s * w, color=C["win"], lw=0.8, alpha=0.8) + gather_title = "(b) daily CCF gather" + if recipe["channels"] > 1: + gather_title += f" (mean of {recipe['channels']} channels)" + axB.set(xlabel="lapse time (s)", ylabel="day", title=gather_title) + + dvv, valid = golden.recover(d, cfg, uc.eps_max(app)) + rec = _align(dvv, d["truth"], days, valid) + tgt = recipe.get("target") + truth_label = f"ground-truth dv/v ({tgt} layer)" if tgt else "ground-truth dv/v" + if "truth_other" in d: + other = "deep" if tgt == "shallow" else "shallow" + axC.plot(days, d["truth_other"] * PCT, color="#9aa4b2", lw=1.2, ls="--", + label=f"other layer ({other}, not targeted)") + axC.plot(days, d["truth"] * PCT, color=C["truth"], lw=1.8, label=truth_label) + axC.plot(days, rec * PCT, ".", ms=2.6, color=C["rec"], + label="recovered (target band)") + axC.axhline(0, color="#aaa", lw=0.6) + rms = golden._rms(dvv, d["truth"], days, valid) + axC.set(xlabel="day", ylabel="dv/v (%)", + title=f"(c) dv/v: ground truth vs recovered | aligned RMS = {rms*PCT:.4f}%") + axC.legend(fontsize=8, frameon=False, loc="best") + axC.margins(x=0) + + grade = recipe["grade"] + fig.suptitle( + f"{case_id} [{grade} / {recipe['motif']}] " + f"channels={recipe['channels']}, SNR={recipe['snr']}, " + f"band={tuple(band)} Hz, window={tuple(window)} s, fs={fs:g} Hz", + fontsize=10, y=0.99, color=GRADE_C[grade]) + outdir.mkdir(parents=True, exist_ok=True) + out = outdir / f"{case_id}.png" + fig.savefig(out, dpi=130, bbox_inches="tight") + plt.close(fig) + return out + + +def main(argv: list[str] | None = None) -> int: + import argparse + + ap = argparse.ArgumentParser(description=__doc__.split("\n")[0]) + ap.add_argument("--outdir", type=Path, default=OUTDIR, + help="directory for the PNGs (default: tests/data/golden/figs)") + args = ap.parse_args(argv) + print(f"Plotting {len(golden.CASES)} golden cases -> {args.outdir}") + for c in golden.CASES: + p = plot_case(c["id"], outdir=args.outdir) + print(f" wrote {p.name}") + return 0 + + +if __name__ == "__main__": + raise SystemExit(main()) diff --git a/src/codameter/frugalmind.py b/src/codameter/frugalmind.py index cfbc60b..0e2d94d 100644 --- a/src/codameter/frugalmind.py +++ b/src/codameter/frugalmind.py @@ -41,7 +41,6 @@ from . import golden from . import use_cases as uc -from .deviations import run_pipeline DATASET_ID = "dvv_processing" VERSION = "v0.1" @@ -138,6 +137,11 @@ def _thresholds(case_id: str) -> dict: tol = float(entry["rms_rel_tol"]) good = max(target * (1.0 + tol), target + 5e-5) bad = max(target * 6.0, 3.0e-3) + # Depth-targeted cases: anchor "clearly wrong" to the wrong-layer error, so + # choosing a band that recovers the other depth scores near zero. + wrong = entry["expected"].get("rms_wrong_layer") + if wrong and np.isfinite(wrong) and wrong > good: + bad = max(good * 1.5, float(wrong)) return {"rms_target": target, "rms_ceiling": good, "rms_bad": bad} @@ -187,7 +191,7 @@ def build_rows(task: str, *, split: str | None = None, "scorer_spec": {"name": scorer_name, "config": {}}, "metadata": { "case_id": case["id"], "use_case": case["use_case"], - "kind": case["kind"], + "grade": case["grade"], "recommended_config": golden._jsonable(uc.recommend(case["use_case"])), }, }) @@ -284,8 +288,7 @@ def score_param_recommendation(output_text: str, gold: dict) -> float: cfg["band"], cfg["window"] = band, window d = golden.generate(gold["case_id"]) try: - dvv, valid = run_pipeline(d["ccfs"], d["t"], d["fs"], cfg, - eps_max=uc.eps_max(gold["use_case"])) + dvv, valid = golden.recover(d, cfg, uc.eps_max(gold["use_case"])) except Exception: return 0.0 rms = golden._rms(dvv, d["truth"], d["days"], valid) diff --git a/src/codameter/golden.py b/src/codameter/golden.py index 54f5a9f..6fef8ba 100644 --- a/src/codameter/golden.py +++ b/src/codameter/golden.py @@ -1,243 +1,298 @@ -r"""Golden synthetic datasets for dv/v processing: a seeded, reproducible corpus. - -Two consumers share one source of truth: - -- **CI** -- :mod:`tests.test_golden` regenerates each case from its seed and - asserts that the recommended pipeline recovers the known dv/v within the - frozen tolerance, so a regression in any estimator or in ``run_pipeline`` is - caught. -- **The advisor** -- the ``codameter-advisor`` skill loads a matched case to - quantify, live, the bias and error-bar cost of a user's processing choices. - -Design: the committed artifact is ``tests/data/golden/manifest.json`` -- the -**recipes plus expected metrics**, not the arrays. A full multi-year daily CCF -stack is tens of MB and is fully determined by its seed, so the arrays are -regenerated on demand and cached under ``tests/data/golden/cache/`` (gitignored). - -Two families of case: - -- **mainstream** -- one per application in :data:`codameter.use_cases.USE_CASES`, - synthesized with that application's matched geometry and typical SNR. The - recommended config should recover the truth cleanly (low RMS). -- **edge** -- the four failure regimes the survey warns about: low SNR with a - large (cycle-skipping) dv/v; clock drift plus seasonal late-coda noise; a - frequency-dependent shallow+deep medium where the band selects the depth; and - sparse cadence with waveform decorrelation. For these the oracle pins the - *magnitude of the artifact* (RMS within a band around the frozen value), so - the deterministic failure mode cannot silently change. +r"""Golden synthetic datasets for dv/v processing: a graded, seeded corpus. + +Three difficulty **grades**, ten cases each (30 total), spanning the monitoring +applications so a grade covers a realistic range rather than one scenario: + +- **easy** -- a pure seasonal signal at high SNR, single channel. The recommended + (best-practice) config should recover it cleanly. +- **medium** -- a transient (coseismic-style) drop with logarithmic, only partial + recovery, plus more measurement noise (lower SNR). +- **hard** -- a *multi-channel*, *depth- and frequency-dependent* problem. A + shallow (high-frequency) layer carries a coseismic drop-and-heal plus a full + hydrological seasonal cycle; a deep (low-frequency) layer carries a long-term + trend. Each case targets one depth, so the measurement band must match it: the + band selects the depth. Low SNR with waveform decorrelation; the channels are + measured independently and aggregated to a network dv/v. + +Every case has an exactly known ground-truth dv/v(t) (imposed by stretching a +band-limited decaying coda in lapse time), so any departure of a recovered series +is an artefact of the processing, not of nature. + +Design: the committed artefact is ``tests/data/golden/manifest.json`` -- the +recipes plus expected metrics, not the arrays. A multi-year (multi-channel) CCF +stack is large and fully determined by its seed, so arrays are regenerated on +demand and cached under ``tests/data/golden/cache/`` (gitignored). + +Consumers: :mod:`tests.test_golden` (regression oracle), the ``codameter-advisor`` +skill (live validation), and the FrugalMind ``dvv_processing`` suites +(:mod:`codameter.frugalmind`). All score through :func:`recover`, so a single +code path handles both single- and multi-channel cases. """ from __future__ import annotations import json +import warnings from pathlib import Path import numpy as np from . import use_cases as uc from .deviations import run_pipeline -from .synthetic_demo import ( - YEAR_D, _days, _seasonal, add_clock_drift, add_seasonal_late_noise, - daily_ccfs, earthquake_truth, groundwater_truth, landslide_truth, - make_coda, volcano_truth, -) +from .synthetic_demo import YEAR_D, _days, _seasonal, daily_ccfs, make_coda DATA_DIR = Path(__file__).resolve().parents[2] / "tests" / "data" / "golden" MANIFEST = DATA_DIR / "manifest.json" CACHE_DIR = DATA_DIR / "cache" -MANIFEST_VERSION = 1 +MANIFEST_VERSION = 2 # --------------------------------------------------------------------------- -# Ground-truth generators, resolved by name from use_cases[...]["dvv"]. -# Volcano/earthquake/landslide/groundwater reuse synthetic_demo; cryosphere and -# geothermal get local builders (adjusted amplitude) until dedicated ones exist. +# Per-application signal amplitudes (fractional dv/v) and timescales. These set +# the physical scale of each motif; they stay within the application's eps_max. # --------------------------------------------------------------------------- -def _cryosphere_truth(days: np.ndarray) -> np.ndarray: - """Sharp seasonal freeze-thaw swing (~ +/-3.5 %), summer velocity minimum.""" - return _seasonal(days, 0.035, 200.0) +AMP = { + "volcano": {"seasonal": 0.0010, "drop": -0.0040, "trend": -0.0015, "tau": 90.0, "phase": 60.0}, + "earthquake_fault": {"seasonal": 0.0006, "drop": -0.0025, "trend": -0.0010, "tau": 160.0, "phase": 30.0}, + "landslide": {"seasonal": 0.0100, "drop": -0.0300, "trend": -0.0050, "tau": 60.0, "phase": 120.0}, + "groundwater": {"seasonal": 0.0015, "drop": -0.0020, "trend": -0.0012, "tau": 120.0, "phase": 250.0}, + "cryosphere": {"seasonal": 0.0300, "drop": -0.0150, "trend": -0.0040, "tau": 45.0, "phase": 200.0}, + "geothermal": {"seasonal": 0.0005, "drop": -0.0060, "trend": -0.0100, "tau": 120.0, "phase": 30.0}, +} + +def _step_heal(days: np.ndarray, app: str, onset_frac: float = 0.5) -> np.ndarray: + """A sharp drop at ~``onset_frac`` of the record with logarithmic partial heal.""" + a = AMP[app] + onset = onset_frac * float(days[-1]) + co = days >= onset + dt = (days[co] - onset).astype(float) + out = np.zeros(len(days), float) + out[co] = a["drop"] * (0.35 + 0.65 * np.exp(-dt / a["tau"])) + return out -def _geothermal_truth(days: np.ndarray) -> np.ndarray: - """Slow injection-driven decline (~ -1 %) with a muted seasonal overlay.""" - ramp = -0.010 * np.clip((days - 0.3 * YEAR_D) / (2.0 * YEAR_D), 0, 1) - return ramp + _seasonal(days, 0.0004, 30.0) +def _motif_seasonal(days: np.ndarray, app: str) -> np.ndarray: + return _seasonal(days, AMP[app]["seasonal"], AMP[app]["phase"]) -TRUTH = { - "volcano": volcano_truth, - "earthquake": earthquake_truth, - "landslide": landslide_truth, - "groundwater_shallow": lambda days: groundwater_truth(days)[0], - "groundwater_deep": lambda days: groundwater_truth(days)[1], - "cryosphere": _cryosphere_truth, - "geothermal": _geothermal_truth, -} + +def _motif_transient(days: np.ndarray, app: str) -> np.ndarray: + # drop + heal, plus a muted seasonal so it is not unrealistically flat. + return _step_heal(days, app) + 0.25 * _motif_seasonal(days, app) + + +def _motif_composite(days: np.ndarray, app: str) -> np.ndarray: + # transient + full hydrological seasonal + a long-term linear trend. + a = AMP[app] + trend = a["trend"] * np.clip((days - 0.15 * days[-1]) / (0.8 * days[-1]), 0, 1) + return _step_heal(days, app) + _motif_seasonal(days, app) + trend + + +MOTIF = {"seasonal": _motif_seasonal, "transient": _motif_transient, + "composite": _motif_composite} + + +# --------------------------------------------------------------------------- +# Depth- and frequency-dependent medium (the hard grade). The shallow layer +# carries the near-surface response (coseismic drop + heal + full hydrological +# seasonal); the deep layer carries the long-term trend with a muted seasonal. +# The two layers live in separated frequency bands, so the measurement *band +# selects the depth*: a shallow (high-frequency) band recovers the shallow truth, +# a deep (low-frequency) band recovers the deep truth. +# --------------------------------------------------------------------------- +def _truth_shallow(days: np.ndarray, app: str) -> np.ndarray: + return _step_heal(days, app) + _motif_seasonal(days, app) + + +def _truth_deep(days: np.ndarray, app: str) -> np.ndarray: + a = AMP[app] + trend = a["trend"] * np.clip((days - 0.15 * days[-1]) / (0.8 * days[-1]), 0, 1) + return trend + 0.3 * _motif_seasonal(days, app) + + +def _depth_bands(app: str) -> tuple[tuple[float, float], tuple[float, float]]: + """Split an application's coda band into separated shallow (high) / deep (low) + sub-bands, each safely inside the generated content.""" + glo, ghi = uc.synth_params(app)["gen_band"] + gm = (glo * ghi) ** 0.5 + shallow = (round(gm * 1.6, 3), round(ghi * 0.9, 3)) + deep = (round(glo * 1.1, 3), round(gm * 0.6, 3)) + return shallow, deep # --------------------------------------------------------------------------- -# Case recipes. Each is a plain dict (no computed metrics); the manifest adds -# the expected metrics at regeneration time. -# -# id : unique case identifier / cache key. -# kind : "mainstream" | "edge". -# use_case : USE_CASES key; sets config, synth geometry, eps_max, truth. -# years, snr, seed : synthesis controls. -# decorr : waveform-decorrelation fraction (daily_ccfs). -# cadence : keep every Nth day (sparse sampling); default 1. -# artifacts : list of injectors applied after daily_ccfs. -# config : optional per-axis overrides on the recommended config. -# probes : optional [{label, config-overrides, truth}] extra measurements -# whose RMS is also frozen (used to prove band-selects-depth). -# rms_rel_tol : allowed fractional drift of RMS around the frozen value. +# Grades and case construction. Each grade cycles through the applications so it +# spans volcano / fault / aquifer / glacier / reservoir at that difficulty. # --------------------------------------------------------------------------- -CASES: list[dict] = [ - # ---- mainstream: one per application ------------------------------- - {"id": "volcano_mainstream", "kind": "mainstream", "use_case": "volcano", - "years": 3.0, "snr": 7.0, "seed": 11, "rms_rel_tol": 0.30}, - {"id": "earthquake_mainstream", "kind": "mainstream", "use_case": "earthquake_fault", - "years": 3.0, "snr": 7.0, "seed": 12, "rms_rel_tol": 0.30}, - {"id": "landslide_mainstream", "kind": "mainstream", "use_case": "landslide", - "years": 3.0, "snr": 8.0, "seed": 13, "rms_rel_tol": 0.30}, - {"id": "groundwater_mainstream", "kind": "mainstream", "use_case": "groundwater", - "years": 3.0, "snr": 8.0, "seed": 14, "rms_rel_tol": 0.30}, - {"id": "cryosphere_mainstream", "kind": "mainstream", "use_case": "cryosphere", - "years": 3.0, "snr": 8.0, "seed": 15, "rms_rel_tol": 0.30}, - {"id": "geothermal_mainstream", "kind": "mainstream", "use_case": "geothermal", - "years": 3.0, "snr": 7.0, "seed": 16, "rms_rel_tol": 0.30}, - - # ---- edge: low SNR + large (cycle-skipping) dv/v ------------------- - {"id": "low_snr_large_dvv", "kind": "edge", "use_case": "landslide", - "years": 3.0, "snr": 2.0, "seed": 21, "rms_rel_tol": 0.35, - "note": "SNR~2 with a several-percent pre-failure drop: the regime that " - "splits stretching from cross-spectral methods."}, - - # ---- edge: clock drift + seasonal late-coda noise ------------------ - {"id": "clock_drift_seasonal", "kind": "edge", "use_case": "volcano", - "years": 3.0, "snr": 7.0, "seed": 22, "rms_rel_tol": 0.35, - "artifacts": [ - {"kind": "clock_drift", "drift_s_per_day": 2.0e-4, "onset_day": 200}, - {"kind": "seasonal_late_noise", "onset_s": 20.0, "dvv_amp": 0.004, - "jitter": 0.06}, - ], - "note": "A growing clock error plus a seasonal late-coda warp inject a " - "spurious dv/v (Zhan 2013 / Daskalakis 2016)."}, - - # ---- edge: frequency-dependent shallow + deep medium --------------- - {"id": "freqdep_shallow_deep", "kind": "edge", "use_case": "groundwater", - "years": 3.0, "snr": 9.0, "seed": 23, "rms_rel_tol": 0.35, - "two_layer": True, - "config": {"band": (4.0, 10.0), "window": (2.0, 8.0)}, - "probes": [ - {"label": "deep_band_recovers_deep", - "config": {"band": (0.2, 1.0), "window": (8.0, 25.0)}, - "truth": "deep"}, - ], - "note": "Shallow (high-freq) and deep (low-freq) layers carry different " - "dv/v; the band selects which one you recover."}, - - # ---- edge: sparse cadence + waveform decorrelation ----------------- - {"id": "sparse_decorr", "kind": "edge", "use_case": "volcano", - "years": 3.0, "snr": 6.0, "seed": 24, "cadence": 3, "decorr": 0.30, - "rms_rel_tol": 0.35, - "note": "Every-third-day sampling with 30 % waveform decorrelation stresses " - "the reference/stacking warm-up."}, -] +GRADES = { + "easy": {"motif": "seasonal", "snr": (8.0, 12.0), "channels": 1, "decorr": 0.00, + "years": 3.0, "split": "validation", "rms_rel_tol": 0.35}, + "medium": {"motif": "transient", "snr": (3.0, 5.0), "channels": 1, "decorr": 0.05, + "years": 3.0, "split": "validation", "rms_rel_tol": 0.45}, + "hard": {"motif": "depth", "snr": (2.0, 4.0), "channels": 4, "decorr": 0.20, + "years": 2.5, "split": "test", "rms_rel_tol": 0.60}, +} + +# 10 application slots per grade (the 6 applications, some repeated). +APP_CYCLE = ["volcano", "earthquake_fault", "landslide", "groundwater", + "cryosphere", "geothermal", "volcano", "groundwater", + "landslide", "earthquake_fault"] + +_SEED_BASE = {"easy": 100, "medium": 200, "hard": 300} + +def _build_cases() -> list[dict]: + cases = [] + for grade, spec in GRADES.items(): + snr_lo, snr_hi = spec["snr"] + for i, app in enumerate(APP_CYCLE): + snr = round(float(np.interp(i, [0, len(APP_CYCLE) - 1], [snr_hi, snr_lo])), 2) + case = { + "id": f"{grade}-{app}-{i + 1:02d}", + "grade": grade, "use_case": app, "motif": spec["motif"], + "snr": snr, "seed": _SEED_BASE[grade] + i, + "channels": spec["channels"], "decorr": spec["decorr"], + "years": spec["years"], "split": spec["split"], + "rms_rel_tol": spec["rms_rel_tol"], + } + if grade == "hard": + # Depth- and frequency-dependent: the case targets one depth, and + # the recommended config's band must match it. Targets alternate. + shallow_band, deep_band = _depth_bands(app) + target = "shallow" if i % 2 == 0 else "deep" + case["two_layer"] = True + case["target"] = target + case["config"] = {"band": shallow_band if target == "shallow" else deep_band} + case["note"] = ( + "The medium is depth-dependent: a shallow near-surface layer " + "(coseismic drop + hydrological seasonal) sits above a deep layer " + "(long-term trend). " + + ("You must resolve the SHALLOW near-surface response." + if target == "shallow" else + "You must resolve the DEEP long-term trend.") + + " The band selects the depth." + ) + cases.append(case) + return cases + + +CASES: list[dict] = _build_cases() CASES_BY_ID = {c["id"]: c for c in CASES} -# Application -> its mainstream golden case, so the advisor can pull a matched -# synthetic for any use case without hardcoding ids. + +def representative_case(use_case: str, grade: str = "easy") -> str: + """A matched case id for an application (default the first easy case). + + Used by the advisor and the sweep to pull a synthetic that matches a user's + application without hardcoding ids. + """ + key = uc.resolve(use_case) + for c in CASES: + if c["use_case"] == key and c["grade"] == grade: + return c["id"] + for c in CASES: # fall back to any grade + if c["use_case"] == key: + return c["id"] + raise KeyError(f"no golden case for use case {use_case!r}") + + +# Back-compat alias: application -> a representative (easy) case. MAINSTREAM_BY_USE_CASE = { - c["use_case"]: c["id"] for c in CASES if c["kind"] == "mainstream" + app: representative_case(app) for app in {c["use_case"] for c in CASES} } def case_split(recipe: dict) -> str: - """Benchmark split for a case: mainstream -> validation, edge -> test. - - (These are the frugalmind splits; a recipe may override with a ``split`` key.) - """ - return recipe.get("split", "validation" if recipe["kind"] == "mainstream" else "test") + return recipe.get("split", "validation") def case_visibility(recipe: dict) -> str: - """Benchmark visibility. Synthetic and seed-reproducible, so public by default.""" + """Synthetic and seed-reproducible, so public by default.""" return recipe.get("visibility", "public") # --------------------------------------------------------------------------- # Synthesis # --------------------------------------------------------------------------- -def _synth_reference(use_case: str): - """One reference coda + its lapse axis, matched to the use case geometry.""" - sp = uc.synth_params(use_case) - t, ref = make_coda(maxlag_s=sp["maxlag_s"], fs=sp["fs"], - band=sp["gen_band"], t_coda_s=sp["t_coda_s"], seed=0) - return t, ref, sp - - def _build(recipe: dict) -> dict: """Regenerate the arrays for one case deterministically from its recipe.""" - use_case = recipe["use_case"] - sp = uc.synth_params(use_case) - fs, gen_band = sp["fs"], sp["gen_band"] + app = recipe["use_case"] + sp = uc.synth_params(app) + fs, gen = sp["fs"], sp["gen_band"] days = _days(recipe["years"]) - seed = recipe["seed"] - snr = recipe["snr"] + seed, snr = recipe["seed"], recipe["snr"] decorr = recipe.get("decorr", 0.0) + nchan = int(recipe.get("channels", 1)) - out: dict = {"fs": fs, "days": days, "use_case": use_case} + t, coda0 = make_coda(maxlag_s=sp["maxlag_s"], fs=fs, band=gen, + t_coda_s=sp["t_coda_s"], seed=0) + out: dict = {"fs": fs, "days": days, "use_case": app, "grade": recipe["grade"], + "t": t} if recipe.get("two_layer"): - # Two band-separated layers so the frequency band selects the depth. - shallow_band, deep_band = (3.0, 11.0), (0.2, 1.1) - t, shallow = make_coda(maxlag_s=sp["maxlag_s"], fs=fs, band=shallow_band, - t_coda_s=sp["t_coda_s"], seed=0) - _, deep = make_coda(maxlag_s=sp["maxlag_s"], fs=fs, band=deep_band, - t_coda_s=sp["t_coda_s"], seed=1) - truth_s = TRUTH["groundwater_shallow"](days) - truth_d = TRUTH["groundwater_deep"](days) - # Noise (and any decorrelation coda) must be broadband across BOTH layer - # bands, or the SNR would be frequency-dependent and confound the - # band-selects-depth test. Span deep_band low to shallow_band high. - two_layer_band = (deep_band[0] * 0.5, shallow_band[1] * 1.1) - ccfs = daily_ccfs(t, [shallow, deep], [truth_s, truth_d], fs=fs, - snr=snr, decorr=decorr, gen_band=two_layer_band, seed=seed) - out.update(t=t, ccfs=ccfs, truth=truth_s, truth_deep=truth_d) + # Depth-dependent medium: shallow (high-freq) and deep (low-freq) layers, + # each carrying its own dv/v. Every channel sums both layers; the band + # (set by the recommended config) selects which depth is recovered. The + # scored truth is the targeted layer. + shallow_band, deep_band = _depth_bands(app) + truth_shallow = _truth_shallow(days, app) + truth_deep = _truth_deep(days, app) + chans = [] + for c in range(nchan): + _, cod_s = make_coda(maxlag_s=sp["maxlag_s"], fs=fs, band=shallow_band, + t_coda_s=sp["t_coda_s"], seed=2 * c) + _, cod_d = make_coda(maxlag_s=sp["maxlag_s"], fs=fs, band=deep_band, + t_coda_s=sp["t_coda_s"], seed=2 * c + 1) + chans.append(daily_ccfs(t, [cod_s, cod_d], [truth_shallow, truth_deep], + fs=fs, snr=snr, decorr=decorr, gen_band=gen, + seed=seed + 7 * c)) + out["channels"] = np.stack(chans) + out["ccfs"] = out["channels"].mean(axis=0) + out["truth"] = truth_shallow if recipe["target"] == "shallow" else truth_deep + out["truth_other"] = truth_deep if recipe["target"] == "shallow" else truth_shallow + return out + + truth = MOTIF[recipe["motif"]](days, app) + out["truth"] = truth + if nchan > 1: + # Independent cross-component channels: distinct coda + distinct noise, + # sharing the medium's truth. Measured per channel, aggregated later. + codas = [coda0] + [make_coda(maxlag_s=sp["maxlag_s"], fs=fs, band=gen, + t_coda_s=sp["t_coda_s"], seed=c)[1] + for c in range(1, nchan)] + chans = [daily_ccfs(t, [cod], [truth], fs=fs, snr=snr, decorr=decorr, + gen_band=gen, seed=seed + 7 * c) + for c, cod in enumerate(codas)] + out["channels"] = np.stack(chans) + out["ccfs"] = out["channels"].mean(axis=0) # a 2D view for plotting else: - t, ref, _ = _synth_reference(use_case) - truth = TRUTH[uc.USE_CASES[use_case]["dvv"]](days) - ccfs = daily_ccfs(t, [ref], [truth], fs=fs, snr=snr, decorr=decorr, - gen_band=gen_band, seed=seed) - out.update(t=t, ccfs=ccfs, truth=truth) - - for art in recipe.get("artifacts", []): - ccfs = out["ccfs"] - if art["kind"] == "clock_drift": - ccfs = add_clock_drift(ccfs, out["t"], - drift_s_per_day=art["drift_s_per_day"], - onset_day=art.get("onset_day", 0)) - elif art["kind"] == "seasonal_late_noise": - ccfs = add_seasonal_late_noise( - ccfs, out["t"], out["days"], fs=fs, onset_s=art["onset_s"], - dvv_amp=art.get("dvv_amp", 0.004), jitter=art.get("jitter", 0.06), - band=gen_band, seed=seed + 100) - else: - raise ValueError(f"unknown artifact {art['kind']!r}") - out["ccfs"] = ccfs - - # Sparse sampling last, so injectors see the full daily record. - cadence = recipe.get("cadence", 1) - if cadence > 1: - idx = np.arange(0, len(out["days"]), cadence) - out["days"] = out["days"][idx] - out["ccfs"] = out["ccfs"][idx] - out["truth"] = out["truth"][idx] - if "truth_deep" in out: - out["truth_deep"] = out["truth_deep"][idx] + out["ccfs"] = daily_ccfs(t, [coda0], [truth], fs=fs, snr=snr, + decorr=decorr, gen_band=gen, seed=seed) return out +def recover(d: dict, cfg: dict, eps_max: float): + """Recover dv/v(t) for a case under ``cfg``; the single scoring entry point. + + Single-channel cases run the pipeline directly. Multi-channel cases run the + pipeline on each channel and aggregate to a network dv/v by averaging the + per-channel series (the "average the per-component dv/v" convention). + """ + if "channels" in d and np.ndim(d["channels"]) == 3: + per = [] + for c in range(d["channels"].shape[0]): + dvv_c, val_c = run_pipeline(d["channels"][c], d["t"], d["fs"], cfg, + eps_max=eps_max) + per.append(np.where(val_c, dvv_c, np.nan)) + with warnings.catch_warnings(): + warnings.simplefilter("ignore", RuntimeWarning) # all-NaN columns + dvv = np.nanmean(np.vstack(per), axis=0) + return dvv, np.isfinite(dvv) + return run_pipeline(d["ccfs"], d["t"], d["fs"], cfg, eps_max=eps_max) + + +# --------------------------------------------------------------------------- +# Cache +# --------------------------------------------------------------------------- def _recipe_hash(recipe: dict) -> str: """Short digest of a recipe, so editing it invalidates the cached arrays.""" import hashlib @@ -247,13 +302,13 @@ def _recipe_hash(recipe: dict) -> str: def generate(case_id: str, *, cache: bool = True) -> dict: - """Return ``{ccfs, t, days, truth[, truth_deep], fs, use_case}`` for a case. + """Return the arrays for a case: ``{ccfs, t, days, truth, fs, use_case, grade}`` + plus ``channels`` (3D) for multi-channel cases. - Deterministic in the recipe seed. When ``cache`` is set the arrays are read - from / written to ``tests/data/golden/cache/-.npz``. The - recipe hash in the filename means a changed recipe misses the stale cache and - rebuilds; a change to the synthesis *code* is not captured by the hash, so - ``regenerate_manifest`` bypasses the cache entirely (``cache=False``). + Deterministic in the seed. Cached to ``cache/-.npz`` (the hash + busts the cache when a recipe changes). ``regenerate_manifest`` uses + ``cache=False`` so a synthesis-*code* change (not captured by the hash) never + scores against stale arrays. """ recipe = CASES_BY_ID[case_id] cache_file = CACHE_DIR / f"{case_id}-{_recipe_hash(recipe)}.npz" @@ -262,15 +317,23 @@ def generate(case_id: str, *, cache: bool = True) -> dict: d = {k: z[k] for k in z.files} d["fs"] = float(d["fs"]) d["use_case"] = str(d["use_case"]) + d["grade"] = str(d["grade"]) return d d = _build(recipe) if cache: CACHE_DIR.mkdir(parents=True, exist_ok=True) - np.savez_compressed( - cache_file, - **{k: np.asarray(v) for k, v in d.items() - if k in ("ccfs", "t", "days", "truth", "truth_deep")}, - fs=np.asarray(d["fs"]), use_case=np.asarray(d["use_case"])) + # Downcast the large CCF arrays to float32 to keep the cache small; the + # lapse/day/truth axes stay float64. + payload = {"t": d["t"], "days": d["days"], "truth": d["truth"], + "ccfs": np.asarray(d["ccfs"], np.float32), + "fs": np.asarray(d["fs"]), + "use_case": np.asarray(d["use_case"]), + "grade": np.asarray(d["grade"])} + if "channels" in d: + payload["channels"] = np.asarray(d["channels"], np.float32) + if "truth_other" in d: + payload["truth_other"] = d["truth_other"] + np.savez_compressed(cache_file, **payload) return d @@ -278,14 +341,12 @@ def generate(case_id: str, *, cache: bool = True) -> dict: # Metrics oracle # --------------------------------------------------------------------------- def _rms(dvv, truth, days, valid, baseline_frac: float = 0.2) -> float: - """Baseline-aligned RMS error of a recovered dv/v series against the truth. - - A dv/v estimate is relative to a reference epoch, so its absolute level (the - DC offset) is not observable: a fixed reference measures change *since the - baseline window*, not since zero. We therefore remove, from both the - recovered and the true series, their mean over the earliest ``baseline_frac`` - of valid epochs -- the reference epoch -- before taking the RMS. Without this - a pure monotonic trend would show a spurious error equal to the trend's mean. + """Baseline-aligned RMS error of a recovered dv/v against the truth. + + A dv/v estimate is relative to a reference epoch, so its DC offset is not + observable. We remove, from both series, their mean over the earliest + ``baseline_frac`` of valid epochs before taking the RMS, or a pure trend would + show a spurious error equal to its mean. """ v = valid & np.isfinite(dvv) if v.sum() < 10: @@ -300,36 +361,29 @@ def _rms(dvv, truth, days, valid, baseline_frac: float = 0.2) -> float: return float(np.sqrt(np.mean((d0 - tr0) ** 2))) +def _jsonable(cfg: dict) -> dict: + """Tuples (band, window) -> lists so the config round-trips through JSON.""" + return {k: (list(v) if isinstance(v, tuple) else v) for k, v in cfg.items()} + + def compute_metrics(case_id: str, data: dict | None = None) -> dict: - """Run the recommended pipeline (+ any probes) and return the RMS oracle.""" + """Recover with the recommended config and return the RMS oracle.""" recipe = CASES_BY_ID[case_id] d = data if data is not None else generate(case_id) - use_case = recipe["use_case"] - cfg = uc.recommend(use_case, **recipe.get("config", {})) - eps = uc.eps_max(use_case) - - dvv, valid = run_pipeline(d["ccfs"], d["t"], d["fs"], cfg, eps_max=eps) + app = recipe["use_case"] + cfg = uc.recommend(app, **recipe.get("config", {})) + eps = uc.eps_max(app) + dvv, valid = recover(d, cfg, eps) res = {"config": _jsonable(cfg), "eps_max": eps, "rms": _rms(dvv, d["truth"], d["days"], valid)} - - probes = [] - for p in recipe.get("probes", []): - pcfg = uc.recommend(use_case, **p.get("config", {})) - truth = d["truth_deep"] if p.get("truth") == "deep" else d["truth"] - pdvv, pvalid = run_pipeline(d["ccfs"], d["t"], d["fs"], pcfg, eps_max=eps) - probes.append({"label": p["label"], "config": _jsonable(pcfg), - "truth": p.get("truth", "shallow"), - "rms": _rms(pdvv, truth, d["days"], pvalid)}) - if probes: - res["probes"] = probes + if "truth_other" in d: + # The error a config would incur by recovering the WRONG depth layer: + # the "clearly wrong" anchor for scoring depth-band selection. + allv = np.ones(len(d["days"]), bool) + res["rms_wrong_layer"] = _rms(d["truth_other"], d["truth"], d["days"], allv) return res -def _jsonable(cfg: dict) -> dict: - """Tuples (band, window) -> lists so the config round-trips through JSON.""" - return {k: (list(v) if isinstance(v, tuple) else v) for k, v in cfg.items()} - - # --------------------------------------------------------------------------- # Manifest read / write # --------------------------------------------------------------------------- @@ -337,31 +391,24 @@ def regenerate_manifest() -> dict: """Recompute every case's expected metrics and rewrite ``manifest.json``.""" cases = [] for c in CASES: - # Always rebuild from current code, never from a possibly stale cache. - d = generate(c["id"], cache=False) + d = generate(c["id"], cache=False) # always from current code m = compute_metrics(c["id"], d) entry = { - "id": c["id"], "kind": c["kind"], "use_case": c["use_case"], - "split": case_split(c), "visibility": case_visibility(c), - "years": c["years"], "snr": c["snr"], "seed": c["seed"], - "cadence": c.get("cadence", 1), "decorr": c.get("decorr", 0.0), + "id": c["id"], "grade": c["grade"], "use_case": c["use_case"], + "motif": c["motif"], "split": case_split(c), + "visibility": case_visibility(c), "years": c["years"], "snr": c["snr"], + "seed": c["seed"], "channels": c["channels"], "decorr": c["decorr"], "rms_rel_tol": c["rms_rel_tol"], "n_days": int(len(d["days"])), "expected": m, } - # Record the remaining recipe fields so the oracle is auditable from the - # manifest alone (not only by reading golden.CASES). if c.get("two_layer"): entry["two_layer"] = True - if "artifacts" in c: - entry["artifacts"] = c["artifacts"] - if "config" in c: - entry["config"] = _jsonable(c["config"]) - if "note" in c: - entry["note"] = c["note"] + entry["target"] = c["target"] cases.append(entry) - print(f" {c['id']:<24} rms={m['rms']:.5f}" - + (f" (+{len(m['probes'])} probe)" if "probes" in m else "")) - manifest = {"version": MANIFEST_VERSION, "cases": cases} + print(f" {c['id']:<26} ch={c['channels']} snr={c['snr']:<4} " + f"rms={m['rms']:.5f}") + manifest = {"version": MANIFEST_VERSION, "grades": list(GRADES), + "cases": cases} DATA_DIR.mkdir(parents=True, exist_ok=True) MANIFEST.write_text(json.dumps(manifest, indent=2) + "\n") return manifest @@ -379,7 +426,8 @@ def expected_metrics(case_id: str) -> dict: def main() -> int: - print(f"Regenerating golden manifest ({len(CASES)} cases) -> {MANIFEST}") + print(f"Regenerating golden manifest ({len(CASES)} cases, " + f"{len(GRADES)} grades) -> {MANIFEST}") regenerate_manifest() print("done.") return 0 diff --git a/tests/data/golden/.gitignore b/tests/data/golden/.gitignore index 07aabdd..5247a43 100644 --- a/tests/data/golden/.gitignore +++ b/tests/data/golden/.gitignore @@ -1,3 +1,6 @@ # Regenerable synthetic CCF arrays (deterministic from manifest seeds). # The committed oracle is manifest.json; the arrays are rebuilt on demand. cache/ + +# Verification figures (regenerate: python scripts/plot_golden.py) +figs/ diff --git a/tests/data/golden/manifest.json b/tests/data/golden/manifest.json index 1d784ad..db02222 100644 --- a/tests/data/golden/manifest.json +++ b/tests/data/golden/manifest.json @@ -1,18 +1,24 @@ { - "version": 1, + "version": 2, + "grades": [ + "easy", + "medium", + "hard" + ], "cases": [ { - "id": "volcano_mainstream", - "kind": "mainstream", + "id": "easy-volcano-01", + "grade": "easy", "use_case": "volcano", + "motif": "seasonal", "split": "validation", "visibility": "public", "years": 3.0, - "snr": 7.0, - "seed": 11, - "cadence": 1, + "snr": 12.0, + "seed": 100, + "channels": 1, "decorr": 0.0, - "rms_rel_tol": 0.3, + "rms_rel_tol": 0.35, "n_days": 1095, "expected": { "config": { @@ -30,21 +36,22 @@ "gate": true }, "eps_max": 0.06, - "rms": 0.00024045031681258644 + "rms": 8.528385265033053e-05 } }, { - "id": "earthquake_mainstream", - "kind": "mainstream", + "id": "easy-earthquake_fault-02", + "grade": "easy", "use_case": "earthquake_fault", + "motif": "seasonal", "split": "validation", "visibility": "public", "years": 3.0, - "snr": 7.0, - "seed": 12, - "cadence": 1, + "snr": 11.56, + "seed": 101, + "channels": 1, "decorr": 0.0, - "rms_rel_tol": 0.3, + "rms_rel_tol": 0.35, "n_days": 1095, "expected": { "config": { @@ -62,21 +69,22 @@ "gate": true }, "eps_max": 0.05, - "rms": 0.00014156185178993663 + "rms": 5.127319326441781e-05 } }, { - "id": "landslide_mainstream", - "kind": "mainstream", + "id": "easy-landslide-03", + "grade": "easy", "use_case": "landslide", + "motif": "seasonal", "split": "validation", "visibility": "public", "years": 3.0, - "snr": 8.0, - "seed": 13, - "cadence": 1, + "snr": 11.11, + "seed": 102, + "channels": 1, "decorr": 0.0, - "rms_rel_tol": 0.3, + "rms_rel_tol": 0.35, "n_days": 1095, "expected": { "config": { @@ -94,21 +102,22 @@ "gate": true }, "eps_max": 0.09, - "rms": 0.0006342923685591622 + "rms": 0.00034448337441832434 } }, { - "id": "groundwater_mainstream", - "kind": "mainstream", + "id": "easy-groundwater-04", + "grade": "easy", "use_case": "groundwater", + "motif": "seasonal", "split": "validation", "visibility": "public", "years": 3.0, - "snr": 8.0, - "seed": 14, - "cadence": 1, + "snr": 10.67, + "seed": 103, + "channels": 1, "decorr": 0.0, - "rms_rel_tol": 0.3, + "rms_rel_tol": 0.35, "n_days": 1095, "expected": { "config": { @@ -126,21 +135,22 @@ "gate": true }, "eps_max": 0.03, - "rms": 0.00011149748859488201 + "rms": 9.91671255794741e-05 } }, { - "id": "cryosphere_mainstream", - "kind": "mainstream", + "id": "easy-cryosphere-05", + "grade": "easy", "use_case": "cryosphere", + "motif": "seasonal", "split": "validation", "visibility": "public", "years": 3.0, - "snr": 8.0, - "seed": 15, - "cadence": 1, + "snr": 10.22, + "seed": 104, + "channels": 1, "decorr": 0.0, - "rms_rel_tol": 0.3, + "rms_rel_tol": 0.35, "n_days": 1095, "expected": { "config": { @@ -158,21 +168,22 @@ "gate": true }, "eps_max": 0.11, - "rms": 0.0013733676670769366 + "rms": 0.0013108903767198416 } }, { - "id": "geothermal_mainstream", - "kind": "mainstream", + "id": "easy-geothermal-06", + "grade": "easy", "use_case": "geothermal", + "motif": "seasonal", "split": "validation", "visibility": "public", "years": 3.0, - "snr": 7.0, - "seed": 16, - "cadence": 1, + "snr": 9.78, + "seed": 105, + "channels": 1, "decorr": 0.0, - "rms_rel_tol": 0.3, + "rms_rel_tol": 0.35, "n_days": 1095, "expected": { "config": { @@ -190,19 +201,86 @@ "gate": true }, "eps_max": 0.03, - "rms": 7.499932473111218e-05 + "rms": 5.4426102327806585e-05 + } + }, + { + "id": "easy-volcano-07", + "grade": "easy", + "use_case": "volcano", + "motif": "seasonal", + "split": "validation", + "visibility": "public", + "years": 3.0, + "snr": 9.33, + "seed": 106, + "channels": 1, + "decorr": 0.0, + "rms_rel_tol": 0.35, + "n_days": 1095, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 0.4, + 1.0 + ], + "window": [ + 10, + 30 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.06, + "rms": 0.0001221091620658341 + } + }, + { + "id": "easy-groundwater-08", + "grade": "easy", + "use_case": "groundwater", + "motif": "seasonal", + "split": "validation", + "visibility": "public", + "years": 3.0, + "snr": 8.89, + "seed": 107, + "channels": 1, + "decorr": 0.0, + "rms_rel_tol": 0.35, + "n_days": 1095, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 2.0, + 4.0 + ], + "window": [ + 2.0, + 8.0 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.03, + "rms": 0.00010207171587653883 } }, { - "id": "low_snr_large_dvv", - "kind": "edge", + "id": "easy-landslide-09", + "grade": "easy", "use_case": "landslide", - "split": "test", + "motif": "seasonal", + "split": "validation", "visibility": "public", "years": 3.0, - "snr": 2.0, - "seed": 21, - "cadence": 1, + "snr": 8.44, + "seed": 108, + "channels": 1, "decorr": 0.0, "rms_rel_tol": 0.35, "n_days": 1095, @@ -222,23 +300,56 @@ "gate": true }, "eps_max": 0.09, - "rms": 0.0008136899296025077 - }, - "note": "SNR~2 with a several-percent pre-failure drop: the regime that splits stretching from cross-spectral methods." + "rms": 0.000346672627405836 + } }, { - "id": "clock_drift_seasonal", - "kind": "edge", - "use_case": "volcano", - "split": "test", + "id": "easy-earthquake_fault-10", + "grade": "easy", + "use_case": "earthquake_fault", + "motif": "seasonal", + "split": "validation", "visibility": "public", "years": 3.0, - "snr": 7.0, - "seed": 22, - "cadence": 1, + "snr": 8.0, + "seed": 109, + "channels": 1, "decorr": 0.0, "rms_rel_tol": 0.35, "n_days": 1095, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 0.5, + 1.5 + ], + "window": [ + 8, + 25 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.05, + "rms": 6.906899962700634e-05 + } + }, + { + "id": "medium-volcano-01", + "grade": "medium", + "use_case": "volcano", + "motif": "transient", + "split": "validation", + "visibility": "public", + "years": 3.0, + "snr": 5.0, + "seed": 200, + "channels": 1, + "decorr": 0.05, + "rms_rel_tol": 0.45, + "n_days": 1095, "expected": { "config": { "estimator": "stretching (TS)", @@ -255,42 +366,434 @@ "gate": true }, "eps_max": 0.06, - "rms": 0.0009051145478551642 - }, - "artifacts": [ - { - "kind": "clock_drift", - "drift_s_per_day": 0.0002, - "onset_day": 200 + "rms": 0.00026258877639392437 + } + }, + { + "id": "medium-earthquake_fault-02", + "grade": "medium", + "use_case": "earthquake_fault", + "motif": "transient", + "split": "validation", + "visibility": "public", + "years": 3.0, + "snr": 4.78, + "seed": 201, + "channels": 1, + "decorr": 0.05, + "rms_rel_tol": 0.45, + "n_days": 1095, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 0.5, + 1.5 + ], + "window": [ + 8, + 25 + ], + "stack": 10, + "reference": "fixed", + "gate": true }, - { - "kind": "seasonal_late_noise", - "onset_s": 20.0, - "dvv_amp": 0.004, - "jitter": 0.06 - } - ], - "note": "A growing clock error plus a seasonal late-coda warp inject a spurious dv/v (Zhan 2013 / Daskalakis 2016)." + "eps_max": 0.05, + "rms": 0.00015733842087016845 + } }, { - "id": "freqdep_shallow_deep", - "kind": "edge", + "id": "medium-landslide-03", + "grade": "medium", + "use_case": "landslide", + "motif": "transient", + "split": "validation", + "visibility": "public", + "years": 3.0, + "snr": 4.56, + "seed": 202, + "channels": 1, + "decorr": 0.05, + "rms_rel_tol": 0.45, + "n_days": 1095, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 4.0, + 12.0 + ], + "window": [ + 0.2, + 1.5 + ], + "stack": 5, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.09, + "rms": 0.001107126031724569 + } + }, + { + "id": "medium-groundwater-04", + "grade": "medium", "use_case": "groundwater", - "split": "test", + "motif": "transient", + "split": "validation", "visibility": "public", "years": 3.0, - "snr": 9.0, - "seed": 23, - "cadence": 1, - "decorr": 0.0, - "rms_rel_tol": 0.35, + "snr": 4.33, + "seed": 203, + "channels": 1, + "decorr": 0.05, + "rms_rel_tol": 0.45, + "n_days": 1095, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 2.0, + 4.0 + ], + "window": [ + 2.0, + 8.0 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.03, + "rms": 0.0001525314605933785 + } + }, + { + "id": "medium-cryosphere-05", + "grade": "medium", + "use_case": "cryosphere", + "motif": "transient", + "split": "validation", + "visibility": "public", + "years": 3.0, + "snr": 4.11, + "seed": 204, + "channels": 1, + "decorr": 0.05, + "rms_rel_tol": 0.45, + "n_days": 1095, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 4.0, + 14.0 + ], + "window": [ + 0.3, + 0.8 + ], + "stack": 5, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.11, + "rms": 0.0007095977946747316 + } + }, + { + "id": "medium-geothermal-06", + "grade": "medium", + "use_case": "geothermal", + "motif": "transient", + "split": "validation", + "visibility": "public", + "years": 3.0, + "snr": 3.89, + "seed": 205, + "channels": 1, + "decorr": 0.05, + "rms_rel_tol": 0.45, + "n_days": 1095, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 0.5, + 2.0 + ], + "window": [ + 10.0, + 25.0 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.03, + "rms": 0.00032151232310283743 + } + }, + { + "id": "medium-volcano-07", + "grade": "medium", + "use_case": "volcano", + "motif": "transient", + "split": "validation", + "visibility": "public", + "years": 3.0, + "snr": 3.67, + "seed": 206, + "channels": 1, + "decorr": 0.05, + "rms_rel_tol": 0.45, + "n_days": 1095, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 0.4, + 1.0 + ], + "window": [ + 10, + 30 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.06, + "rms": 0.00032005473057429297 + } + }, + { + "id": "medium-groundwater-08", + "grade": "medium", + "use_case": "groundwater", + "motif": "transient", + "split": "validation", + "visibility": "public", + "years": 3.0, + "snr": 3.44, + "seed": 207, + "channels": 1, + "decorr": 0.05, + "rms_rel_tol": 0.45, + "n_days": 1095, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 2.0, + 4.0 + ], + "window": [ + 2.0, + 8.0 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.03, + "rms": 0.000168247967817991 + } + }, + { + "id": "medium-landslide-09", + "grade": "medium", + "use_case": "landslide", + "motif": "transient", + "split": "validation", + "visibility": "public", + "years": 3.0, + "snr": 3.22, + "seed": 208, + "channels": 1, + "decorr": 0.05, + "rms_rel_tol": 0.45, "n_days": 1095, "expected": { "config": { "estimator": "stretching (TS)", "band": [ 4.0, - 10.0 + 12.0 + ], + "window": [ + 0.2, + 1.5 + ], + "stack": 5, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.09, + "rms": 0.0011341164599683257 + } + }, + { + "id": "medium-earthquake_fault-10", + "grade": "medium", + "use_case": "earthquake_fault", + "motif": "transient", + "split": "validation", + "visibility": "public", + "years": 3.0, + "snr": 3.0, + "seed": 209, + "channels": 1, + "decorr": 0.05, + "rms_rel_tol": 0.45, + "n_days": 1095, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 0.5, + 1.5 + ], + "window": [ + 8, + 25 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.05, + "rms": 0.00020088776612501967 + } + }, + { + "id": "hard-volcano-01", + "grade": "hard", + "use_case": "volcano", + "motif": "depth", + "split": "test", + "visibility": "public", + "years": 2.5, + "snr": 4.0, + "seed": 300, + "channels": 4, + "decorr": 0.2, + "rms_rel_tol": 0.6, + "n_days": 913, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 1.131, + 9.0 + ], + "window": [ + 10, + 30 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.06, + "rms": 0.00025212395764457387, + "rms_wrong_layer": 0.000930391397212034 + }, + "two_layer": true, + "target": "shallow" + }, + { + "id": "hard-earthquake_fault-02", + "grade": "hard", + "use_case": "earthquake_fault", + "motif": "depth", + "split": "test", + "visibility": "public", + "years": 2.5, + "snr": 3.78, + "seed": 301, + "channels": 4, + "decorr": 0.2, + "rms_rel_tol": 0.6, + "n_days": 913, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 0.055, + 0.424 + ], + "window": [ + 8, + 25 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.05, + "rms": 0.00012877251474255628, + "rms_wrong_layer": 0.0007607270475765937 + }, + "two_layer": true, + "target": "deep" + }, + { + "id": "hard-landslide-03", + "grade": "hard", + "use_case": "landslide", + "motif": "depth", + "split": "test", + "visibility": "public", + "years": 2.5, + "snr": 3.56, + "seed": 302, + "channels": 4, + "decorr": 0.2, + "rms_rel_tol": 0.6, + "n_days": 913, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 7.155, + 18.0 + ], + "window": [ + 0.2, + 1.5 + ], + "stack": 5, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.09, + "rms": 0.0012734011077726465, + "rms_wrong_layer": 0.007530894565057449 + }, + "two_layer": true, + "target": "shallow" + }, + { + "id": "hard-groundwater-04", + "grade": "hard", + "use_case": "groundwater", + "motif": "depth", + "split": "test", + "visibility": "public", + "years": 2.5, + "snr": 3.33, + "seed": 303, + "channels": 4, + "decorr": 0.2, + "rms_rel_tol": 0.6, + "n_days": 913, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 0.55, + 1.2 ], "window": [ 2.0, @@ -301,61 +804,104 @@ "gate": true }, "eps_max": 0.03, - "rms": 0.00010868054317771904, - "probes": [ - { - "label": "deep_band_recovers_deep", - "config": { - "estimator": "stretching (TS)", - "band": [ - 0.2, - 1.0 - ], - "window": [ - 8.0, - 25.0 - ], - "stack": 10, - "reference": "fixed", - "gate": true - }, - "truth": "deep", - "rms": 0.0002856315920472662 - } - ] + "rms": 0.00022902404646362183, + "rms_wrong_layer": 0.0010007890743546153 }, "two_layer": true, - "config": { - "band": [ - 4.0, - 10.0 - ], - "window": [ - 2.0, - 8.0 - ] + "target": "deep" + }, + { + "id": "hard-cryosphere-05", + "grade": "hard", + "use_case": "cryosphere", + "motif": "depth", + "split": "test", + "visibility": "public", + "years": 2.5, + "snr": 3.11, + "seed": 304, + "channels": 4, + "decorr": 0.2, + "rms_rel_tol": 0.6, + "n_days": 913, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 11.314, + 22.5 + ], + "window": [ + 0.3, + 0.8 + ], + "stack": 5, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.11, + "rms": 0.0010750047910110425, + "rms_wrong_layer": 0.017753287172458793 }, - "note": "Shallow (high-freq) and deep (low-freq) layers carry different dv/v; the band selects which one you recover." + "two_layer": true, + "target": "shallow" }, { - "id": "sparse_decorr", - "kind": "edge", + "id": "hard-geothermal-06", + "grade": "hard", + "use_case": "geothermal", + "motif": "depth", + "split": "test", + "visibility": "public", + "years": 2.5, + "snr": 2.89, + "seed": 305, + "channels": 4, + "decorr": 0.2, + "rms_rel_tol": 0.6, + "n_days": 913, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 0.11, + 0.537 + ], + "window": [ + 10.0, + 25.0 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.03, + "rms": 0.00035681825048443534, + "rms_wrong_layer": 0.0037482773313763365 + }, + "two_layer": true, + "target": "deep" + }, + { + "id": "hard-volcano-07", + "grade": "hard", "use_case": "volcano", + "motif": "depth", "split": "test", "visibility": "public", - "years": 3.0, - "snr": 6.0, - "seed": 24, - "cadence": 3, - "decorr": 0.3, - "rms_rel_tol": 0.35, - "n_days": 365, + "years": 2.5, + "snr": 2.67, + "seed": 306, + "channels": 4, + "decorr": 0.2, + "rms_rel_tol": 0.6, + "n_days": 913, "expected": { "config": { "estimator": "stretching (TS)", "band": [ - 0.4, - 1.0 + 1.131, + 9.0 ], "window": [ 10, @@ -366,9 +912,119 @@ "gate": true }, "eps_max": 0.06, - "rms": 0.0004359819565566644 + "rms": 0.00025180197315354186, + "rms_wrong_layer": 0.000930391397212034 }, - "note": "Every-third-day sampling with 30 % waveform decorrelation stresses the reference/stacking warm-up." + "two_layer": true, + "target": "shallow" + }, + { + "id": "hard-groundwater-08", + "grade": "hard", + "use_case": "groundwater", + "motif": "depth", + "split": "test", + "visibility": "public", + "years": 2.5, + "snr": 2.44, + "seed": 307, + "channels": 4, + "decorr": 0.2, + "rms_rel_tol": 0.6, + "n_days": 913, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 0.55, + 1.2 + ], + "window": [ + 2.0, + 8.0 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.03, + "rms": 0.00031478238532623977, + "rms_wrong_layer": 0.0010007890743546153 + }, + "two_layer": true, + "target": "deep" + }, + { + "id": "hard-landslide-09", + "grade": "hard", + "use_case": "landslide", + "motif": "depth", + "split": "test", + "visibility": "public", + "years": 2.5, + "snr": 2.22, + "seed": 308, + "channels": 4, + "decorr": 0.2, + "rms_rel_tol": 0.6, + "n_days": 913, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 7.155, + 18.0 + ], + "window": [ + 0.2, + 1.5 + ], + "stack": 5, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.09, + "rms": 0.0012432270713038142, + "rms_wrong_layer": 0.007530894565057449 + }, + "two_layer": true, + "target": "shallow" + }, + { + "id": "hard-earthquake_fault-10", + "grade": "hard", + "use_case": "earthquake_fault", + "motif": "depth", + "split": "test", + "visibility": "public", + "years": 2.5, + "snr": 2.0, + "seed": 309, + "channels": 4, + "decorr": 0.2, + "rms_rel_tol": 0.6, + "n_days": 913, + "expected": { + "config": { + "estimator": "stretching (TS)", + "band": [ + 0.055, + 0.424 + ], + "window": [ + 8, + 25 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.05, + "rms": 0.00021111898027523828, + "rms_wrong_layer": 0.0007607270475765937 + }, + "two_layer": true, + "target": "deep" } ] } diff --git a/tests/test_frugalmind_export.py b/tests/test_frugalmind_export.py index fb5f64d..935242f 100644 --- a/tests/test_frugalmind_export.py +++ b/tests/test_frugalmind_export.py @@ -43,7 +43,7 @@ def test_split_filter(): def test_param_scorer_rewards_recovery_and_punishes_wrong_choice(): scorer = fm.make_scorer_from_spec({"name": "dvv_recovery"}) - case = golden.CASES_BY_ID["landslide_mainstream"] + case = golden.CASES_BY_ID["easy-landslide-03"] gold = fm._gold(case, "param_recommendation") good = json.dumps(golden._jsonable(uc.recommend("landslide"))) @@ -61,7 +61,7 @@ def test_param_scorer_rewards_recovery_and_punishes_wrong_choice(): def test_param_scorer_accepts_partial_config_filling_noncore_axes(): scorer = fm.make_scorer_from_spec({"name": "dvv_recovery"}) - case = golden.CASES_BY_ID["volcano_mainstream"] + case = golden.CASES_BY_ID["easy-volcano-01"] gold = fm._gold(case, "param_recommendation") # Only the three scientific choices; stack/reference/gate fall back to the # use-case default and the pipeline still recovers. @@ -72,9 +72,9 @@ def test_param_scorer_accepts_partial_config_filling_noncore_axes(): def test_series_scorer_truth_vs_null(): scorer = fm.make_scorer_from_spec({"name": "dvv_series_regression"}) - case = golden.CASES_BY_ID["volcano_mainstream"] + case = golden.CASES_BY_ID["easy-volcano-01"] gold = fm._gold(case, "dvv_series") - d = golden.generate("volcano_mainstream") + d = golden.generate("easy-volcano-01") truth_txt = json.dumps(list(map(float, d["truth"]))) zeros_txt = json.dumps([0.0] * int(gold["n_days"])) diff --git a/tests/test_golden.py b/tests/test_golden.py index 7dd8c1b..846864a 100644 --- a/tests/test_golden.py +++ b/tests/test_golden.py @@ -1,9 +1,9 @@ -"""Golden-dataset regression tests. +"""Golden-dataset regression tests (graded benchmark). Each case in ``tests/data/golden/manifest.json`` is regenerated from its seed and -the recommended pipeline is re-run; the recovered RMS must match the frozen value -within the case's ``rms_rel_tol``. This locks both the estimators and the -use-case recommendation: a change in either shifts the RMS and fails here. +the recommended pipeline is re-run through :func:`codameter.golden.recover`; the +recovered RMS must match the frozen value within the case's ``rms_rel_tol``. This +locks the estimators, the aggregation, and the use-case recommendation. """ from __future__ import annotations @@ -27,8 +27,18 @@ def _rel_close(got: float, want: float, rel_tol: float) -> bool: def test_manifest_is_current(): assert MANIFEST["version"] == golden.MANIFEST_VERSION - # The manifest and the code's case list must not drift apart. assert CASE_IDS == [c["id"] for c in golden.CASES] + assert MANIFEST["grades"] == list(golden.GRADES) + + +def test_thirty_cases_ten_per_grade(): + assert len(golden.CASES) == 30 + for grade in golden.GRADES: + n = sum(1 for c in golden.CASES if c["grade"] == grade) + assert n == 10, f"{grade}: {n}" + # Every application appears in every grade span. + apps = {c["use_case"] for c in golden.CASES} + assert apps == set(golden.AMP) @pytest.mark.parametrize("case_id", CASE_IDS) @@ -36,53 +46,73 @@ def test_case_recovers_within_tolerance(case_id): entry = CASES_BY_ID[case_id] data = golden.generate(case_id) got = golden.compute_metrics(case_id, data) - - want = entry["expected"]["rms"] - tol = entry["rms_rel_tol"] + want, tol = entry["expected"]["rms"], entry["rms_rel_tol"] assert _rel_close(got["rms"], want, tol), ( f"{case_id}: rms {got['rms']:.5f} not within {tol:.0%} of frozen {want:.5f}") - # Probes (e.g. band-selects-depth) must also match their frozen RMS. - for got_p, want_p in zip(got.get("probes", []), - entry["expected"].get("probes", []), strict=True): - assert _rel_close(got_p["rms"], want_p["rms"], tol), ( - f"{case_id}/{got_p['label']}: rms {got_p['rms']:.5f} " - f"not within {tol:.0%} of frozen {want_p['rms']:.5f}") + +def test_easy_cases_recover_cleanly(): + # Best-practice recovery on the easy grade should be well under 0.2 % RMS. + for entry in MANIFEST["cases"]: + if entry["grade"] == "easy": + assert entry["expected"]["rms"] < 2e-3, entry["id"] -def test_mainstream_cases_recover_cleanly(): - # Every mainstream case should recover its truth to well under 0.2 % RMS. +def test_hard_cases_are_multichannel(): for entry in MANIFEST["cases"]: - if entry["kind"] != "mainstream": - continue - assert entry["expected"]["rms"] < 2e-3, entry["id"] + if entry["grade"] == "hard": + assert entry["channels"] > 1, entry["id"] + +def test_recover_handles_single_and_multichannel(): + from codameter import use_cases as uc + + # single channel (easy) + d1 = golden.generate("easy-volcano-01") + assert "channels" not in d1 or np.ndim(d1["channels"]) != 3 + dvv1, val1 = golden.recover(d1, uc.recommend("volcano"), uc.eps_max("volcano")) + assert val1.sum() > 10 and golden._rms(dvv1, d1["truth"], d1["days"], val1) < 2e-3 + + # multi channel (hard): channels present, aggregate recovers the composite + d2 = golden.generate("hard-earthquake_fault-02") + assert np.ndim(d2["channels"]) == 3 and d2["channels"].shape[0] == 4 + dvv2, val2 = golden.recover(d2, uc.recommend("earthquake_fault"), + uc.eps_max("earthquake_fault")) + assert val2.sum() > 10 and np.isfinite(golden._rms(dvv2, d2["truth"], d2["days"], val2)) -def test_freqdep_band_selects_depth(): - # The core edge-case claim: the shallow-band config recovers the shallow - # layer, the deep-band config recovers the deep layer, and using the shallow - # band to read the deep truth is clearly worse. + +def test_hard_grade_band_selects_depth(): + # The hard grade is depth-dependent: the target band recovers the targeted + # layer, a wrong band recovers the other layer and scores clearly worse. from codameter import use_cases as uc - from codameter.deviations import run_pipeline - d = golden.generate("freqdep_shallow_deep") - eps = uc.eps_max("groundwater") - shallow_cfg = uc.recommend("groundwater", band=(4.0, 10.0), window=(2.0, 8.0)) - deep_cfg = uc.recommend("groundwater", band=(0.2, 1.0), window=(8.0, 25.0)) + hard = next(c for c in golden.CASES if c["grade"] == "hard") + assert hard["two_layer"] and hard["target"] in ("shallow", "deep") + d = golden.generate(hard["id"]) + app = hard["use_case"] + eps = uc.eps_max(app) + shallow, deep = golden._depth_bands(app) + wrong_band = deep if hard["target"] == "shallow" else shallow + + dvv_t, vt = golden.recover(d, uc.recommend(app, **hard["config"]), eps) + dvv_w, vw = golden.recover(d, uc.recommend(app, band=wrong_band), eps) + rms_target = golden._rms(dvv_t, d["truth"], d["days"], vt) + rms_wrong = golden._rms(dvv_w, d["truth"], d["days"], vw) + rms_target_on_other = golden._rms(dvv_t, d["truth_other"], d["days"], vt) - dvv_s, val_s = run_pipeline(d["ccfs"], d["t"], d["fs"], shallow_cfg, eps_max=eps) - rms_shallow_on_shallow = golden._rms(dvv_s, d["truth"], d["days"], val_s) - rms_shallow_on_deep = golden._rms(dvv_s, d["truth_deep"], d["days"], val_s) + assert rms_wrong > 2 * rms_target # wrong band is clearly worse + assert rms_target < rms_target_on_other # target band tracks the target layer - dvv_d, val_d = run_pipeline(d["ccfs"], d["t"], d["fs"], deep_cfg, eps_max=eps) - rms_deep_on_deep = golden._rms(dvv_d, d["truth_deep"], d["days"], val_d) - assert rms_shallow_on_shallow < rms_deep_on_deep * 3 # both recover their own layer - assert rms_shallow_on_deep > 2 * rms_shallow_on_shallow # band matters +def test_hard_manifest_records_depth_target(): + for entry in MANIFEST["cases"]: + if entry["grade"] == "hard": + assert entry.get("two_layer") and entry["target"] in ("shallow", "deep") + assert "rms_wrong_layer" in entry["expected"] def test_generate_is_deterministic(): - a = golden.generate("volcano_mainstream", cache=False) - b = golden.generate("volcano_mainstream", cache=False) + a = golden.generate("easy-volcano-01", cache=False) + b = golden.generate("easy-volcano-01", cache=False) assert np.array_equal(a["ccfs"], b["ccfs"]) assert np.array_equal(a["truth"], b["truth"])