From 678f00ff993705ba1d3eed344fc1938c0450484b Mon Sep 17 00:00:00 2001 From: "Simon D.A. Thomas" Date: Sun, 16 Aug 2026 18:54:15 +0100 Subject: [PATCH 01/17] move comp to eval --- README.md | 8 +-- REPRODUCE.md | 18 +++---- {comp => adforce/eval}/README.md | 20 ++++---- {comp => adforce/eval}/__init__.py | 2 +- {comp => adforce/eval}/annual_max.py | 20 ++++---- {comp => adforce/eval}/constants.py | 12 ++--- {comp => adforce/eval}/coops.py | 2 +- {comp => adforce/eval}/detide_sensitivity.py | 12 ++--- {comp => adforce/eval}/nulltest.py | 6 +-- {comp => adforce/eval}/sensitivity.py | 4 +- {comp => adforce/eval}/validate.py | 12 ++--- docs/comp.rst | 53 -------------------- docs/conf.py | 2 +- docs/eval.rst | 53 ++++++++++++++++++++ docs/index.rst | 2 +- setup.py | 29 +++++++---- tests/{test_comp.py => test_eval.py} | 16 +++--- tests/test_gauge_evt.py | 8 +-- worst/gauge_fit.py | 4 +- worst/gauge_fit_ns.py | 2 +- 20 files changed, 146 insertions(+), 139 deletions(-) rename {comp => adforce/eval}/README.md (86%) rename {comp => adforce/eval}/__init__.py (90%) rename {comp => adforce/eval}/annual_max.py (94%) rename {comp => adforce/eval}/constants.py (94%) rename {comp => adforce/eval}/coops.py (98%) rename {comp => adforce/eval}/detide_sensitivity.py (96%) rename {comp => adforce/eval}/nulltest.py (97%) rename {comp => adforce/eval}/sensitivity.py (97%) rename {comp => adforce/eval}/validate.py (98%) delete mode 100644 docs/comp.rst create mode 100644 docs/eval.rst rename tests/{test_comp.py => test_eval.py} (97%) diff --git a/README.md b/README.md index 0fd6aa9c..14125225 100755 --- a/README.md +++ b/README.md @@ -80,9 +80,9 @@ Key scripts: python -m adbo.exp_1d --test True --exp_name test ``` -### `comp` +### `adforce.eval` -Validates the historical ADCIRC surge simulations (the SurgeNet training set, published on Hugging Face) against de-tided NOAA CO-OPS tide-gauge observations, scoring peak skill, time-series skill, and peak timing, with permutation/negative-control null tests and sensitivity checks. See [`comp/README.md`](https://github.com/sdat2/worstsurge/blob/main/comp/README.md) for the validation methodology, negative controls, and results. +Validates the historical ADCIRC surge simulations (the SurgeNet training set, published on Hugging Face) against de-tided NOAA CO-OPS tide-gauge observations, scoring peak skill, time-series skill, and peak timing, with permutation/negative-control null tests and sensitivity checks. See [`adforce/eval/README.md`](https://github.com/sdat2/worstsurge/blob/main/adforce/eval/README.md) for the validation methodology, negative controls, and results. (Formerly the top-level `comp/` package; caches remain under `data/comp/`.) ### `worst` @@ -104,8 +104,8 @@ A motivating toy example (`toy_example.py`): a single small Keras ReLU MLP is tr ``` adbo/ Bayesian optimization loop (trieste) adforce/ ADCIRC wrapper and forcing utilities +adforce/eval/ Historical surge validation against NOAA tide gauges (see adforce/eval/README.md) cle15/ CLE15 wind profile implementations (pure-Python & Numba) -comp/ Historical surge validation against NOAA tide gauges (see comp/README.md) data/ Key input data (fort.22.nc, IBTrACS, ERA5, CMIP6, etc.) docs/ ReadTheDocs source img/ Key figures @@ -152,7 +152,7 @@ dependency groups are available as extras: pip install -e .[bo] # trieste/TensorFlow Bayesian optimization (adbo, worst) pip install -e .[cmip] # intake/xESMF CMIP6 download + regridding (tcpips) pip install -e .[mpi] # dask_mpi/mpi4py HPC parallelism (tcpips) -pip install -e .[comp] # utide/huggingface_hub tide-gauge validation (comp) +pip install -e .[eval] # utide/huggingface_hub tide-gauge validation (adforce.eval; [comp] is a deprecated alias) pip install -e .[all] # everything ``` diff --git a/REPRODUCE.md b/REPRODUCE.md index 47e1b96a..c4e564ed 100644 --- a/REPRODUCE.md +++ b/REPRODUCE.md @@ -78,7 +78,7 @@ Bayesian optimization (`adbo.exp*`) → EVT with an upper bound (`worst`)** | 6. ADCIRC forcing | `python -m adforce.wrap` | ARCHER2 only (needs `padcirc`/`adcprep` binaries, see `adforce/config/files/`) | `wrap.slurm`, `n01_wrap.slurm` | | 7. Bayesian optimization | `python -m adbo.exp_1d/_2d/_3d ...` | ARCHER2 only (drives ADCIRC via `adforce`) | `[bo]` extra; `bo_*.slurm`, `exp.slurm`, per-city scripts, `n01_bo_first..fourth` | | 8. EVT with upper bound | `python -m worst.vary_samples_ns`, `worst.vary_noise`, `worst.vary_nonstationary`, `worst.sigma_robustness`, `worst.ns_evt_figs` | local (CPU; `[bo]` extra for TensorFlow fits) | synthetic experiments — no ADCIRC needed | -| Validation (side chain) | `python -m comp.validate` / `comp.nulltest` / `comp.detide_sensitivity` / `comp.sensitivity` | local (needs internet + `[comp]` extra) | Hugging Face archived runs vs NOAA gauges | +| Validation (side chain) | `python -m adforce.eval.validate` / `adforce.eval.nulltest` / `adforce.eval.detide_sensitivity` / `adforce.eval.sensitivity` | local (needs internet + `[eval]` extra) | Hugging Face archived runs vs NOAA gauges | | Observations (side chain, thesis) | `python -m tcpips.ibtracs` | local or ARCHER2 (`ibtracs.slurm`) | IBTrACS vs ERA5-derived PI/PS | Stages 6–7 require ARCHER2 (or another HPC machine with a ported @@ -107,17 +107,17 @@ Some plotting entry points are driven from commented-out calls in the module | **Fig. 10** — AEP with known vs unknown upper bound, GEV resampling (Nr=600) (`worst:fig:evt`) | `python -m worst.vary_samples_ns` (hydra config `worst/config/`) | synthetic GEV draws (TensorFlow fit) | `img/worst/evt_fig_tens_*.pdf` | — (local) | | **Fig. 11** — effect of upper-bound uncertainty σ_ẑ* on return-value estimates (`worst:fig:vary_z_star_sigma`) | `python -m worst.vary_noise` | synthetic GEV draws; cache `data/worst/vary_z_star_*.nc` | `img/worst/vary_z_star_sigma.pdf` | — (local) | -### Appendix A — historical surge validation (`comp`; see [`comp/README.md`](comp/README.md)) +### Appendix A — historical surge validation (`adforce.eval`; see [`adforce/eval/README.md`](adforce/eval/README.md)) | Paper item | Command | Key inputs | Output | SLURM job | | --- | --- | --- | --- | --- | -| **Table 3** — de-tiding-method sensitivity (`tab:detide-robustness`) | `python -m comp.detide_sensitivity` | de-tided gauge cache (`data/comp/ts_cache/`) | `/paper/comp_detide_table.tex` (`\input` by the appendix) | — (local) | -| **Table 4** — per-storm skill vs de-tided NOAA gauges (`tab:gauge-validation`) | `python -m comp.validate` | Hugging Face [`sdat2/surgenet-train`](https://huggingface.co/datasets/sdat2/surgenet-train); NOAA CO-OPS gauges (de-tided with `utide`) | `/paper/comp_val_table.tex` | — (local) | -| **Fig. 12** — peak-skill scatter (`fig:gauge-validation-scatter`) | `python -m comp.validate` (full 14-storm sweep) | as above | `img/comp/val_scatter.png` · `/img/comp_val_scatter.pdf` | — (local) | -| **Fig. 13** — example de-tided surge time series (`fig:gauge-validation-examples`) | `python -m comp.validate` (or `--examples-only` from cache) | as above (cache: `data/comp/ts_cache/`) | `img/comp/val_examples.png` · `/img/comp_val_examples.pdf` | — (local) | -| **Fig. 14** — negative-control null tests (`fig:gauge-validation-null`) | `python -m comp.nulltest` | `data/comp/out/val_summary.csv` (+ netCDFs for the lag null) | `/img/comp_val_nulltests.pdf` · `img/comp/val_nulltests.png` | — (local) | -| Summary metrics (support; pinned by `tests/test_comp.py`) | `python -m comp.validate` | as above | `data/comp/out/val_summary.csv` | — (local) | -| Threshold/node-selection sensitivity (text numbers) | `python -m comp.sensitivity` | `val_summary.csv` | console report (headline r stability) | — (local) | +| **Table 3** — de-tiding-method sensitivity (`tab:detide-robustness`) | `python -m adforce.eval.detide_sensitivity` | de-tided gauge cache (`data/comp/ts_cache/`) | `/paper/comp_detide_table.tex` (`\input` by the appendix) | — (local) | +| **Table 4** — per-storm skill vs de-tided NOAA gauges (`tab:gauge-validation`) | `python -m adforce.eval.validate` | Hugging Face [`sdat2/surgenet-train`](https://huggingface.co/datasets/sdat2/surgenet-train); NOAA CO-OPS gauges (de-tided with `utide`) | `/paper/comp_val_table.tex` | — (local) | +| **Fig. 12** — peak-skill scatter (`fig:gauge-validation-scatter`) | `python -m adforce.eval.validate` (full 14-storm sweep) | as above | `img/comp/val_scatter.png` · `/img/comp_val_scatter.pdf` | — (local) | +| **Fig. 13** — example de-tided surge time series (`fig:gauge-validation-examples`) | `python -m adforce.eval.validate` (or `--examples-only` from cache) | as above (cache: `data/comp/ts_cache/`) | `img/comp/val_examples.png` · `/img/comp_val_examples.pdf` | — (local) | +| **Fig. 14** — negative-control null tests (`fig:gauge-validation-null`) | `python -m adforce.eval.nulltest` | `data/comp/out/val_summary.csv` (+ netCDFs for the lag null) | `/img/comp_val_nulltests.pdf` · `img/comp/val_nulltests.png` | — (local) | +| Summary metrics (support; pinned by `tests/test_eval.py`) | `python -m adforce.eval.validate` | as above | `data/comp/out/val_summary.csv` | — (local) | +| Threshold/node-selection sensitivity (text numbers) | `python -m adforce.eval.sensitivity` | `val_summary.csv` | console report (headline r stability) | — (local) | ### Appendix B — non-stationary EVT diff --git a/comp/README.md b/adforce/eval/README.md similarity index 86% rename from comp/README.md rename to adforce/eval/README.md index 7ae2f62c..0fc64a92 100644 --- a/comp/README.md +++ b/adforce/eval/README.md @@ -35,10 +35,10 @@ For each storm: ## Run ```bash -python -m comp.validate # full 19-storm sweep (Gulf + Florida) -python -m comp.validate --storms "Ida 2021" # one storm (skips example panels + table) -python -m comp.validate --examples-only # just the example figure, from cache (fast) -python -m comp.validate --examples-only --refresh-cache # recompute the cached series first +python -m adforce.eval.validate # full 19-storm sweep (Gulf + Florida) +python -m adforce.eval.validate --storms "Ida 2021" # one storm (skips example panels + table) +python -m adforce.eval.validate --examples-only # just the example figure, from cache (fast) +python -m adforce.eval.validate --examples-only --refresh-cache # recompute the cached series first ``` A full sweep regenerates **everything the paper uses**, in one step, so the figures, the @@ -65,17 +65,17 @@ The thesis tree is located by searching for `paper/appendix.tex`; override with `val_summary.csv` is enough to falsify the result without re-running ADCIRC: ```bash -python -m comp.nulltest # permutation, cross-storm, and temporal-lag nulls -python -m comp.nulltest --no-lag # peak-level nulls only (no netCDF) -python -m comp.sensitivity # threshold + node-selection robustness -python -m comp.sensitivity --no-node +python -m adforce.eval.nulltest # permutation, cross-storm, and temporal-lag nulls +python -m adforce.eval.nulltest --no-lag # peak-level nulls only (no netCDF) +python -m adforce.eval.sensitivity # threshold + node-selection robustness +python -m adforce.eval.sensitivity --no-node ``` -`comp.nulltest` writes `img/comp_val_nulltests.pdf` and reports: a label-permutation null +`adforce.eval.nulltest` writes `img/comp_val_nulltests.pdf` and reports: a label-permutation null (observed r=0.89 vs null max 0.31 over 5000 shuffles, p<1e-3); a within-storm permutation (observed spatial r=0.86 vs null max 0.61) showing real *spatial* skill; a cross-storm same-gauge null (collapses to r~0.1–0.2); and a temporal-lag curve (time-series r peaks -sharply at lag 0 and decays to zero within two days). `comp.sensitivity` shows the headline +sharply at lag 0 and decays to zero within two days). `adforce.eval.sensitivity` shows the headline is stable across the clean-filter cut-offs (r 0.81–0.89) and the node-selection knobs (r changes <0.001 with the wet-depth threshold). diff --git a/comp/__init__.py b/adforce/eval/__init__.py similarity index 90% rename from comp/__init__.py rename to adforce/eval/__init__.py index a0c5cfae..f5f44034 100644 --- a/comp/__init__.py +++ b/adforce/eval/__init__.py @@ -5,7 +5,7 @@ Entry point:: - python -m comp.validate + python -m adforce.eval.validate """ from .validate import run, validate_storm, metrics diff --git a/comp/annual_max.py b/adforce/eval/annual_max.py similarity index 94% rename from comp/annual_max.py rename to adforce/eval/annual_max.py index 09aaf527..4581a649 100644 --- a/comp/annual_max.py +++ b/adforce/eval/annual_max.py @@ -4,20 +4,20 @@ (one calendar-year request at a time -- the CO-OPS ``hourly_height`` product accepts up to a year per request), de-tides each calendar year separately with a robust ``utide`` harmonic fit (mean + trend included, following -:func:`comp.coops.observed_residual`), and records the ANNUAL MAXIMUM of the +:func:`adforce.eval.coops.observed_residual`), and records the ANNUAL MAXIMUM of the surge residual for every year with adequate data. Missing-data policy: a year is skipped unless (a) it has at least ``UTIDE_MIN_SAMPLES`` hourly samples (stable harmonic fit, same threshold as -``comp.coops``) and (b) at least ``MIN_YEAR_COVERAGE`` (80%) of the year's +``adforce.eval.coops``) and (b) at least ``MIN_YEAR_COVERAGE`` (80%) of the year's hours are present, since with large gaps the true annual maximum may fall in a gap and the recorded maximum would be biased low. Even in accepted years the maximum can be truncated if the gauge failed *during* the peak (e.g. Grand -Isle in Ida 2021, see ``comp.constants.KNOWN_FAILED``); such years are flagged +Isle in Ida 2021, see ``adforce.eval.constants.KNOWN_FAILED``); such years are flagged with ``max_at_gap_edge`` rather than dropped. Everything is cached under ``data/comp/``: the raw CO-OPS responses land in -``COOPS_CACHE`` (via :mod:`comp.coops`), the per-year de-tided residuals and +``COOPS_CACHE`` (via :mod:`adforce.eval.coops`), the per-year de-tided residuals and the annual-maxima table are Parquet files under ``ANNUAL_MAX_CACHE``, keyed by the de-tiding parameters so the cache self-invalidates if those change. Reruns are therefore free. @@ -28,8 +28,8 @@ Run:: - python -m comp.annual_max --station 8761724 # Grand Isle, LA - python -m comp.annual_max --station 8735180 --start 1980 # Dauphin Island, AL + python -m adforce.eval.annual_max --station 8761724 # Grand Isle, LA + python -m adforce.eval.annual_max --station 8735180 --start 1980 # Dauphin Island, AL """ from __future__ import annotations @@ -68,7 +68,7 @@ def _hours_in_year(year: int) -> int: def station_meta(station: str) -> Optional[dict]: """Name/lat/lon for a CO-OPS station from the cached station list. - Reuses the same ``stations.json`` cache as :func:`comp.coops.gulf_gauges` + Reuses the same ``stations.json`` cache as :func:`adforce.eval.coops.gulf_gauges` (all water-level stations, no bounding box), fetching it once if absent. Args: @@ -96,7 +96,7 @@ def station_meta(station: str) -> Optional[dict]: # --------------------------------------------------------------------------- # -# Fetching. comp.coops._coops caches every response, including transport +# Fetching. adforce.eval.coops._coops caches every response, including transport # failures (as zero-byte files). A genuine "no data" year is a NON-empty API # error message and is accepted at once; a zero-byte cache entry is deleted and # the year re-requested a few times, so a one-off network blip cannot @@ -127,7 +127,7 @@ def fetch_year_wl( ) -> pd.Series: """Hourly water level (MSL, GMT) for one calendar year, cached + polite. - Delegates to :func:`comp.coops.fetch_year` (verified ``hourly_height``, + Delegates to :func:`adforce.eval.coops.fetch_year` (verified ``hourly_height``, falling back to preliminary ``water_level``); on an empty result caused by a cached transport failure the bad cache entry is cleared and the request retried with a growing sleep. Tests monkeypatch this function. @@ -164,7 +164,7 @@ def fetch_year_wl( def detide_year(wl: pd.Series, lat: float, method: str = "robust") -> pd.Series: """De-tide one calendar year of hourly water levels with ``utide``. - Mirrors :func:`comp.coops.observed_residual` (harmonic fit with mean and + Mirrors :func:`adforce.eval.coops.observed_residual` (harmonic fit with mean and linear trend over the full year, residual = observed - reconstruction), with one deliberate difference: the DatetimeIndex is passed to ``utide`` directly instead of via ``matplotlib.dates.date2num``. With matplotlib's diff --git a/comp/constants.py b/adforce/eval/constants.py similarity index 94% rename from comp/constants.py rename to adforce/eval/constants.py index 0f7f7abc..bd333065 100644 --- a/comp/constants.py +++ b/adforce/eval/constants.py @@ -1,8 +1,8 @@ -"""Constants and configuration for the ``comp`` (observational comparison) module. +"""Constants and configuration for the ``adforce.eval`` (observational comparison) module. This module validates the historical ADCIRC storm-surge simulations (the SurgeNet training set, published on Hugging Face) against de-tided NOAA CO-OPS tide-gauge -observations. See :mod:`comp.validate`. +observations. See :mod:`adforce.eval.validate`. """ import os @@ -10,7 +10,7 @@ # --- paths ----------------------------------------------------------------- SRC_PATH = Path(__file__).parent -PROJ_PATH = SRC_PATH.parent +PROJ_PATH = SRC_PATH.parent.parent # adforce/eval/ -> repo root DATA_PATH = os.path.join(PROJ_PATH, "data") COMP_DATA_PATH = os.path.join(DATA_PATH, "comp") # caches (git-ignored) HF_STORM_CACHE = os.path.join(COMP_DATA_PATH, "hf_storms") # downloaded storm netCDFs @@ -26,7 +26,7 @@ # (thesis/worstsurge -> worstsurge), so the module's own parent is NOT the thesis # root: the Environmental Data Science paper builds from a separate that # holds paper/appendix.tex and img/. We write the final PDFs and the generated -# LaTeX table straight there, so ``python -m comp.validate`` reproduces exactly +# LaTeX table straight there, so ``python -m adforce.eval.validate`` reproduces exactly # what the paper \includegraphics/\inputs -- no manual copy/convert step that # could silently drift. The root is located by looking for paper/appendix.tex # (override with WORSTSURGE_PAPER_ROOT); absent it, we fall back to module dirs. @@ -77,7 +77,7 @@ def _find_paper_root(): # moved from the single New Orleans site to the three-city comparison # (New Orleans / Galveston / Miami). A separate box (rather than widening # GAUGE_BOX) keeps the Gulf storms' time-series cache valid: the box is part of -# the cache tag (see comp.validate._ts_cache_tag). +# the cache tag (see adforce.eval.validate._ts_cache_tag). FLORIDA_BOX = dict(lon=(-82.3, -79.7), lat=(24.4, 30.8)) # Storms scored against FLORIDA_BOX instead of the (default Gulf) GAUGE_BOX. @@ -171,7 +171,7 @@ def box_for(storm: str) -> dict: # Minimum hourly samples in a calendar year for a stable utide harmonic fit. UTIDE_MIN_SAMPLES = 2000 -# --- annual-maximum pipeline (comp.annual_max) ------------------------------ +# --- annual-maximum pipeline (adforce.eval.annual_max) ------------------------------ ANNUAL_MAX_CACHE = os.path.join( COMP_DATA_PATH, "annual_max" ) # residuals + annual maxima diff --git a/comp/coops.py b/adforce/eval/coops.py similarity index 98% rename from comp/coops.py rename to adforce/eval/coops.py index db74277b..def0cd28 100644 --- a/comp/coops.py +++ b/adforce/eval/coops.py @@ -118,7 +118,7 @@ def observed_residual( # selects ZERO constituents (residual keeps the full tide). # Verified on Grand Isle 2020: date2num path -> empty constituent # list, tide std 0.03 m; DatetimeIndex -> SA/K1/O1 amplitudes - # 0.16/0.12/0.12 m, tide std 0.17 m. Same fix as comp.annual_max. + # 0.16/0.12/0.12 m, tide std 0.17 m. Same fix as adforce.eval.annual_max. coef = utide.solve( wl.index, wl.values, diff --git a/comp/detide_sensitivity.py b/adforce/eval/detide_sensitivity.py similarity index 96% rename from comp/detide_sensitivity.py rename to adforce/eval/detide_sensitivity.py index a5fc459a..209ceb34 100644 --- a/comp/detide_sensitivity.py +++ b/adforce/eval/detide_sensitivity.py @@ -12,14 +12,14 @@ A full run also writes the appendix table ``/paper/comp_detide_table.tex`` (the ``tabular`` that ``tab:detide-robustness`` ``\\input``s), so the paper table regenerates with -the data and cannot drift -- like the per-storm table from :func:`comp.validate.latex_table`. +the data and cannot drift -- like the per-storm table from :func:`adforce.eval.validate.latex_table`. Run:: - python -m comp.detide_sensitivity # all methods + skew, all pairs, writes the table - python -m comp.detide_sensitivity --methods godin_lowpass noaa_predictions - python -m comp.detide_sensitivity --limit 40 # quick subset (no table written) - python -m comp.detide_sensitivity --skew # instantaneous vs skew-surge metric only + python -m adforce.eval.detide_sensitivity # all methods + skew, all pairs, writes the table + python -m adforce.eval.detide_sensitivity --methods godin_lowpass noaa_predictions + python -m adforce.eval.detide_sensitivity --limit 40 # quick subset (no table written) + python -m adforce.eval.detide_sensitivity --skew # instantaneous vs skew-surge metric only """ from __future__ import annotations @@ -243,7 +243,7 @@ def _latex_table(pooled: Dict[str, dict], path: str) -> None: """Emit the de-tiding-robustness ``tabular`` the appendix \\inputs. Only the ``tabular`` (the ``table`` float, caption and label live in - ``paper/appendix.tex``), matching :func:`comp.validate.latex_table`, so the prose stays + ``paper/appendix.tex``), matching :func:`adforce.eval.validate.latex_table`, so the prose stays hand-edited while every number is generated -- the table cannot drift from the data. """ lines = [ diff --git a/comp/nulltest.py b/adforce/eval/nulltest.py similarity index 97% rename from comp/nulltest.py rename to adforce/eval/nulltest.py index e00d34c5..2decee43 100644 --- a/comp/nulltest.py +++ b/adforce/eval/nulltest.py @@ -21,8 +21,8 @@ Run:: - python -m comp.nulltest # peak-level nulls + lag curve (+ figure) - python -m comp.nulltest --no-lag # skip the netCDF-heavy lag test + python -m adforce.eval.nulltest # peak-level nulls + lag curve (+ figure) + python -m adforce.eval.nulltest --no-lag # skip the netCDF-heavy lag test """ from __future__ import annotations @@ -141,7 +141,7 @@ def lag_curve( per_lag: Dict[float, List[float]] = {L: [] for L in lags_days} for storm, fname in items.items(): try: - # per-storm gauge region (Gulf vs Florida), matching comp.validate + # per-storm gauge region (Gulf vs Florida), matching adforce.eval.validate _, series = validate_storm(storm, fname, gulf_gauges(C.box_for(storm))) except Exception as e: # pragma: no cover print(f"!! {storm}: {e}") diff --git a/comp/sensitivity.py b/adforce/eval/sensitivity.py similarity index 97% rename from comp/sensitivity.py rename to adforce/eval/sensitivity.py index 8ca28d73..63335ce2 100644 --- a/comp/sensitivity.py +++ b/adforce/eval/sensitivity.py @@ -13,8 +13,8 @@ Run:: - python -m comp.sensitivity # cheap filter sweep + node-selection sweep - python -m comp.sensitivity --no-node # cheap filter sweep only (no netCDF) + python -m adforce.eval.sensitivity # cheap filter sweep + node-selection sweep + python -m adforce.eval.sensitivity --no-node # cheap filter sweep only (no netCDF) """ from __future__ import annotations diff --git a/comp/validate.py b/adforce/eval/validate.py similarity index 98% rename from comp/validate.py rename to adforce/eval/validate.py index 5d8780ce..7483d544 100644 --- a/comp/validate.py +++ b/adforce/eval/validate.py @@ -4,7 +4,7 @@ 1. download the storm's netCDF from Hugging Face (``HF_REPO``); 2. extract the simulated surge (SSH = WD + DEM) at the nearest *wet* mesh element centroid (the archived dual-graph node) to each NOAA CO-OPS gauge in the box; - 3. fetch + de-tide the gauge record (:func:`comp.coops.observed_residual`); + 3. fetch + de-tide the gauge record (:func:`adforce.eval.coops.observed_residual`); 4. score peak surge (bias/RMSE/correlation, with bootstrap CIs and a within-storm spatial correlation), the full hydrograph (:func:`timeseries_skill`), and peak timing; tag "clean" pairs and regenerate the paper figures + LaTeX table. @@ -20,10 +20,10 @@ Run:: - python -m comp.validate # full sweep, all STORMS (populates the cache) - python -m comp.validate --storms "Ida 2021" "Katrina 2005" - python -m comp.validate --examples-only # just the example figure, from cache (fast) - python -m comp.validate --examples-only --refresh-cache # recompute the series first + python -m adforce.eval.validate # full sweep, all STORMS (populates the cache) + python -m adforce.eval.validate --storms "Ida 2021" "Katrina 2005" + python -m adforce.eval.validate --examples-only # just the example figure, from cache (fast) + python -m adforce.eval.validate --examples-only --refresh-cache # recompute the series first """ from __future__ import annotations @@ -600,7 +600,7 @@ def plot_failures(n_panels: int = 6, refresh: bool = False) -> None: """ csv = os.path.join(C.OUT_PATH, "val_summary.csv") if not os.path.exists(csv): - raise SystemExit(f"{csv} not found: run `python -m comp.validate` first") + raise SystemExit(f"{csv} not found: run `python -m adforce.eval.validate` first") df = pd.read_csv(csv) df["sid"] = df["sid"].astype(str) # gauge coordinates over both region boxes diff --git a/docs/comp.rst b/docs/comp.rst deleted file mode 100644 index 27e7071b..00000000 --- a/docs/comp.rst +++ /dev/null @@ -1,53 +0,0 @@ -comp package -============ - -Submodules ----------- - -comp.constants module ---------------------- - -.. automodule:: comp.constants - :members: - :undoc-members: - :show-inheritance: - -comp.coops module ------------------ - -.. automodule:: comp.coops - :members: - :undoc-members: - :show-inheritance: - -comp.nulltest module --------------------- - -.. automodule:: comp.nulltest - :members: - :undoc-members: - :show-inheritance: - -comp.sensitivity module ------------------------ - -.. automodule:: comp.sensitivity - :members: - :undoc-members: - :show-inheritance: - -comp.validate module --------------------- - -.. automodule:: comp.validate - :members: - :undoc-members: - :show-inheritance: - -Module contents ---------------- - -.. automodule:: comp - :members: - :undoc-members: - :show-inheritance: diff --git a/docs/conf.py b/docs/conf.py index 58b2aa52..b100ce69 100755 --- a/docs/conf.py +++ b/docs/conf.py @@ -71,7 +71,7 @@ # numba - cle15.cle15n # (not in install_requires; remove from this list if added there) # Not mocked because they are only imported inside functions and so cannot -# break the docs build: utide, huggingface_hub (comp), mpi4py, oct2py. +# break the docs build: utide, huggingface_hub (adforce.eval), mpi4py, oct2py. # Not mocked because they are core install_requires: imageio, datatree # (xarray-datatree), slurmpy, tcpyPI (tcpypi), hydra/omegaconf. autodoc_mock_imports = [ diff --git a/docs/eval.rst b/docs/eval.rst new file mode 100644 index 00000000..76d5449f --- /dev/null +++ b/docs/eval.rst @@ -0,0 +1,53 @@ +adforce.eval package +==================== + +Submodules +---------- + +adforce.eval.constants module +----------------------------- + +.. automodule:: adforce.eval.constants + :members: + :undoc-members: + :show-inheritance: + +adforce.eval.coops module +------------------------- + +.. automodule:: adforce.eval.coops + :members: + :undoc-members: + :show-inheritance: + +adforce.eval.nulltest module +---------------------------- + +.. automodule:: adforce.eval.nulltest + :members: + :undoc-members: + :show-inheritance: + +adforce.eval.sensitivity module +------------------------------- + +.. automodule:: adforce.eval.sensitivity + :members: + :undoc-members: + :show-inheritance: + +adforce.eval.validate module +---------------------------- + +.. automodule:: adforce.eval.validate + :members: + :undoc-members: + :show-inheritance: + +Module contents +--------------- + +.. automodule:: adforce.eval + :members: + :undoc-members: + :show-inheritance: diff --git a/docs/index.rst b/docs/index.rst index 3b6c007a..17f92fa4 100755 --- a/docs/index.rst +++ b/docs/index.rst @@ -23,7 +23,7 @@ at the University of Cambridge. cle15 adforce adbo - comp + eval worst surgenet about diff --git a/setup.py b/setup.py index 3a638515..68f06203 100755 --- a/setup.py +++ b/setup.py @@ -2,7 +2,7 @@ from typing import Dict, List # Core dependencies: what the pure-python, paper-relevant modules -# (cle15, worst, comp core, w22) need from a bare `pip install .`. +# (cle15, worst, adforce.eval core, w22) need from a bare `pip install .`. # Heavy stacks (Bayesian optimisation, CMIP6 ingestion, MPI, surge # comparison data) are split into the extras_require lists below. REQUIRED: List[str] = [ @@ -20,7 +20,7 @@ # turn imported at module level by w22.ps_runs/w22.stats2 (paper-relevant); # xarray[complete] would pull it in anyway. "dask[complete]", # to process netCDF4 files lazily (adforce, tcpips). - "scipy", # nearest-node KD-tree for gauge matching (comp) + "scipy", # nearest-node KD-tree for gauge matching (adforce.eval) "uncertainties", # common utility for linear error propagation # "sithom @ git+https://github.com/sdat2/sithom", "sithom >= 0.1.1", # personal common utilities for timing, plotting, and fitting @@ -69,20 +69,21 @@ "dask_jobqueue", # SLURMCluster in tcpips/run_dask_calculation.py ] -# "comp": historical surge validation against tide gauges (comp). -COMP: List[str] = [ - "huggingface_hub", # to download historical surge datasets (comp) - "utide", # to de-tide tide-gauge records for surge validation (comp) - "pyarrow", # Parquet cache of de-tided gauge time series (comp) +# "eval": historical surge validation against tide gauges (adforce.eval). +EVAL: List[str] = [ + "huggingface_hub", # to download historical surge datasets (adforce.eval) + "utide", # to de-tide tide-gauge records for surge validation (adforce.eval) + "pyarrow", # Parquet cache of de-tided gauge time series (adforce.eval) ] EXTRAS: Dict[str, List[str]] = { "bo": BO, "cmip": CMIP, "mpi": MPI, - "comp": COMP, + "eval": EVAL, + "comp": EVAL, # deprecated alias (comp/ moved to adforce/eval/) } -EXTRAS["all"] = BO + CMIP + MPI + COMP # union of all extras +EXTRAS["all"] = BO + CMIP + MPI + EVAL # union of all extras setup( @@ -119,8 +120,14 @@ "adbo", "worst", "cle15", - "comp", + "adforce.eval", ], + # include_package_data=True only ships files named in a MANIFEST.in (absent); + # package_data is what actually gets the hydra YAML trees into a wheel. + package_data={ + "adforce": ["config/*.yaml", "config/*/*.yaml"], + "adforce.eval": ["config/*.yaml", "config/*/*.yaml", "README.md"], + }, package_dir={ "tcpips": "tcpips", # Calculate potential intensity and prerequisites for potential size "w22": "w22", # Calculate the Chavas, Lin and Emanuel (2015) profile using matlab (octave), calculate potential size @@ -129,6 +136,6 @@ "adbo": "adbo", # All of the tensorflow/trieste Bayesian optimization stuff "worst": "worst", # Extreme value theory using the upper bound limit using tensorflow for fitting. "cle15": "cle15", # Chavas, Lin & Emanuel (2015) TC wind profile implementations - "comp": "comp", # Compare historical ADCIRC surge against de-tided NOAA tide gauges + "adforce.eval": "adforce/eval", # Compare ADCIRC surge configurations against de-tided NOAA tide gauges }, ) diff --git a/tests/test_comp.py b/tests/test_eval.py similarity index 97% rename from tests/test_comp.py rename to tests/test_eval.py index fcfa88d5..15f93979 100644 --- a/tests/test_comp.py +++ b/tests/test_eval.py @@ -1,6 +1,6 @@ -"""Unit + regression tests for the comp (tide-gauge validation) module. +"""Unit + regression tests for the adforce.eval (tide-gauge validation) module. -Covers the pure, network-free functions in comp.validate (skill metrics, time-series +Covers the pure, network-free functions in adforce.eval.validate (skill metrics, time-series alignment, the valid/clean gating, the LaTeX-table generator) and pins the headline numbers from a completed sweep so they cannot silently drift. The regression test is skipped when the summary CSV is absent (e.g. a fresh checkout with no cached data), so @@ -19,8 +19,8 @@ if str(REPO_ROOT) not in sys.path: sys.path.insert(0, str(REPO_ROOT)) -from comp import constants as C -from comp.validate import ( +from adforce.eval import constants as C +from adforce.eval.validate import ( add_flags, bootstrap_ci, classify_setting, @@ -313,7 +313,7 @@ def test_regression_per_storm_counts(): reason="val_summary.csv not present (no cached sweep)", ) def test_global_permutation_separates_signal_from_null(): - from comp.nulltest import perm_global, load_clean + from adforce.eval.nulltest import perm_global, load_clean res = perm_global(load_clean(), n=1000, seed=0) # the real correlation must sit far above every shuffled (null) correlation @@ -328,7 +328,7 @@ def test_global_permutation_separates_signal_from_null(): reason="val_summary.csv not present (no cached sweep)", ) def test_within_storm_permutation_shows_real_spatial_skill(): - from comp.nulltest import perm_within_storm, load_clean + from adforce.eval.nulltest import perm_within_storm, load_clean res = perm_within_storm(load_clean(), n=1000, seed=0) # spatial skill (storm means removed) beats the within-storm-shuffled null @@ -341,7 +341,7 @@ def test_within_storm_permutation_shows_real_spatial_skill(): reason="val_summary.csv not present (no cached sweep)", ) def test_cross_storm_pairing_collapses_skill(): - from comp.nulltest import cross_storm, load_clean + from adforce.eval.nulltest import cross_storm, load_clean res = cross_storm(load_clean(), seed=0) # pairing a gauge's sim with a DIFFERENT storm's obs must lose most of the skill @@ -488,7 +488,7 @@ def test_regression_population_and_methods(): # Network-free -- uses synthetic ragged (sim, obs) series in a tmp cache dir. # --------------------------------------------------------------------------- # def test_series_cache_roundtrip_and_staleness(tmp_path, monkeypatch): - import comp.validate as cv + import adforce.eval.validate as cv monkeypatch.setattr(C, "TS_CACHE", str(tmp_path)) # isolate from the real cache h = pd.date_range("2021-08-27", periods=48, freq="h") diff --git a/tests/test_gauge_evt.py b/tests/test_gauge_evt.py index 81bb7744..2f99b457 100644 --- a/tests/test_gauge_evt.py +++ b/tests/test_gauge_evt.py @@ -2,7 +2,7 @@ Two layers are covered without touching the CO-OPS API: -1. ``comp.annual_max``: a synthetic hourly water-level series (known harmonic +1. ``adforce.eval.annual_max``: a synthetic hourly water-level series (known harmonic tide + trend + one dominant storm-surge bump per year) is pushed through the de-tide + annual-maximum pipeline with the fetch monkeypatched. The recovered annual maxima must match the injected storm amplitudes, the @@ -30,8 +30,8 @@ if str(REPO_ROOT) not in sys.path: sys.path.insert(0, str(REPO_ROOT)) -import comp.annual_max as am -from comp import constants as C +import adforce.eval.annual_max as am +from adforce.eval import constants as C from worst.gauge_fit import ( POTENTIAL_HEIGHT_M, bootstrap_bands, @@ -68,7 +68,7 @@ def _synth_year(year: int, storm_peak: float, seed: int) -> pd.Series: @pytest.fixture() def synth_pipeline(monkeypatch, tmp_path): - """Route comp.annual_max at the synthetic gauge and an isolated cache.""" + """Route adforce.eval.annual_max at the synthetic gauge and an isolated cache.""" series = {y: _synth_year(y, p, seed=y) for y, p in STORM_PEAKS.items()} calls = {"fetch": 0} diff --git a/worst/gauge_fit.py b/worst/gauge_fit.py index 5b0c2b03..6435f948 100644 --- a/worst/gauge_fit.py +++ b/worst/gauge_fit.py @@ -2,7 +2,7 @@ Closes the loop between the paper's Bayesian-optimization potential-height results and its EVT theory: the annual maxima of de-tided NOAA tide-gauge -surge residuals (from :mod:`comp.annual_max`) are fitted with +surge residuals (from :mod:`adforce.eval.annual_max`) are fitted with * case I -- Weibull-class GEV with the upper bound FIXED at the site's potential height ``z*`` (``worst.tens.fit_gev_upper_bound_known``), and @@ -349,7 +349,7 @@ def analyze_station( Returns: dict: The results dictionary that is also written to json. """ - from comp.annual_max import annual_maxima + from adforce.eval.annual_max import annual_maxima if station not in POTENTIAL_HEIGHT_M: raise ValueError( diff --git a/worst/gauge_fit_ns.py b/worst/gauge_fit_ns.py index 60069ec8..cead2107 100644 --- a/worst/gauge_fit_ns.py +++ b/worst/gauge_fit_ns.py @@ -142,7 +142,7 @@ def analyze_station_ns( method: str = "robust", make_figure: bool = True, ) -> dict: - from comp.annual_max import annual_maxima + from adforce.eval.annual_max import annual_maxima df = annual_maxima(station, start, end, method=method) if df.empty: From 54d9d1ae56e9825d10e834037cd1aa72c523b1f8 Mon Sep 17 00:00:00 2001 From: "Simon D.A. Thomas" Date: Sun, 16 Aug 2026 19:18:20 +0100 Subject: [PATCH 02/17] appendix --- REPRODUCE.md | 2 +- adforce/eval/README.md | 14 +++---- adforce/eval/annual_max.py | 51 ++++++++++++----------- adforce/eval/constants.py | 36 ++++++++++------ adforce/eval/coops.py | 3 ++ adforce/eval/detide_sensitivity.py | 54 ++++++++++++------------ adforce/eval/nulltest.py | 26 +++++++----- adforce/eval/sensitivity.py | 24 +++++++---- adforce/eval/validate.py | 66 ++++++++++++++---------------- tests/test_eval.py | 40 ++++++++++++++++++ 10 files changed, 190 insertions(+), 126 deletions(-) diff --git a/REPRODUCE.md b/REPRODUCE.md index c4e564ed..548a1328 100644 --- a/REPRODUCE.md +++ b/REPRODUCE.md @@ -114,7 +114,7 @@ Some plotting entry points are driven from commented-out calls in the module | **Table 3** — de-tiding-method sensitivity (`tab:detide-robustness`) | `python -m adforce.eval.detide_sensitivity` | de-tided gauge cache (`data/comp/ts_cache/`) | `/paper/comp_detide_table.tex` (`\input` by the appendix) | — (local) | | **Table 4** — per-storm skill vs de-tided NOAA gauges (`tab:gauge-validation`) | `python -m adforce.eval.validate` | Hugging Face [`sdat2/surgenet-train`](https://huggingface.co/datasets/sdat2/surgenet-train); NOAA CO-OPS gauges (de-tided with `utide`) | `/paper/comp_val_table.tex` | — (local) | | **Fig. 12** — peak-skill scatter (`fig:gauge-validation-scatter`) | `python -m adforce.eval.validate` (full 14-storm sweep) | as above | `img/comp/val_scatter.png` · `/img/comp_val_scatter.pdf` | — (local) | -| **Fig. 13** — example de-tided surge time series (`fig:gauge-validation-examples`) | `python -m adforce.eval.validate` (or `--examples-only` from cache) | as above (cache: `data/comp/ts_cache/`) | `img/comp/val_examples.png` · `/img/comp_val_examples.pdf` | — (local) | +| **Fig. 13** — example de-tided surge time series (`fig:gauge-validation-examples`) | `python -m adforce.eval.validate` (or `validate.examples_only=true` from cache) | as above (cache: `data/comp/ts_cache/`) | `img/comp/val_examples.png` · `/img/comp_val_examples.pdf` | — (local) | | **Fig. 14** — negative-control null tests (`fig:gauge-validation-null`) | `python -m adforce.eval.nulltest` | `data/comp/out/val_summary.csv` (+ netCDFs for the lag null) | `/img/comp_val_nulltests.pdf` · `img/comp/val_nulltests.png` | — (local) | | Summary metrics (support; pinned by `tests/test_eval.py`) | `python -m adforce.eval.validate` | as above | `data/comp/out/val_summary.csv` | — (local) | | Threshold/node-selection sensitivity (text numbers) | `python -m adforce.eval.sensitivity` | `val_summary.csv` | console report (headline r stability) | — (local) | diff --git a/adforce/eval/README.md b/adforce/eval/README.md index 0fc64a92..60e2df72 100644 --- a/adforce/eval/README.md +++ b/adforce/eval/README.md @@ -36,9 +36,9 @@ For each storm: ```bash python -m adforce.eval.validate # full 19-storm sweep (Gulf + Florida) -python -m adforce.eval.validate --storms "Ida 2021" # one storm (skips example panels + table) -python -m adforce.eval.validate --examples-only # just the example figure, from cache (fast) -python -m adforce.eval.validate --examples-only --refresh-cache # recompute the cached series first +python -m adforce.eval.validate 'storms=["Ida 2021"]' # one storm (skips example panels + table) +python -m adforce.eval.validate validate.examples_only=true # example figure, from cache (fast) +python -m adforce.eval.validate validate.examples_only=true validate.refresh=true # recompute first ``` A full sweep regenerates **everything the paper uses**, in one step, so the figures, the @@ -54,8 +54,8 @@ table and the prose cannot drift apart: The slow step is the per-gauge `utide` de-tiding. A full sweep **caches** each storm's de-tided `(sim, obs)` series as Parquet under `data/comp/ts_cache/` (write-through), keyed by the node-selection + de-tiding parameters so the cache self-invalidates if any of those change. -`--examples-only` then re-renders `comp_val_examples.pdf` from that cache in seconds (vs minutes) -— use it to iterate on the figure's layout without re-detiding; `--refresh-cache` forces a recompute. +`validate.examples_only=true` then re-renders `comp_val_examples.pdf` from that cache in seconds (vs minutes) +— use it to iterate on the figure's layout without re-detiding; `validate.refresh=true` forces a recompute. The thesis tree is located by searching for `paper/appendix.tex`; override with the `WORSTSURGE_PAPER_ROOT` env var. Downloads/caches live under `data/comp/` (git-ignored). @@ -66,9 +66,9 @@ The thesis tree is located by searching for `paper/appendix.tex`; override with ```bash python -m adforce.eval.nulltest # permutation, cross-storm, and temporal-lag nulls -python -m adforce.eval.nulltest --no-lag # peak-level nulls only (no netCDF) +python -m adforce.eval.nulltest nulltest.lag=false # peak-level nulls only (no netCDF) python -m adforce.eval.sensitivity # threshold + node-selection robustness -python -m adforce.eval.sensitivity --no-node +python -m adforce.eval.sensitivity sensitivity.node=false ``` `adforce.eval.nulltest` writes `img/comp_val_nulltests.pdf` and reports: a label-permutation null diff --git a/adforce/eval/annual_max.py b/adforce/eval/annual_max.py index 4581a649..c04796d4 100644 --- a/adforce/eval/annual_max.py +++ b/adforce/eval/annual_max.py @@ -26,20 +26,22 @@ insensitive to the datum because the harmonic fit absorbs the mean (and a linear trend, which also removes most local sea-level rise within a year). -Run:: +Run (hydra overrides; config root adforce/eval/config/annual_max_config.yaml):: - python -m adforce.eval.annual_max --station 8761724 # Grand Isle, LA - python -m adforce.eval.annual_max --station 8735180 --start 1980 # Dauphin Island, AL + python -m adforce.eval.annual_max station=8761724 # Grand Isle, LA + python -m adforce.eval.annual_max station=8735180 start=1980 # Dauphin Island, AL """ from __future__ import annotations -import argparse import json import os import time from typing import Optional +import hydra +from omegaconf import DictConfig + import numpy as np import pandas as pd @@ -260,6 +262,7 @@ def year_residual( Returns: Optional[pd.Series]: Residual series, or None if the year is skipped. """ + C.ensure_dirs() # cache dirs are created lazily, not at import fp = _resid_path(station, year, method) if not refresh and os.path.exists(fp): df = pd.read_parquet(fp) @@ -305,6 +308,7 @@ def annual_maxima( ``station, name, year, ann_max_m, t_max, n_obs, coverage, max_at_gap_edge``. """ + C.ensure_dirs() # cache dirs are created lazily, not at import tp = _table_path(station, start, end, method) if not refresh and os.path.exists(tp): df = pd.read_parquet(tp) @@ -351,31 +355,27 @@ def annual_maxima( return df -def main() -> None: - ap = argparse.ArgumentParser(description=__doc__) - ap.add_argument( - "--station", - default="8761724", - help="CO-OPS station id (default 8761724, Grand Isle LA)", - ) - ap.add_argument("--start", type=int, default=C.AM_START_YEAR) - ap.add_argument("--end", type=int, default=C.AM_END_YEAR) - ap.add_argument( - "--method", - default="robust", - choices=["robust", "ols"], - help="utide.solve method (default robust)", - ) - ap.add_argument( - "--refresh", action="store_true", help="recompute instead of reading the caches" +_LEGACY_FLAGS = { + "--station": "station=8761724", + "--start": "start=1980", + "--end": "end=2025", + "--method": "method=robust", + "--refresh": "refresh=true", +} + + +@hydra.main(version_base=None, config_path="config", config_name="annual_max_config") +def main(cfg: DictConfig) -> None: + if cfg.method not in ("robust", "ols"): + raise SystemExit(f"method must be robust|ols, got {cfg.method!r}") + df = annual_maxima( + str(cfg.station), cfg.start, cfg.end, method=cfg.method, refresh=cfg.refresh ) - a = ap.parse_args() - df = annual_maxima(a.station, a.start, a.end, method=a.method, refresh=a.refresh) if df.empty: print("no usable years") return print( - f"\n{len(df)} usable years of {a.end - a.start + 1}; " + f"\n{len(df)} usable years of {cfg.end - cfg.start + 1}; " f"max residual {df.ann_max_m.max():.2f} m in " f"{int(df.loc[df.ann_max_m.idxmax(), 'year'])}; " f"{int(df.max_at_gap_edge.sum())} year(s) flagged max-near-gap" @@ -383,4 +383,7 @@ def main() -> None: if __name__ == "__main__": + from ._cli import reject_legacy_flags + + reject_legacy_flags(_LEGACY_FLAGS, "adforce.eval.annual_max") main() diff --git a/adforce/eval/constants.py b/adforce/eval/constants.py index bd333065..89276f50 100644 --- a/adforce/eval/constants.py +++ b/adforce/eval/constants.py @@ -47,17 +47,28 @@ def _find_paper_root(): PAPER_IMG_PATH = str(PAPER_ROOT / "img") if PAPER_ROOT else FIGURE_PATH PAPER_TEX_PATH = str(PAPER_ROOT / "paper") if PAPER_ROOT else OUT_PATH -for _p in ( - COMP_DATA_PATH, - HF_STORM_CACHE, - COOPS_CACHE, - TS_CACHE, - FIGURE_PATH, - OUT_PATH, - PAPER_IMG_PATH, - PAPER_TEX_PATH, -): - os.makedirs(_p, exist_ok=True) + +def ensure_dirs() -> None: + """Create every cache/output directory (idempotent, cheap). + + Called lazily by the entry points and the low-level writers rather than at + import time: importing adforce.eval must stay side-effect-free, because + ``pytest --doctest-modules`` imports it on clean CI machines and the + launch/orchestration path will import it on HPC nodes (where creating + thesis-tree directories would be wrong). + """ + for _p in ( + COMP_DATA_PATH, + HF_STORM_CACHE, + COOPS_CACHE, + TS_CACHE, + FIGURE_PATH, + OUT_PATH, + PAPER_IMG_PATH, + PAPER_TEX_PATH, + ANNUAL_MAX_CACHE, + ): + os.makedirs(_p, exist_ok=True) # --- data sources ---------------------------------------------------------- # Historical ADCIRC simulations (228 IBTrACS N-Atlantic landfalling TCs, 1980-2024, @@ -174,8 +185,7 @@ def box_for(storm: str) -> dict: # --- annual-maximum pipeline (adforce.eval.annual_max) ------------------------------ ANNUAL_MAX_CACHE = os.path.join( COMP_DATA_PATH, "annual_max" -) # residuals + annual maxima -os.makedirs(ANNUAL_MAX_CACHE, exist_ok=True) +) # residuals + annual maxima (created lazily via ensure_dirs) AM_START_YEAR = 1980 # default year range for annual-maximum extraction AM_END_YEAR = 2025 # A year only contributes an annual maximum if at least this fraction of its diff --git a/adforce/eval/coops.py b/adforce/eval/coops.py index def0cd28..6942496c 100644 --- a/adforce/eval/coops.py +++ b/adforce/eval/coops.py @@ -24,6 +24,7 @@ COOPS_MDAPI, GAUGE_BOX, UTIDE_MIN_SAMPLES, + ensure_dirs, ) Gauge = Tuple[str, str, float, float] # (station_id, name, lat, lon) @@ -31,6 +32,7 @@ def gulf_gauges(box: dict = GAUGE_BOX) -> List[Gauge]: """Return CO-OPS water-level stations within ``box`` (cached metadata).""" + ensure_dirs() # cache dirs are created lazily, not at import cache = os.path.join(COOPS_CACHE, "stations.json") if os.path.exists(cache): stations = json.load(open(cache)) @@ -51,6 +53,7 @@ def _coops( station: str, begin: str, end: str, product: str, datum: str = "MSL" ) -> pd.Series: """Fetch one hourly CO-OPS series (GMT, metric), cached on disk by request.""" + ensure_dirs() # cache dirs are created lazily, not at import key = f"{station}_{product}_{begin}_{end}_{datum}.csv" fp = os.path.join(COOPS_CACHE, key) if os.path.exists(fp): diff --git a/adforce/eval/detide_sensitivity.py b/adforce/eval/detide_sensitivity.py index 209ceb34..0d371241 100644 --- a/adforce/eval/detide_sensitivity.py +++ b/adforce/eval/detide_sensitivity.py @@ -14,20 +14,22 @@ ``tabular`` that ``tab:detide-robustness`` ``\\input``s), so the paper table regenerates with the data and cannot drift -- like the per-storm table from :func:`adforce.eval.validate.latex_table`. -Run:: +Run (hydra overrides; config root adforce/eval/config/eval_config.yaml):: - python -m adforce.eval.detide_sensitivity # all methods + skew, all pairs, writes the table - python -m adforce.eval.detide_sensitivity --methods godin_lowpass noaa_predictions - python -m adforce.eval.detide_sensitivity --limit 40 # quick subset (no table written) - python -m adforce.eval.detide_sensitivity --skew # instantaneous vs skew-surge metric only + python -m adforce.eval.detide_sensitivity # all methods + skew, all pairs, writes the table + python -m adforce.eval.detide_sensitivity 'detide_sensitivity.methods=[godin_lowpass,noaa_predictions]' + python -m adforce.eval.detide_sensitivity detide_sensitivity.limit=40 # quick subset (no table) + python -m adforce.eval.detide_sensitivity detide_sensitivity.skew=true # skew-surge metric only """ from __future__ import annotations -import argparse import os from typing import Dict, Optional, Tuple +import hydra +from omegaconf import DictConfig + import numpy as np import pandas as pd import xarray as xr @@ -400,30 +402,28 @@ def write_table_from_csv(path: Optional[str] = None) -> None: ) -def main() -> None: - ap = argparse.ArgumentParser(description=__doc__) - ap.add_argument("--methods", nargs="*", default=None, help="subset of method ids") - ap.add_argument( - "--limit", type=int, default=None, help="only the first N clean pairs (quick)" - ) - ap.add_argument( - "--skew", - action="store_true", - help="instantaneous-residual vs skew-surge peak metric (phase-insensitive)", - ) - ap.add_argument( - "--table-only", - action="store_true", - help="regenerate the appendix table from the cached CSV (no re-detiding)", - ) - a = ap.parse_args() - if a.table_only: +_LEGACY_FLAGS = { + "--methods": "'detide_sensitivity.methods=[godin_lowpass]'", + "--limit": "detide_sensitivity.limit=40", + "--skew": "detide_sensitivity.skew=true", + "--table-only": "detide_sensitivity.table_only=true", +} + + +@hydra.main(version_base=None, config_path="config", config_name="eval_config") +def main(cfg: DictConfig) -> None: + C.ensure_dirs() + d = cfg.detide_sensitivity + if d.table_only: write_table_from_csv() - elif a.skew: - run_skew(a.limit) + elif d.skew: + run_skew(d.limit) else: - run(a.methods, a.limit) + run(list(d.methods) if d.methods else None, d.limit) if __name__ == "__main__": + from ._cli import reject_legacy_flags + + reject_legacy_flags(_LEGACY_FLAGS, "adforce.eval.detide_sensitivity") main() diff --git a/adforce/eval/nulltest.py b/adforce/eval/nulltest.py index 2decee43..5538b36f 100644 --- a/adforce/eval/nulltest.py +++ b/adforce/eval/nulltest.py @@ -19,19 +19,21 @@ recompute the time-series correlation; real skill peaks sharply at lag 0. -Run:: +Run (hydra overrides; config root adforce/eval/config/eval_config.yaml):: - python -m adforce.eval.nulltest # peak-level nulls + lag curve (+ figure) - python -m adforce.eval.nulltest --no-lag # skip the netCDF-heavy lag test + python -m adforce.eval.nulltest # peak-level nulls + lag curve (+ figure) + python -m adforce.eval.nulltest nulltest.lag=false # skip the netCDF-heavy lag test """ from __future__ import annotations -import argparse import os import warnings from typing import Dict, List, Optional, Tuple +import hydra +from omegaconf import DictConfig + import numpy as np import pandas as pd @@ -215,6 +217,7 @@ def plot(g: Dict, w: Dict, lag: Optional[pd.DataFrame], paths: List[str]) -> Non def run(do_lag: bool = True) -> None: + C.ensure_dirs() c = load_clean() print(f"clean pairs: n={len(c)}\n") g = perm_global(c) @@ -249,13 +252,16 @@ def run(do_lag: bool = True) -> None: ) -def main() -> None: - ap = argparse.ArgumentParser(description=__doc__) - ap.add_argument( - "--no-lag", action="store_true", help="skip the netCDF-heavy lag test" - ) - run(do_lag=not ap.parse_args().no_lag) +_LEGACY_FLAGS = {"--no-lag": "nulltest.lag=false"} + + +@hydra.main(version_base=None, config_path="config", config_name="eval_config") +def main(cfg: DictConfig) -> None: + run(do_lag=cfg.nulltest.lag) if __name__ == "__main__": + from ._cli import reject_legacy_flags + + reject_legacy_flags(_LEGACY_FLAGS, "adforce.eval.nulltest") main() diff --git a/adforce/eval/sensitivity.py b/adforce/eval/sensitivity.py index 63335ce2..1b76a8af 100644 --- a/adforce/eval/sensitivity.py +++ b/adforce/eval/sensitivity.py @@ -11,19 +11,21 @@ sampled and therefore the simulated peak, so they require the cached storm netCDFs. Each storm is loaded ONCE and re-sampled for every combination. -Run:: +Run (hydra overrides; config root adforce/eval/config/eval_config.yaml):: - python -m adforce.eval.sensitivity # cheap filter sweep + node-selection sweep - python -m adforce.eval.sensitivity --no-node # cheap filter sweep only (no netCDF) + python -m adforce.eval.sensitivity # filter sweep + node-selection sweep + python -m adforce.eval.sensitivity sensitivity.node=false # cheap filter sweep only (no netCDF) """ from __future__ import annotations -import argparse import os import warnings from typing import Dict, List, Optional, Tuple +import hydra +from omegaconf import DictConfig + import numpy as np import pandas as pd import xarray as xr @@ -166,6 +168,7 @@ def node_sweep( def run(do_node: bool = True) -> None: + C.ensure_dirs() fs = filter_sweep() print("=== clean-filter sweep (from val_summary.csv) ===") print(fs.to_string(index=False, float_format=lambda v: f"{v:.3f}")) @@ -177,11 +180,16 @@ def run(do_node: bool = True) -> None: ns.to_csv(os.path.join(C.OUT_PATH, "sensitivity_node.csv"), index=False) -def main() -> None: - ap = argparse.ArgumentParser(description=__doc__) - ap.add_argument("--no-node", action="store_true", help="skip the netCDF node sweep") - run(do_node=not ap.parse_args().no_node) +_LEGACY_FLAGS = {"--no-node": "sensitivity.node=false"} + + +@hydra.main(version_base=None, config_path="config", config_name="eval_config") +def main(cfg: DictConfig) -> None: + run(do_node=cfg.sensitivity.node) if __name__ == "__main__": + from ._cli import reject_legacy_flags + + reject_legacy_flags(_LEGACY_FLAGS, "adforce.eval.sensitivity") main() diff --git a/adforce/eval/validate.py b/adforce/eval/validate.py index 7483d544..2efbc5ba 100644 --- a/adforce/eval/validate.py +++ b/adforce/eval/validate.py @@ -18,20 +18,22 @@ then regenerates the example-panel figure from that cache without re-running the (slow) utide de-tiding. -Run:: +Run (hydra overrides; config root adforce/eval/config/eval_config.yaml):: python -m adforce.eval.validate # full sweep, all STORMS (populates the cache) - python -m adforce.eval.validate --storms "Ida 2021" "Katrina 2005" - python -m adforce.eval.validate --examples-only # just the example figure, from cache (fast) - python -m adforce.eval.validate --examples-only --refresh-cache # recompute the series first + python -m adforce.eval.validate 'storms=["Ida 2021","Katrina 2005"]' + python -m adforce.eval.validate validate.examples_only=true # example figure, from cache (fast) + python -m adforce.eval.validate validate.examples_only=true validate.refresh=true """ from __future__ import annotations -import argparse import os from typing import Dict, List, Optional, Tuple +import hydra +from omegaconf import DictConfig + import numpy as np import pandas as pd import xarray as xr @@ -146,6 +148,7 @@ def _series_cache_path(storm: str) -> str: def _save_series_cache(storm: str, series: Dict[str, tuple]) -> None: try: + C.ensure_dirs() # cache dirs are created lazily, not at import frames = [] for name, (sim, obs) in series.items(): for kind, s in (("sim", sim), ("obs", obs)): @@ -531,6 +534,7 @@ def latex_table(df: pd.DataFrame, path: str) -> None: def run(storms: Optional[List[str]] = None) -> pd.DataFrame: + C.ensure_dirs() items = {k: C.STORMS[k] for k in (storms or C.STORMS)} print( f"{len(gulf_gauges())} candidate gauges in Gulf box, " @@ -658,47 +662,37 @@ def plot_failures(n_panels: int = 6, refresh: bool = False) -> None: ) -def main() -> None: - ap = argparse.ArgumentParser(description=__doc__) - ap.add_argument( - "--storms", nargs="*", default=None, help="subset of storm names (default: all)" - ) - ap.add_argument( - "--examples-only", - action="store_true", - help="regenerate only the example-panel figure from the cached " - "time series (no sweep) -- instant once the cache exists", - ) - ap.add_argument( - "--refresh-cache", - action="store_true", - help="force-recompute the cached time series instead of reading it", - ) - ap.add_argument( - "--failures", - action="store_true", - help="per-city worst-case panels (largest |sim-obs| peak difference) " - "from the cached sweep -- no re-run", - ) - ap.add_argument( - "--n-failures", type=int, default=6, help="panels per city for --failures" - ) - a = ap.parse_args() - if a.failures: - plot_failures(n_panels=a.n_failures, refresh=a.refresh_cache) +_LEGACY_FLAGS = { + "--storms": "'storms=[\"Ida 2021\"]'", + "--examples-only": "validate.examples_only=true", + "--refresh-cache": "validate.refresh=true", + "--failures": "validate.failures=true", + "--n-failures": "validate.n_failures=6", +} + + +@hydra.main(version_base=None, config_path="config", config_name="eval_config") +def main(cfg: DictConfig) -> None: + C.ensure_dirs() + v = cfg.validate + if v.failures: + plot_failures(n_panels=v.n_failures, refresh=v.refresh) return - if a.examples_only: + if v.examples_only: plot_examples( C.EXAMPLE_PANELS, [ os.path.join(C.FIGURE_PATH, "val_examples.png"), os.path.join(C.PAPER_IMG_PATH, "comp_val_examples.pdf"), ], - refresh=a.refresh_cache, + refresh=v.refresh, ) return - run(a.storms) + run(list(cfg.storms) if cfg.storms else None) if __name__ == "__main__": + from ._cli import reject_legacy_flags + + reject_legacy_flags(_LEGACY_FLAGS, "adforce.eval.validate") main() diff --git a/tests/test_eval.py b/tests/test_eval.py index 15f93979..03a1409d 100644 --- a/tests/test_eval.py +++ b/tests/test_eval.py @@ -516,3 +516,43 @@ def test_series_cache_roundtrip_and_staleness(tmp_path, monkeypatch): monkeypatch.setattr(cv, "_ts_cache_tag", lambda: "DIFFERENT-PARAMS") assert cv._load_series_cache("Test 2021") is None # stale tag -> recompute assert cv._load_series_cache("Never Cached 1999") is None # missing -> recompute + + +# --------------------------------------------------------------------------- # +# Hydra config vs constants anti-drift, and cache-tag byte-stability. +# The scoring YAML mirrors constants.py; the ts-cache tag string keys the +# utide-expensive Parquet caches, so its byte-stability IS the cache validity. +# --------------------------------------------------------------------------- # +def test_scoring_config_matches_constants(): + from hydra import compose, initialize + + with initialize(version_base=None, config_path="../adforce/eval/config"): + cfg = compose(config_name="eval_config") + am = compose(config_name="annual_max_config") + s = cfg.scoring + assert s.max_node_deg == C.MAX_NODE_DEG + assert s.wet_min_m == C.WET_MIN_M + assert s.knn == C.KNN + assert s.max_timing_hr == C.MAX_TIMING_HR + assert s.min_obs_peak_m == C.MIN_OBS_PEAK_M + assert s.ts_min_overlap == C.TS_MIN_OVERLAP + assert s.n_bootstrap == C.N_BOOTSTRAP + assert s.bootstrap_seed == C.BOOTSTRAP_SEED + assert s.utide_min_samples == C.UTIDE_MIN_SAMPLES + assert am.start == C.AM_START_YEAR + assert am.end == C.AM_END_YEAR + + +def test_cache_tag_matches_legacy_literal(): + """The tag must stay byte-identical to the comp/-era string: a formatting + slip would silently invalidate every cached utide fit and re-hammer CO-OPS.""" + from adforce.eval.validate import _ts_cache_tag + + assert _ts_cache_tag("Katrina 2005") == ( + "v1|deg0.12|wet0.3|knn60|ut2000" + "|box{'lon': (-97.6, -84.0), 'lat': (27.3, 30.9)}" + ) + assert _ts_cache_tag("Irma 2017") == ( # Florida-box storm + "v1|deg0.12|wet0.3|knn60|ut2000" + "|box{'lon': (-82.3, -79.7), 'lat': (24.4, 30.8)}" + ) From 9386b05f7e1f90873bed845a316617937bc937f6 Mon Sep 17 00:00:00 2001 From: "Simon D.A. Thomas" Date: Sun, 16 Aug 2026 19:23:34 +0100 Subject: [PATCH 03/17] hydra-ize adforce.eval: config tree, 5 mains, lazy dirs, drift tests - adforce/eval/config/: eval_config.yaml (shared root: validate/nulltest/ sensitivity/detide_sensitivity), scoring/default.yaml (mirrors constants, drift-tested), annual_max_config.yaml; hydra output_subdir null + run.dir . - five argparse mains -> @hydra.main; _cli.reject_legacy_flags prints the exact override translation for stale --flag invocations - constants: import-time makedirs loop -> lazy ensure_dirs() (pytest --doctest-modules and future HPC imports stay side-effect-free); called by mains and low-level cache writers (coops, ts-cache, annual_max) - tests: test_scoring_config_matches_constants, test_cache_tag_matches_ legacy_literal (pins the v1 ts-cache tag byte-for-byte) - docs/REPRODUCE: hydra override syntax Gate: 47 passed 1 skipped; doctest collection clean; wrap tests 19 passed; full sweep via hydra -> val_summary.csv byte-identical (md5 f689004e); Katrina parquet deterministic (63650a83). Co-Authored-By: Claude Fable 5 --- adforce/eval/_cli.py | 31 +++++++++++++++++++ adforce/eval/config/annual_max_config.yaml | 12 ++++++++ adforce/eval/config/eval_config.yaml | 36 ++++++++++++++++++++++ adforce/eval/config/scoring/default.yaml | 18 +++++++++++ 4 files changed, 97 insertions(+) create mode 100644 adforce/eval/_cli.py create mode 100644 adforce/eval/config/annual_max_config.yaml create mode 100644 adforce/eval/config/eval_config.yaml create mode 100644 adforce/eval/config/scoring/default.yaml diff --git a/adforce/eval/_cli.py b/adforce/eval/_cli.py new file mode 100644 index 00000000..e375a93c --- /dev/null +++ b/adforce/eval/_cli.py @@ -0,0 +1,31 @@ +"""Shared CLI helpers for the adforce.eval hydra entry points. + +The comp/-era argparse flags were replaced by hydra overrides when the module +moved to adforce/eval/. Stale invocations (old docs, shell history, thesis +Makefiles) should fail loudly with the exact translation rather than be +half-parsed by hydra -- so each entry point calls :func:`reject_legacy_flags` +with its flag map before handing argv to hydra. Only known legacy flags are +rejected; hydra's own ``--help``/``--cfg``/etc. pass through untouched. +""" + +import sys +from typing import Dict + + +def reject_legacy_flags(mapping: Dict[str, str], module: str) -> None: + """Exit with the hydra translation if an old argparse flag is on argv. + + Args: + mapping: legacy flag -> replacement hydra override (shown verbatim). + module: dotted module path for the usage line, e.g. "adforce.eval.validate". + """ + hits = [a for a in sys.argv[1:] if a.split("=")[0] in mapping] + if not hits: + return + lines = [ + f"{module}: argparse flags were replaced by hydra overrides " + "when comp/ moved to adforce/eval/. Translation:" + ] + lines += [f" {old:<18} -> {new}" for old, new in mapping.items()] + lines.append(f"e.g. python -m {module} " + mapping[hits[0].split("=")[0]]) + raise SystemExit("\n".join(lines)) diff --git a/adforce/eval/config/annual_max_config.yaml b/adforce/eval/config/annual_max_config.yaml new file mode 100644 index 00000000..00c78232 --- /dev/null +++ b/adforce/eval/config/annual_max_config.yaml @@ -0,0 +1,12 @@ +# Root config for python -m adforce.eval.annual_max (obs-side annual-maximum +# de-tided residual pipeline; feeds the gauge-EVT fits in worst/gauge_fit*.py). +hydra: + output_subdir: null + run: + dir: . + +station: "8761724" # CO-OPS station id (default: Grand Isle, LA) +start: 1980 # first year; = C.AM_START_YEAR (drift-tested) +end: 2025 # last year; = C.AM_END_YEAR (drift-tested) +method: robust # utide.solve method: robust | ols +refresh: false # recompute instead of reading the Parquet caches diff --git a/adforce/eval/config/eval_config.yaml b/adforce/eval/config/eval_config.yaml new file mode 100644 index 00000000..5aee5726 --- /dev/null +++ b/adforce/eval/config/eval_config.yaml @@ -0,0 +1,36 @@ +# Root config for the adforce.eval analysis entry points +# (validate, nulltest, sensitivity, detide_sensitivity). +# +# House style: flat keys + comments (infrastructure/analysis config, like +# adforce/config/slurm|files); the {value, description, units} nesting is +# reserved for physical parameters in the wrap tree. +defaults: + - scoring: default + - _self_ + +# hydra hygiene: every output goes to an absolute path from +# adforce.eval.constants; keep the CWD and write no .hydra/ run dirs. +hydra: + output_subdir: null + run: + dir: . + +storms: null # subset of storm display names, e.g. 'storms=["Ida 2021","Katrina 2005"]' (null = all) + +validate: + refresh: false # force-recompute the cached time series (was --refresh-cache) + examples_only: false # only the example-panel figure, from cache (was --examples-only) + failures: false # per-city worst-case panels from the cached sweep (was --failures) + n_failures: 6 # panels per city for validate.failures=true (was --n-failures) + +nulltest: + lag: true # temporal-lag null (netCDF-heavy); --no-lag -> nulltest.lag=false + +sensitivity: + node: true # netCDF node-selection sweep; --no-node -> sensitivity.node=false + +detide_sensitivity: + methods: null # subset of method ids, e.g. '["godin_lowpass","noaa_predictions"]' (null = all) + limit: null # only the first N clean pairs (quick look; suppresses the table) + skew: false # instantaneous-residual vs skew-surge peak metric (was --skew) + table_only: false # regenerate the appendix table from the cached CSV (was --table-only) diff --git a/adforce/eval/config/scoring/default.yaml b/adforce/eval/config/scoring/default.yaml new file mode 100644 index 00000000..ed7a986b --- /dev/null +++ b/adforce/eval/config/scoring/default.yaml @@ -0,0 +1,18 @@ +# Node-selection, clean-pair, and uncertainty gates for the gauge validation. +# +# Mirrors adforce/eval/constants.py EXACTLY: constants.py remains the single +# source of truth consumed by library code (and by the ts-cache tag, whose +# byte-stability keeps the utide-expensive caches valid), while this group +# makes the values visible/overridable from the CLI for future sweep studies. +# tests/test_eval.py::test_scoring_config_matches_constants pins the equality +# so the two cannot drift (the worstsurge audit flagged constants drift as a +# live failure mode). +max_node_deg: 0.12 # = C.MAX_NODE_DEG +wet_min_m: 0.3 # = C.WET_MIN_M +knn: 60 # = C.KNN +max_timing_hr: 6.0 # = C.MAX_TIMING_HR +min_obs_peak_m: 0.4 # = C.MIN_OBS_PEAK_M +ts_min_overlap: 24 # = C.TS_MIN_OVERLAP +n_bootstrap: 2000 # = C.N_BOOTSTRAP +bootstrap_seed: 0 # = C.BOOTSTRAP_SEED +utide_min_samples: 2000 # = C.UTIDE_MIN_SAMPLES From 7051fbf746d773feecdb0018a846dcdc324ef996 Mon Sep 17 00:00:00 2001 From: "Simon D.A. Thomas" Date: Sun, 16 Aug 2026 19:54:31 +0100 Subject: [PATCH 04/17] continued refactor --- adforce/eval/cells.py | 103 +++++++++++++ adforce/eval/config/extract_config.yaml | 15 ++ adforce/eval/extract.py | 181 ++++++++++++++++++++++ adforce/eval/sources.py | 184 +++++++++++++++++++++++ adforce/eval/status.py | 75 +++++++++ adforce/eval/validate.py | 24 +-- adforce/fort61.py | 48 +++++- tests/test_eval.py | 192 ++++++++++++++++++++++++ 8 files changed, 812 insertions(+), 10 deletions(-) create mode 100644 adforce/eval/cells.py create mode 100644 adforce/eval/config/extract_config.yaml create mode 100644 adforce/eval/extract.py create mode 100644 adforce/eval/sources.py create mode 100644 adforce/eval/status.py diff --git a/adforce/eval/cells.py b/adforce/eval/cells.py new file mode 100644 index 00000000..b5e4f464 --- /dev/null +++ b/adforce/eval/cells.py @@ -0,0 +1,103 @@ +"""Configuration-cell identity for adforce run comparison. + +A *cell* is one point in the comparison matrix (resolution x tide x swan x +forcing x physics tag). Its canonical string (:func:`cell_id`) names run +directories, cache files, and model-vs-model join keys; provenance is always +established from the ``/config.yaml`` that ``adforce.config.save_config`` +writes (never from a directory name alone). + +Storm identity is kept separate from the cell: the same cell is run once per +storm, in ``////``. +""" + +from __future__ import annotations + +import re +from dataclasses import dataclass +from typing import Optional + +# Wrap-config paths for the physical axes (the one place the wrap tree's +# {value, ...} nesting appears on the eval side). +AXIS_PATHS = { + "resolution": "adcirc.resolution.value", + "tide": "adcirc.tide.value", + "swan": "adcirc.swan.value", +} + + +@dataclass(frozen=True, order=True) +class ConfigCell: + """One comparison-matrix cell (storm-independent).""" + + resolution: str = "mid" # low | mid | high + tide: bool = False # tidal forcing in fort.15 + swan: bool = False # SWAN wave coupling (padcswan) + forcing: str = "storm" # storm | tide | both (tide = wind-off control run) + physics_tag: str = "default" # opaque label for future physics variants + + +def _onoff(b: bool) -> str: + return "on" if b else "off" + + +def cell_id(c: ConfigCell) -> str: + """Deterministic cell name, e.g. ``res-mid_tide-off_swan-off_f-storm``. + + The ``physics_tag`` is appended only when non-default, so today's standard + cells stay short and stable. + """ + parts = [ + f"res-{c.resolution}", + f"tide-{_onoff(c.tide)}", + f"swan-{_onoff(c.swan)}", + f"f-{c.forcing}", + ] + if c.physics_tag != "default": + parts.append(f"p-{c.physics_tag}") + return "_".join(parts) + + +def dir_to_storm(d: str) -> Optional[str]: + """``'22_MICHAEL_2018'`` -> ``'Michael 2018'`` (val_summary storm naming). + + Ported from ``rerun/adcirc/score_resolution.py`` so legacy GCP/ARCHER2 + sweep directories (``__``) can be ingested read-only. + Returns None when the basename does not match the convention. + """ + m = re.match(r"\d+_(.+)_(\d{4})$", d) + if m is None: + return None + name = m.group(1).replace("_", " ").title() + return f"{name} {m.group(2)}" + + +def storm_to_slug(storm: str, fname: str) -> str: + """``('Katrina 2005', '152_KATRINA_2005.nc')`` -> ``'152_KATRINA_2005'``. + + The slug matches both the HF archive filename stem and the run-dir + basename produced by the training driver, so one name keys both sources. + """ + return fname[:-3] if fname.endswith(".nc") else fname + + +def provenance_match(run_cfg, cell: ConfigCell) -> bool: + """Does a run's dumped ``config.yaml`` match this cell's physical axes? + + Args: + run_cfg: OmegaConf ``DictConfig`` loaded from ``/config.yaml`` + (via :func:`adforce.config.load_config`). + cell (ConfigCell): The cell the directory is expected to hold. + + Returns: + bool: True when resolution/tide/swan all agree. Missing keys count as + a mismatch (FOREIGN) -- a run without provenance is never trusted. + """ + try: + adcirc = run_cfg.adcirc + return ( + str(adcirc.resolution.value) == cell.resolution + and bool(adcirc.tide.value) == cell.tide + and bool(adcirc.swan.value) == cell.swan + ) + except Exception: + return False diff --git a/adforce/eval/config/extract_config.yaml b/adforce/eval/config/extract_config.yaml new file mode 100644 index 00000000..13026ab3 --- /dev/null +++ b/adforce/eval/config/extract_config.yaml @@ -0,0 +1,15 @@ +# Root config for python -m adforce.eval.extract (run dir -> gauge_ts.parquet +# reducer; runs remote where the fort.63.nc files are). +defaults: + - scoring: default + - _self_ + +hydra: + output_subdir: null + run: + dir: . + +runs_root: ??? # parent of the per-storm run directories (required) +gauges_csv: null # sid,name,lat,lon CSV; null -> CO-OPS panel over both boxes (needs cache/net) +out: null # optional combined parquet (legacy data/comp/{lowres,midres} convention) +per_run: true # write /gauge_ts.parquet beside each fort.63.nc (harvest convention) diff --git a/adforce/eval/extract.py b/adforce/eval/extract.py new file mode 100644 index 00000000..7964cce5 --- /dev/null +++ b/adforce/eval/extract.py @@ -0,0 +1,181 @@ +"""Reduce adforce run directories to per-gauge water-level series. + +The remote-side reducer of the eval pipeline: a tidal mid-res ``fort.63.nc`` +is 5-8 GB, so each run is reduced to a small long-format parquet (columns +``storm, sid, gauge, time, zeta`` -- the exact schema of the legacy +``data/comp/{lowres,midres}`` sweep extracts) right where it ran; only the +reduced artifact travels to the laptop. + +Ported from ``rerun/adcirc/extract_gauge_series.py`` (which produced the +published resolution-comparison numbers); the node selection mirrors +``adforce.eval.validate._nearest_wet`` with the fort.63-specific guard that a +node must never dry (NaN zeta) during the run. + +Run (hydra; config root adforce/eval/config/extract_config.yaml):: + + python -m adforce.eval.extract runs_root=/work/lowres_runs + python -m adforce.eval.extract runs_root=... out=combined.parquet per_run=false +""" + +from __future__ import annotations + +import glob +import os +from typing import Optional + +import hydra +import numpy as np +import pandas as pd +from omegaconf import DictConfig + +from . import constants as C + + +def load_fort63(path: str): + """Read ``fort.63.nc`` -> ``(x, y, depth, zeta, t)``; zeta NaN when dry.""" + import netCDF4 as nc + + ds = nc.Dataset(path) + x = np.asarray(ds.variables["x"][:]) + y = np.asarray(ds.variables["y"][:]) + depth = np.asarray(ds.variables["depth"][:]) + zeta = ds.variables["zeta"][:] # (time, node), masked/fill when dry + zeta = np.ma.filled(zeta, np.nan) if np.ma.isMaskedArray(zeta) else np.asarray(zeta) + t = nc.num2date( + ds.variables["time"][:], + ds.variables["time"].units, + only_use_cftime_datetimes=False, + only_use_python_datetimes=True, + ) + ds.close() + return x, y, depth, zeta, pd.DatetimeIndex(t) + + +def gauge_frame(boxes: Optional[list] = None) -> pd.DataFrame: + """CO-OPS gauge panel as a DataFrame (``sid, name, lat, lon``). + + Defaults to the union of the Gulf and Florida boxes -- the panel the + validation sweep scores against (the remote side ships this as a CSV so + extraction needs no CO-OPS access). + """ + from .coops import gulf_gauges + + boxes = boxes if boxes is not None else [C.GAUGE_BOX, C.FLORIDA_BOX] + rows, seen = [], set() + for box in boxes: + for sid, name, lat, lon in gulf_gauges(box): + if sid not in seen: + seen.add(sid) + rows.append(dict(sid=str(sid), name=name, lat=lat, lon=lon)) + return pd.DataFrame(rows) + + +def extract_run( + run_dir: str, + gauges: pd.DataFrame, + max_deg: float = C.MAX_NODE_DEG, + wet_min: float = C.WET_MIN_M, + knn: int = C.KNN, +) -> pd.DataFrame: + """Extract every matchable gauge's zeta series from one run directory. + + The ``storm`` column is the run-dir basename (legacy ``__`` + convention decodes via :func:`adforce.eval.cells.dir_to_storm`). + + A node qualifies for a gauge when it is within ``max_deg`` of it and its + total depth ``zeta + depth`` stays finite (never dries) and above + ``wet_min`` for the whole run -- exactly the published behaviour of + ``rerun/adcirc/extract_gauge_series.py``. + """ + from scipy.spatial import cKDTree + + storm = os.path.basename(os.path.normpath(run_dir)) + x, y, depth, zeta, t = load_fort63(os.path.join(run_dir, "fort.63.nc")) + total = zeta + depth[None, :] + tree = cKDTree(np.c_[x, y]) + frames = [] + for g in gauges.itertuples(): + d, idx = tree.query([g.lon, g.lat], k=min(knn, len(x))) + pick = None + for di, ii in zip(np.atleast_1d(d), np.atleast_1d(idx)): + if di > max_deg: + break + col = total[:, ii] + if np.isfinite(col).all() and np.nanmin(col) > wet_min: + pick = int(ii) + break + if pick is None: + continue + frames.append( + pd.DataFrame( + dict(storm=storm, sid=str(g.sid), gauge=g.name, time=t, zeta=zeta[:, pick]) + ) + ) + if not frames: + return pd.DataFrame(columns=["storm", "sid", "gauge", "time", "zeta"]) + return pd.concat(frames, ignore_index=True) + + +def extract_runs( + runs_root: str, + gauges: pd.DataFrame, + out: Optional[str] = None, + per_run: bool = True, + max_deg: float = C.MAX_NODE_DEG, + wet_min: float = C.WET_MIN_M, + knn: int = C.KNN, +) -> pd.DataFrame: + """Extract every run dir under ``runs_root`` that has a ``fort.63.nc``. + + Args: + runs_root (str): Parent of the per-storm run directories. + gauges (pd.DataFrame): ``sid, name, lat, lon`` rows. + out (Optional[str]): Combined-parquet path (legacy sweep convention). + per_run (bool): Also write ``/gauge_ts.parquet`` beside each + ``fort.63.nc`` (the minimal-artifact convention the laptop + harvests). + + Returns: + pd.DataFrame: The combined long-format table. + """ + frames = [] + for run in sorted(p for p in glob.glob(os.path.join(runs_root, "*")) if os.path.isdir(p)): + if not os.path.exists(os.path.join(run, "fort.63.nc")): + continue + df = extract_run(run, gauges, max_deg=max_deg, wet_min=wet_min, knn=knn) + n_gauges = df.sid.nunique() if not df.empty else 0 + print(f"{os.path.basename(run)}: {n_gauges}/{len(gauges)} gauges matched") + if per_run and not df.empty: + df.to_parquet(os.path.join(run, "gauge_ts.parquet"), index=False) + frames.append(df) + combined = ( + pd.concat(frames, ignore_index=True) + if frames + else pd.DataFrame(columns=["storm", "sid", "gauge", "time", "zeta"]) + ) + if out: + combined.to_parquet(out, index=False) + print(f"wrote {out} ({len(combined)} rows)") + return combined + + +@hydra.main(version_base=None, config_path="config", config_name="extract_config") +def main(cfg: DictConfig) -> None: + gauges = ( + pd.read_csv(cfg.gauges_csv, dtype={"sid": str}) + if cfg.gauges_csv + else gauge_frame() + ) + extract_runs( + cfg.runs_root, + gauges, + out=cfg.out, + per_run=cfg.per_run, + max_deg=cfg.scoring.max_node_deg, + wet_min=cfg.scoring.wet_min_m, + knn=cfg.scoring.knn, + ) + + +if __name__ == "__main__": + main() diff --git a/adforce/eval/sources.py b/adforce/eval/sources.py new file mode 100644 index 00000000..b7d976a8 --- /dev/null +++ b/adforce/eval/sources.py @@ -0,0 +1,184 @@ +"""Simulation sources: where a storm's simulated water-level field comes from. + +Two granularities: + +* :class:`HFArchiveSource` / :class:`Fort63RunSource` provide the full field + ``(x, y, wet, elev, t)`` that :func:`adforce.eval.validate.validate_storm` + samples with its nearest-wet-node lookup. ``wet`` is the array the + wet-criterion runs on (total water depth), ``elev`` the elevation series + that is scored (surge for the tide-off HF archive, raw zeta -- total water + level under tidal forcing -- for run dirs). +* :class:`RunDirSource` provides per-gauge series directly from a completed + adforce run directory, preferring the cheapest artifact present: + ``gauge_ts.parquet`` (remote-extracted) > ``fort.61.nc`` (station output) + > ``fort.63.nc`` (full field, nearest-wet sampling). + +Geometry note: the HF archive stores element-centroid dual-graph fields with +a gap-free ``WD`` (so ``np.nanmin`` is a sufficient wet criterion), while +``fort.63.nc`` is node-based with NaN when a node dries -- Fort63RunSource +maps dry samples to ``-inf`` so the shared nearest-wet criterion rejects any +node that ever dries, reproducing the published behaviour of +``rerun/adcirc/extract_gauge_series.py``. +""" + +from __future__ import annotations + +import os +from typing import Dict, Optional, Protocol, Tuple + +import numpy as np +import pandas as pd + +from . import constants as C + +#: per-gauge result: (gauge display name, series [m], node/station distance +#: in degrees (NaN when unknown), artifact the series came from) +GaugeSeries = Tuple[str, pd.Series, float, str] + + +class FieldSource(Protocol): + """A storm's full simulated field, in validate_storm's expected shape.""" + + def load( + self, storm: str, fname: Optional[str] = None + ) -> Tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray, pd.DatetimeIndex]: + """Return ``(x, y, wet, elev, t)`` for one storm.""" + ... + + +class HFArchiveSource: + """The published SurgeNet HF archive (element centroids, tide-off). + + A thin delegation around the extraction previously inlined in + ``validate_storm``: ``ssh = WD + DEM`` is the surge elevation, ``WD`` the + wet criterion. Behaviour (dtypes, values) is identical by construction -- + the byte-identical ``val_summary.csv`` regression gate pins it. + """ + + def load(self, storm: str, fname: Optional[str] = None): + import xarray as xr + + from .validate import download_storm + + if fname is None: + fname = C.STORMS[storm] + ds = xr.open_dataset(download_storm(fname)) + x, y, DEM, WD = ds.x.values, ds.y.values, ds.DEM.values, ds.WD.values + ssh = WD + DEM[None, :] + t = pd.to_datetime(ds.time.values) + return x, y, WD, ssh, t + + +class Fort63RunSource: + """Full-field source for one adforce run directory's ``fort.63.nc``. + + ``elev`` is raw zeta: total water level when the run had tidal forcing, + surge-only otherwise (read the run's ``config.yaml`` to know which). + """ + + def __init__(self, run_dir: str): + self.run_dir = run_dir + + def load(self, storm: str, fname: Optional[str] = None): + from .extract import load_fort63 + + x, y, depth, zeta, t = load_fort63(os.path.join(self.run_dir, "fort.63.nc")) + total = zeta + depth[None, :] + # NaN (dry) -> -inf so nanmin-based wet criteria reject drying nodes. + wet = np.where(np.isfinite(total), total, -np.inf) + return x, y, wet, zeta, t + + +class RunDirSource: + """Per-gauge series from a completed run dir, cheapest artifact first. + + Args: + run_dir (str): The per-storm run directory. + prefer (tuple): Artifact preference order among + ``("gauge_ts", "fort61", "fort63")``. Gauges the preferred + artifact cannot supply fall through to the next one. + max_station_deg (float): fort.61 stations are matched to gauges by + nearest Euclidean distance in degrees within this radius (station + lists are configured locations, not guaranteed to be the CO-OPS + panel). + """ + + def __init__( + self, + run_dir: str, + prefer: Tuple[str, ...] = ("gauge_ts", "fort61", "fort63"), + max_station_deg: float = 0.02, + ): + self.run_dir = run_dir + self.prefer = prefer + self.max_station_deg = max_station_deg + + # -- one loader per artifact kind; each returns {sid: GaugeSeries} ------ # + def _from_gauge_ts(self, gauges: pd.DataFrame) -> Dict[str, GaugeSeries]: + fp = os.path.join(self.run_dir, "gauge_ts.parquet") + if not os.path.exists(fp): + return {} + df = pd.read_parquet(fp) + out: Dict[str, GaugeSeries] = {} + wanted = set(gauges.sid.astype(str)) + for sid, grp in df.groupby("sid", sort=False): + sid = str(sid) + if sid not in wanted: + continue + grp = grp.sort_values("time") + s = pd.Series(grp.zeta.values, index=pd.DatetimeIndex(grp.time)) + out[sid] = (str(grp.gauge.iloc[0]), s, float("nan"), "gauge_ts") + return out + + def _from_fort61(self, gauges: pd.DataFrame) -> Dict[str, GaugeSeries]: + fp = os.path.join(self.run_dir, "fort.61.nc") + if not os.path.exists(fp): + return {} + from adforce.fort61 import read_fort61 + + df = read_fort61(fp) + stations = df.drop_duplicates("station")[["station", "x", "y"]] + out: Dict[str, GaugeSeries] = {} + for g in gauges.itertuples(): + d = np.hypot(stations.x - g.lon, stations.y - g.lat) + i = int(d.idxmin()) + if float(d.loc[i]) > self.max_station_deg: + continue + grp = df[df.station == stations.station.loc[i]].sort_values("time") + s = pd.Series(grp.zeta.values, index=pd.DatetimeIndex(grp.time)) + out[str(g.sid)] = (str(g.name), s, float(d.loc[i]), "fort61") + return out + + def _from_fort63(self, gauges: pd.DataFrame) -> Dict[str, GaugeSeries]: + fp = os.path.join(self.run_dir, "fort.63.nc") + if not os.path.exists(fp): + return {} + from scipy.spatial import cKDTree + + from .validate import _nearest_wet + + x, y, wet, elev, t = Fort63RunSource(self.run_dir).load("") + tree = cKDTree(np.column_stack([x, y])) + out: Dict[str, GaugeSeries] = {} + for g in gauges.itertuples(): + idx, dist = _nearest_wet(tree, wet, g.lon, g.lat) + if idx is None: + continue + s = pd.Series(elev[:, idx], index=t) + out[str(g.sid)] = (str(g.name), s, dist, "fort63") + return out + + def sim_series(self, gauges: pd.DataFrame) -> Dict[str, GaugeSeries]: + """Series for each gauge (``sid, name, lat, lon`` DataFrame rows).""" + loaders = { + "gauge_ts": self._from_gauge_ts, + "fort61": self._from_fort61, + "fort63": self._from_fort63, + } + out: Dict[str, GaugeSeries] = {} + for kind in self.prefer: + missing = gauges[~gauges.sid.astype(str).isin(out)] + if missing.empty: + break + out.update(loaders[kind](missing)) + return out diff --git a/adforce/eval/status.py b/adforce/eval/status.py new file mode 100644 index 00000000..071ca0cd --- /dev/null +++ b/adforce/eval/status.py @@ -0,0 +1,75 @@ +"""Run-directory status for eval orchestration and discovery. + +Generalizes the receipt logic of ``adforce.check_training_runs.status_list`` +(the ``'Job completed successfully.'`` line that the wrap/slurm path prints +into ``/slurm.out``) with provenance checking against the cell a +directory is supposed to hold. + +Deliberately light on imports: the laptop scoring path must work without the +launch stack, so this module never imports ``adforce.wrap`` or slurmpy. +Slurm-queue polling (RUNNING) and the maxele-based fallback for subprocess +runs without a ``slurm.out`` belong to the launch backend (Commit E), not +here. +""" + +from __future__ import annotations + +import os +from enum import Enum +from typing import Optional + +from .cells import ConfigCell, provenance_match + +#: Written by the wrap slurm/subprocess paths on success (see +#: adforce/check_training_runs.py, which scans for the same line). +SUCCESS_MARKER = "Job completed successfully." + + +class RunStatus(str, Enum): + MISSING = "missing" # no directory or no config.yaml provenance + FOREIGN = "foreign" # config.yaml disagrees with the expected cell + EXTRACTED = "extracted" # success + gauge_ts.parquet already reduced + SUCCESS = "success" # success receipt in slurm.out + FAILED = "failed" # provenance present but no success receipt + + +def _has_success_receipt(run_dir: str) -> bool: + out = os.path.join(run_dir, "slurm.out") + if not os.path.exists(out): + return False + try: + with open(out, errors="ignore") as f: + return any(SUCCESS_MARKER in line for line in f) + except OSError: + return False + + +def run_status(run_dir: str, cell: Optional[ConfigCell] = None) -> RunStatus: + """Classify one run directory. + + Args: + run_dir (str): The per-storm run directory. + cell (Optional[ConfigCell]): Expected cell; when given, the dumped + ``config.yaml`` must match its physical axes or the directory is + FOREIGN (never scored, never overwritten without ``overwrite``). + + Returns: + RunStatus: See the enum docstrings. + """ + cfg_path = os.path.join(run_dir, "config.yaml") + if not os.path.isdir(run_dir) or not os.path.exists(cfg_path): + return RunStatus.MISSING + if cell is not None: + from adforce.config import load_config # omegaconf only, no wrap import + + try: + run_cfg = load_config(cfg_path) + except Exception: + return RunStatus.MISSING + if not provenance_match(run_cfg, cell): + return RunStatus.FOREIGN + if not _has_success_receipt(run_dir): + return RunStatus.FAILED + if os.path.exists(os.path.join(run_dir, "gauge_ts.parquet")): + return RunStatus.EXTRACTED + return RunStatus.SUCCESS diff --git a/adforce/eval/validate.py b/adforce/eval/validate.py index 2efbc5ba..f6d73b24 100644 --- a/adforce/eval/validate.py +++ b/adforce/eval/validate.py @@ -36,7 +36,6 @@ import numpy as np import pandas as pd -import xarray as xr from scipy.spatial import cKDTree from . import constants as C @@ -210,27 +209,34 @@ def load_storm_series( def validate_storm( - storm: str, fname: str, gauges: List[Gauge] + storm: str, fname: str, gauges: List[Gauge], source=None ) -> Tuple[List[dict], Dict[str, tuple]]: - """Return (rows, series) for one storm. ``series[name] = (sim, obs)``.""" + """Return (rows, series) for one storm. ``series[name] = (sim, obs)``. + + ``source`` is a :class:`adforce.eval.sources.FieldSource` supplying + ``(x, y, wet, elev, t)``; the default HF-archive source reproduces the + original inline extraction exactly (``wet = WD``, ``elev = WD + DEM``). + The nearest-wet sampling and scoring below are source-agnostic. + """ year = int(storm.split()[-1]) - ds = xr.open_dataset(download_storm(fname)) - x, y, DEM, WD = ds.x.values, ds.y.values, ds.DEM.values, ds.WD.values - ssh = WD + DEM[None, :] - t = pd.to_datetime(ds.time.values) + if source is None: + from .sources import HFArchiveSource + + source = HFArchiveSource() + x, y, wet, elev, t = source.load(storm, fname) s0, s1 = t[0], t[-1] tree = cKDTree(np.column_stack([x, y])) rows: List[dict] = [] series: Dict[str, tuple] = {} for sid, name, lat, lon in gauges: - idx, dist = _nearest_wet(tree, WD, lon, lat) + idx, dist = _nearest_wet(tree, wet, lon, lat) if idx is None: continue obs, method = observed_residual(sid, lat, year, s0, s1) if obs.empty or obs.size < 12: continue - sim = pd.Series(ssh[:, idx], index=t) + sim = pd.Series(elev[:, idx], index=t) series[name] = (sim, obs) timing = (sim.idxmax() - obs.idxmax()).total_seconds() / 3600.0 ts_r, ts_rmse, ts_n = timeseries_skill(sim, obs) diff --git a/adforce/fort61.py b/adforce/fort61.py index 7480037a..9148a272 100755 --- a/adforce/fort61.py +++ b/adforce/fort61.py @@ -1,10 +1,56 @@ -"""Fort.61""" +"""Fort.61 (ADCIRC station-elevation output) reading and plotting.""" import os + import matplotlib.pyplot as plt +import numpy as np +import pandas as pd + from .mesh import xr_loader +def read_fort61(path: str) -> pd.DataFrame: + """Read an ADCIRC ``fort.61.nc`` into a tidy long-format DataFrame. + + Columns: ``station`` (0-based index), ``name`` (trimmed station name, or + the index as a string when the file has no ``station_name`` variable), + ``x``/``y`` (station lon/lat), ``time``, ``zeta`` (elevation in metres; + NaN where the file uses a fill/dry value). + + Args: + path (str): Path to a ``fort.61.nc`` file (or a run directory + containing one). + + Returns: + pd.DataFrame: One row per (station, time step). + """ + if os.path.isdir(path): + path = os.path.join(path, "fort.61.nc") + ds = xr_loader(path) + zeta = np.asarray(ds.zeta.values, dtype=float) # (time, station) + t = pd.DatetimeIndex(pd.to_datetime(ds.time.values)) + x = np.asarray(ds.x.values, dtype=float) + y = np.asarray(ds.y.values, dtype=float) + n_station = zeta.shape[1] + if "station_name" in ds: + raw = ds.station_name.values + names = [ + (b.tobytes().decode(errors="ignore") if hasattr(b, "tobytes") else str(b)) + .strip("\x00 ") + .strip() + for b in raw + ] + else: + names = [str(i) for i in range(n_station)] + frames = [ + pd.DataFrame( + dict(station=i, name=names[i], x=x[i], y=y[i], time=t, zeta=zeta[:, i]) + ) + for i in range(n_station) + ] + return pd.concat(frames, ignore_index=True) + + def plot_model_tgauges(folder: str = "../../kat.nws13/") -> None: """ Plot the zeta timeseries from the fort.61 file. diff --git a/tests/test_eval.py b/tests/test_eval.py index 03a1409d..5357be2b 100644 --- a/tests/test_eval.py +++ b/tests/test_eval.py @@ -556,3 +556,195 @@ def test_cache_tag_matches_legacy_literal(): "v1|deg0.12|wet0.3|knn60|ut2000" "|box{'lon': (-82.3, -79.7), 'lat': (24.4, 30.8)}" ) + + +# --------------------------------------------------------------------------- # +# Cell identity, run-dir status, sources, fort.61/63 readers (Commit C). +# All synthetic / tmp_path -- no network, no real run dirs. +# --------------------------------------------------------------------------- # +def test_cell_id_and_dir_to_storm(): + from adforce.eval.cells import ConfigCell, cell_id, dir_to_storm, storm_to_slug + + assert cell_id(ConfigCell()) == "res-mid_tide-off_swan-off_f-storm" + assert ( + cell_id(ConfigCell(resolution="low", tide=True, forcing="both")) + == "res-low_tide-on_swan-off_f-both" + ) + assert cell_id(ConfigCell(physics_tag="rh08")).endswith("_p-rh08") + assert dir_to_storm("22_MICHAEL_2018") == "Michael 2018" + assert dir_to_storm("152_KATRINA_2005") == "Katrina 2005" + assert dir_to_storm("not-a-run-dir") is None + assert storm_to_slug("Katrina 2005", "152_KATRINA_2005.nc") == "152_KATRINA_2005" + + +def test_run_status_lifecycle(tmp_path): + from omegaconf import OmegaConf + + from adforce.eval.cells import ConfigCell + from adforce.eval.status import SUCCESS_MARKER, RunStatus, run_status + + cell = ConfigCell(resolution="mid", tide=False, swan=False) + run = tmp_path / "152_KATRINA_2005" + assert run_status(str(run), cell) is RunStatus.MISSING # no dir + run.mkdir() + assert run_status(str(run), cell) is RunStatus.MISSING # no config.yaml + cfg = OmegaConf.create( + {"adcirc": {"resolution": {"value": "mid"}, "tide": {"value": False}, "swan": {"value": False}}} + ) + OmegaConf.save(cfg, str(run / "config.yaml")) + assert run_status(str(run), cell) is RunStatus.FAILED # no receipt yet + (run / "slurm.out").write_text(f"stuff\n{SUCCESS_MARKER}\n") + assert run_status(str(run), cell) is RunStatus.SUCCESS + (run / "gauge_ts.parquet").write_bytes(b"") + assert run_status(str(run), cell) is RunStatus.EXTRACTED + foreign = ConfigCell(resolution="low") + assert run_status(str(run), foreign) is RunStatus.FOREIGN + assert run_status(str(run), None) is RunStatus.EXTRACTED # no provenance check + + +def _write_fort63(path, x, y, depth, zeta, times): + """Synthetic node-based fort.63.nc via xarray.""" + import xarray as xr + + ds = xr.Dataset( + dict( + x=("node", np.asarray(x, dtype=float)), + y=("node", np.asarray(y, dtype=float)), + depth=("node", np.asarray(depth, dtype=float)), + zeta=(("time", "node"), np.asarray(zeta, dtype=float)), + ), + coords=dict(time=times), + ) + ds.to_netcdf(path) + + +def test_extract_run_node_selection(tmp_path): + """Nearest node that never dries wins; drying (NaN) nodes are rejected.""" + from adforce.eval.extract import extract_run + + run = tmp_path / "22_MICHAEL_2018" + run.mkdir() + t = pd.date_range("2018-10-09", periods=6, freq="h") + # node 0: nearest to the gauge but dries (NaN); node 1: wet throughout; + # node 2: far away (> max_deg). + zeta = np.array( + [ + [np.nan, 0.5, 0.1], + [0.2, 0.6, 0.1], + [0.3, 0.9, 0.1], + [0.2, 0.7, 0.1], + [np.nan, 0.5, 0.1], + [0.1, 0.4, 0.1], + ] + ) + _write_fort63( + run / "fort.63.nc", + x=[-90.00, -90.02, -91.5], + y=[29.00, 29.02, 29.5], + depth=[5.0, 5.0, 5.0], + zeta=zeta, + times=t, + ) + gauges = pd.DataFrame( + [dict(sid="8761724", name="Grand Isle", lat=29.0, lon=-90.0)] + ) + df = extract_run(str(run), gauges, max_deg=0.12, wet_min=0.3, knn=3) + assert list(df.columns) == ["storm", "sid", "gauge", "time", "zeta"] + assert df.storm.unique().tolist() == ["22_MICHAEL_2018"] + assert df.sid.unique().tolist() == ["8761724"] + np.testing.assert_allclose(df.zeta.values, zeta[:, 1]) # picked the wet node + + # An all-drying mesh yields no rows. + run2 = tmp_path / "all_dry" + run2.mkdir() + _write_fort63( + run2 / "fort.63.nc", + x=[-90.0], + y=[29.0], + depth=[5.0], + zeta=np.full((6, 1), np.nan), + times=t, + ) + assert extract_run(str(run2), gauges).empty + + +def test_read_fort61_roundtrip(tmp_path): + import xarray as xr + + from adforce.fort61 import read_fort61 + + t = pd.date_range("2018-10-09", periods=4, freq="h") + zeta = np.array([[0.1, 1.0], [0.2, 1.1], [0.3, 1.2], [0.2, 1.3]]) + ds = xr.Dataset( + dict( + x=("station", [-90.0, -89.5]), + y=("station", [29.0, 29.5]), + zeta=(("time", "station"), zeta), + ), + coords=dict(time=t), + ) + ds.to_netcdf(tmp_path / "fort.61.nc") + df = read_fort61(str(tmp_path / "fort.61.nc")) + assert set(df.columns) == {"station", "name", "x", "y", "time", "zeta"} + assert df.station.nunique() == 2 + np.testing.assert_allclose( + df[df.station == 1].sort_values("time").zeta.values, zeta[:, 1] + ) + # directory form resolves fort.61.nc inside + assert len(read_fort61(str(tmp_path))) == len(df) + + +def test_rundir_source_prefers_cheapest_artifact(tmp_path): + from adforce.eval.sources import RunDirSource + + run = tmp_path / "run" + run.mkdir() + t = pd.date_range("2018-10-09", periods=4, freq="h") + gauges = pd.DataFrame( + [dict(sid="1", name="A", lat=29.0, lon=-90.0)] + ) + # only fort.63 present -> fort63 sampling + _write_fort63( + run / "fort.63.nc", + x=[-90.0], + y=[29.0], + depth=[5.0], + zeta=np.array([[0.5], [0.6], [0.7], [0.6]]), + times=t, + ) + out = RunDirSource(str(run)).sim_series(gauges) + assert out["1"][3] == "fort63" + # gauge_ts.parquet appears -> preferred over fort63 + pd.DataFrame( + dict(storm="run", sid="1", gauge="A", time=t, zeta=[1.0, 2.0, 3.0, 2.0]) + ).to_parquet(run / "gauge_ts.parquet", index=False) + out = RunDirSource(str(run)).sim_series(gauges) + name, series, deg, src = out["1"] + assert src == "gauge_ts" and series.max() == 3.0 + + +def test_validate_storm_source_seam(tmp_path, monkeypatch): + """validate_storm scores an injected FieldSource without touching HF/CO-OPS.""" + import adforce.eval.validate as cv + + monkeypatch.setattr(C, "TS_CACHE", str(tmp_path)) + t = pd.date_range("2021-08-28", periods=48, freq="h") + surge = np.concatenate([np.linspace(0, 2.0, 24), np.linspace(2.0, 0, 24)]) + + class FakeSource: + def load(self, storm, fname=None): + x = np.array([-90.0]) + y = np.array([29.0]) + wet = np.full((48, 1), 5.0) # always wet + elev = surge[:, None] + return x, y, wet, elev, pd.DatetimeIndex(t) + + obs = pd.Series(surge * 0.9, index=t) + monkeypatch.setattr(cv, "observed_residual", lambda *a, **k: (obs, "utide")) + rows, series = cv.validate_storm( + "Fake 2021", "fake.nc", [("42", "Fake Gauge", 29.0, -90.0)], source=FakeSource() + ) + assert len(rows) == 1 and rows[0]["sid"] == "42" + assert rows[0]["sim_peak"] == 2.0 + assert abs(rows[0]["peak_dt_hr"]) < 1e-9 # aligned peaks + assert "Fake Gauge" in series From ff5cccf693d17ea78539551563493e02441d2de0 Mon Sep 17 00:00:00 2001 From: "Simon D.A. Thomas" Date: Sun, 16 Aug 2026 20:23:04 +0100 Subject: [PATCH 05/17] eval pairs + matrix grammar: model-vs-model, published-artifact pins - cells.expand_matrix(): axes cartesian product - exclude + include/cells, name_keys aliases (res-low_tide-off style), duplicate/baseline validation - config/matrix/{res_x_tide,archive_default}.yaml, config/source/ {hf_archive,run_dirs}.yaml, config/pairs_config.yaml - pairs.py: resolution_bias_table (port of rerun/adcirc/score_resolution.py), tide_surge_interaction (canonical home for the forcing-triple analysis whose producer script was never committed; series + HF-archive-val_summary fallback with storm_src provenance), compare_cells (general pairwise), archive_peak_table; hydra main (action=resolution_bias|interaction) - SUPERSEDED headers on rerun/adcirc/{extract_gauge_series,score_resolution}.py Acceptance (pinned as skipif-gated tests): resolution_bias.csv reproduced exactly (shape+columns+allclose); tide_surge_interaction.csv low rows 683/683 exact incl. provenance labels; mid rows verified on the 182 locally-supportable rows (all 78 committed-archive rows allclose) -- the remaining mid rows need remote extracts whose local copies are truncated (mid_storm_runs_gs_*.parquet fail parquet magic-byte check). Co-Authored-By: Claude Fable 5 --- adforce/eval/cells.py | 46 ++++ .../eval/config/matrix/archive_default.yaml | 14 + adforce/eval/config/matrix/res_x_tide.yaml | 13 + adforce/eval/config/pairs_config.yaml | 20 ++ adforce/eval/config/source/hf_archive.yaml | 5 + adforce/eval/config/source/run_dirs.yaml | 6 + adforce/eval/pairs.py | 250 ++++++++++++++++++ rerun/adcirc/extract_gauge_series.py | 4 + rerun/adcirc/score_resolution.py | 4 + tests/test_eval.py | 127 +++++++++ 10 files changed, 489 insertions(+) create mode 100644 adforce/eval/config/matrix/archive_default.yaml create mode 100644 adforce/eval/config/matrix/res_x_tide.yaml create mode 100644 adforce/eval/config/pairs_config.yaml create mode 100644 adforce/eval/config/source/hf_archive.yaml create mode 100644 adforce/eval/config/source/run_dirs.yaml create mode 100644 adforce/eval/pairs.py diff --git a/adforce/eval/cells.py b/adforce/eval/cells.py index b5e4f464..88165092 100644 --- a/adforce/eval/cells.py +++ b/adforce/eval/cells.py @@ -80,6 +80,52 @@ def storm_to_slug(storm: str, fname: str) -> str: return fname[:-3] if fname.endswith(".nc") else fname +def expand_matrix(matrix_cfg) -> list: + """Expand the matrix-YAML grammar into ``[(cell_name, overrides), ...]``. + + Grammar (see ``adforce/eval/config/matrix/``): ``axes`` maps WRAP config + paths (``adcirc.resolution.value`` etc.) to value lists whose cartesian + product, minus partial-match ``exclude`` dicts, forms the grid; + ``cells``/``include`` append fully explicit ``{name, overrides}`` entries. + Cell names come from ``name_keys`` aliases in axes order (bools render + ``on``/``off``), e.g. ``res-low_tide-off``. + + Raises: + ValueError: On duplicate cell names or a ``baseline`` naming no cell. + """ + from itertools import product + + axes = dict(matrix_cfg.get("axes") or {}) + name_keys = dict(matrix_cfg.get("name_keys") or {}) + excludes = [dict(e) for e in (matrix_cfg.get("exclude") or [])] + + def render(path, val) -> str: + alias = name_keys.get(path, path.split(".")[-1]) + v = ("on" if val else "off") if isinstance(val, bool) else str(val) + return f"{alias}-{v}" + + cells = [] + if axes: + paths = list(axes) + for combo in product(*(list(axes[p]) for p in paths)): + ov = dict(zip(paths, combo)) + if any(all(ov.get(k) == v for k, v in ex.items()) for ex in excludes): + continue + cells.append(("_".join(render(p, ov[p]) for p in paths), ov)) + for group in ("cells", "include"): + for item in matrix_cfg.get(group) or []: + cells.append((str(item["name"]), dict(item["overrides"]))) + + names = [n for n, _ in cells] + dups = sorted({n for n in names if names.count(n) > 1}) + if dups: + raise ValueError(f"duplicate cell names in matrix: {dups}") + baseline = matrix_cfg.get("baseline") + if baseline and baseline not in names: + raise ValueError(f"baseline {baseline!r} is not one of the cells {names}") + return cells + + def provenance_match(run_cfg, cell: ConfigCell) -> bool: """Does a run's dumped ``config.yaml`` match this cell's physical axes? diff --git a/adforce/eval/config/matrix/archive_default.yaml b/adforce/eval/config/matrix/archive_default.yaml new file mode 100644 index 00000000..8c7e0900 --- /dev/null +++ b/adforce/eval/config/matrix/archive_default.yaml @@ -0,0 +1,14 @@ +# The published HF EC95d no-tide archive as a degenerate one-cell matrix: +# scoring with this matrix reproduces the paper's validation numbers inside +# the general framework. +axes: {} +name_keys: {} +exclude: [] +include: [] +cells: + - name: ec95d-archive + overrides: + adcirc.resolution.value: mid + adcirc.tide.value: false + adcirc.swan.value: false +baseline: ec95d-archive diff --git a/adforce/eval/config/matrix/res_x_tide.yaml b/adforce/eval/config/matrix/res_x_tide.yaml new file mode 100644 index 00000000..3fc5bcd2 --- /dev/null +++ b/adforce/eval/config/matrix/res_x_tide.yaml @@ -0,0 +1,13 @@ +# low/mid resolution x tide on/off comparison matrix. +# Axis keys are WRAP config paths (hence the .value nesting -- the wrap +# tree's convention; eval's own keys stay flat). +axes: + adcirc.resolution.value: [low, mid] + adcirc.tide.value: [false, true] +name_keys: + adcirc.resolution.value: res + adcirc.tide.value: tide +exclude: [] # partial-match dicts over axis paths, e.g. {adcirc.resolution.value: low} +include: [] # extra fully-explicit {name, overrides} cells appended to the grid +cells: [] # alternative: fully explicit named cells instead of axes +baseline: res-mid_tide-off # model-vs-model reference cell diff --git a/adforce/eval/config/pairs_config.yaml b/adforce/eval/config/pairs_config.yaml new file mode 100644 index 00000000..e8014c1f --- /dev/null +++ b/adforce/eval/config/pairs_config.yaml @@ -0,0 +1,20 @@ +# Root config for python -m adforce.eval.pairs (model-vs-model comparison). +hydra: + output_subdir: null + run: + dir: . + +action: resolution_bias # resolution_bias | interaction + +# resolution_bias inputs +series: ??? # candidate-cell gauge-series parquet (long format) +val_summary: null # baseline val_summary.csv (null -> data/comp/out/val_summary.csv) +out: null # CSV to write +min_obs_peak: 1.0 # pooled-ratio report threshold on observed peak [m] + +# interaction inputs +res: low # resolution label for the output rows +storm_series: null # storm-only gauge-series parquet (null -> archive fallback only) +tide_series: null # tide-only gauge-series parquet +both_series: null # storm+tide gauge-series parquet +use_archive: true # fall back to HF-archive val_summary peaks for storm-only diff --git a/adforce/eval/config/source/hf_archive.yaml b/adforce/eval/config/source/hf_archive.yaml new file mode 100644 index 00000000..2b9c50da --- /dev/null +++ b/adforce/eval/config/source/hf_archive.yaml @@ -0,0 +1,5 @@ +# The published SurgeNet HF archive (sdat2/surgenet-train): element-centroid +# dual-graph fields, ssh = WD + DEM, tides excluded (pure surge). +# Consumed by adforce.eval.sources.HFArchiveSource. +kind: hf_archive +repo: sdat2/surgenet-train # = constants.HF_REPO (drift would break download_storm) diff --git a/adforce/eval/config/source/run_dirs.yaml b/adforce/eval/config/source/run_dirs.yaml new file mode 100644 index 00000000..e6f85746 --- /dev/null +++ b/adforce/eval/config/source/run_dirs.yaml @@ -0,0 +1,6 @@ +# Completed adforce run directories as the simulation source. +# Consumed by adforce.eval.sources.RunDirSource. +kind: run_dirs +runs_root: ??? # parent of the per-storm run directories (required) +prefer: [gauge_ts, fort61, fort63] # artifact preference; gauges fall through to the next +verify_config: true # provenance check vs /config.yaml (see eval.status) diff --git a/adforce/eval/pairs.py b/adforce/eval/pairs.py new file mode 100644 index 00000000..ba01c319 --- /dev/null +++ b/adforce/eval/pairs.py @@ -0,0 +1,250 @@ +"""Model-vs-model comparison: resolution bias and tide-surge interaction. + +Consolidates and generalizes ``rerun/adcirc/score_resolution.py`` and the +(uncommitted) script behind ``rerun/results/tide_surge_interaction.csv``: + +* :func:`resolution_bias_table` -- join a candidate cell's gauge-series peaks + onto ``val_summary.csv`` (baseline peaks + de-tided observed peaks + clean + flags), keyed by (storm, sid). +* :func:`tide_surge_interaction` -- the forcing triple: per (storm, gauge) + ``interaction = peak(storm+tide) - peak(storm) - peak(tide)``. Storm-only + peaks fall back to the HF-archive ``val_summary`` peaks where no storm-only + run exists (``storm_src`` records ``series`` vs ``archive``). +* :func:`compare_cells` -- general pairwise peaks of any two cells' series. + +Inputs are the long-format gauge-series tables (``storm, sid, gauge, time, +zeta``) that :mod:`adforce.eval.extract` writes (and the legacy +``data/comp/{lowres,midres}`` sweep parquets already hold). + +Run (hydra; config root adforce/eval/config/pairs_config.yaml):: + + python -m adforce.eval.pairs action=resolution_bias series=data/comp/lowres/low_storm_gauge_series.parquet +""" + +from __future__ import annotations + +import os +import re +from typing import Dict, Optional + +import hydra +import numpy as np +import pandas as pd +from omegaconf import DictConfig + +from . import constants as C + + +def storm_key(storm: str) -> str: + """``'140_LAURA_2020'`` or ``'Laura 2020'`` -> ``'LAURA_2020'``.""" + s = re.sub(r"^\d+_", "", str(storm)) + return s.replace(" ", "_").upper() + + +def peak_table(series: pd.DataFrame) -> pd.DataFrame: + """Per-(storm, gauge) peak of a long-format series table. + + Returns columns ``key, sid, gauge, peak`` (key = ``NAME_YEAR``). + """ + df = series.copy() + df["key"] = df.storm.map(storm_key) + df["sid"] = df.sid.astype(str) + return df.groupby(["key", "sid"], as_index=False).agg( + gauge=("gauge", "first"), peak=("zeta", "max") + ) + + +def archive_peak_table(val_summary: Optional[str] = None) -> pd.DataFrame: + """HF-archive storm-only peaks from ``val_summary.csv`` (columns as + :func:`peak_table`); the fallback for storms without a storm-only run.""" + path = val_summary or os.path.join(C.OUT_PATH, "val_summary.csv") + vs = pd.read_csv(path) + return pd.DataFrame( + dict( + key=vs.storm.map(storm_key), + sid=vs.sid.astype(str), + gauge=vs.name, + peak=vs.sim_peak, + ) + ) + + +def resolution_bias_table( + series_parquet: str, + val_summary: Optional[str] = None, + out: Optional[str] = None, + min_obs_peak: float = 1.0, +) -> pd.DataFrame: + """Join candidate-cell peaks onto the validation summary (exact port of + ``rerun/adcirc/score_resolution.py``; writes the same CSV shape). + + Args: + series_parquet (str): Long-format gauge series of the candidate cell + (e.g. the low-res sweep extract). + val_summary (Optional[str]): Baseline ``val_summary.csv`` (default: + the standard output path). + out (Optional[str]): CSV to write (default: none). + min_obs_peak (float): Pooled-ratio report threshold on observed peak. + + Returns: + pd.DataFrame: ``val_summary`` columns + ``low_peak`` per (storm, sid). + """ + from .cells import dir_to_storm + + ser = pd.read_parquet(series_parquet) + ser["storm_name"] = ser.storm.map(dir_to_storm) + low = ( + ser.groupby(["storm_name", "sid"], as_index=False) + .zeta.max() + .rename(columns={"zeta": "low_peak", "storm_name": "storm"}) + ) + low["sid"] = low.sid.astype(str) + + val = pd.read_csv(val_summary or os.path.join(C.OUT_PATH, "val_summary.csv")) + val["sid"] = val.sid.astype(str) + df = val.merge(low, on=["storm", "sid"], how="left") + clean = df[ + df.clean.astype(bool) & df.low_peak.notna() & (df.obs_peak >= min_obs_peak) + ].copy() + clean["mid_ratio"] = clean.sim_peak / clean.obs_peak + clean["low_ratio"] = clean.low_peak / clean.obs_peak + + print(f"clean pairs with obs>={min_obs_peak} m and low coverage: {len(clean)}") + for name, r in (("mid/obs", clean.mid_ratio), ("low/obs", clean.low_ratio)): + print(f" {name}: mean {r.mean():.2f} median {np.median(r):.2f}") + print(f" low/mid peak ratio: mean {(clean.low_peak / clean.sim_peak).mean():.2f}") + + if out: + df.to_csv(out, index=False) + print(f"wrote {out}") + return df + + +def tide_surge_interaction( + storm_series: Optional[pd.DataFrame], + tide_series: pd.DataFrame, + both_series: pd.DataFrame, + res: str, + archive_peaks: Optional[pd.DataFrame] = None, +) -> pd.DataFrame: + """Peak-based nonlinear tide-surge interaction for one resolution. + + ``interaction = peak(storm+tide) - peak(storm) - peak(tide)`` per + (storm, gauge); a (tide, both) pair is required, the storm-only peak + comes from ``storm_series`` where available (``storm_src='series'``) and + otherwise from ``archive_peaks`` (``storm_src='archive'``, the HF-archive + ``val_summary`` convention of the published mid-res rows). + + Returns: + pd.DataFrame: ``res, key, sid, gauge, zeta_storm, zeta_tide, + zeta_both, interaction, storm_src`` (the exact schema of + ``rerun/results/tide_surge_interaction.csv``). + """ + tide = peak_table(tide_series).rename(columns={"peak": "zeta_tide"}) + both = peak_table(both_series).rename(columns={"peak": "zeta_both"}) + df = tide.merge(both.drop(columns="gauge"), on=["key", "sid"], how="inner") + + parts = [] + if storm_series is not None and len(storm_series): + s = peak_table(storm_series).drop(columns="gauge") + s = s.rename(columns={"peak": "zeta_storm"}) + s["storm_src"] = "series" + parts.append(s) + if archive_peaks is not None and len(archive_peaks): + a = archive_peaks.drop(columns="gauge").rename(columns={"peak": "zeta_storm"}) + a = a.copy() + a["storm_src"] = "archive" + parts.append(a) + if not parts: + raise ValueError("need storm_series and/or archive_peaks for storm-only peaks") + storm = pd.concat(parts, ignore_index=True).drop_duplicates( + ["key", "sid"], keep="first" # series wins over archive + ) + + df = df.merge(storm, on=["key", "sid"], how="inner") + df["interaction"] = df.zeta_both - df.zeta_storm - df.zeta_tide + df["res"] = res + return df[ + [ + "res", + "key", + "sid", + "gauge", + "zeta_storm", + "zeta_tide", + "zeta_both", + "interaction", + "storm_src", + ] + ].sort_values(["key", "sid"], ignore_index=True) + + +def compare_cells( + frames: Dict[str, pd.DataFrame], baseline: str, out: Optional[str] = None +) -> pd.DataFrame: + """Pairwise peak comparison of every cell against a baseline cell. + + Args: + frames: ``{cell_id: long-format gauge-series table}``. + baseline: Key in ``frames`` to compare everything against. + out: Optional CSV path. + + Returns: + pd.DataFrame: one row per (cell, storm, gauge) with + ``cell, key, sid, gauge, peak, base_peak, peak_diff``. + """ + if baseline not in frames: + raise ValueError(f"baseline {baseline!r} not in frames {sorted(frames)}") + base = peak_table(frames[baseline]).rename(columns={"peak": "base_peak"}) + rows = [] + for cell, series in frames.items(): + if cell == baseline: + continue + p = peak_table(series).drop(columns="gauge") + m = p.merge(base, on=["key", "sid"], how="inner") + m["cell"] = cell + m["peak_diff"] = m.peak - m.base_peak + rows.append(m) + df = ( + pd.concat(rows, ignore_index=True)[ + ["cell", "key", "sid", "gauge", "peak", "base_peak", "peak_diff"] + ] + if rows + else pd.DataFrame( + columns=["cell", "key", "sid", "gauge", "peak", "base_peak", "peak_diff"] + ) + ) + if out: + df.to_csv(out, index=False) + return df + + +@hydra.main(version_base=None, config_path="config", config_name="pairs_config") +def main(cfg: DictConfig) -> None: + C.ensure_dirs() + if cfg.action == "resolution_bias": + resolution_bias_table( + cfg.series, + val_summary=cfg.val_summary, + out=cfg.out, + min_obs_peak=cfg.min_obs_peak, + ) + elif cfg.action == "interaction": + storm = pd.read_parquet(cfg.storm_series) if cfg.storm_series else None + archive = archive_peak_table(cfg.val_summary) if cfg.use_archive else None + df = tide_surge_interaction( + storm, + pd.read_parquet(cfg.tide_series), + pd.read_parquet(cfg.both_series), + res=cfg.res, + archive_peaks=archive, + ) + if cfg.out: + df.to_csv(cfg.out, index=False) + print(f"wrote {cfg.out} ({len(df)} rows)") + else: + raise SystemExit(f"unknown action {cfg.action!r} (resolution_bias|interaction)") + + +if __name__ == "__main__": + main() diff --git a/rerun/adcirc/extract_gauge_series.py b/rerun/adcirc/extract_gauge_series.py index d1dd6d23..6e0a1ded 100644 --- a/rerun/adcirc/extract_gauge_series.py +++ b/rerun/adcirc/extract_gauge_series.py @@ -1,5 +1,9 @@ """Extract simulated water-level series at NOAA gauges from ADCIRC runs. +SUPERSEDED by ``adforce.eval.extract`` (hydra-driven port, same node +selection and output schema); kept as the frozen record of the GCP sweeps +that produced ``data/comp/{lowres,midres}`` and ``rerun/results/*``. + For every run directory ``/__/`` containing a ``fort.63.nc``, find the nearest mesh node to each gauge in ``--gauges`` (CSV: sid,name,lat,lon) that stays wet (min total depth zeta+depth > diff --git a/rerun/adcirc/score_resolution.py b/rerun/adcirc/score_resolution.py index 8fcdcb4f..66d79e86 100644 --- a/rerun/adcirc/score_resolution.py +++ b/rerun/adcirc/score_resolution.py @@ -1,5 +1,9 @@ """Score the low-resolution historical sweep against the mid (EC95d) archive. +SUPERSEDED by ``adforce.eval.pairs.resolution_bias_table`` (which reproduces +``rerun/results/resolution_bias.csv`` exactly -- pinned by +``tests/test_eval.py``); kept as the frozen record of the published run. + Joins the low-res gauge series (rerun/adcirc/extract_gauge_series.py output) with comp.validate's val_summary.csv (mid-resolution sim peaks + de-tided observed peaks + clean flags), keyed by (storm "Name YYYY", gauge sid). The diff --git a/tests/test_eval.py b/tests/test_eval.py index 5357be2b..88cfe7f5 100644 --- a/tests/test_eval.py +++ b/tests/test_eval.py @@ -723,6 +723,133 @@ def test_rundir_source_prefers_cheapest_artifact(tmp_path): assert src == "gauge_ts" and series.max() == 3.0 +# --------------------------------------------------------------------------- # +# Matrix grammar + model-vs-model (Commit D). The reproduction tests pin the +# ported pairs.py against the published rerun/results artifacts and are +# skipped when the cached sweep data is absent. +# --------------------------------------------------------------------------- # +def test_expand_matrix_grammar(): + from adforce.eval.cells import expand_matrix + + cells = expand_matrix( + dict( + axes={ + "adcirc.resolution.value": ["low", "mid"], + "adcirc.tide.value": [False, True], + }, + name_keys={ + "adcirc.resolution.value": "res", + "adcirc.tide.value": "tide", + }, + exclude=[{"adcirc.resolution.value": "low", "adcirc.tide.value": True}], + include=[ + {"name": "swan", "overrides": {"adcirc.swan.value": True}} + ], + baseline="res-mid_tide-off", + ) + ) + names = [n for n, _ in cells] + assert names == ["res-low_tide-off", "res-mid_tide-off", "res-mid_tide-on", "swan"] + assert dict(cells)["res-mid_tide-on"] == { + "adcirc.resolution.value": "mid", + "adcirc.tide.value": True, + } + with pytest.raises(ValueError, match="baseline"): + expand_matrix(dict(axes={}, cells=[], baseline="nope")) + with pytest.raises(ValueError, match="duplicate"): + expand_matrix( + dict( + cells=[ + {"name": "a", "overrides": {}}, + {"name": "a", "overrides": {}}, + ] + ) + ) + + +def test_compare_cells_synthetic(): + from adforce.eval.pairs import compare_cells + + t = pd.date_range("2020-01-01", periods=3, freq="h") + + def frame(peak): + return pd.DataFrame( + dict(storm="1_TEST_2020", sid="42", gauge="G", time=t, zeta=[0.0, peak, 0.1]) + ) + + df = compare_cells({"base": frame(1.0), "cand": frame(1.5)}, baseline="base") + assert len(df) == 1 + r = df.iloc[0] + assert (r.cell, r.key, r.peak, r.base_peak, r.peak_diff) == ( + "cand", + "TEST_2020", + 1.5, + 1.0, + 0.5, + ) + with pytest.raises(ValueError, match="baseline"): + compare_cells({"a": frame(1.0)}, baseline="missing") + + +_LOWRES = os.path.join( + str(REPO_ROOT), "data", "comp", "lowres", "low_storm_gauge_series.parquet" +) +_RERUN_RB = os.path.join(str(REPO_ROOT), "rerun", "results", "resolution_bias.csv") +_RERUN_TSI = os.path.join( + str(REPO_ROOT), "rerun", "results", "tide_surge_interaction.csv" +) + + +@pytest.mark.skipif( + not (os.path.exists(_LOWRES) and os.path.exists(_RERUN_RB) and os.path.exists(SUMMARY_CSV)), + reason="cached low-res sweep / published artifact not present", +) +def test_resolution_bias_reproduces_published(): + """pairs.resolution_bias_table must reproduce rerun/results/resolution_bias.csv + (the published low-vs-mid resolution comparison) from the cached extract.""" + import io + + from adforce.eval.pairs import resolution_bias_table + + new = resolution_bias_table(_LOWRES) + buf = io.StringIO() + new.to_csv(buf, index=False) + buf.seek(0) + new_rt = pd.read_csv(buf) + old = pd.read_csv(_RERUN_RB) + assert list(new_rt.columns) == list(old.columns) + assert new_rt.shape == old.shape + np.testing.assert_allclose( + new_rt.select_dtypes("number").fillna(-999).values, + old.select_dtypes("number").fillna(-999).values, + ) + + +@pytest.mark.skipif( + not (os.path.exists(_LOWRES) and os.path.exists(_RERUN_TSI)), + reason="cached low-res sweep / published artifact not present", +) +def test_tide_surge_interaction_low_reproduces_published(): + """The low-res forcing triple must reproduce every published row exactly + (storm/tide/both peaks, interaction, and series-vs-archive provenance).""" + from adforce.eval.pairs import tide_surge_interaction + + lowres = os.path.dirname(_LOWRES) + new = tide_surge_interaction( + pd.read_parquet(_LOWRES), + pd.read_parquet(os.path.join(lowres, "low_tide_runs_gauge_series.parquet")), + pd.read_parquet(os.path.join(lowres, "low_both_full.parquet")), + "low", + ) + old = pd.read_csv(_RERUN_TSI, dtype={"sid": str}) + old = old[old.res == "low"] + m = new.merge(old, on=["res", "key", "sid"], suffixes=("_n", "_o")) + assert len(m) == len(old) == len(new) == 683 + for col in ("zeta_storm", "zeta_tide", "zeta_both", "interaction"): + np.testing.assert_allclose(m[f"{col}_n"], m[f"{col}_o"]) + assert (m.storm_src_n == m.storm_src_o).all() + + def test_validate_storm_source_seam(tmp_path, monkeypatch): """validate_storm scores an injected FieldSource without touching HF/CO-OPS.""" import adforce.eval.validate as cv From 88d2d152f878c06a3ea327328525badf8277fc62 Mon Sep 17 00:00:00 2001 From: "Simon D.A. Thomas" Date: Sun, 16 Aug 2026 20:33:59 +0100 Subject: [PATCH 06/17] eval launch orchestration: drive_storm seam, plan/dry-run, harvest - training/driver.py: extract drive_storm(storm, storm_ds, run_directory, cfg, resolution, mode, spinup_days, ...) + is_run_successful() from drive_all_adcirc (whose CLI behaviour is unchanged); the seam is the ONLY sanctioned path for tide-on historical runs (per-storm fort.15 windows vs wrap.stage_input_files' Katrina-pinned static decks) - eval/launch.py + config/launch_config.yaml: matrix -> per-(cell,storm) plan table (run/skip/blocked + reason); dry_run=true default; guards: swan cells blocked (driver has no SWAN staging), missing resolution decks blocked (fort.14.high is ARCHER2-only), FOREIGN dirs blocked unless overwrite; controls=true auto-adds tide-only cells for interaction triples; sequential subprocess execution with extract-after-run + strip_after_extract (hist_sweep disk lesson); eval_manifest.json with full-resolved-config sha256 (provenance-drift detection); slurm backend deliberately deferred with loud documentation - eval/harvest.py + config/harvest_config.yaml: rsync of the minimal artifact set (config.yaml, gauge_ts.parquet, fort.61.nc, maxele.63.nc, slurm.out, manifest) remote -> laptop - tcpips/constants.py: import-time makedirs -> _ensure_dir tolerant of the dangling data/era5 symlink (unmounted volume broke `import tcpips`, and transitively the training driver, on any machine without the drive) Gate: 52 passed 1 skipped (driver-seam test skips where adcircpy/stormevents absent -- same reason pytest.ini ignores training modules); dry-run smoke: 2x2 matrix x 2 storms -> 8-row plan, correct run dirs, nothing created. Co-Authored-By: Claude Fable 5 --- adforce/eval/config/harvest_config.yaml | 10 + adforce/eval/config/launch_config.yaml | 26 +++ adforce/eval/harvest.py | 72 +++++++ adforce/eval/launch.py | 262 ++++++++++++++++++++++++ adforce/training/driver.py | 128 ++++++++---- tcpips/constants.py | 58 ++++-- tests/test_eval.py | 112 ++++++++++ 7 files changed, 603 insertions(+), 65 deletions(-) create mode 100644 adforce/eval/config/harvest_config.yaml create mode 100644 adforce/eval/config/launch_config.yaml create mode 100644 adforce/eval/harvest.py create mode 100644 adforce/eval/launch.py diff --git a/adforce/eval/config/harvest_config.yaml b/adforce/eval/config/harvest_config.yaml new file mode 100644 index 00000000..46b1ffc5 --- /dev/null +++ b/adforce/eval/config/harvest_config.yaml @@ -0,0 +1,10 @@ +# Root config for python -m adforce.eval.harvest (remote -> laptop rsync). +hydra: + output_subdir: null + run: + dir: . + +remote: ??? # host:/path/to/exp/eval (parent of study dirs), e.g. archer2:/work/n02/n02/sithom/exp/eval +study: ??? # study subtree to pull +dest: null # local mirror root (null -> data/comp/runs/) +dry_run: true # print the rsync command without running it diff --git a/adforce/eval/config/launch_config.yaml b/adforce/eval/config/launch_config.yaml new file mode 100644 index 00000000..cd34087d --- /dev/null +++ b/adforce/eval/config/launch_config.yaml @@ -0,0 +1,26 @@ +# Root config for python -m adforce.eval.launch (sweep orchestrator). +defaults: + - matrix: res_x_tide + - _self_ + +hydra: + output_subdir: null + run: + dir: . + +study: default # namespaces /// and the manifest +runs_root: ${oc.env:WORSTSURGE_RUNS_ROOT,runs/eval} # relative -> under the repo root +storms: null # subset of storm display names (null = all 19) +controls: false # auto-add tide-only control cells (interaction triples) +spinup_days: auto # auto: 6.0 for tidal cells, else 0.0 (hist_sweep SPINUP) +recommended_dt: 5.0 # ADCIRC timestep [s] passed to input generation +dry_run: true # DEFAULT: print the plan table, launch nothing +skip_completed: true # resume semantics: skip runs with a success receipt +overwrite: false # refuse FOREIGN run dirs (config.yaml mismatch) unless true +extract_after_run: true # reduce each run to gauge_ts.parquet right after it finishes +strip_after_extract: # then delete the bulky outputs (hist_sweep disk lesson) + - fort.63.nc + - fort.64.nc + - fort.73.nc + - fort.74.nc + - "PE*" diff --git a/adforce/eval/harvest.py b/adforce/eval/harvest.py new file mode 100644 index 00000000..8af52a00 --- /dev/null +++ b/adforce/eval/harvest.py @@ -0,0 +1,72 @@ +"""Harvest the minimal artifact set of a remote sweep to the laptop. + +Runs execute on ARCHER2/GCP; scoring runs locally. Per run only the small +artifacts travel (~MBs instead of the 5-8 GB ``fort.63.nc``): +``config.yaml`` (provenance), ``gauge_ts.parquet`` (the science), +``fort.61.nc`` (station cross-check), ``maxele.63.nc`` (run health), +``slurm.out`` (receipts), and the study's ``eval_manifest.json``. + +Run (hydra; config root adforce/eval/config/harvest_config.yaml):: + + python -m adforce.eval.harvest remote=archer2:/work/n02/n02/sithom/exp/eval study=kat-ida + python -m adforce.eval.harvest remote=... study=... dry_run=false +""" + +from __future__ import annotations + +import os +import subprocess +from typing import List + +import hydra +from omegaconf import DictConfig + +from . import constants as C + +#: files worth moving; everything else is excluded +INCLUDE = [ + "config.yaml", + "gauge_ts.parquet", + "fort.61.nc", + "maxele.63.nc", + "slurm.out", + "eval_manifest.json", +] + + +def rsync_command(remote: str, study: str, dest: str) -> List[str]: + """Build the rsync invocation (pure; unit-testable).""" + cmd = ["rsync", "-av", "--prune-empty-dirs", "--include=*/"] + cmd += [f"--include={name}" for name in INCLUDE] + cmd += ["--exclude=*", f"{remote.rstrip('/')}/{study}/", f"{dest.rstrip('/')}/"] + return cmd + + +def harvest(remote: str, study: str, dest: str = None, dry_run: bool = True) -> List[str]: + """Rsync one study's minimal artifacts into the local runs mirror. + + Args: + remote (str): ``host:/path/to/exp/eval`` (parent of study dirs). + study (str): Study name (the ``/`` subtree to pull). + dest (str): Local mirror root (default ``data/comp/runs/``). + dry_run (bool): Print the command instead of running it. + + Returns: + List[str]: The rsync argv (executed unless ``dry_run``). + """ + dest = dest or os.path.join(C.COMP_DATA_PATH, "runs", study) + cmd = rsync_command(remote, study, dest) + print(" ".join(cmd)) + if not dry_run: + os.makedirs(dest, exist_ok=True) + subprocess.run(cmd, check=True) + return cmd + + +@hydra.main(version_base=None, config_path="config", config_name="harvest_config") +def main(cfg: DictConfig) -> None: + harvest(str(cfg.remote), str(cfg.study), dest=cfg.dest, dry_run=cfg.dry_run) + + +if __name__ == "__main__": + main() diff --git a/adforce/eval/launch.py b/adforce/eval/launch.py new file mode 100644 index 00000000..9d9e98e9 --- /dev/null +++ b/adforce/eval/launch.py @@ -0,0 +1,262 @@ +"""Sweep orchestrator: launch the comparison matrix of adforce runs. + +Replaces ``rerun/adcirc/hist_sweep.sh``'s env-var matrix with the eval matrix +YAML (``config/matrix/``). Every historical cell routes through +``adforce.training.driver.drive_storm`` -- the per-storm input generation +that builds each storm's fort.15 with the correct cold-start/tidal window. +``wrap.stage_input_files``'s static decks (Katrina-pinned tidal window) are +never used here, which is what makes tide-on cells for the other 18 storms +legal at all. + +``dry_run: true`` is the default: ``python -m adforce.eval.launch`` prints +the plan table and exits without firing anything. + +Deliberate limits (fail loudly, not silently): +* ``swan`` cells are BLOCKED -- the training driver's input generation does + not stage SWAN decks (idealized-TC swan runs go through ``adforce.wrap``). +* the slurm backend is not implemented yet -- run the subprocess backend + inside one allocation (as ``hist_sweep.sh`` did on GCP), or sbatch a shell + loop; a throttled slurmpy backend is a known follow-up. + +All runner imports (wrap, training.driver, tcpips) stay function-local so the +laptop scoring path never needs the launch stack. + +Run (hydra; config root adforce/eval/config/launch_config.yaml):: + + python -m adforce.eval.launch study=kat-ida matrix=res_x_tide \ + 'storms=["Katrina 2005","Ida 2021"]' # plan only (dry_run) + python -m adforce.eval.launch study=kat-ida matrix=res_x_tide \ + 'storms=["Katrina 2005","Ida 2021"]' dry_run=false +""" + +from __future__ import annotations + +import hashlib +import json +import os +from typing import Optional + +import hydra +import pandas as pd +from omegaconf import DictConfig, OmegaConf + +from . import constants as C +from .cells import AXIS_PATHS, ConfigCell, cell_id, expand_matrix, storm_to_slug +from .status import RunStatus, run_status + + +def _cell_from_overrides(overrides: dict, forcing: str = None) -> ConfigCell: + """Matrix overrides (wrap-tree paths) -> ConfigCell.""" + get = lambda axis, default: overrides.get(AXIS_PATHS[axis], default) + tide = bool(get("tide", False)) + return ConfigCell( + resolution=str(get("resolution", "mid")), + tide=tide, + swan=bool(get("swan", False)), + forcing=forcing or ("both" if tide else "storm"), + ) + + +def _deck_exists(resolution: str) -> bool: + from adforce.constants import SETUP_PATH # light: paths only + + return os.path.exists(os.path.join(SETUP_PATH, f"fort.14.{resolution}")) + + +def _spinup(cfg, cell: ConfigCell) -> float: + if str(cfg.spinup_days) != "auto": + return float(cfg.spinup_days) + return 6.0 if (cell.tide or cell.forcing in ("tide", "both")) else 0.0 + + +def _runs_root(cfg) -> str: + root = str(cfg.runs_root) + return root if os.path.isabs(root) else os.path.join(C.PROJ_PATH, root) + + +def _expand_cells(cfg: DictConfig) -> list: + """Matrix cells (+ auto tide-only controls) as ``[(name, ConfigCell)]``. + + Shared by :func:`plan` and :func:`launch` so both see identical cells -- + including the control cells, which exist only here, not in the matrix. + """ + cells = [ + (name, _cell_from_overrides(dict(ov))) + for name, ov in expand_matrix(cfg.matrix) + ] + if cfg.controls: + # tide-only control run per unique (resolution, swan) with a tide-on + # cell -- the third leg of the tide-surge-interaction triple. + seen = {(c.resolution, c.swan) for _, c in cells if c.tide} + for res, swan in sorted(seen): + ctrl = ConfigCell(resolution=res, tide=True, swan=swan, forcing="tide") + cells.append((cell_id(ctrl), ctrl)) + return cells + + +def plan(cfg: DictConfig) -> pd.DataFrame: + """Expand the matrix into per-(cell, storm) rows with status and action. + + Pure (no launching, no heavy imports): usable from the laptop to preview + exactly what a remote launch would do. + + Returns: + pd.DataFrame: columns ``cell, storm, slug, run_dir, status, action, + reason`` -- action in {run, skip, blocked}. + """ + cells = _expand_cells(cfg) + storms = list(cfg.storms) if cfg.storms else list(C.STORMS) + rows = [] + for name, cell in cells: + for storm in storms: + slug = storm_to_slug(storm, C.STORMS[storm]) + run_dir = os.path.join(_runs_root(cfg), str(cfg.study), name, slug) + status = run_status(run_dir, cell) + action, reason = "run", "" + if cell.swan: + action, reason = "blocked", "driver input generation has no SWAN staging" + elif not _deck_exists(cell.resolution): + action, reason = ( + "blocked", + f"no fort.14.{cell.resolution} deck in adforce/setup (ARCHER2-only?)", + ) + elif status is RunStatus.FOREIGN and not cfg.overwrite: + action, reason = "blocked", "run dir holds a different config (FOREIGN)" + elif status in (RunStatus.SUCCESS, RunStatus.EXTRACTED) and cfg.skip_completed: + action, reason = "skip", "already successful" + rows.append( + dict( + cell=name, + storm=storm, + slug=slug, + run_dir=run_dir, + status=status.value, + action=action, + reason=reason, + ) + ) + return pd.DataFrame(rows) + + +def config_hash(resolved_cfg) -> str: + """sha256 of the full resolved run config: the provenance fingerprint + stored in the manifest so model-vs-model can detect drift in + un-overridden defaults between launches (see plan gap G5).""" + return hashlib.sha256( + OmegaConf.to_yaml(resolved_cfg, resolve=False).encode() + ).hexdigest() + + +def _write_manifest(runs_root: str, study: str, entries: dict) -> str: + path = os.path.join(runs_root, study, "eval_manifest.json") + os.makedirs(os.path.dirname(path), exist_ok=True) + existing = {} + if os.path.exists(path): + try: + existing = json.load(open(path)) + except Exception: + pass + existing.update(entries) + with open(path, "w") as f: + json.dump(existing, f, indent=2, sort_keys=True) + return path + + +def launch(cfg: DictConfig) -> pd.DataFrame: + """Execute the plan's ``run`` rows sequentially via the driver seam. + + Sequential on purpose (the hist_sweep disk-pressure lesson): each run is + reduced to ``gauge_ts.parquet`` and stripped before the next starts when + ``extract_after_run`` is set. + """ + import glob as _glob + + from adforce.wrap import get_default_config # runner imports stay local + + table = plan(cfg) + todo = table[table.action == "run"] + print(table.to_string(index=False)) + if todo.empty: + print("nothing to launch") + return table + + from tcpips.ibtracs import na_landing_tcs + + from adforce.training.driver import drive_storm + from adforce.training.storms import Storm + + target = na_landing_tcs() + base_cfg = get_default_config() + cells = dict(_expand_cells(cfg)) + gauges = None + + for row in todo.itertuples(): + cell = cells[row.cell] + i = int(row.slug.split("_")[0]) # archive index == IBTrACS enumeration + raw_times = pd.to_datetime(target.time[i].values, errors="coerce") + storm = Storm( + sid=target.sid[i].item(), + name=target.name[i].item(), + time=raw_times.dropna().to_pydatetime().tolist(), + ) + wrap_cfg = base_cfg.copy() + for axis, path in AXIS_PATHS.items(): + OmegaConf.update(wrap_cfg, path, getattr(cell, axis), merge=False) + os.makedirs(row.run_dir, exist_ok=True) + OmegaConf.update(wrap_cfg.files, "run_folder", row.run_dir) + OmegaConf.save(wrap_cfg, os.path.join(row.run_dir, "config.yaml")) + entry = { + f"{row.cell}/{row.slug}": dict( + status="launched", config_hash=config_hash(wrap_cfg) + ) + } + _write_manifest(_runs_root(cfg), str(cfg.study), entry) + + try: + drive_storm( + storm, + target.isel(storm=i), + row.run_dir, + wrap_cfg, + resolution=cell.resolution, + mode=cell.forcing, + spinup_days=_spinup(cfg, cell), + recommended_dt=cfg.recommended_dt, + ) + status = run_status(row.run_dir, cell) + if cfg.extract_after_run and status is RunStatus.SUCCESS: + from .extract import extract_run, gauge_frame + + gauges = gauge_frame() if gauges is None else gauges + df = extract_run(row.run_dir, gauges) + if not df.empty: + df.to_parquet( + os.path.join(row.run_dir, "gauge_ts.parquet"), index=False + ) + for pat in cfg.strip_after_extract: + for f in _glob.glob(os.path.join(row.run_dir, str(pat))): + os.remove(f) + status = run_status(row.run_dir, cell) + except Exception as e: # per-storm failures never stop the sweep + print(f"!!! FAILED {row.cell}/{row.slug}: {e}") + status = RunStatus.FAILED + entry[f"{row.cell}/{row.slug}"]["status"] = ( + status.value if isinstance(status, RunStatus) else str(status) + ) + _write_manifest(_runs_root(cfg), str(cfg.study), entry) + return plan(cfg) + + +@hydra.main(version_base=None, config_path="config", config_name="launch_config") +def main(cfg: DictConfig) -> None: + if cfg.dry_run: + table = plan(cfg) + print(table.to_string(index=False)) + n = (table.action == "run").sum() + print(f"\ndry_run=true: {n} runs would launch; pass dry_run=false to execute") + return + launch(cfg) + + +if __name__ == "__main__": + main() diff --git a/adforce/training/driver.py b/adforce/training/driver.py index 77d04426..5decd6d3 100644 --- a/adforce/training/driver.py +++ b/adforce/training/driver.py @@ -39,6 +39,84 @@ RUNS_PARENT_DIR = os.path.join(PROJ_PATH, "runs") os.makedirs(RUNS_PARENT_DIR, exist_ok=True) +#: Success line written to /slurm.out (same marker check_training_runs +#: and adforce.eval.status key on). +SUCCESS_MARKER = "Job completed successfully.\n" + + +def is_run_successful(run_directory: str) -> bool: + """True when ``/slurm.out`` carries the success marker.""" + slurm_path = os.path.join(run_directory, "slurm.out") + if not os.path.exists(slurm_path): + return False + try: + with open(slurm_path, "r") as slurm_out_file: + return any(SUCCESS_MARKER in line for line in slurm_out_file) + except (IOError, FileNotFoundError) as e: + print(f"Warning: Could not read {slurm_path}. Will attempt to rerun. Error: {e}") + return False + + +def drive_storm( + storm, + storm_ds, + run_directory: str, + cfg, + resolution: str = "mid", + mode: str = "storm", + spinup_days: float = 0.0, + recommended_dt: Optional[float] = None, + no_subprocess: bool = False, +) -> str: + """Generate inputs for ONE historical storm and (optionally) run it. + + The per-storm seam extracted from :func:`drive_all_adcirc` so callers + (``adforce.eval.launch``) control the run-directory naming and get the + path back. This is the ONLY sanctioned path for tide-on historical runs: + ``generate_adcirc_inputs`` builds each storm's fort.15 with the correct + cold-start/tidal window, unlike ``wrap.stage_input_files``'s static decks + (which hard-pin the Katrina 2005 window). + + Args: + storm (Storm): The storm metadata object. + storm_ds: The storm's IBTrACS row (``target_storms_ds.isel(storm=i)``). + run_directory (str): Absolute run directory to create/populate. + cfg: Base wrap config (``adforce.wrap.get_default_config()``). + resolution (str): Mesh resolution deck. + mode (str): ``storm`` | ``tide`` | ``both`` forcing. + spinup_days (float): Tidal spinup (6.0 for tidal modes in the sweeps). + recommended_dt (Optional[float]): ADCIRC timestep [s]. + no_subprocess (bool): Only generate inputs; do not run ADCIRC. + + Returns: + str: ``run_directory``. + """ + storm_name_safe = storm.name.upper().replace(" ", "_") + generate_adcirc_inputs( + storm, + storm_ds, + run_directory, + recommended_dt=recommended_dt, + resolution=resolution, + wind=mode in ("storm", "both"), + tides=mode in ("tide", "both"), + spinup_days=spinup_days, + ) + + storm_cfg = cfg.copy() + OmegaConf.update(storm_cfg.files, "run_folder", run_directory) + OmegaConf.update(storm_cfg, "name", f"{storm_name_safe}_{storm.year}") + # ASWIP converts pre_aswip_fort.22 -> fort.22 (NWS=20); tide-only runs + # have no met forcing -> no ASWIP. + OmegaConf.update(storm_cfg, "use_aswip", mode in ("storm", "both")) + # each storm is a local subprocess even when the sweep itself is one + # SLURM job + OmegaConf.update(storm_cfg, "use_slurm", False) + + if not no_subprocess: + setoff_subprocess_job_and_wait(run_directory, storm_cfg) + return run_directory + def drive_all_adcirc( test_single=False, @@ -156,22 +234,7 @@ def _match(st) -> bool: ) # --- Check for existing successful run (implements TODO #1) --- - slurm_path = os.path.join(run_directory, "slurm.out") - is_successful = False - if os.path.exists(slurm_path): - try: - with open(slurm_path, "r") as slurm_out_file: - for line in slurm_out_file: - if "Job completed successfully.\n" in line: - is_successful = True - break - except (IOError, FileNotFoundError) as e: - print( - f"Warning: Could not read {slurm_path}. Will attempt to rerun. Error: {e}" - ) - is_successful = False - - if is_successful: + if is_run_successful(run_directory): print(f"Run {run_directory} already completed successfully. Skipping.") continue # Skip to the next storm elif os.path.exists(run_directory): @@ -184,42 +247,19 @@ def _match(st) -> bool: print( f"\n--- Processing Storm {i+1}/{len(target_storms)}: {storm.name} {storm.year} ---" ) - - # 1. Generate ADCIRC input files (fort.15, pre_aswip_fort.22, fort.13) print(f"Generating inputs in: {run_directory}") - generate_adcirc_inputs( + drive_storm( storm, target_storms_ds.isel(storm=i), run_directory, - recommended_dt=recommended_dt, + cfg, resolution=resolution, - wind=mode in ("storm", "both"), - tides=mode in ("tide", "both"), + mode=mode, spinup_days=spinup_days, + recommended_dt=recommended_dt, + no_subprocess=test_nosubprocess, ) - - # 2. Create a storm-specific config to pass to subprocess - # Start with the default config - storm_cfg = cfg.copy() - - # Set the specific run folder and name for logging - OmegaConf.update(storm_cfg.files, "run_folder", run_directory) - OmegaConf.update(storm_cfg, "name", f"{storm_name_safe}_{storm.year}") - - # Tell the subprocess runner to execute ASWIP - # This converts 'pre_aswip_fort.22' to 'fort.22' with the NWS=20 - # format. Tide-only runs have no met forcing -> no ASWIP. - OmegaConf.update(storm_cfg, "use_aswip", mode in ("storm", "both")) - - # Ensure we're NOT using SLURM for this subprocess - # (The main script is one SLURM job, but each storm is a local subprocess) - OmegaConf.update(storm_cfg, "use_slurm", False) - - # 3. Run the simulation (ASWIP, adcprep, padcirc) - print(f"Running ADCIRC simulation for {storm.name} via subprocess...") - # This function will chdir into run_directory if not test_nosubprocess: - setoff_subprocess_job_and_wait(run_directory, storm_cfg) print(f"✅ Successfully completed run for {storm.name} {storm.year}") else: print(f"✅ Successfully made inputs for {storm.name} {storm.year}") diff --git a/tcpips/constants.py b/tcpips/constants.py index e108abfc..a9c095f6 100755 --- a/tcpips/constants.py +++ b/tcpips/constants.py @@ -25,6 +25,22 @@ CONFIG_PATH: str = os.path.join(SRC_PATH, "config") DATA_PATH: str = os.path.join(PROJECT_PATH, "data") + +def _ensure_dir(path) -> None: + """``mkdir -p`` that tolerates dangling symlinks. + + ``data/era5`` is a symlink onto an external volume; when the volume is + unmounted, ``os.makedirs(..., exist_ok=True)`` on it (or any child) still + raises FileExistsError, which used to make ``import tcpips`` fail on any + machine without the drive attached (breaking transitive importers like + ``adforce.training.driver``). Directory creation is best-effort here -- + the code that actually writes will surface a real error at write time. + """ + try: + os.makedirs(path, exist_ok=True) + except OSError: + pass + # General data from e.g. paper or cmip etc. GOM: Tuple[float] = (25.443701, -90.013120) # Centre of Gulf of Mexico, lat, lon MONTHS: List[str] = [ # 3 letter month names @@ -58,54 +74,54 @@ # Directories for ERA5 data processing steps: make them ahead of time. ERA5_PATH: str = os.path.join(DATA_PATH, "era5") # main data folder -os.makedirs(ERA5_PATH, exist_ok=True) +_ensure_dir(ERA5_PATH) ERA5_RAW_PATH: str = os.path.join(ERA5_PATH, "raw") # download data here -os.makedirs(ERA5_RAW_PATH, exist_ok=True) +_ensure_dir(ERA5_RAW_PATH) ERA5_REGRIDDED_PATH: str = os.path.join(ERA5_PATH, "regrid") # regridded data here -os.makedirs(ERA5_REGRIDDED_PATH, exist_ok=True) +_ensure_dir(ERA5_REGRIDDED_PATH) ERA5_PI_OG_PATH: str = os.path.join(ERA5_PATH, "pi_og") # pi no bias correction before -os.makedirs(ERA5_PI_OG_PATH, exist_ok=True) # potential intensity on original grid +_ensure_dir(ERA5_PI_OG_PATH) # potential intensity on original grid ERA5_PI_PATH: str = os.path.join(ERA5_PATH, "pi") # pi on new grid -os.makedirs(ERA5_PI_PATH, exist_ok=True) # potential intensity on new grid +_ensure_dir(ERA5_PI_PATH) # potential intensity on new grid ERA5_PS_OG_PATH: str = os.path.join( ERA5_PATH, "ps_og" ) # potential size on original grid -os.makedirs(ERA5_PS_OG_PATH, exist_ok=True) # potential size on +_ensure_dir(ERA5_PS_OG_PATH) # potential size on ERA5_PS_PATH = os.path.join(ERA5_PATH, "ps") # potential size -os.makedirs(ERA5_PS_PATH, exist_ok=True) # potential size on new grid +_ensure_dir(ERA5_PS_PATH) # potential size on new grid ERA5_PRODUCTS_PATH = os.path.join(ERA5_PATH, "products") # products from ERA5 -os.makedirs(ERA5_PRODUCTS_PATH, exist_ok=True) # products from ERA5 +_ensure_dir(ERA5_PRODUCTS_PATH) # products from ERA5 ERA5_FIGURE_PATH: str = os.path.join(FIGURE_PATH, "era5") # figures from ERA5 -os.makedirs(ERA5_FIGURE_PATH, exist_ok=True) # figures from ERA5 +_ensure_dir(ERA5_FIGURE_PATH) # figures from ERA5 # Directories for CMIP6 data processing steps: make them ahead of time. CMIP6_PATH: str = os.path.join(DATA_PATH, "cmip6") # main data folder -os.makedirs(CMIP6_PATH, exist_ok=True) +_ensure_dir(CMIP6_PATH) RAW_PATH: str = os.path.join(CMIP6_PATH, "raw") # download data here -os.makedirs(RAW_PATH, exist_ok=True) +_ensure_dir(RAW_PATH) REGRIDDED_PATH: str = os.path.join(CMIP6_PATH, "regridded") # regridded data here -os.makedirs(REGRIDDED_PATH, exist_ok=True) +_ensure_dir(REGRIDDED_PATH) CDO_PATH: str = os.path.join(CMIP6_PATH, "regrid") -os.makedirs(CDO_PATH, exist_ok=True) +_ensure_dir(CDO_PATH) BIAS_CORRECTED_PATH = os.path.join( CMIP6_PATH, "bias_corrected" ) # bias corrected data here -os.makedirs(BIAS_CORRECTED_PATH, exist_ok=True) +_ensure_dir(BIAS_CORRECTED_PATH) PI_PATH: str = os.path.join(CMIP6_PATH, "pi") # pi no bias correction before -os.makedirs(PI_PATH, exist_ok=True) +_ensure_dir(PI_PATH) PI2_PATH: str = os.path.join(CMIP6_PATH, "pi2") # pi after cdo regridding -os.makedirs(PI2_PATH, exist_ok=True) +_ensure_dir(PI2_PATH) PI3_PATH: str = os.path.join(CMIP6_PATH, "pi3") # pi after temp profile fix -os.makedirs(PI3_PATH, exist_ok=True) +_ensure_dir(PI3_PATH) PI4_PATH: str = os.path.join(CMIP6_PATH, "pi4") # pi with different temp profile fix -os.makedirs(PI4_PATH, exist_ok=True) +_ensure_dir(PI4_PATH) BC_PI_PATH: str = os.path.join(CMIP6_PATH, "bc_pi") # pi after bias correction -os.makedirs(BC_PI_PATH, exist_ok=True) +_ensure_dir(BC_PI_PATH) PS_PATH = os.path.join(CMIP6_PATH, "ps") # potential size -os.makedirs(PS_PATH, exist_ok=True) # potential size +_ensure_dir(PS_PATH) # potential size BIAS_PATH = os.path.join(CMIP6_PATH, "bias") # bias data -os.makedirs(BIAS_PATH, exist_ok=True) # bias data +_ensure_dir(BIAS_PATH) # bias data # Constants for converting CMIP6 variables to PI input variables CONVERSION_NAMES: Dict[str, str] = { diff --git a/tests/test_eval.py b/tests/test_eval.py index 88cfe7f5..9f7e8f04 100644 --- a/tests/test_eval.py +++ b/tests/test_eval.py @@ -723,6 +723,118 @@ def test_rundir_source_prefers_cheapest_artifact(tmp_path): assert src == "gauge_ts" and series.max() == 3.0 +# --------------------------------------------------------------------------- # +# Launch planning + harvest (Commit E). plan() is pure -- no ADCIRC, no +# network; everything runs against tmp_path run dirs. +# --------------------------------------------------------------------------- # +def _launch_cfg(tmp_path, **over): + from hydra import compose, initialize + + overrides = [f"runs_root={tmp_path}", "study=t", 'storms=["Katrina 2005","Ida 2021"]'] + overrides += [f"{k}={v}" for k, v in over.items()] + with initialize(version_base=None, config_path="../adforce/eval/config"): + return compose(config_name="launch_config", overrides=overrides) + + +def test_launch_plan_statuses_and_actions(tmp_path): + from omegaconf import OmegaConf + + from adforce.eval.launch import plan + from adforce.eval.status import SUCCESS_MARKER + + cfg = _launch_cfg(tmp_path) + table = plan(cfg) + # res_x_tide: 4 cells x 2 storms; low/mid decks ship in adforce/setup + assert len(table) == 8 + assert set(table.cell) == { + "res-low_tide-off", + "res-low_tide-on", + "res-mid_tide-off", + "res-mid_tide-on", + } + assert (table.status == "missing").all() and (table.action == "run").all() + + # fabricate a successful run -> skipped on the next plan + run = tmp_path / "t" / "res-mid_tide-off" / "152_KATRINA_2005" + run.mkdir(parents=True) + OmegaConf.save( + OmegaConf.create( + {"adcirc": {"resolution": {"value": "mid"}, "tide": {"value": False}, "swan": {"value": False}}} + ), + str(run / "config.yaml"), + ) + (run / "slurm.out").write_text(SUCCESS_MARKER + "\n") + table = plan(cfg) + row = table[(table.cell == "res-mid_tide-off") & (table.slug == "152_KATRINA_2005")] + assert row.action.item() == "skip" and row.status.item() == "success" + + # a FOREIGN dir (wrong resolution inside) blocks unless overwrite + OmegaConf.save( + OmegaConf.create( + {"adcirc": {"resolution": {"value": "low"}, "tide": {"value": False}, "swan": {"value": False}}} + ), + str(run / "config.yaml"), + ) + table = plan(cfg) + row = table[(table.cell == "res-mid_tide-off") & (table.slug == "152_KATRINA_2005")] + assert row.action.item() == "blocked" and row.status.item() == "foreign" + + +def test_launch_plan_controls_and_blocked_cells(tmp_path): + from adforce.eval.launch import plan + + # controls add one tide-only cell per resolution with a tide-on cell + table = plan(_launch_cfg(tmp_path, controls=True)) + assert "res-low_tide-on_swan-off_f-tide" in set(table.cell) + assert "res-mid_tide-on_swan-off_f-tide" in set(table.cell) + assert len(table) == (4 + 2) * 2 + + # high resolution (no local deck) and swan cells are blocked, loudly + cfg = _launch_cfg(tmp_path, matrix="archive_default") + cfg.matrix.cells = [ + dict(name="hi", overrides={"adcirc.resolution.value": "high"}), + dict(name="sw", overrides={"adcirc.swan.value": True}), + ] + cfg.matrix.baseline = "hi" + table = plan(cfg) + assert (table.action == "blocked").all() + reasons = " ".join(table.reason) + assert "fort.14.high" in reasons and "SWAN" in reasons + + +def test_config_hash_and_harvest_command(tmp_path): + from omegaconf import OmegaConf + + from adforce.eval.harvest import rsync_command + from adforce.eval.launch import config_hash + + a = OmegaConf.create({"x": 1, "y": {"z": "s"}}) + assert config_hash(a) == config_hash(OmegaConf.create({"x": 1, "y": {"z": "s"}})) + assert config_hash(a) != config_hash(OmegaConf.create({"x": 2, "y": {"z": "s"}})) + + cmd = rsync_command("host:/work/exp/eval", "study1", str(tmp_path)) + assert cmd[0] == "rsync" and cmd[-2] == "host:/work/exp/eval/study1/" + assert "--include=gauge_ts.parquet" in cmd and "--exclude=*" in cmd + + +def test_driver_seam_importable(): + """The extracted per-storm seam exists with the expected signature. + + Skipped where the training driver's HPC-only deps (adcircpy/stormevents) + are not installed -- the same reason pytest.ini --ignores the module.""" + import inspect + + try: + from adforce.training.driver import drive_storm, is_run_successful + except ModuleNotFoundError as e: + pytest.skip(f"training-driver dependency absent locally: {e.name}") + + params = list(inspect.signature(drive_storm).parameters) + assert params[:4] == ["storm", "storm_ds", "run_directory", "cfg"] + assert "mode" in params and "spinup_days" in params + assert is_run_successful("/nonexistent") is False + + # --------------------------------------------------------------------------- # # Matrix grammar + model-vs-model (Commit D). The reproduction tests pin the # ported pairs.py against the published rerun/results artifacts and are From 6f0f49a65b6f9bdc4dadcc6a8d7823a5ffa5585d Mon Sep 17 00:00:00 2001 From: "Simon D.A. Thomas" Date: Sun, 16 Aug 2026 20:37:47 +0100 Subject: [PATCH 07/17] finished eval move --- adforce/eval/README.md | 63 +++++++++++++++++++++++++++++++++++++++--- 1 file changed, 59 insertions(+), 4 deletions(-) diff --git a/adforce/eval/README.md b/adforce/eval/README.md index 60e2df72..4b2cc44e 100644 --- a/adforce/eval/README.md +++ b/adforce/eval/README.md @@ -1,4 +1,8 @@ -# comp — historical surge validation against tide gauges +# adforce.eval — ADCIRC configuration evaluation against tide gauges + +(Formerly the top-level `comp/` package; on-disk caches stay under +`data/comp/` because the raw CO-OPS responses are not reproducibly +re-downloadable and the utide fits are expensive.) Validates the historical ADCIRC storm-surge simulations (the SurgeNet training set, 228 IBTrACS North-Atlantic landfalling TCs on the EC95d mesh, published on Hugging @@ -19,7 +23,7 @@ For each storm: NOAA CO-OPS water-level gauge in the storm's region box — NW Gulf (Texas → Florida panhandle; the New Orleans and Galveston study regions) for Gulf storms, or the Florida peninsula (Key West → Fernandina Beach; the Miami study region) for the - Atlantic/Florida storms (`FLORIDA_STORMS` in `comp/constants.py`); + Atlantic/Florida storms (`FLORIDA_STORMS` in [`constants.py`](constants.py)); 3. de-tide the gauge record with a robust [`utide`](https://github.com/wesleybowman/UTide) harmonic fit on the storm's calendar year (falls back to CO-OPS `predictions`); 4. score skill on three axes: @@ -79,12 +83,12 @@ sharply at lag 0 and decays to zero within two days). `adforce.eval.sensitivity` is stable across the clean-filter cut-offs (r 0.81–0.89) and the node-selection knobs (r changes <0.001 with the wet-depth threshold). -Unit + regression tests live in [`../tests/test_comp.py`](../tests/test_comp.py) (skill +Unit + regression tests live in [`../../tests/test_eval.py`](../../tests/test_eval.py) (skill metrics, time-series alignment, the valid/clean gating, the LaTeX table, and pinned headline numbers + negative-control separations): ```bash -python -m pytest tests/test_comp.py -o addopts="" +python -m pytest tests/test_eval.py -o addopts="" ``` ## Configuration @@ -112,6 +116,57 @@ The slight low bias is consistent with the omitted wave setup/runup and medium m resolution (datum and node-sampling effects checked and found minor); over-predictions concentrate at shallow semi-enclosed bay/pass gauges during direct landfalls. +## Comparing adforce configurations (general framework) + +Beyond the historical HF-archive validation above, the submodule compares +**any** adforce configurations (resolution low/mid/high, tides on/off, SWAN) +against gauges and against each other, driven entirely by YAML under +[`config/`](config/): + +- **matrix** ([`config/matrix/`](config/matrix/)) — the comparison grid: + cartesian `axes` over wrap-config paths (`adcirc.resolution.value`, + `adcirc.tide.value`, ...) minus `exclude`, plus explicit `cells`, with a + designated `baseline` for model-vs-model. +- **launch** (`python -m adforce.eval.launch`) — expands the matrix into + per-(cell, storm) runs under `////`, routes + every historical run through the training driver's per-storm input + generation (correct tidal windows — never the Katrina-pinned static + decks), skips completed runs, refuses FOREIGN directories, and records a + full-config sha256 per run in `eval_manifest.json`. `dry_run=true` is the + default: it prints the plan table and fires nothing. `controls=true` adds + the tide-only runs of the tide-surge-interaction triple. +- **extract** (`python -m adforce.eval.extract`) — remote-side reducer: + `fort.63.nc` (5–8 GB) → `gauge_ts.parquet` per run (long format + `storm, sid, gauge, time, zeta`), same node selection as the validation. +- **harvest** (`python -m adforce.eval.harvest`) — rsyncs the minimal + artifact set (`config.yaml`, `gauge_ts.parquet`, `fort.61.nc`, + `maxele.63.nc`, `slurm.out`, manifest) to the laptop. +- **pairs** (`python -m adforce.eval.pairs`) — model-vs-model: + `resolution_bias_table` (reproduces `rerun/results/resolution_bias.csv` + exactly) and `tide_surge_interaction` + (`peak(both) − peak(storm) − peak(tide)`; reproduces the published + low-res rows of `rerun/results/tide_surge_interaction.csv` 683/683). + Both reproductions are pinned by tests. + +Example — low vs mid × tide on/off for two storms: + +```bash +# remote (ARCHER2/GCP): preview, then launch +python -m adforce.eval.launch study=kat-ida matrix=res_x_tide \ + 'storms=["Katrina 2005","Ida 2021"]' # dry run (default) +python -m adforce.eval.launch study=kat-ida matrix=res_x_tide \ + 'storms=["Katrina 2005","Ida 2021"]' dry_run=false controls=true +# laptop: +python -m adforce.eval.harvest remote=archer2:/work/.../exp/eval study=kat-ida dry_run=false +python -m adforce.eval.pairs action=interaction res=mid \ + tide_series=... both_series=... out=data/comp/out/kat-ida_tsi.csv +``` + +Caveat (tide-on scoring): comparing a tide-on run against raw gauge water +level needs the datum/steric offset handled (pre-storm-mean alignment) and +skew surge as the headline metric; that obs-side path is staged follow-up +work — do not publish tide-on skill numbers before it lands. + ## Dependencies `huggingface_hub`, `utide`, `xarray`, `scipy`, `pandas`, `pyarrow`, `requests`, `matplotlib`, From e1e0eeee9bc804b23fafb0c7f222672a24eccdae Mon Sep 17 00:00:00 2001 From: "Simon D.A. Thomas" Date: Sun, 16 Aug 2026 21:16:53 +0100 Subject: [PATCH 08/17] eval Phase 2a: tide-on obs machinery, tide validation, TWL/skew scoring - eval/detide.py: windowed CO-OPS fetches (polite on cache miss), align_pair (pre-forcing datum alignment, offset reported not hidden), skew_surge_peak promoted from detide_sensitivity (re-exported there; ragged-index guard) - eval/tidecheck.py + config: tide-only runs vs CO-OPS predictions -- datum_offset_m, amp_ratio, sub-hourly cross-correlation lag, r, rmse - eval/twl.py + config: tide-on scoring -- storm+tide TWL vs raw gauge WL (aligned peaks, hydrograph skill) with SKEW SURGE as headline (model skew vs its own tide-only run, obs skew vs the prediction, inner-joined) - tests: synthetic recovery of amp/lag/offset, alignment, phase-insensitive skew, cached-predictions parse, gappy-obs TWL end-to-end (57 passed) First results (data/comp/out/tidecheck_{low,mid}.csv, twl_low.csv): * low tides: amp_ratio median 1.45, lag -66 min, r 0.82 (1019 pairs) * mid tides: 1.41, -42 min, r 0.855 (969 pairs) -> the ~40% coastal tide over-amplification is RESOLUTION-INDEPENDENT (forcing/friction, not mesh); TX semi-enclosed bays are structurally broken at both (amp ~5.8x, r~0.1) * low TWL: aligned peak bias -0.02 m but skew sim 0.07 vs obs 0.36 m -- tide error compensates surge deficit; skew exposes it => gate verdict: fix EC95d tidal forcing before any tide-on skill claims. Co-Authored-By: Claude Fable 5 --- adforce/eval/config/tidecheck_config.yaml | 10 ++ adforce/eval/config/twl_config.yaml | 12 ++ adforce/eval/detide.py | 113 +++++++++++++ adforce/eval/detide_sensitivity.py | 25 +-- adforce/eval/tidecheck.py | 154 ++++++++++++++++++ adforce/eval/twl.py | 184 ++++++++++++++++++++++ tests/test_eval.py | 104 ++++++++++++ 7 files changed, 580 insertions(+), 22 deletions(-) create mode 100644 adforce/eval/config/tidecheck_config.yaml create mode 100644 adforce/eval/config/twl_config.yaml create mode 100644 adforce/eval/detide.py create mode 100644 adforce/eval/tidecheck.py create mode 100644 adforce/eval/twl.py diff --git a/adforce/eval/config/tidecheck_config.yaml b/adforce/eval/config/tidecheck_config.yaml new file mode 100644 index 00000000..7fbd5641 --- /dev/null +++ b/adforce/eval/config/tidecheck_config.yaml @@ -0,0 +1,10 @@ +# Root config for python -m adforce.eval.tidecheck (tide-only run validation +# vs CO-OPS astronomical predictions -- the gate before tide-on surge skill). +hydra: + output_subdir: null + run: + dir: . + +series: ??? # gauge-series parquet path or glob (tide-only runs) +label: tide # names the default output data/comp/out/tidecheck_