From bd573b077a45c38ad1b501ef683505e79e1df287 Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Thu, 9 Jul 2026 06:37:21 -0700 Subject: [PATCH 1/2] Graded 30-case golden benchmark (easy/medium/hard, multi-channel hard) Replace the 10 mainstream/edge golden cases with a graded benchmark: three difficulty grades, ten cases each, spanning the applications. - easy (validation): pure seasonal signal, high SNR (8-12), single channel. Best-practice recovery is well under 0.2 % RMS. - medium (validation): a transient coseismic-style drop with logarithmic partial healing, plus more measurement noise (SNR 3-5). - hard (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 to a network dv/v. golden.py rewrite: - motif builders (seasonal / transient / composite) with per-application amplitudes and timescales; grades cycle through the 6 applications. - golden.recover(d, cfg, eps): single scoring entry point handling single- and multi-channel cases (per-channel run_pipeline, then average the per-channel dv/v). frugalmind and the manifest oracle both go through it. - generate() emits `channels` (3D) for multi-channel cases and downcasts the large CCF arrays to float32 in the cache. representative_case()/ MAINSTREAM_BY_USE_CASE kept for the advisor and sweep. - manifest v2: grade/motif/channels/snr/decorr per case. Verified: the benchmark discriminates on the choices that matter (a moving reference scores 0.0 on a strong-trend hard case, ~0.19 on a drop case), while being band-insensitive for a homogeneous medium (documented). frugalmind rows now 30 per suite; tests updated to grade ids. Full suite: 214 passed, 1 skipped. Follow-ups (separate branches): bench (#13) must switch score_cell to golden.recover and drop the old ids; the live-site golden gallery regenerates for the 30 cases. Co-Authored-By: Claude Opus 4.8 (1M context) --- .../references/golden_datasets.md | 45 +- .../references/validation_loop.md | 14 +- scripts/plot_golden.py | 120 +++ src/codameter/frugalmind.py | 6 +- src/codameter/golden.py | 476 +++++----- tests/data/golden/.gitignore | 3 + tests/data/golden/manifest.json | 886 +++++++++++++++--- tests/test_frugalmind_export.py | 8 +- tests/test_golden.py | 78 +- 9 files changed, 1176 insertions(+), 460 deletions(-) create mode 100644 scripts/plot_golden.py 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..62e4703 --- /dev/null +++ b/scripts/plot_golden.py @@ -0,0 +1,120 @@ +#!/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) + 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) + axC.plot(days, d["truth"] * PCT, color=C["truth"], lw=1.8, label="ground-truth dv/v") + axC.plot(days, rec * PCT, ".", ms=2.6, color=C["rec"], + label="recovered (recommended config)") + 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..0962d13 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" @@ -187,7 +186,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 +283,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..79fdf39 100644 --- a/src/codameter/golden.py +++ b/src/codameter/golden.py @@ -1,243 +1,226 @@ -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* cross-correlation problem whose truth combines a + transient drop-and-heal, a full hydrological seasonal cycle, and a long-term + trend, at low SNR with waveform decorrelation. The channels are measured + independently and aggregated to a network dv/v (per-channel measure, then mean). + +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 _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 _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 -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_seasonal(days: np.ndarray, app: str) -> np.ndarray: + return _seasonal(days, AMP[app]["seasonal"], AMP[app]["phase"]) + + +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} # --------------------------------------------------------------------------- -# 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": "composite", "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) + cases.append({ + "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"], + }) + 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) - - out: dict = {"fs": fs, "days": days, "use_case": use_case} - - 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) + nchan = int(recipe.get("channels", 1)) + truth = MOTIF[recipe["motif"]](days, app) + + 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, "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 +230,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 +245,21 @@ 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) + np.savez_compressed(cache_file, **payload) return d @@ -278,14 +267,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 +287,23 @@ def _rms(dvv, truth, days, valid, baseline_frac: float = 0.2) -> float: return float(np.sqrt(np.mean((d0 - tr0) ** 2))) -def compute_metrics(case_id: str, data: dict | None = None) -> dict: - """Run the recommended pipeline (+ any probes) 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) - 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 - 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()} +def compute_metrics(case_id: str, data: dict | None = None) -> dict: + """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) + app = recipe["use_case"] + cfg = uc.recommend(app, **recipe.get("config", {})) + eps = uc.eps_max(app) + dvv, valid = recover(d, cfg, eps) + return {"config": _jsonable(cfg), "eps_max": eps, + "rms": _rms(dvv, d["truth"], d["days"], valid)} + + # --------------------------------------------------------------------------- # Manifest read / write # --------------------------------------------------------------------------- @@ -337,31 +311,20 @@ 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), + cases.append({ + "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"] - 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 +342,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..d7b688b 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,425 @@ "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, + 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": "composite", + "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": [ + 0.4, + 1.0 + ], + "window": [ + 10, + 30 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.06, + "rms": 0.00026356217881596736 + } + }, + { + "id": "hard-earthquake_fault-02", + "grade": "hard", + "use_case": "earthquake_fault", + "motif": "composite", + "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.5, + 1.5 + ], + "window": [ + 8, + 25 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.05, + "rms": 0.00015763246578731757 + } + }, + { + "id": "hard-landslide-03", + "grade": "hard", + "use_case": "landslide", + "motif": "composite", + "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": [ 4.0, - 10.0 + 12.0 + ], + "window": [ + 0.2, + 1.5 + ], + "stack": 5, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.09, + "rms": 0.001207971655195018 + } + }, + { + "id": "hard-groundwater-04", + "grade": "hard", + "use_case": "groundwater", + "motif": "composite", + "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": [ + 2.0, + 4.0 ], "window": [ 2.0, @@ -301,55 +795,89 @@ "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 - } - ] - }, - "two_layer": true, - "config": { - "band": [ - 4.0, - 10.0 - ], - "window": [ - 2.0, - 8.0 - ] - }, - "note": "Shallow (high-freq) and deep (low-freq) layers carry different dv/v; the band selects which one you recover." + "rms": 0.00019152716187270316 + } + }, + { + "id": "hard-cryosphere-05", + "grade": "hard", + "use_case": "cryosphere", + "motif": "composite", + "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": [ + 4.0, + 14.0 + ], + "window": [ + 0.3, + 0.8 + ], + "stack": 5, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.11, + "rms": 0.0014250912049133585 + } + }, + { + "id": "hard-geothermal-06", + "grade": "hard", + "use_case": "geothermal", + "motif": "composite", + "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.5, + 2.0 + ], + "window": [ + 10.0, + 25.0 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.03, + "rms": 0.0003501958538763985 + } }, { - "id": "sparse_decorr", - "kind": "edge", + "id": "hard-volcano-07", + "grade": "hard", "use_case": "volcano", + "motif": "composite", "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)", @@ -366,9 +894,107 @@ "gate": true }, "eps_max": 0.06, - "rms": 0.0004359819565566644 - }, - "note": "Every-third-day sampling with 30 % waveform decorrelation stresses the reference/stacking warm-up." + "rms": 0.00032773731442452806 + } + }, + { + "id": "hard-groundwater-08", + "grade": "hard", + "use_case": "groundwater", + "motif": "composite", + "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": [ + 2.0, + 4.0 + ], + "window": [ + 2.0, + 8.0 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.03, + "rms": 0.00020391435312739175 + } + }, + { + "id": "hard-landslide-09", + "grade": "hard", + "use_case": "landslide", + "motif": "composite", + "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": [ + 4.0, + 12.0 + ], + "window": [ + 0.2, + 1.5 + ], + "stack": 5, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.09, + "rms": 0.0012465392609767101 + } + }, + { + "id": "hard-earthquake_fault-10", + "grade": "hard", + "use_case": "earthquake_fault", + "motif": "composite", + "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.5, + 1.5 + ], + "window": [ + 8, + 25 + ], + "stack": 10, + "reference": "fixed", + "gate": true + }, + "eps_max": 0.05, + "rms": 0.00020047248208209103 + } } ] } 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..791ddc1 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,43 @@ 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_mainstream_cases_recover_cleanly(): - # Every mainstream case should recover its truth to well under 0.2 % RMS. +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["kind"] != "mainstream": - continue - assert entry["expected"]["rms"] < 2e-3, entry["id"] + if entry["grade"] == "easy": + assert entry["expected"]["rms"] < 2e-3, entry["id"] -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. - from codameter import use_cases as uc - from codameter.deviations import run_pipeline +def test_hard_cases_are_multichannel(): + for entry in MANIFEST["cases"]: + if entry["grade"] == "hard": + assert entry["channels"] > 1, entry["id"] - 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)) - 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) +def test_recover_handles_single_and_multichannel(): + from codameter import use_cases as uc - 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) + # 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 - 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 + # 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_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"]) From 0930f327ed3551fd426488108d37099f22f78386 Mon Sep 17 00:00:00 2001 From: Marine Denolle Date: Thu, 9 Jul 2026 07:59:09 -0700 Subject: [PATCH 2/2] Make the hard grade depth- and frequency-dependent On top of multi-channel + composite, the hard grade is now a depth-dependent medium: a shallow (high-frequency) layer carries the coseismic drop-and-heal plus the full hydrological seasonal cycle; a deep (low-frequency) layer carries the long-term trend. The two layers live in separated bands, so the measurement band selects the depth. Each hard case targets one depth (alternating), and the recommended config's band must match it. - golden.py: _truth_shallow / _truth_deep / _depth_bands; the hard _build path sums both layers per channel; d["truth"] is the targeted layer (d["truth_other"] kept for plotting). compute_metrics adds rms_wrong_layer (the error from recovering the other depth) as the scoring anchor; manifest records two_layer/target. - frugalmind.py: anchor rms_bad to rms_wrong_layer so a band that recovers the wrong depth scores ~0 (verified: target band 1.0, wrong band 0.0). - plot_golden.py: use the target band and overlay the non-targeted layer. - tests: band-selects-depth (target band beats wrong band > 2x and tracks the target layer) and manifest depth-target checks. Verified: measuring in the target band recovers the target layer; a wrong band recovers the other layer and is ~4-6x worse in RMS, 0.0 in score. Full suite: 216 passed, 1 skipped. Co-Authored-By: Claude Opus 4.8 (1M context) --- scripts/plot_golden.py | 12 ++- src/codameter/frugalmind.py | 5 ++ src/codameter/golden.py | 110 +++++++++++++++++++++++---- tests/data/golden/manifest.json | 130 ++++++++++++++++++++------------ tests/test_golden.py | 30 ++++++++ 5 files changed, 221 insertions(+), 66 deletions(-) diff --git a/scripts/plot_golden.py b/scripts/plot_golden.py index 62e4703..43b76d2 100644 --- a/scripts/plot_golden.py +++ b/scripts/plot_golden.py @@ -51,7 +51,7 @@ def plot_case(case_id: str, outdir: Path = OUTDIR) -> Path: 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) + 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)) @@ -79,9 +79,15 @@ def plot_case(case_id: str, outdir: Path = OUTDIR) -> Path: dvv, valid = golden.recover(d, cfg, uc.eps_max(app)) rec = _align(dvv, d["truth"], days, valid) - axC.plot(days, d["truth"] * PCT, color=C["truth"], lw=1.8, label="ground-truth dv/v") + 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 (recommended config)") + 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 (%)", diff --git a/src/codameter/frugalmind.py b/src/codameter/frugalmind.py index 0962d13..0e2d94d 100644 --- a/src/codameter/frugalmind.py +++ b/src/codameter/frugalmind.py @@ -137,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} diff --git a/src/codameter/golden.py b/src/codameter/golden.py index 79fdf39..6fef8ba 100644 --- a/src/codameter/golden.py +++ b/src/codameter/golden.py @@ -7,10 +7,12 @@ (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* cross-correlation problem whose truth combines a - transient drop-and-heal, a full hydrological seasonal cycle, and a long-term - trend, at low SNR with waveform decorrelation. The channels are measured - independently and aggregated to a network dv/v (per-channel measure, then mean). +- **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 @@ -89,6 +91,34 @@ def _motif_composite(days: np.ndarray, app: str) -> np.ndarray: "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 + + # --------------------------------------------------------------------------- # Grades and case construction. Each grade cycles through the applications so it # spans volcano / fault / aquifer / glacier / reservoir at that difficulty. @@ -98,7 +128,7 @@ def _motif_composite(days: np.ndarray, app: str) -> np.ndarray: "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": "composite", "snr": (2.0, 4.0), "channels": 4, "decorr": 0.20, + "hard": {"motif": "depth", "snr": (2.0, 4.0), "channels": 4, "decorr": 0.20, "years": 2.5, "split": "test", "rms_rel_tol": 0.60}, } @@ -116,14 +146,32 @@ def _build_cases() -> list[dict]: 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) - cases.append({ + 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 @@ -174,13 +222,37 @@ def _build(recipe: dict) -> dict: seed, snr = recipe["seed"], recipe["snr"] decorr = recipe.get("decorr", 0.0) nchan = int(recipe.get("channels", 1)) - truth = MOTIF[recipe["motif"]](days, app) 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, "truth": truth} + "t": t} + + if recipe.get("two_layer"): + # 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. @@ -259,6 +331,8 @@ def generate(case_id: str, *, cache: bool = True) -> dict: "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 @@ -300,8 +374,14 @@ def compute_metrics(case_id: str, data: dict | None = None) -> dict: cfg = uc.recommend(app, **recipe.get("config", {})) eps = uc.eps_max(app) dvv, valid = recover(d, cfg, eps) - return {"config": _jsonable(cfg), "eps_max": eps, - "rms": _rms(dvv, d["truth"], d["days"], valid)} + res = {"config": _jsonable(cfg), "eps_max": eps, + "rms": _rms(dvv, d["truth"], d["days"], valid)} + 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 # --------------------------------------------------------------------------- @@ -313,14 +393,18 @@ def regenerate_manifest() -> dict: for c in CASES: d = generate(c["id"], cache=False) # always from current code m = compute_metrics(c["id"], d) - cases.append({ + entry = { "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, - }) + } + if c.get("two_layer"): + entry["two_layer"] = True + entry["target"] = c["target"] + cases.append(entry) print(f" {c['id']:<26} ch={c['channels']} snr={c['snr']:<4} " f"rms={m['rms']:.5f}") manifest = {"version": MANIFEST_VERSION, "grades": list(GRADES), diff --git a/tests/data/golden/manifest.json b/tests/data/golden/manifest.json index d7b688b..db02222 100644 --- a/tests/data/golden/manifest.json +++ b/tests/data/golden/manifest.json @@ -670,7 +670,7 @@ "id": "hard-volcano-01", "grade": "hard", "use_case": "volcano", - "motif": "composite", + "motif": "depth", "split": "test", "visibility": "public", "years": 2.5, @@ -684,8 +684,8 @@ "config": { "estimator": "stretching (TS)", "band": [ - 0.4, - 1.0 + 1.131, + 9.0 ], "window": [ 10, @@ -696,14 +696,17 @@ "gate": true }, "eps_max": 0.06, - "rms": 0.00026356217881596736 - } + "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": "composite", + "motif": "depth", "split": "test", "visibility": "public", "years": 2.5, @@ -717,8 +720,8 @@ "config": { "estimator": "stretching (TS)", "band": [ - 0.5, - 1.5 + 0.055, + 0.424 ], "window": [ 8, @@ -729,14 +732,17 @@ "gate": true }, "eps_max": 0.05, - "rms": 0.00015763246578731757 - } + "rms": 0.00012877251474255628, + "rms_wrong_layer": 0.0007607270475765937 + }, + "two_layer": true, + "target": "deep" }, { "id": "hard-landslide-03", "grade": "hard", "use_case": "landslide", - "motif": "composite", + "motif": "depth", "split": "test", "visibility": "public", "years": 2.5, @@ -750,8 +756,8 @@ "config": { "estimator": "stretching (TS)", "band": [ - 4.0, - 12.0 + 7.155, + 18.0 ], "window": [ 0.2, @@ -762,14 +768,17 @@ "gate": true }, "eps_max": 0.09, - "rms": 0.001207971655195018 - } + "rms": 0.0012734011077726465, + "rms_wrong_layer": 0.007530894565057449 + }, + "two_layer": true, + "target": "shallow" }, { "id": "hard-groundwater-04", "grade": "hard", "use_case": "groundwater", - "motif": "composite", + "motif": "depth", "split": "test", "visibility": "public", "years": 2.5, @@ -783,8 +792,8 @@ "config": { "estimator": "stretching (TS)", "band": [ - 2.0, - 4.0 + 0.55, + 1.2 ], "window": [ 2.0, @@ -795,14 +804,17 @@ "gate": true }, "eps_max": 0.03, - "rms": 0.00019152716187270316 - } + "rms": 0.00022902404646362183, + "rms_wrong_layer": 0.0010007890743546153 + }, + "two_layer": true, + "target": "deep" }, { "id": "hard-cryosphere-05", "grade": "hard", "use_case": "cryosphere", - "motif": "composite", + "motif": "depth", "split": "test", "visibility": "public", "years": 2.5, @@ -816,8 +828,8 @@ "config": { "estimator": "stretching (TS)", "band": [ - 4.0, - 14.0 + 11.314, + 22.5 ], "window": [ 0.3, @@ -828,14 +840,17 @@ "gate": true }, "eps_max": 0.11, - "rms": 0.0014250912049133585 - } + "rms": 0.0010750047910110425, + "rms_wrong_layer": 0.017753287172458793 + }, + "two_layer": true, + "target": "shallow" }, { "id": "hard-geothermal-06", "grade": "hard", "use_case": "geothermal", - "motif": "composite", + "motif": "depth", "split": "test", "visibility": "public", "years": 2.5, @@ -849,8 +864,8 @@ "config": { "estimator": "stretching (TS)", "band": [ - 0.5, - 2.0 + 0.11, + 0.537 ], "window": [ 10.0, @@ -861,14 +876,17 @@ "gate": true }, "eps_max": 0.03, - "rms": 0.0003501958538763985 - } + "rms": 0.00035681825048443534, + "rms_wrong_layer": 0.0037482773313763365 + }, + "two_layer": true, + "target": "deep" }, { "id": "hard-volcano-07", "grade": "hard", "use_case": "volcano", - "motif": "composite", + "motif": "depth", "split": "test", "visibility": "public", "years": 2.5, @@ -882,8 +900,8 @@ "config": { "estimator": "stretching (TS)", "band": [ - 0.4, - 1.0 + 1.131, + 9.0 ], "window": [ 10, @@ -894,14 +912,17 @@ "gate": true }, "eps_max": 0.06, - "rms": 0.00032773731442452806 - } + "rms": 0.00025180197315354186, + "rms_wrong_layer": 0.000930391397212034 + }, + "two_layer": true, + "target": "shallow" }, { "id": "hard-groundwater-08", "grade": "hard", "use_case": "groundwater", - "motif": "composite", + "motif": "depth", "split": "test", "visibility": "public", "years": 2.5, @@ -915,8 +936,8 @@ "config": { "estimator": "stretching (TS)", "band": [ - 2.0, - 4.0 + 0.55, + 1.2 ], "window": [ 2.0, @@ -927,14 +948,17 @@ "gate": true }, "eps_max": 0.03, - "rms": 0.00020391435312739175 - } + "rms": 0.00031478238532623977, + "rms_wrong_layer": 0.0010007890743546153 + }, + "two_layer": true, + "target": "deep" }, { "id": "hard-landslide-09", "grade": "hard", "use_case": "landslide", - "motif": "composite", + "motif": "depth", "split": "test", "visibility": "public", "years": 2.5, @@ -948,8 +972,8 @@ "config": { "estimator": "stretching (TS)", "band": [ - 4.0, - 12.0 + 7.155, + 18.0 ], "window": [ 0.2, @@ -960,14 +984,17 @@ "gate": true }, "eps_max": 0.09, - "rms": 0.0012465392609767101 - } + "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": "composite", + "motif": "depth", "split": "test", "visibility": "public", "years": 2.5, @@ -981,8 +1008,8 @@ "config": { "estimator": "stretching (TS)", "band": [ - 0.5, - 1.5 + 0.055, + 0.424 ], "window": [ 8, @@ -993,8 +1020,11 @@ "gate": true }, "eps_max": 0.05, - "rms": 0.00020047248208209103 - } + "rms": 0.00021111898027523828, + "rms_wrong_layer": 0.0007607270475765937 + }, + "two_layer": true, + "target": "deep" } ] } diff --git a/tests/test_golden.py b/tests/test_golden.py index 791ddc1..846864a 100644 --- a/tests/test_golden.py +++ b/tests/test_golden.py @@ -81,6 +81,36 @@ def test_recover_handles_single_and_multichannel(): assert val2.sum() > 10 and np.isfinite(golden._rms(dvv2, d2["truth"], d2["days"], val2)) +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 + + 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) + + assert rms_wrong > 2 * rms_target # wrong band is clearly worse + assert rms_target < rms_target_on_other # target band tracks the target layer + + +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("easy-volcano-01", cache=False) b = golden.generate("easy-volcano-01", cache=False)