diff --git a/.claude/skills/ddrs-dev/SKILL.md b/.claude/skills/ddrs-dev/SKILL.md index f725d0bf..ffc42e0a 100644 --- a/.claude/skills/ddrs-dev/SKILL.md +++ b/.claude/skills/ddrs-dev/SKILL.md @@ -27,6 +27,7 @@ config reference, test-authoring patterns, and the current research status. | Diagnose a failed, hung, or wrong-result run | `references/traps.md` §Symptom → trap | | Check whether a question is already settled | `references/research-status.md` | | Build or change the training gauge CSV / population | `references/gauge-population.md` | +| Launch / watch / resume / audit a training run, or use any `ddrs` subcommand | skill `ddrs-run` | | Plot or interpret eval output | skill `ddrs-eval-plots` | ## Contents @@ -114,6 +115,7 @@ Full commands and rationale in `references/testing.md`. | `config/**/*.yaml` only | **D** | `ddrs plan --config … --workspace …` exits 0, no drift | | `.github/workflows/ci.yml`, `.githooks/**` | — | push to a branch and let the PR run vet it; there is no local gate for workflow changes | | `examples/juniata/**` (bundle, config, README) | — | `cargo test --test juniata_bundle` (never skips — bundle is committed), `cargo test --release --test juniata_acceptance` (end-to-end metric floors), and `ddrs --config examples/juniata/ddrs.yaml plan` exits 0 from the repo root | +| `src/adjacency/gridded.rs`, `resolve_or_build_gridded`, `examples/juniata_gridded/**`, the `Qr` axis-order sniff in `src/data/store/icechunk.rs` | **C** | `cargo test --test gridded_bundle` (never skips), `cargo test --test hourly_streamflow` (time-major fixture), `cargo test --release --test gridded_acceptance`, then Tier C | | Plotting / analysis scripts only | — | no gate | | `epochs`, `learning_rate`, `batch_size`, loss weights within documented ranges | — | no gate | @@ -209,9 +211,49 @@ Two deviations from DDR's bundle: `data/statistics/*.json` is **committed** `ddrs.yaml` ignore. Regenerate the bundle in the ddr repo (`extract_bundle.py`), then re-copy `data/` plus the generated statistics JSON. +## Gridded (ISIMIP DDM30) routing — `examples/juniata_gridded/` (2026-09-09) + +ddrs consumes DDR PR #194's gridded products as a **data source**; it ports +none of the raster/regridding engine. `data_sources.gridded_network` points at +DDR's sub-reach adjacency zarr (`order`, `parent_cell`, `indices_0/1`, +`length_m`, `slope`) and `ddrs plan` relabels it into the existing subdivided +store layout (parent = cell, pieces = the cell's contiguous sub-reaches, gauge +read at the last piece) via `adjacency::gridded::GriddedNetwork` + +`cache::resolve_or_build_gridded`. Attributes and Q′ are keyed on the cell id, +so `parent_order` lookups and `pieces_per_row_divisor` (Q′/k) need no change. +The gauge CSV's `cell` column is a serde alias of `COMID`. Guards: exclusive +with `geospatial_fabric` and the explicit zarr pair; `params.subdivision.enabled` +is rejected with it. Spec: +`docs/superpowers/specs/2026-09-08-ddrs-gridded-routing-design.md`. + +```bash +target/release/ddrs --config examples/juniata_gridded/ddrs.yaml plan +target/release/ddrs --config examples/juniata_gridded/ddrs.yaml run --workflow train-and-test --backend cpu +``` + +Verified 2026-09-09 (4 cells / 27 sub-reaches, 30 epochs, CPU, ~15 s): baseline +NSE 0.594 / KGE 0.672 (DDR 0.594 / 0.671); routed NSE 0.751 / KGE 0.730 at seed +42, seeds {42, 7, 123} span 0.744–0.758 / 0.728–0.732 (DDR 0.795 / 0.737 — RNG +windows plus ddrs's per-cell KAN head vs DDR's per-sub-reach `log10_uparea`). +CONUS: `config/sources/conus-gridded.yaml` + `config/experiments/gridded_conus.yaml` +(620 gauges from `scripts/snap_gridded_gauges.py`, DDR's population); `plan` +builds 5,198 cells → 91,867 sub-reaches in 0.6 s, baseline median NSE 0.354 / +KGE 0.505 over 1995-10-01..1997-09-30 (DDR's 0.390 includes its `da_ratio` +output correction, which ddrs does not apply). + +Two things this work found, both fixed: (1) `GageSubgraph::upstream_comids` +returned one entry per piece on subdivided stores, so the summed-Q′ baseline +counted each divide `m` times — now deduplicated by id (also affected MERIT +`params.subdivision` baselines); (2) DDR's gridded Q′ stores are +`Qr(time, divide_id)`, and zarrs returns fill values for an out-of-range subset +instead of erroring, so the old fixed-axis reader silently read NaN — the +reader now sniffs the axis order (traps.md T11). + ## Maintenance -This skill and `ddrs-eval-plots` are the only two skills in this repo. When a run, +Three skills live in this repo: this one (build / code / configure / test / debug), +`ddrs-run` (launch, watch, resume, and audit runs; the full `ddrs` command +reference), and `ddrs-eval-plots` (visualize and interpret run output). When a run, eval, or export completes, update the relevant section here in the same session that produced the knowledge — do not leave it only in a findings doc. If a rule here is superseded, correct it in place with the new nuance rather than deleting it. diff --git a/.claude/skills/ddrs-dev/references/build-and-env.md b/.claude/skills/ddrs-dev/references/build-and-env.md index 7f4e3263..2acd2b72 100644 --- a/.claude/skills/ddrs-dev/references/build-and-env.md +++ b/.claude/skills/ddrs-dev/references/build-and-env.md @@ -25,6 +25,13 @@ Verified against `Cargo.toml` and source on 2026-07-30. If a system HDF5 leaks into the static build: `unset HDF5_DIR; unset NETCDF_DIR`. +- **CUDA newer than `cudarc`'s table** (13.3.1 on this host since 2026-09 vs + cudarc 0.19.7's max 13.2): the build script panics with `Unsupported cuda + toolkit version`. Set `CUDARC_CUDA_VERSION=13020` — a gitignored + `.cargo/config.toml` with `[env] CUDARC_CUDA_VERSION = "13020"` is the + per-checkout fix (traps.md T12). Cached `target/` dirs mask it until the + build script reruns. + ## Fork pins All **13** `burn-*` crates must resolve from `github.com/taddyb/burn` branch diff --git a/.claude/skills/ddrs-dev/references/config.md b/.claude/skills/ddrs-dev/references/config.md index 44f344a0..5f0db599 100644 --- a/.claude/skills/ddrs-dev/references/config.md +++ b/.claude/skills/ddrs-dev/references/config.md @@ -4,9 +4,13 @@ Struct: `src/config.rs::Config`. Loaded via `Config::from_yaml_file_with_mode(path, ConfigMode::Training|Testing)`. Six top-level sections. Verified against source 2026-07-30. -**No `deny_unknown_fields` except `DisaggregationSection`** (added 2026-08-03) — -everywhere else a typo'd key silently takes its default instead of erroring. This -is the single most common cause of "my config change did nothing". +**Every section sets `deny_unknown_fields`** as of 2026-09-09 (`ConfigRaw`, +`DataSources`, `Experiment`, `KanHeadConfigSection`, `LossConfig`, `ParamsRaw`, +`TestingOverridesRaw`, plus `DisaggregationSection` and `Subdivision`, which had it +already), so a typo'd key is a load error naming the key. Historically only the last +two had it and everywhere else a typo silently took its default, which was the +single most common cause of "my config change did nothing". On a binary older than +2026-09-09, that silence is still the first thing to suspect. ## Contents @@ -37,10 +41,14 @@ to §`params:` if you are looking for `tau` (it is not a routing sub-step count) (feature-concatenated on COMID, NaN-filled, `deserialize_one_or_many_paths`). An empty list is a hard error. -Adjacency rule: provide **either** both `conus_adjacency` + `gages_adjacency`, -**or** `geospatial_fabric` (managed build into `.ddrs/adjacency//`). -Exactly one of the pair ⇒ error; neither source ⇒ `"adjacency sources are missing"`. -For multi-layer gpkg set `geospatial_fabric_layer` (participates in the cache key). +Adjacency rule: provide **exactly one** of (a) both `conus_adjacency` + +`gages_adjacency`, (b) `geospatial_fabric` (MERIT managed build into +`.ddrs/adjacency//`), (c) `gridded_network` (DDR's DDM30 sub-reach +adjacency zarr, managed build via `resolve_or_build_gridded`; added 2026-09-09). +Exactly one of the pair ⇒ error; neither source ⇒ `"adjacency sources are missing"`; +`gridded_network` with (a) or (b) ⇒ error. For multi-layer gpkg set +`geospatial_fabric_layer` (participates in the cache key). `params.subdivision.enabled` +is rejected alongside `gridded_network` (the store is already split by DDR). `aorc_precip` is required whenever `kan_head.disaggregation:` is present — see below. @@ -177,6 +185,10 @@ Four validators run at `Config::from_yaml_file`, plus one at dataset open. | `validate_data_sources` | one of the adjacency pair | `` "`gages_adjacency` is missing" `` | | | neither adjacency nor fabric | `"adjacency sources are missing"` | | | `geospatial_fabric_layer` on a non-gpkg | `"geospatial_fabric_layer"` + `".gpkg"` | +| | `gridded_network` + `geospatial_fabric` | `"gridded_network"` + `"geospatial_fabric"` | +| | `gridded_network` + explicit adjacency pair | `"gridded_network"` + `"conus_adjacency"` | +| `validate_subdivision` | `subdivision.enabled: true` + `gridded_network` | `"params.subdivision"` + `"gridded_network"` | +| `validate_geodataset` | `geodataset:` contradicting the adjacency source (`ddm30` with `geospatial_fabric`, `merit` with `gridded_network`) | `"geodataset"` + the source key. Absent ⇒ inferred; explicit adjacency paths ⇒ any label allowed | | `validate_leakance` | `use_leakance` + `use_cuda_graphs` | both key names | | `validate_ddr_match` | `use_cuda_graphs: true` without the deprecated `ddr_match: true` | `"use_cuda_graphs: true` requires the DEPRECATED `ddr_match: true"` | | `validate_disagg_pretrained` | `freeze: true` without `pretrained_checkpoint` | `"freeze: true requires pretrained_checkpoint"` | diff --git a/.claude/skills/ddrs-dev/references/testing.md b/.claude/skills/ddrs-dev/references/testing.md index e1835bd3..15a21517 100644 --- a/.claude/skills/ddrs-dev/references/testing.md +++ b/.claude/skills/ddrs-dev/references/testing.md @@ -48,6 +48,15 @@ debug_assertions** (a debug train takes minutes), so it only runs with readers when you want end-to-end confirmation. Verified 2026-09-02: NSE 0.7903 / KGE 0.8810 / baseline 0.6947 in 18.3 s. +```bash +cargo test --release --test gridded_acceptance -- --nocapture # ~30 s after build +``` + +The gridded (DDM30) twin on `examples/juniata_gridded/`: floors NSE ≥ 0.70 / +KGE ≥ 0.69 (seed-42 reference 0.751 / 0.730; seeds {42, 7, 123} span +0.744–0.758 / 0.728–0.732), baseline NSE ∈ [0.57, 0.62] (deterministic 0.594, +matches DDR). The baseline band is the one that catches a wrong `Qr` axis order. + ### Tier B — KAN head `src/nn/`, `Cargo.toml` rskan tag @@ -77,12 +86,12 @@ Exit 0, no drift warnings. `.github/workflows/ci.yml`, on every PR and master push (no path filters: required checks must never be skipped): job `test` = debug -`cargo test --features fixtures --no-fail-fast` (13 min warm, 12–17 min -cold); job `acceptance` = release `compare_ddr_sandbox` + -`juniata_acceptance` (17 min warm, ~25 min cold). Warm is only modestly -faster because `rust-cache` keeps dependency artifacts only: the ddrs -crate and every test binary rebuild on each run (measured 2026-09-04, -runs 33817851139 cold and 33830503587 warm). The acceptance job +`cargo test --features fixtures --no-fail-fast` (~12–17 min cold; warm is +shorter); job `acceptance` = release `compare_ddr_sandbox` + +`juniata_acceptance` + `gridded_acceptance` (~25 min cold; warm is shorter; +both acceptance tests self-skip under `debug_assertions`, so the debug `test` +job never runs them — the release job is their only CI coverage). The +acceptance job builds with `CARGO_PROFILE_RELEASE_LTO=false` (env override in the workflow only): the thin-LTO link of the test binary was 30 of 40 minutes on the 2-core runner while the 30 training epochs took ~13 s, so the training is @@ -97,8 +106,10 @@ local tier gates. Both jobs install a CUDA toolkit only because compilation requires it (build-and-env.md). Local hook: `git config core.hooksPath .githooks` -enables `.githooks/pre-push` (runs `ddr_sandbox_match`; bypass with -`git push --no-verify`). +enables `.githooks/pre-push`, which runs `ddr_sandbox_match` + +`gridded_bundle` in one cargo invocation (~0.4 s warm; both never skip, so a +broken checkout fails locally rather than in CI). Bypass with +`git push --no-verify`. ## What covers what @@ -106,7 +117,8 @@ enables `.githooks/pre-push` (runs `ddr_sandbox_match`; bypass with |---|---| | Routing | `mmc`, `routing_utils`, `geometry` | | DDR parity (invariant 1) | `ddr_sandbox_match` (in-process, plain `cargo test`) + `compare_ddr_sandbox` example (diagnostics, exits 1 on mismatch) | -| End-to-end acceptance | `juniata_acceptance` (release-only; metric floors + beats-baseline on the committed bundle) | +| End-to-end acceptance | `juniata_acceptance`, `gridded_acceptance` (release-only; metric floors + beats-baseline on the committed bundles) | +| Gridded (DDM30) ingestion | `gridded_bundle` (sub-reach store → subdivided layout, gauge at the cell's last piece, cache hit, dataset opens), `cargo test --lib adjacency::gridded` (synthetic-store validation), `hourly_streamflow::time_major_store_reads_identically_to_divide_major`, `cargo test --lib zarr::tests::upstream_comids_names_each_subdivided_parent_once` | | Sparse / autograd | `sparse_gradcheck`, `sp8_gradcheck` | | KAN head | the 4 `kan_head_*` fixture tests (need `--features fixtures`) | | Leakance | `leakance_gradcheck`, `leakance_off_parity`, `zeta_accum` | diff --git a/.claude/skills/ddrs-dev/references/traps.md b/.claude/skills/ddrs-dev/references/traps.md index 9786b70c..7c154e8c 100644 --- a/.claude/skills/ddrs-dev/references/traps.md +++ b/.claude/skills/ddrs-dev/references/traps.md @@ -18,6 +18,8 @@ is what a future session needs, not the narrative. | Mid-eval `object not found` from icechunk | T8 | | `.ddrs/` appeared somewhere unexpected | T9 | | Resumed run trains zero batches | T10 | +| Baseline predicts ~0 (FHV −100 %) on a new Q′ store, no error | T11 | +| Fresh checkout/worktree fails in `cudarc`'s build script: `Unsupported cuda toolkit version` | T12 | --- @@ -195,6 +197,31 @@ Resume additionally requires `experiment.epochs > E` or the resumed run trains z batches. Stored weights are f16, so a resumed trajectory drifts slowly — expected, not a bug. +## T11 — Time-major `Qr(time, divide_id)` read as fill values (fixed 2026-09-09) + +DDR's gridded (DDM30) `build_gridded_qprime.py` writes `Qr(time, divide_id)`; +the MERIT stores and the contract are `Qr(divide_id, time)`. zarrs answers an +out-of-range `retrieve_array_subset` with **fill values, not an error**, so the +old fixed-axis `read_slab` returned NaN for every cell and the summed-Q′ +baseline scored NSE −0.655 with FHV −100 % — no message anywhere. `open` now +decides the axis order from the array's `dimension_names`, falling back to its +shape against the coordinate lengths, and refuses a shape matching neither +(`detect_time_major`, `src/data/store/icechunk.rs`). Fixture: +`tests/fixtures/qr_daily_time_major.ic`. A baseline that is *deterministic and +matches DDR* is the check that catches this class of bug; that is why both +acceptance tests band the baseline tightly. + +## T12 — CUDA toolkit newer than `cudarc` knows (2026-09-09) + +`cudarc 0.19.7` (via `cuda-version-from-build-system`) parses `nvcc --version` +and panics on a version outside its table (max 13.2). This host moved to CUDA +13.3.1 in 2026-09; the main tree's `target/` still carried the 13.2 build-script +output, so only a fresh worktree/checkout hit it. Fix: `CUDARC_CUDA_VERSION=13020` +(gitignored `.cargo/config.toml` `[env]` block, or export it). 13.3 is +driver/ABI-compatible with 13.2 for the dynamic-loading path ddrs uses. The +main tree will hit the same panic the next time its cudarc build script reruns +(any `CUDA_*` env change or a `cargo clean`). + ## Exit codes `src/cli/types.rs`: 0 Success · 1 Generic · 2 ConfigInvalid · 3 DataSourceMissing · diff --git a/.claude/skills/ddrs-eval-plots/SKILL.md b/.claude/skills/ddrs-eval-plots/SKILL.md index 0699ac6a..09348897 100644 --- a/.claude/skills/ddrs-eval-plots/SKILL.md +++ b/.claude/skills/ddrs-eval-plots/SKILL.md @@ -200,7 +200,11 @@ the first dual win on the dHBV2-UH store. Previous best was 0.7152 / 0.7106 +0.037 NSE and −0.007 KGE). Full table and the 2026-07-30 KGE qualification: `ddrs-dev` → `references/research-status.md`. -`ddrs show | grep -E "nse|kge|loss"` is the cheap pre-plot triage. +`ddrs show ` is the cheap pre-plot triage: since 2026-09-09 it prints the +`metrics` block and, for a `train-and-test` run, the summed-Q' baseline median with +the routed value and delta beside it, so "did it beat the baseline" is answered +without opening a notebook. On a binary older than that, text mode printed no +metrics at all and you need `--json`. See skill `ddrs-run`. ### Have the parameters converged? diff --git a/.claude/skills/ddrs-run/SKILL.md b/.claude/skills/ddrs-run/SKILL.md new file mode 100644 index 00000000..f9979138 --- /dev/null +++ b/.claude/skills/ddrs-run/SKILL.md @@ -0,0 +1,259 @@ +--- +name: ddrs-run +description: Use when launching, watching, resuming, or auditing a ddrs training or evaluation job: starting a train / train-and-test run, switching datasets between MERIT and gridded DDM30, finding a past run's metrics or git SHA, resuming after a crash, pruning old runs, or when a run used the wrong code or the wrong data. Trigger on "train a model", "run ddrs", "which run was that", "resume from epoch N", "did it beat the baseline", "ddrs status/show/gc/sources", "my config change did nothing". +--- + +# Running and tracking ddrs models + +The `ddrs` binary is a terraform-style lifecycle: **`plan`** prepares and validates, +**`run`** executes a workflow into `.ddrs/runs//`, **`show`/`status`** audit what +happened, **`gc`** prunes. Every run writes a self-describing manifest (config +snapshot, input fingerprints, git SHA), so a finished run can be reproduced without +guessing what produced it. + +This skill is the runbook. For building/testing ddrs itself use `ddrs-dev`; for +plotting a finished run use `ddrs-eval-plots`. + +## Preflight: do these three before EVERY run + +Skipping any one of them produces a run that looks fine and is wrong. + +**1. Refresh the binary if `src/` changed since the last install.** + +```bash +cargo install --path . # canonical; updates ~/.cargo/bin/ddrs +``` + +`cargo build` does **not** update the `ddrs` on your PATH. A stale binary does not +say it is stale. It fails as if your *config* were wrong. Reproduced 2026-09-09: +a pre-gridded `~/.cargo/bin/ddrs` given a valid `gridded_network:` config answered + +``` +error: config invalid: data_sources: adjacency sources are missing — +either set both `conus_adjacency` and `gages_adjacency`, or set `geospatial_fabric` +``` + +The key was right there in the file. Alternatives to `cargo install`: run +`cargo run --release --bin ddrs -- `, or use an **absolute** path to +`target/release/ddrs` (a relative one resolves to the main tree from a worktree). + +`cargo install` writes the machine-global `~/.cargo/bin/ddrs`, shared by every +checkout and worktree on the box. Installing from a worktree silently changes what +plain `ddrs` means everywhere until someone installs from a different checkout. When +several branches are in play, prefer the absolute `target/release/ddrs` of the tree +you mean. + +**2. Pass `--workspace` whenever you pass `--config`.** + +The workspace defaults to `.ddrs/` *beside the config*, so +`--config config/experiments/x.yaml` silently creates +`config/experiments/.ddrs/`, a fresh workspace with no run history and a cold +adjacency/baseline cache. Always: + +```bash +ddrs --config config/experiments/x.yaml --workspace /path/to/repo/.ddrs plan +``` + +The one sanctioned exception is `examples/juniata*/ddrs.yaml`, whose workspace is +meant to land beside it. + +**3. `plan` before `run`, and read the baseline it prints.** + +`run` re-plans internally, but running `plan` first is where you see the summed-Q' +baseline, source drift, and the resolved adjacency cache **before** burning hours. +If the trained median NSE does not beat that baseline number, the routing is not +earning its keep. + +## Start from a ready-made config, not from scratch + +`ddrs.yaml` is gitignored and per-workspace; the committed configs are the +starting points. Copy one in rather than hand-building a config: + +| Want | Copy | +|---|---| +| MERIT CONUS | `config/merit_training.yaml` | +| Gridded ISIMIP DDM30 CONUS | `config/experiments/gridded_conus.yaml` | +| Single-basin MERIT sample | `examples/juniata/ddrs.yaml` (run in place) | +| Single-basin gridded sample | `examples/juniata_gridded/ddrs.yaml` (run in place) | + +```bash +cp config/experiments/gridded_conus.yaml ddrs.yaml # gitignored; safe to overwrite +``` + +**`ddrs plan`'s interactive bootstrap does not offer these.** With no `ddrs.yaml` +present it offers exactly two choices: the last successful run's config snapshot, +or the clean template, and the clean template is hardcoded to +`config/merit_training.yaml` (`src/cli/plan_bootstrap.rs`). Bootstrapping and then +expecting a gridded config gets you a MERIT one. Copy the file you want first. + +## Launch a training run + +```bash +# from the repo root, with ./ddrs.yaml as the config +ddrs plan # validate + baseline (cached) +ddrs run --workflow train-and-test # train, eval, compare vs baseline +ddrs run --workflow train # train only, no eval +ddrs run --workflow train --max-mini-batches 2 # mechanics smoke, minutes not hours +ddrs run --workflow train-and-test --backend cpu # CPU (diagnostics; slower) +``` + +Always smoke-test a new config with `--max-mini-batches 2` first. It exercises the +whole data → forward → backward → checkpoint path in minutes. + +`--workflow` must agree with the config's `mode:` (`training` ↔ `train` / +`train-and-test`, `testing` ↔ `eval`); `plan` rejects a contradiction at load. + +## Switch datasets without editing YAML + +Data-source groups are named save files under `config/sources/`. Switching is one +command and re-locks the workspace: + +```bash +ddrs sources list # '*' marks the group matching the current config +ddrs sources use conus-gridded # gridded ISIMIP DDM30 network +ddrs sources use conus # MERIT CONUS +ddrs sources save my-group # snapshot the current data_sources block +``` + +**`ddrs sources` resolves groups at `/config/sources/`.** Run it from +the repo root against `./ddrs.yaml`. Pointing `--config` at +`config/experiments/foo.yaml` makes it look in `config/experiments/config/sources/` +and report *"no groups"*. The groups are fine, the lookup was wrong. + +Switching is textual: the whole `data_sources:` block is replaced, so keys that +belong to only one network (`gridded_network`, `geospatial_fabric`, `aorc_precip`) +appear and disappear correctly. Everything outside the block is untouched, so +`geodataset:` and `params.subdivision` do **not** follow the switch. + +Fix `geodataset:` after switching, or just delete the key. It is a provenance label +(written into managed-build zarr attrs and the run's `config.yaml` snapshot) and +nothing at runtime reads it, but since 2026-09-09 a value contradicting the +adjacency source is a **load error**, and an absent one is **inferred** from the +source. So omitting it is always right, and a stale one now stops the run instead of +silently mislabelling its audit trail. Only the two managed-build sources are +decidable; with explicit `conus_adjacency`/`gages_adjacency` any label is accepted, +because a pre-built store could be either network. + +## Track a run + +**While it runs**, every run tees a timestamped log: + +`run` prints `run output → ` to **stderr as it starts**, before training +begins, so take the path from there rather than guessing. If you already lost it: + +```bash +tail -f "$(ls -dt .ddrs/runs/*/ | head -1)run.log" +``` + +It captures stdout+stderr at fd level (`dup2`), so CUDA driver messages and any +child process land there too, and it is opened `O_APPEND` and flushed per line, so +the log survives a crash. That makes `run.log` the first thing to read after a run +dies. Per-epoch lines carry loss and median `n`; a flat loss across epochs is a real +signal, not noise (see `ddrs-dev` traps). + +**Afterwards:** + +```bash +ddrs status # workspace summary: last run, lock state, disk usage +ddrs show # status, git SHA, adjacency, metrics, baseline delta +ddrs show --json # the full manifest, machine-readable +``` + +`ddrs show` prints a `metrics` block and, for a `train-and-test` run, the +summed-Q' baseline reduced to a median per metric with the routed value and the +delta beside it, flagging the case where routed does not beat baseline: + +``` +metrics + median_nse_finite 0.7511 + median_kge_finite 0.7299 +baseline (summed Q', median over gauges) + NSE 0.5943 routed 0.7511 Δ +0.1568 + KGE 0.6715 routed 0.7299 Δ +0.0583 +``` + +(Before 2026-09-09 text mode printed no metrics at all, so `ddrs show | grep +nse` came back empty on a healthy run. If you are on an older binary and see that, +it is the binary, not the run.) + +`.ddrs/runs//` holds `manifest.json`, `config.yaml` (the exact config that ran), +`run.log`, `checkpoints/`, and `Cargo.lock`. The manifest is the audit record: +`git` (sha, **dirty**, branch), `sources` (blake3 fingerprint per input file), +`source_lock.drift`, `resolved_adjacency`, and `metrics`. + +**`metrics` differs by workflow.** A `train` run reports only +`epochs_completed`, `final_mini_batch`, `phase1_seconds`. There are no NSE/KGE +numbers, because nothing was evaluated. Skill numbers +(`median_nse_finite`, `median_kge_finite`, `n_gauges_total`, `phase2_seconds`) come +from `train-and-test`. Comparing a `train` run's metrics to a baseline is a category +error; re-run as `train-and-test`. + +**If you read the baseline file yourself, reduce it first.** `train-and-test` copies +the baseline to `/baseline/manifest.json`, and its `metrics` are +**per-gauge arrays** (`nse`, `kge`, `rmse`, `bias`, `fhv`, `flv`, one entry per +gauge in `gage_ids` order), not pre-reduced medians. Comparing that raw array to the +run manifest's `median_nse_finite` compares a scalar against a 620-element list. +`ddrs show` does the reduction for you; `--json` does not, and `run` never re-prints +the baseline table that `plan` shows. + +`git.dirty: true` means the working tree had uncommitted edits, so the SHA alone +does not reproduce that run. + +## Resume after a crash + +A checkpoint is a **directory** `.ddrs/runs//checkpoints/epoch_E_mb_M/` holding +`head.mpk`, `optim.mpk`, `state.json` (epoch, next mini-batch, RNG, sampler +permutation). To resume, set in the config: + +```yaml +experiment: + checkpoint: .ddrs/runs//checkpoints/epoch_6_mb_0 + epochs: 20 # MUST exceed the checkpoint's epoch or zero batches train +``` + +Resume restores all three files, so the run draws the same gauge batches and +windows the original would have. Weights are stored f16, so a resumed trajectory +drifts slowly from an uninterrupted one. That is expected, not a bug. + +**A resume produces a NEW run id; it does not continue the old directory.** The run +id is built from the invocation's own start timestamp, so every `ddrs run` mints a +fresh `.ddrs/runs//`. The interrupted run's manifest, log, and checkpoints stay +frozen at the crash, and the resumed run starts an empty log of its own. Expect two +run ids for one training curve, and keep both: the first holds the epochs the second +does not. + +Flat `epoch_E_mb_M.mpk` **files** instead of directories mean a stale pre-resume +binary wrote them. Go back to preflight step 1. + +## Prune + +```bash +ddrs gc --keep 5 --keep-successful --dry-run # always dry-run first +ddrs gc --keep 5 --keep-successful +ddrs gc --older-than 30d +``` + +`gc` prunes only `.ddrs/runs/`. Adjacency and baseline caches are kept, because they are +content-addressed and expensive to rebuild. + +## What does not work + +| You might try | Reality | +|---|---| +| `ddrs run --workflow eval` | Fails: *"standalone --workflow eval needs a --from-run flag"*, and `--from-run` is unimplemented. Use `train-and-test`, or the legacy `eval` binary against a checkpoint. | +| `ddrs init` | Removed. Prints *"ddrs init has been merged into ddrs plan"*. Use `ddrs plan`. | +| `ddrs experiment ` | **Not on master.** The subcommand and `src/experiment/` live on the `experiment-adjoint` branch. On master it exits with *"unrecognized subcommand"*. | +| `ddrs run --epochs 10` | No such flag. `experiment.epochs` is YAML-only; edit the config. | +| `ddrs run --checkpoint ` | No such flag. Resume goes through `experiment.checkpoint:` in the config. | +| A typo'd config key | Now a load error naming the key. Every config section sets `deny_unknown_fields` as of 2026-09-09; before that a typo silently took its default, which was the most common cause of "my config change did nothing". On an older binary, that silence is still the first thing to suspect. | + +## Red flags: stop and re-check preflight + +- You edited `src/` and then typed `ddrs …` without `cargo install --path .` +- You passed `--config` without `--workspace` +- A config error names a key you did not set, or denies a key you did set +- `ddrs sources list` reports "no groups" +- Checkpoints are flat `.mpk` files rather than directories +- You are comparing a `train` run's metrics against a baseline + +Full flag-by-flag reference: `references/commands.md`. diff --git a/.claude/skills/ddrs-run/references/commands.md b/.claude/skills/ddrs-run/references/commands.md new file mode 100644 index 00000000..4c983347 --- /dev/null +++ b/.claude/skills/ddrs-run/references/commands.md @@ -0,0 +1,178 @@ +# `ddrs` command reference + +Every subcommand and flag, captured from `ddrs --help` on master, 2026-09-09. +Global flags `--config` and `--workspace` are accepted by **every** subcommand. + +## Global + +| Flag | Default | Notes | +|---|---|---| +| `--config ` | discover `ddrs.yaml` upward from cwd, stopping at the first `.git` ancestor | | +| `--workspace ` | `.ddrs/` **beside the config** | Pass it explicitly whenever you pass `--config` | + +Config discovery walking up to `.git` means a bare `ddrs plan` from a subdirectory +finds the repo-root `ddrs.yaml`. It does **not** search sibling directories. + +## `ddrs plan` + +Prepare + preview. Idempotent, but **not side-effect free**: on first run it opens +the forcing store to compute the summed-Q' baseline and builds the managed adjacency +caches. Subsequent plans on the same inputs are cache hits and instant. + +| Flag | Notes | +|---|---| +| `--workflow ` | Override the config's `workflow:` for this invocation | +| `--json` | Plan result as JSON | +| `--force` | Re-run the GPU smoke test even if a cached verdict exists | +| `--min-free-gpu-gb ` | Warn below this free GPU memory at probe time (default 8) | + +What it does, in order: workspace skeleton + GPU probe + cached smoke test → +locate/bootstrap `ddrs.yaml` → load and validate the config → fingerprint data +sources and diff against `.ddrs/sources.lock` → resolve adjacency (explicit paths, +or a managed build from `geospatial_fabric` / `gridded_network`) → summed-Q' +baseline. + +Lock semantics: `plan` reports drift, then **refreshes** the lock ("sources as of my +last plan"). Use `run --strict` to abort on drift *before* relocking, preserving the +evidence. + +If `ddrs.yaml` is missing, `plan` bootstraps one interactively, offering the last +successful run's `config.yaml` snapshot or the bundled template. **Non-TTY callers +must pass `--config`**, otherwise it errors with +*"no ddrs.yaml found; pass --config or run interactively"*. + +## `ddrs run` + +Executes a workflow: re-plans internally, then trains and/or evaluates into +`.ddrs/runs//`. + +| Flag | Notes | +|---|---| +| `--workflow ` | Overrides the config's `workflow:` | +| `--plot` | After success, dump per-COMID KAN parameters to `plot/kan_parameters.nc` | +| `--strict` | Exit 4 if sources drifted since the last plan, instead of warning + relocking | +| `--max-mini-batches ` | Stop each epoch after N mini-batches (smoke testing) | +| `--batch-order-from ` | Replay a captured mini-batch order from JSON; overrides the per-epoch shuffle. Schema: `[{"epoch": int, "mb": int, "staids": [str, …]}]` | +| `--backend ` | **Default `cuda`.** `cpu` is NdArray, deterministic, and forces `sparse_solver` to cpu | +| `--json` | Run result as JSON | + +Diagnostics belong on `--backend cpu`: it is deterministic, and the default is +CUDA, so a diagnostic run that omits the flag is not the run you think you made. +Never mix backends within one comparison set. + +Run id format: `-[-]`, where the optional group +segment is the active data-source group when the config's `data_sources` block +matches a saved `config/sources/.yaml` (compared structurally, so comments +and key order do not break the match). + +Consequence of the group lookup path: a run launched with +`--config config/experiments/foo.yaml` **never** gets a group segment, because +`active_group` looks in `config/experiments/config/sources/`, which does not +exist. Launch from the repo root against `./ddrs.yaml` if you want the dataset +named in the run id. + +## `ddrs show ` + +Prints a past run's manifest. `--json` for the raw record. + +Human output: run id, status, workflow, started/finished, `git (dirty)`, +drift, and the resolved adjacency store paths + cache key with hit/built. + +Manifest keys: `run_id`, `ddrs_version`, `git` {sha, dirty, branch}, `workflow`, +`config_path`, `started_at`, `finished_at`, `status`, `exit_reason`, `system`, +`sources` (per input: path, mtime, size, `fp: blake3:…`), `resolved_adjacency`, +`source_lock` {lockfile, matched, drift}, `outputs` (checkpoint list), +`metrics`, `max_mini_batches`. + +## `ddrs status` + +Workspace summary: workspace path, lockfile present/absent, last run id, and disk +usage of `.ddrs/runs/` and `.ddrs/adjacency/`. `--json` supported. Use it to list +run ids for `ddrs show`. + +## `ddrs gc` + +Deletes run directories under `.ddrs/runs/`. + +| Flag | Notes | +|---|---| +| `--keep ` | Keep the N most recent runs | +| `--keep-successful` | Never delete successful runs | +| `--older-than ` | Only delete runs older than e.g. `30d`, `12h` | +| `--dry-run` | List what would be deleted, delete nothing | + +Prints `note: .ddrs/adjacency/ caches are kept (not pruned by gc in v1)`. + +## `ddrs sources` + +Named data-source groups stored at `/config/sources/.yaml`. For +the standard repo-root `ddrs.yaml` that is `config/sources/`. + +| Subcommand | Notes | +|---|---| +| `list` | Lists groups; `*` marks the one matching the current config | +| `use ` | Splice the group into the config's `data_sources:` and refresh `sources.lock` | +| `save ` | Snapshot the current `data_sources` block as a group (`--force` to overwrite) | + +`save`/`use` are textual: comments inside the block travel with the group, +everything outside it is untouched, and `use` validates that the spliced config +parses before committing. + +Groups shipped in-repo: `conus`, `conus-gridded`, `conus-hourly`, +`conus-experimental`, `conus-hydrodl2`, `global`, `daily-lstm`, `hourly-lstm`, +`aorc_dhbv_distributed`. + +## `ddrs import ` + +Validates a Q' store against the DDR store contract +(`docs/nh-qprime-store-contract.md`) and registers it as a source group. + +| Flag | Notes | +|---|---| +| `--name ` | Group to register under `config/sources/` | +| `--dry-run` | Validate + report coverage only | +| `--force` | Overwrite an existing group of that name | + +The reader sniffs daily vs hourly-native resolution from the CF time units. Hourly +stores are sliced natively. An hourly store combined with +`kan_head.disaggregation` is a config error. + +## Workspace layout + +| Path | Written by | Holds | +|---|---|---| +| `ddrs.yaml` | `plan` (via `$EDITOR`) | Workflow + experiment config (gitignored) | +| `.ddrs/system.json` | `plan` | GPU/driver/smoke-test record | +| `.ddrs/sources.lock` | `plan` | Fingerprints of `data_sources` paths | +| `.ddrs/adjacency//` | `plan`, first time | Cached adjacency stores, content-addressed | +| `.ddrs/baselines//` | `plan`, first time | Cached summed-Q' baseline | +| `.ddrs/runs//manifest.json` | `run` | Per-run audit record | +| `.ddrs/runs//config.yaml` | `run` | Snapshot of the config that produced the run | +| `.ddrs/runs//run.log` | `run` | fd-level timestamped tee of stdout+stderr | +| `.ddrs/runs//Cargo.lock` | `run` | Dependency snapshot | +| `.ddrs/runs//checkpoints/epoch_E_mb_M/` | `run`, train phase | `head.mpk`, `optim.mpk`, `state.json` | +| `.ddrs/runs//baseline/` | `run`, train-and-test | `manifest.json` + raw `predictions.f32` / `observations.f32` | +| `.ddrs/runs//eval/predictions.zarr` | `run`, eval phase | Routed predictions over the eval window | + +## Exit codes + +`src/cli/types.rs`: 0 Success · 1 Generic · 2 ConfigInvalid · 3 DataSourceMissing · +4 LockDrift · 5 RuntimeFailure · 6 WorkspaceNotInitialized. + +## Legacy binaries (deprecated, removed in 0.4) + +`train`, `eval`, `train_and_test` still build and each prints a deprecation warning +pointing at the `ddrs run` equivalent. The legacy `eval` binary remains the only way +to evaluate an existing checkpoint without retraining, since +`ddrs run --workflow eval` is unimplemented: + +```bash +cargo build --release --bin eval +target/release/eval --config \ + --checkpoint .ddrs/runs//checkpoints/epoch_E_mb_M \ + --output /tmp/out.zarr +``` + +`--checkpoint` means different things per binary: `eval` and `probe_zeta_gradient` +take the **directory**; `dump_parameters` takes the **head base** +(`epoch_E_mb_M/head`, no `.mpk`). diff --git a/.githooks/pre-push b/.githooks/pre-push index ea658bff..1ca888a8 100755 --- a/.githooks/pre-push +++ b/.githooks/pre-push @@ -1,7 +1,12 @@ #!/usr/bin/env bash -# Fast local gate before any push: the DDR sandbox parity test (invariant 1). +# Fast local gate before any push: +# ddr_sandbox_match — DDR parity (invariant 1) +# gridded_bundle — gridded (DDM30) ingestion contract against the +# committed examples/juniata_gridded/ bundle +# Both are debug tests that never skip, so a broken checkout fails here +# rather than in CI. One cargo invocation, so the build is shared. # Enable per-clone: git config core.hooksPath .githooks # Bypass once: git push --no-verify set -euo pipefail cd "$(git rev-parse --show-toplevel)" -cargo test --test ddr_sandbox_match +cargo test --test ddr_sandbox_match --test gridded_bundle diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml index 85713b6d..b2bed1a1 100644 --- a/.github/workflows/ci.yml +++ b/.github/workflows/ci.yml @@ -139,3 +139,8 @@ jobs: run: mkdir -p output && cargo run --release --example compare_ddr_sandbox - name: Juniata end-to-end acceptance run: cargo test --release --test juniata_acceptance -- --nocapture + # Same chain on the gridded (ISIMIP DDM30) network: both acceptance + # tests self-skip under debug_assertions, so the debug `test` job never + # runs them. Shares this job's release build; the test itself is ~10 s. + - name: Gridded (DDM30) end-to-end acceptance + run: cargo test --release --test gridded_acceptance -- --nocapture diff --git a/.gitignore b/.gitignore index a1d4fbb2..56620905 100644 --- a/.gitignore +++ b/.gitignore @@ -62,6 +62,8 @@ ddrs.yaml # …but the Juniata sample's config is a committed example, not a personal # workspace config (its .ddrs/ workspace is still ignored by the rule above). !examples/juniata/ddrs.yaml +# …same exception for the gridded Juniata example's committed config. +!examples/juniata_gridded/ddrs.yaml # Personal research journal (local only) JOURNAL.md diff --git a/CLAUDE.md b/CLAUDE.md index 5079a0e8..8c0ab35b 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -72,6 +72,8 @@ cargo test --test ddr_sandbox_match # assert routed NSE/KGE floors + beats the summed-Q' baseline. Release-only # (silently skips in debug builds — a debug train takes minutes): cargo test --release --test juniata_acceptance +# Same, on the gridded (ISIMIP DDM30) Juniata bundle (examples/juniata_gridded/): +cargo test --release --test gridded_acceptance # Hydrograph plot for visual sanity-check: cargo run --release --example benchmark_hydrograph @@ -85,12 +87,14 @@ cd ~/projects/ddr && uv run python ~/projects/ddrs/scripts/export_ddr_sandbox.py `.github/workflows/ci.yml` runs on every PR and every master push: job `test` (debug `cargo test --features fixtures`) and job `acceptance` -(release `compare_ddr_sandbox` + `juniata_acceptance`). Branch protection +(release `compare_ddr_sandbox` + `juniata_acceptance` + +`gridded_acceptance`). Branch protection requires both for PR merges; direct pushes stay possible and are vetted by the push-triggered run. CI has no GPU and no real data: green CI does NOT cover `--features cuda` tests or the data-dependent suite; run those locally per the tier gates. Local pre-push hook (opt-in): -`git config core.hooksPath .githooks`. +`git config core.hooksPath .githooks` — runs `ddr_sandbox_match` + +`gridded_bundle` (both never skip; ~0.4 s warm). ### `ddrs` CLI (preferred entrypoint) @@ -322,6 +326,27 @@ math. Read it before touching `src/routing/` or `src/sparse/`. | Global streamflow forcing | `/gpfs/hjj5218/data/dmc_forcing/streamflow/zarr/8km/merit_global_v2.7` (multi-zone zarr v2, `(time, COMID)` f64 per pfaf-2 zone) | `zarrs` | | Global gage metadata | `/gpfs/hjj5218/data/dmc_forcing/gage_information/formatted_gage_csvs/v3.1/8km/` (57 per-zone `_all.csv`) | `csv` | +**Gridded (ISIMIP DDM30, 0.5°) routing** (2026-09-09, mirrors DDR PR #194): +ddrs has no `engine/` module and ports none of DDR's raster aggregation or Q′ +regridding — it consumes DDR's products as a data source. `data_sources. +gridded_network` names DDR's **sub-reach** adjacency zarr (`order`, +`parent_cell`, `indices_0/1`, `length_m`, `slope`); `ddrs plan` relabels it into +the existing subdivided store layout (parent = cell, pieces = sub-reaches, gauge +read at the cell's last piece) and cuts gauge subgraphs from it +(`src/adjacency/gridded.rs`, `cache::resolve_or_build_gridded`). Attributes +(`COMID` = cell id) and Q′ (`divide_id` = cell id, split across pieces by +`pieces_per_row_divisor`) need no changes; the gauge CSV's `cell` column is an +alias of `COMID`. It is exclusive with `geospatial_fabric` and the explicit +zarr pair, and `params.subdivision.enabled` is rejected with it. Committed +sample: `examples/juniata_gridded/` (routed NSE 0.751 / KGE 0.730 vs baseline +0.594 / 0.672, which matches DDR's baseline to rounding); CONUS: +`config/sources/conus-gridded.yaml` (620 gauges via +`scripts/snap_gridded_gauges.py`). Known deviations from DDR (per-cell KAN head, +no `da_ratio` output correction) are in +`docs/superpowers/specs/2026-09-08-ddrs-gridded-routing-design.md` §7. +DDR's gridded Q′ stores are `Qr(time, divide_id)`; the icechunk reader sniffs +the axis order (`detect_time_major`) — see the ddrs-dev skill's trap T11. + The `observations` and `streamflow` data sources auto-detect their format (`ObservationsStore::open` / `StreamflowSource::open`, `src/data/store/mod.rs`); the icechunk CONUS path is the fallback when the diff --git a/config/experiments/gridded_conus.yaml b/config/experiments/gridded_conus.yaml new file mode 100644 index 00000000..5f70179c --- /dev/null +++ b/config/experiments/gridded_conus.yaml @@ -0,0 +1,94 @@ +# Gridded (ISIMIP DDM30, 0.5°) CONUS training — the ddrs counterpart of DDR's +# examples/juniata_gridded/train_conus.py (DeepGroundwater/ddr PR #194). +# +# Data: config/sources/conus-gridded.yaml (inlined below). 620 gages_3000 +# gauges above 2,000 km² snapped to cells (scripts/snap_gridded_gauges.py), +# routed over DDR's 91,867-sub-reach Courant≈1 network. +# +# Run from the repo root — pass --workspace, or .ddrs lands beside this file: +# +# ddrs --config config/experiments/gridded_conus.yaml --workspace .ddrs plan +# ddrs --config config/experiments/gridded_conus.yaml --workspace .ddrs run --workflow train-and-test +# +# Deviations from train_conus.py (see the 2026-09-08 gridded design spec §7): +# the KAN runs per cell (not per sub-reach), and predictions are NOT divided +# by the gauge's drainage-area ratio (`da_ratio` in the gauge CSV is +# informational). DDR's reference: 45 epochs, median NSE 0.542 vs summed 0.390. + +mode: training +workflow: train-and-test +geodataset: ddm30 +device: 0 +seed: 42 +np_seed: 42 + +data_sources: + gridded_network: /home/tbindas/projects/ddr/data/ddm30/ddm30_subreach_adjacency.zarr + attributes: /mnt/ssd1/data/icechunk/ddm30_conus_attributes.nc + streamflow: /mnt/ssd1/data/icechunk/ddm30_conus_uh_retrospective_regridded.ic + observations: /mnt/ssd1/data/icechunk/usgs_daily_observations + gages: /home/tbindas/projects/ddr/references/gage_info/ddm30_conus_gauges.csv + +experiment: + # train_conus.py defaults: batch 16 gauges, rho 90, warmup 5, lr 1e-3 + # halved at epoch max(2, epochs // 2). + batch_size: 16 + start_time: 1981/10/01 + end_time: 1995/09/30 + epochs: 10 + rho: 90 + shuffle: true + warmup: 5 + learning_rate: + 1: 0.001 + 5: 0.0005 + grad_clip_max_norm: 1.0 + +kan_head: + hidden_size: 21 + num_hidden_layers: 2 + grid: 50 + k: 2 + input_var_names: + - SoilGrids1km_clay + - aridity + - meanelevation + - meanP + - NDVI + - meanslope + - log10_uparea + - SoilGrids1km_sand + - ETPOT_Hargr + - Porosity + learnable_parameters: + - n + - q_spatial + - p_spatial + +params: + parameter_ranges: + n: [0.02, 0.2] + q_spatial: [0.0, 1.0] + p_spatial: [1.0, 200.0] + attribute_minimums: + discharge: 1.0e-4 + slope: 1.0e-3 + velocity: 0.01 + depth: 0.01 + bottom_width: 0.01 + defaults: + p_spatial: 21.0 + log_space_parameters: + - p_spatial + tau: 9 + subdivision: + enabled: false + sparse_solver: cuda + use_cuda_graphs: false + +# train_conus.py's default test window (two water years). +testing: + start_time: 1995/10/01 + end_time: 1997/09/30 + batch_size: 15 + rho: null diff --git a/config/sources/conus-gridded.yaml b/config/sources/conus-gridded.yaml new file mode 100644 index 00000000..3a7060fd --- /dev/null +++ b/config/sources/conus-gridded.yaml @@ -0,0 +1,17 @@ +data_sources: + # Gridded (ISIMIP DDM30, 0.5°) CONUS routing on DDR PR #194's products. + # The sub-reach adjacency is DDR's `scripts/build_subdivided_adjacency.py` + # output (91,867 Courant≈1 sub-reaches over 5,198 cells); ddrs relabels it + # into its subdivided store layout on the first `ddrs plan`. + gridded_network: /home/tbindas/projects/ddr/data/ddm30/ddm30_subreach_adjacency.zarr + # 25 cell attributes, `COMID` = cell id (statistics sidecar: + # scripts/build_attribute_statistics.py writes it beside the file). + attributes: /mnt/ssd1/data/icechunk/ddm30_conus_attributes.nc + # dHBV2 UH retrospective Q' area-weighted onto cells (`Qr(divide_id, time)`, + # 4,085 of 5,198 cells; the rest read as the 0.001 m³/s floor). + streamflow: /mnt/ssd1/data/icechunk/ddm30_conus_uh_retrospective_regridded.ic + observations: /mnt/ssd1/data/icechunk/usgs_daily_observations + # gages_3000 gauges above 2,000 km² snapped to cells by drainage-area match + # (scripts/snap_gridded_gauges.py; `COMID` = snapped cell, 620 gauges as in + # DDR's train_conus.py). + gages: /home/tbindas/projects/ddr/references/gage_info/ddm30_conus_gauges.csv diff --git a/docs/nh-qprime-store-contract.md b/docs/nh-qprime-store-contract.md index abdf40eb..58dea1e5 100644 --- a/docs/nh-qprime-store-contract.md +++ b/docs/nh-qprime-store-contract.md @@ -15,7 +15,14 @@ lives here. - An **icechunk repository** (`main` branch, local filesystem), root group. - One data variable **`Qr(divide_id, time)`**, dtype **float32**, attr - `units: m^3/s`. The dtype is not checked at `open` — `open` only + `units: m^3/s`. The transposed layout **`Qr(time, divide_id)`** — what + DDR's gridded (DDM30) `build_gridded_qprime.py` writes — is also read: + `open` decides the axis order from the array's dimension names, falling + back to matching its shape against the coordinate lengths, and refuses a + shape that fits neither (`detect_time_major`). The check exists because + zarrs answers an out-of-range subset with fill values rather than an + error, so a wrong axis guess would read NaN silently (2026-09-08). + The dtype is not checked at `open` — `open` only resolves the array handle (`src/data/store/icechunk.rs:239-241`). It is enforced at the first *read*, where `retrieve_array_subset::>` errors on a non-f32 array. `ddrs import` triggers that read itself via diff --git a/docs/superpowers/specs/2026-09-08-ddrs-gridded-routing-design.md b/docs/superpowers/specs/2026-09-08-ddrs-gridded-routing-design.md new file mode 100644 index 00000000..cb763977 --- /dev/null +++ b/docs/superpowers/specs/2026-09-08-ddrs-gridded-routing-design.md @@ -0,0 +1,240 @@ +# Gridded (ISIMIP DDM30) routing in ddrs — design + +**Date:** 2026-09-08 +**Source:** DeepGroundwater/ddr PR #194 ("Gridded (ISIMIP DDM30) routing: +network, attributes, forcing, and training", merged to DDR master). +**Status:** designed autonomously from the user's directive; assumptions are +listed in §7 for review. + +## 1. What DDR #194 actually is + +DDR #194 is a **data-preparation** feature, not a solver change. Its PR body +says so: "The routing engine, `mmc.py`, `merit.py` and the MERIT/Lynker paths +are untouched. Everything here is additive." It adds: + +- `ddr_engine.gridded` — builds a river network from the ISIMIP DDM30 0.5° + D8 grid, aggregates 25 attributes from native rasters onto cells, regrids + per-catchment Q′ onto cells (mass-conserving area weights), splits cells + into Courant≈1 sub-reaches, and validates the products. +- Data products keyed on the **cell id** (`row·720 + col`, row 0 at 55.75°S): + an attributes NetCDF in the `merit_global_attributes_v2` layout (`COMID` = + cell id), a Q′ icechunk store (`Qr(divide_id, time)`, `divide_id` = cell id), + and a **sub-reach adjacency zarr** carrying `order` (node id = + `cell·1000 + sub`), `parent_cell`, `indices_0/1`, `values`, `length_m`, + `slope`, `lat`, `lon`. +- Two trainers (`examples/juniata_gridded/train_gridded.py`, `train_conus.py`) + that feed those products to the unchanged `dmc` router, plus a 1.6 MB + self-contained Juniata bundle. + +ddrs has no `engine/` module and the user asked to keep it that way: **ddrs +consumes DDR's gridded products as another data source; it does not port the +raster/regridding engine.** + +## 2. The key observation: ddrs already has the abstraction + +ddrs's reach subdivision (`params.subdivision`, `.claude/REACH-SUBDIVISION.md`) +already models exactly the structure the gridded sub-reach network has: + +| Gridded sub-reach zarr (DDR) | ddrs subdivided store | +|---|---| +| `parent_cell` per row (contiguous, upstream-first) | `order` (parent COMID repeated per piece), `parent_order`, `parent_offset` | +| attributes looked up by parent cell | `AttributesStore::open(…, &conus.parent_order)` | +| Q′ per cell ÷ `k_per_node` | `pieces_per_row_divisor` in `MuskingumCunge::setup_inputs` | +| gauge outlet = cell's last sub-reach | `position_lookup` = COMID → **last** row; `outlet(p) = parent_offset[p+1]−1` | +| `length_m`, `slope` in the adjacency zarr | `length_m`, `slope` in the adjacency zarr | +| summed-Q′ baseline sums each cell once | `upstream_comids` dedups via `parent_order` | + +So the whole routing/training/eval/baseline path is untouched. What is new is +**one builder input**: a reader that turns the DDM30 sub-reach zarr into a +`SubdividedAdjacency`, and the plumbing that lets `ddrs plan` build the managed +stores from it. + +``` + data_sources.gridded_network data_sources.gages (CSV, `cell` col) + (DDM30 sub-reach zarr from DDR) │ + │ │ + ▼ │ + adjacency::gridded::GriddedNetwork::open │ + order/parent_cell/indices/length_m/slope │ + validate: lower-tri, contiguous parents │ + │ │ + ▼ to_subdivided() │ + SubdividedAdjacency { order = parent_cell, │ + parent_order, parent_offset, rows, cols, │ + length_m, slope } │ + │ │ + ├──► write_conus_store_subdivided ──► _adjacency.zarr + │ │ + └──► build_gauge_subgraphs(expanded_view, csv) ◄──┘ + │ (BFS upstream from the cell's OUTLET piece) + ▼ + write_gauges_store ──► _gages_adjacency.zarr + + .ddrs/adjacency// … then MeritGagesDataset, baseline, train, eval + run exactly as they do for a subdivided MERIT network. +``` + +## 3. Changes (blast radius) + +| Area | Change | Tier | +|---|---|---| +| `src/config.rs` | `data_sources.gridded_network: Option`; guards (§4) | C | +| `src/data/store/gage_csv.rs` | `COMID` column accepts `cell` as a serde alias | C | +| `src/adjacency/gridded.rs` (new) | `GriddedNetwork::open`, validation, `to_subdivided` | C | +| `src/adjacency/cache.rs` | `resolve_or_build_gridded`; zarr-dir fingerprint; manifest `network_kind` | C | +| `src/adjacency/zarr_write.rs` | writers take a `geodataset` attr value (`"ddm30"`) | C | +| `src/cli/plan.rs`, `src/cli/sources.rs` | third `resolve_adjacency` branch; lock `gridded_network` | C | +| `examples/juniata_gridded/` | DDR's bundle byte-identical + `ddrs.yaml` + statistics JSON + README | — | +| `scripts/build_attribute_statistics.py` | writes the statistics sidecar for any attributes NetCDF (DDR venv) | — | +| `scripts/snap_gridded_gauges.py` | CONUS gauge→cell snapping (DDR venv) → gauge CSV with `COMID` = cell | — | +| `config/sources/conus-gridded.yaml` | source group for the on-disk CONUS gridded stores | D | +| `tests/gridded_bundle.rs`, `tests/gridded_acceptance.rs` | contract + end-to-end floors | — | + +Not touched: `src/routing/`, `src/sparse/`, `src/nn/`, `src/training/`, +`src/data/dataset.rs`, `src/baseline/`. Invariants 1–7 are unaffected; the +Tier A gates still run because `cargo test` includes `ddr_sandbox_match`. + +## 4. Config contract + +```yaml +data_sources: + gridded_network: examples/juniata_gridded/data/juniata_subreach_adjacency.zarr + attributes: examples/juniata_gridded/data/ddm30_conus_attributes.nc + streamflow: examples/juniata_gridded/data/juniata_qprime.ic + observations: examples/juniata_gridded/data/juniata_obs.ic + gages: examples/juniata_gridded/data/juniata_gage.csv +``` + +Load-time guards (all in `validate_data_sources` / `validate_subdivision`): + +- `gridded_network` is a **third mutually-exclusive adjacency source**: + with `geospatial_fabric` → error; with `conus_adjacency`/`gages_adjacency` + → error; with neither of the other two → managed build. +- `geospatial_fabric_layer` with `gridded_network` → error (existing rule: + the layer key needs a `.gpkg` fabric). +- `params.subdivision.enabled: true` with `gridded_network` → error. The + gridded graph is already subdivided upstream (DDR sizes sub-reaches from + celerity at mean flow); ddrs's MERIT subdivision does not apply to it and + the flag would be silently inert. + +The gauge CSV column that names the outlet is `COMID` **or** `cell` (DDR's +gridded CSVs write `cell`). Its value is the DDM30 cell id; the gauge is read +at that cell's most downstream sub-reach. + +## 5. Cache key and manifest + +`key = blake3(gridded_fp ∥ gages_fp ∥ "gridded" ∥ BUILDER_VERSION)[..16]` +where `gridded_fp` hashes every file under the zarr store directory in sorted +relative-path order (name ∥ NUL ∥ bytes). `params.subdivision` does not enter +the key because it is rejected with `gridded_network`. The manifest reuses +`CacheManifest` with `fabric_path`/`fabric_resolved_path` = the store dir and a +new `network_kind: "ddm30"` field (`#[serde(default)]`, so MERIT manifests +stay readable). Store names use the zarr stem +(`juniata_subreach_adjacency_adjacency.zarr` is ugly; the stem's trailing +`_adjacency` is stripped first → `juniata_subreach_adjacency.zarr` / +`juniata_subreach_gages_adjacency.zarr`). + +## 6. Example bundle and tests + +`examples/juniata_gridded/data/` is DDR's bundle **byte-identical** (same +readers on both sides is the point, as with `examples/juniata`), plus +`statistics/merit_attribute_statistics_ddm30_conus_attributes.nc.json` +generated by `scripts/build_attribute_statistics.py` (same six statistics DDR's +`set_statistics` writes; `merit_` prefix kept because ddrs derives the sidecar +name from a fixed template). `ddrs.yaml` mirrors `train_gridded.py`: +`n ∈ [0.02, 0.2]`, `q ∈ [0, 1]`, `p ∈ [1, 200]` log-space, `p` default 21, +lr 1e-3 constant, 30 epochs, rho 90, warmup 5, seed 42, CPU solver. + +Tests: + +- `tests/gridded_bundle.rs` (never skips): bundle opens through + `GriddedNetwork`; 27 rows / 26 edges / 4 parents with pieces 9,6,6,6; the + managed build into a tempdir yields a gauges store whose `01567000` subgraph + covers all 27 rows with `gage_idx` = the last row of cell 138445; the + statistics cover the ten KAN inputs; `MeritGagesDataset::open` succeeds. +- `tests/gridded_acceptance.rs` (release-only, like `juniata_acceptance`): + full `train-and-test`, summed-Q′ baseline in a tight band around DDR's + 0.594 NSE (deterministic, a cross-implementation reader check), routed NSE/KGE + above floors set from a measured ddrs run, routed beats baseline. +- Unit tests in `gridded.rs` (synthetic store), `config.rs` (guards), + `gage_csv.rs` (`cell` alias), `cache.rs` (key sensitivity). + +## 7. Assumptions and known deviations from DDR + +1. **KAN head runs in parent (cell) space.** ddrs's subdivision design feeds + the head `(n_parent, F)` attributes and gathers its outputs onto the pieces. + DDR's gridded trainer runs the head per sub-reach with `log10_uparea` + recomputed per piece by accumulating `cell_area/k`. In ddrs every piece of + a cell shares the cell's parameters and the cell's `log10_uparea`. + *Why accepted:* mirroring it would change `training/forward.rs` and the + attribute path for every network (Tier C on the training core) to add + within-cell variation that DDR's own README calls minor ("with only four + cells the KAN has four distinct attribute rows"). Revisit only if a CONUS + comparison shows a material gap. +2. **No `da_ratio` output correction.** `train_conus.py` divides predictions + by the cell/gauge drainage-area ratio; `train_gridded.py` (the Juniata + reference) does not. ddrs's `FLOW_SCALE` scales the gauge reach's own + lateral inflow, a different quantity, so it is *not* used to emulate it. + The CONUS gauge CSV carries `da_ratio` as an informational column only. +3. **Normalisation statistics are CONUS-wide** (all 5,198 cells), matching + DDR's `conus.describe()`; the script uses `nanstd` (ddof 0) like DDR's + MERIT sidecar rather than pandas' ddof 1 — a 1e-4 relative difference. +4. **Cells with no Q′** (21 % of CONUS sub-reach parents) get the icechunk + reader's 0.001 m³/s fill, DDR's gridded trainer uses 0. Irrelevant for the + Juniata bundle (all four cells present). +5. **The sub-reach zarr is trusted as topologically ordered and + lower-triangular**; `GriddedNetwork::open` verifies `rows > cols` and + parent contiguity and refuses otherwise, rather than re-sorting. +6. **The `gridded_network` key is the dispatch**, not `geodataset:`. The + example config sets `geodataset: ddm30` for documentation; the field is + still unread (as it is today for `merit`). + +## 8. Concerns + +- **Routed skill will not reproduce DDR's 0.795/0.737 exactly**: different + RNG windows (as with `examples/juniata`) plus deviation 1. The acceptance + floors are set from ddrs's own measured run, and the baseline band is the + DDR cross-check. +- **Only the sub-reach form is supported.** DDR also ships an unsplit cell + adjacency (`ddm30_adjacency.zarr`, no `parent_cell`). It is out of scope: + DDR's trainers use the sub-reach store, and its README calls the cell + network "available for cheap iteration" only. A store without + `parent_cell` is rejected with a message saying so. +- **`ddrs sources use conus-gridded` needs a gauge CSV that does not exist in + DDR's references yet**; `scripts/snap_gridded_gauges.py` produces it under + DDR's venv, following the `build_gages_2000_area_balanced.py` precedent. + +## 9. Outcome (2026-09-09) + +Implemented on branch `worktree-gridded-routing` as designed, with two +findings the design did not anticipate: + +1. **`GageSubgraph::upstream_comids` was per-piece on subdivided stores.** It + deduplicated on row position, so a cell with `k` sub-reaches appeared `k` + times and the summed-Q′ baseline counted its Q′ `k` times. Fixed to + deduplicate on the id (a no-op for un-subdivided stores). The same bug + applied to MERIT `params.subdivision` baselines. +2. **DDR's gridded Q′ stores are `Qr(time, divide_id)`**, the transpose of the + contract, and zarrs returns fill values for an out-of-range subset rather + than an error — the fixed-axis reader produced an all-NaN baseline with no + message. `StreamflowStore::open` now detects the axis order from + `dimension_names` (shape fallback, refuses neither) and `read_slab` addresses + the raw buffer accordingly; `tests/fixtures/qr_daily_time_major.ic` pins it. + +Results on the Juniata bundle (CPU, 30 epochs): + +| | NSE | KGE | +|---|---|---| +| summed-Q′ baseline, ddrs / DDR | 0.594 / 0.594 | 0.672 / 0.671 | +| routed, seed 42 | 0.751 | 0.730 | +| routed, seeds {42, 7, 123} | 0.744–0.758 | 0.728–0.732 | +| routed, DDR Python | 0.795 | 0.737 | + +CONUS (`config/experiments/gridded_conus.yaml`, 620 gauges, test window +1995-10-01..1997-09-30): `plan` builds 5,198 cells → 91,867 sub-reaches in +0.6 s; baseline median NSE 0.354 / KGE 0.505. DDR reports 0.390 for its summed +baseline, which includes the `da_ratio` output correction ddrs omits +(deviation 2 in §7). Training beyond a mechanics smoke was not run. + +Environment note: the worktree build needed `CUDARC_CUDA_VERSION=13020` +(host CUDA 13.3.1 vs cudarc 0.19.7's table); recorded as trap T12. diff --git a/docs/usage/inputs-reading.md b/docs/usage/inputs-reading.md index d176bb43..9d236ff9 100644 --- a/docs/usage/inputs-reading.md +++ b/docs/usage/inputs-reading.md @@ -73,7 +73,19 @@ Adjacency is an either/or, enforced at config load by `.ddrs/adjacency//{conus,gages}`. Subsequent plans are cache hits. - **Pre-built zarr** — set *both* `conus_adjacency` and `gages_adjacency` and omit `geospatial_fabric`. -- Exactly one adjacency key, or neither key and no fabric, is a load-time +- **Gridded (ISIMIP DDM30) network** — set `gridded_network` to DDR's + sub-reach adjacency zarr (DeepGroundwater/ddr PR #194, + `scripts/build_subdivided_adjacency.py`: `order`, `parent_cell`, + `indices_0/1`, `length_m`, `slope`) and omit the other two. `ddrs plan` + relabels it into the subdivided store layout below — parent = grid cell, + pieces = the cell's contiguous sub-reaches, gauge read at the last piece — + and cuts per-gauge subgraphs from it + (`adjacency::gridded::GriddedNetwork`, `cache::resolve_or_build_gridded`). + Attributes and Q′ stay keyed on the cell id, so nothing downstream changes; + the gauge CSV may name the outlet under `cell` (DDR's column) instead of + `COMID`. See `examples/juniata_gridded/` and the 2026-09-08 design spec. +- Exactly one adjacency key, neither key and no fabric or gridded network, + or `gridded_network` combined with either of the others, is a load-time error. The cache key is `blake3(fabric_fingerprint ∥ gages_fingerprint ∥ diff --git a/examples/juniata_gridded/README.md b/examples/juniata_gridded/README.md new file mode 100644 index 00000000..3f421ad7 --- /dev/null +++ b/examples/juniata_gridded/README.md @@ -0,0 +1,87 @@ +# Gridded (ISIMIP DDM30) Juniata sample + +The gridded counterpart of [`examples/juniata`](../juniata/README.md), and +the ddrs mirror of DDR's `examples/juniata_gridded` (DeepGroundwater/ddr +PR #194). Same basin and gauge (USGS 01567000, 8,657 km²), but the river +network is the ISIMIP DDM30 0.5° drainage-direction grid instead of MERIT +flowlines: four grid cells, split by DDR into 27 sub-reaches of about 7 km. + +Everything needed is in `data/` (1.6 MB, byte-identical to DDR's bundle plus +the statistics sidecar). No HPC, S3, CUDA, or external stores. + +## Quickstart + +From the **repo root** (the config's data paths are repo-root-relative): + + cargo build --release --bin ddrs + target/release/ddrs --config examples/juniata_gridded/ddrs.yaml plan + target/release/ddrs --config examples/juniata_gridded/ddrs.yaml run --workflow train-and-test --backend cpu + +`plan` relabels DDR's sub-reach store into ddrs's subdivided adjacency +layout, cuts the gauge's upstream subgraph, caches both under +`examples/juniata_gridded/.ddrs/adjacency//` (gitignored), and prints +the summed-Q' baseline. `run` trains the KAN head for 30 optimizer steps +(one random 90-day window per epoch on a single gauge), then evaluates over +water years 1996–2010. The whole workflow takes well under a minute on CPU. + +Reference results (CPU, 30 epochs, test 1995-10-01 – 2010-09-30): + +| | NSE | KGE | +|---|---|---| +| ddrs routed (seed 42) | 0.751 | 0.730 | +| ddrs routed, seeds {42, 7, 123} | 0.744–0.758 | 0.728–0.732 | +| DDR (Python) routed | 0.795 | 0.737 | +| summed-Q' baseline (ddrs) | 0.594 | 0.672 | +| summed-Q' baseline (DDR) | 0.594 | 0.671 | + +The baseline is deterministic and matches DDR to rounding: the two +implementations read the same bundle through independent readers, so it is +a live cross-implementation check (it is what exposed that DDR's gridded Q' +store is written `Qr(time, divide_id)`, the transpose of the MERIT layout; +ddrs now detects both). The routed gap to DDR has two known sources: the 30 +random training windows come from different RNG streams (as in +`examples/juniata`), and ddrs's KAN head runs **per cell** — every +sub-reach of a cell shares the cell's parameters and the cell's +`log10_uparea` — where DDR's gridded trainer runs it per sub-reach with an +accumulated `log10_uparea`. With four cells the field is nearly uniform +either way; DDR's own README calls the sample a demonstration that the +gridded network trains and routes, not a converged model. + +## How the gridded inputs map onto ddrs + +ddrs has no `engine/` module and ports none of DDR's raster aggregation or +Q' regridding. It consumes DDR's products as a data source: the new +`data_sources.gridded_network` key points at DDR's sub-reach adjacency +zarr, and everything downstream is the existing subdivided-network path: + +| DDR gridded product | ddrs reads it as | +|---|---| +| `juniata_subreach_adjacency.zarr` — `order` (node id = cell·1000 + sub), `parent_cell`, `indices_0/1`, `length_m`, `slope` | a subdivided CONUS store: parent = cell, pieces = the cell's contiguous sub-reaches, outlet = last piece | +| `ddm30_conus_attributes.nc` — 25 attributes, `COMID` = cell id (all 5,198 CONUS cells, so normalisation matches a CONUS run) | the attributes table, looked up per cell; sidecar `data/statistics/merit_attribute_statistics_ddm30_conus_attributes.nc.json` from `scripts/build_attribute_statistics.py` | +| `juniata_qprime.ic` — `Qr(time, divide_id)` daily m³/s per cell, 1980–2010 | daily Q' per cell, split evenly across the cell's sub-reaches at routing time (`pieces_per_row_divisor`) | +| `juniata_obs.ic` — `streamflow(gage_id, time)`, USGS 01567000 | observations, unchanged | +| `juniata_gage.csv` — one row with the snapped `cell` and `da_ratio` | the gauge table; `cell` is read as the outlet id, routed at that cell's last sub-reach. `da_ratio` is informational (no output correction, as in DDR's `train_gridded.py`) | + +Design and the full list of deviations: +`docs/superpowers/specs/2026-09-08-ddrs-gridded-routing-design.md`. + +## Tests + + cargo test --test gridded_bundle # bundle contract, never skips + cargo test --release --test gridded_acceptance # end-to-end metric floors + +## CONUS + +`config/sources/conus-gridded.yaml` names the on-disk CONUS products +(91,867 sub-reaches, 620 gauges from `scripts/snap_gridded_gauges.py`), and +`config/experiments/gridded_conus.yaml` mirrors DDR's `train_conus.py` +defaults. Both need the workstation stores under `/mnt/ssd1` and +`~/projects/ddr/data/ddm30`. + +## Regenerating the bundle + +Maintainer only. Regenerate in the ddr repo +(`examples/juniata_gridded/extract_bundle.py`), copy `data/` here, then + + ~/projects/ddr/.venv/bin/python scripts/build_attribute_statistics.py \ + examples/juniata_gridded/data/ddm30_conus_attributes.nc diff --git a/examples/juniata_gridded/data/ddm30_conus_attributes.nc b/examples/juniata_gridded/data/ddm30_conus_attributes.nc new file mode 100644 index 00000000..b0097706 Binary files /dev/null and b/examples/juniata_gridded/data/ddm30_conus_attributes.nc differ diff --git a/examples/juniata_gridded/data/juniata_gage.csv b/examples/juniata_gridded/data/juniata_gage.csv new file mode 100644 index 00000000..fac7774b --- /dev/null +++ b/examples/juniata_gridded/data/juniata_gage.csv @@ -0,0 +1,2 @@ +STAID,STANAME,DRAIN_SQKM,LAT_GAGE,LNG_GAGE,cell,da_ratio +01567000,"JUNIATA RIVER AT NEWPORT, PA",8657.0,40.4785,-77.1294,138445,1.0880212544761465 diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/0JAYFEV4XTV7Z0NMHZFG b/examples/juniata_gridded/data/juniata_obs.ic/chunks/0JAYFEV4XTV7Z0NMHZFG new file mode 100644 index 00000000..ef154772 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/0JAYFEV4XTV7Z0NMHZFG differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/20P4Z9TTV5JCVRTP5NFG b/examples/juniata_gridded/data/juniata_obs.ic/chunks/20P4Z9TTV5JCVRTP5NFG new file mode 100644 index 00000000..83a10b4b Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/20P4Z9TTV5JCVRTP5NFG differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/2MGHGS8R7AN9RNHYFAS0 b/examples/juniata_gridded/data/juniata_obs.ic/chunks/2MGHGS8R7AN9RNHYFAS0 new file mode 100644 index 00000000..02fd74e9 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/2MGHGS8R7AN9RNHYFAS0 differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/4R696W5CMNH287T1N560 b/examples/juniata_gridded/data/juniata_obs.ic/chunks/4R696W5CMNH287T1N560 new file mode 100644 index 00000000..8f16acb0 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/4R696W5CMNH287T1N560 differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/CBTKMXW2S2ZJ22F808DG b/examples/juniata_gridded/data/juniata_obs.ic/chunks/CBTKMXW2S2ZJ22F808DG new file mode 100644 index 00000000..4891a774 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/CBTKMXW2S2ZJ22F808DG differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/D8V36GXME8HHBW236MZG b/examples/juniata_gridded/data/juniata_obs.ic/chunks/D8V36GXME8HHBW236MZG new file mode 100644 index 00000000..26d7b139 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/D8V36GXME8HHBW236MZG differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/DGM5YQ1HNSHJFWTNX7R0 b/examples/juniata_gridded/data/juniata_obs.ic/chunks/DGM5YQ1HNSHJFWTNX7R0 new file mode 100644 index 00000000..8418a901 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/DGM5YQ1HNSHJFWTNX7R0 differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/E6X05RWV1MZRX90AHJVG b/examples/juniata_gridded/data/juniata_obs.ic/chunks/E6X05RWV1MZRX90AHJVG new file mode 100644 index 00000000..9494796f Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/E6X05RWV1MZRX90AHJVG differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/PHM8GS06CPVYYQ5XEQTG b/examples/juniata_gridded/data/juniata_obs.ic/chunks/PHM8GS06CPVYYQ5XEQTG new file mode 100644 index 00000000..4ffa498c Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/PHM8GS06CPVYYQ5XEQTG differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/PWE693HB7A0MC7TTDJ1G b/examples/juniata_gridded/data/juniata_obs.ic/chunks/PWE693HB7A0MC7TTDJ1G new file mode 100644 index 00000000..466aa454 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/PWE693HB7A0MC7TTDJ1G differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/Q8MEWCPV4R28DEDGWZJG b/examples/juniata_gridded/data/juniata_obs.ic/chunks/Q8MEWCPV4R28DEDGWZJG new file mode 100644 index 00000000..c49c76e7 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/Q8MEWCPV4R28DEDGWZJG differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/RNY757368CZE2AKZBFFG b/examples/juniata_gridded/data/juniata_obs.ic/chunks/RNY757368CZE2AKZBFFG new file mode 100644 index 00000000..c3b4f58b Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/RNY757368CZE2AKZBFFG differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/V8J21HGT7H3ZYF8F3C0G b/examples/juniata_gridded/data/juniata_obs.ic/chunks/V8J21HGT7H3ZYF8F3C0G new file mode 100644 index 00000000..5248d40a Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/V8J21HGT7H3ZYF8F3C0G differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/chunks/XE382X12ETH0RA2Z8XC0 b/examples/juniata_gridded/data/juniata_obs.ic/chunks/XE382X12ETH0RA2Z8XC0 new file mode 100644 index 00000000..d58eeebd Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/chunks/XE382X12ETH0RA2Z8XC0 differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/manifests/2Y5BXPWC61SW9G9FCW20 b/examples/juniata_gridded/data/juniata_obs.ic/manifests/2Y5BXPWC61SW9G9FCW20 new file mode 100644 index 00000000..f48689d5 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/manifests/2Y5BXPWC61SW9G9FCW20 differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/manifests/J994M2111PCCVC67SJ20 b/examples/juniata_gridded/data/juniata_obs.ic/manifests/J994M2111PCCVC67SJ20 new file mode 100644 index 00000000..ffd295fc Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/manifests/J994M2111PCCVC67SJ20 differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/manifests/XY7CDE93NPGS33T5JSB0 b/examples/juniata_gridded/data/juniata_obs.ic/manifests/XY7CDE93NPGS33T5JSB0 new file mode 100644 index 00000000..a6f648ed Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/manifests/XY7CDE93NPGS33T5JSB0 differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/overwritten/repo.30714761398729.4286S7ZF82C2YVDWRQ00 b/examples/juniata_gridded/data/juniata_obs.ic/overwritten/repo.30714761398729.4286S7ZF82C2YVDWRQ00 new file mode 100644 index 00000000..6d1cbae5 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/overwritten/repo.30714761398729.4286S7ZF82C2YVDWRQ00 differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/repo b/examples/juniata_gridded/data/juniata_obs.ic/repo new file mode 100644 index 00000000..51586d7e Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/repo differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/snapshots/1CECHNKREP0F1RSTCMT0 b/examples/juniata_gridded/data/juniata_obs.ic/snapshots/1CECHNKREP0F1RSTCMT0 new file mode 100644 index 00000000..666e08aa Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/snapshots/1CECHNKREP0F1RSTCMT0 differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/snapshots/FTT0PBZ068Q649YQZPRG b/examples/juniata_gridded/data/juniata_obs.ic/snapshots/FTT0PBZ068Q649YQZPRG new file mode 100644 index 00000000..42fdbcb9 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/snapshots/FTT0PBZ068Q649YQZPRG differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/transactions/1CECHNKREP0F1RSTCMT0 b/examples/juniata_gridded/data/juniata_obs.ic/transactions/1CECHNKREP0F1RSTCMT0 new file mode 100644 index 00000000..2a1bf888 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/transactions/1CECHNKREP0F1RSTCMT0 differ diff --git a/examples/juniata_gridded/data/juniata_obs.ic/transactions/FTT0PBZ068Q649YQZPRG b/examples/juniata_gridded/data/juniata_obs.ic/transactions/FTT0PBZ068Q649YQZPRG new file mode 100644 index 00000000..cd67eea2 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_obs.ic/transactions/FTT0PBZ068Q649YQZPRG differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/08PA56AQTKRSGG58PFTG b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/08PA56AQTKRSGG58PFTG new file mode 100644 index 00000000..132132e9 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/08PA56AQTKRSGG58PFTG differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/1N1AXXYYGWMRC887BAVG b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/1N1AXXYYGWMRC887BAVG new file mode 100644 index 00000000..2d7d6cf4 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/1N1AXXYYGWMRC887BAVG differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/6YZMPPEHXRRYRYBM02P0 b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/6YZMPPEHXRRYRYBM02P0 new file mode 100644 index 00000000..95994c40 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/6YZMPPEHXRRYRYBM02P0 differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/935T0S2BG6F1GPR03AC0 b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/935T0S2BG6F1GPR03AC0 new file mode 100644 index 00000000..a6bcb4dd Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/935T0S2BG6F1GPR03AC0 differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/D4RJKVDD0VZ29MTAKKY0 b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/D4RJKVDD0VZ29MTAKKY0 new file mode 100644 index 00000000..6be37d21 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/D4RJKVDD0VZ29MTAKKY0 differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/E9Y90M6R0V6Q9VE67BC0 b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/E9Y90M6R0V6Q9VE67BC0 new file mode 100644 index 00000000..b87880db Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/E9Y90M6R0V6Q9VE67BC0 differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/ERDTKF1N4QMKCMEFF53G b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/ERDTKF1N4QMKCMEFF53G new file mode 100644 index 00000000..516f2aa0 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/ERDTKF1N4QMKCMEFF53G differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/FXDT100W5GJN994AX400 b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/FXDT100W5GJN994AX400 new file mode 100644 index 00000000..6cefb692 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/FXDT100W5GJN994AX400 differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/J7E3V2DC8JK18XY7BCWG b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/J7E3V2DC8JK18XY7BCWG new file mode 100644 index 00000000..c4828e24 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/J7E3V2DC8JK18XY7BCWG differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/JYHXZ2CGJ380PJ94XZ00 b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/JYHXZ2CGJ380PJ94XZ00 new file mode 100644 index 00000000..69d6cd45 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/JYHXZ2CGJ380PJ94XZ00 differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/K765M9ZZ0ARX5RM3VY1G b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/K765M9ZZ0ARX5RM3VY1G new file mode 100644 index 00000000..50eb7c19 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/K765M9ZZ0ARX5RM3VY1G differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/KMK1SW4GSD4C1Y1HW7Y0 b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/KMK1SW4GSD4C1Y1HW7Y0 new file mode 100644 index 00000000..55111660 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/KMK1SW4GSD4C1Y1HW7Y0 differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/MK2Z5AJ49K7ZFNJMB74G b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/MK2Z5AJ49K7ZFNJMB74G new file mode 100644 index 00000000..74623b02 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/MK2Z5AJ49K7ZFNJMB74G differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/chunks/ZQBR3098145BTY3TSGMG b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/ZQBR3098145BTY3TSGMG new file mode 100644 index 00000000..9e702a65 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/chunks/ZQBR3098145BTY3TSGMG differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/manifests/VX301M7GS451GPW7GVK0 b/examples/juniata_gridded/data/juniata_qprime.ic/manifests/VX301M7GS451GPW7GVK0 new file mode 100644 index 00000000..5e0e559b Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/manifests/VX301M7GS451GPW7GVK0 differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/manifests/W2Y2A8QP3PQ8G5RH6B4G b/examples/juniata_gridded/data/juniata_qprime.ic/manifests/W2Y2A8QP3PQ8G5RH6B4G new file mode 100644 index 00000000..35bc515f Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/manifests/W2Y2A8QP3PQ8G5RH6B4G differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/manifests/ZJNAC88KHQTKWQT1Q2EG b/examples/juniata_gridded/data/juniata_qprime.ic/manifests/ZJNAC88KHQTKWQT1Q2EG new file mode 100644 index 00000000..850b061b Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/manifests/ZJNAC88KHQTKWQT1Q2EG differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/overwritten/repo.30714761398770.VP16GAW6RSMQV1DVQTW0 b/examples/juniata_gridded/data/juniata_qprime.ic/overwritten/repo.30714761398770.VP16GAW6RSMQV1DVQTW0 new file mode 100644 index 00000000..55a3dbfc Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/overwritten/repo.30714761398770.VP16GAW6RSMQV1DVQTW0 differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/repo b/examples/juniata_gridded/data/juniata_qprime.ic/repo new file mode 100644 index 00000000..ca76533c Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/repo differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/snapshots/1CECHNKREP0F1RSTCMT0 b/examples/juniata_gridded/data/juniata_qprime.ic/snapshots/1CECHNKREP0F1RSTCMT0 new file mode 100644 index 00000000..c8b984db Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/snapshots/1CECHNKREP0F1RSTCMT0 differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/snapshots/95YRZS8ZRBY80YW9D5T0 b/examples/juniata_gridded/data/juniata_qprime.ic/snapshots/95YRZS8ZRBY80YW9D5T0 new file mode 100644 index 00000000..add55ce7 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/snapshots/95YRZS8ZRBY80YW9D5T0 differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/transactions/1CECHNKREP0F1RSTCMT0 b/examples/juniata_gridded/data/juniata_qprime.ic/transactions/1CECHNKREP0F1RSTCMT0 new file mode 100644 index 00000000..2a1bf888 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/transactions/1CECHNKREP0F1RSTCMT0 differ diff --git a/examples/juniata_gridded/data/juniata_qprime.ic/transactions/95YRZS8ZRBY80YW9D5T0 b/examples/juniata_gridded/data/juniata_qprime.ic/transactions/95YRZS8ZRBY80YW9D5T0 new file mode 100644 index 00000000..8786ac9d Binary files /dev/null and b/examples/juniata_gridded/data/juniata_qprime.ic/transactions/95YRZS8ZRBY80YW9D5T0 differ diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/indices_0/c/0 b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/indices_0/c/0 new file mode 100644 index 00000000..880fd36d Binary files /dev/null and b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/indices_0/c/0 differ diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/indices_0/zarr.json b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/indices_0/zarr.json new file mode 100644 index 00000000..1539ea0f --- /dev/null +++ b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/indices_0/zarr.json @@ -0,0 +1,40 @@ +{ + "shape": [ + 26 + ], + "data_type": "int32", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 26 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": 0, + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "zstd", + "configuration": { + "level": 0, + "checksum": false + } + } + ], + "attributes": {}, + "zarr_format": 3, + "node_type": "array", + "storage_transformers": [] +} diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/indices_1/c/0 b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/indices_1/c/0 new file mode 100644 index 00000000..a05f07a5 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/indices_1/c/0 differ diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/indices_1/zarr.json b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/indices_1/zarr.json new file mode 100644 index 00000000..1539ea0f --- /dev/null +++ b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/indices_1/zarr.json @@ -0,0 +1,40 @@ +{ + "shape": [ + 26 + ], + "data_type": "int32", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 26 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": 0, + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "zstd", + "configuration": { + "level": 0, + "checksum": false + } + } + ], + "attributes": {}, + "zarr_format": 3, + "node_type": "array", + "storage_transformers": [] +} diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/lat/c/0 b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/lat/c/0 new file mode 100644 index 00000000..fd833d9b Binary files /dev/null and b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/lat/c/0 differ diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/lat/zarr.json b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/lat/zarr.json new file mode 100644 index 00000000..0455e05a --- /dev/null +++ b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/lat/zarr.json @@ -0,0 +1,40 @@ +{ + "shape": [ + 27 + ], + "data_type": "float32", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 27 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": 0.0, + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "zstd", + "configuration": { + "level": 0, + "checksum": false + } + } + ], + "attributes": {}, + "zarr_format": 3, + "node_type": "array", + "storage_transformers": [] +} diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/length_m/c/0 b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/length_m/c/0 new file mode 100644 index 00000000..2564f5f3 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/length_m/c/0 differ diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/length_m/zarr.json b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/length_m/zarr.json new file mode 100644 index 00000000..0455e05a --- /dev/null +++ b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/length_m/zarr.json @@ -0,0 +1,40 @@ +{ + "shape": [ + 27 + ], + "data_type": "float32", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 27 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": 0.0, + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "zstd", + "configuration": { + "level": 0, + "checksum": false + } + } + ], + "attributes": {}, + "zarr_format": 3, + "node_type": "array", + "storage_transformers": [] +} diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/lon/c/0 b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/lon/c/0 new file mode 100644 index 00000000..25160223 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/lon/c/0 differ diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/lon/zarr.json b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/lon/zarr.json new file mode 100644 index 00000000..0455e05a --- /dev/null +++ b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/lon/zarr.json @@ -0,0 +1,40 @@ +{ + "shape": [ + 27 + ], + "data_type": "float32", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 27 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": 0.0, + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "zstd", + "configuration": { + "level": 0, + "checksum": false + } + } + ], + "attributes": {}, + "zarr_format": 3, + "node_type": "array", + "storage_transformers": [] +} diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/order/c/0 b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/order/c/0 new file mode 100644 index 00000000..60ff0fb1 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/order/c/0 differ diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/order/zarr.json b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/order/zarr.json new file mode 100644 index 00000000..c63fe6b6 --- /dev/null +++ b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/order/zarr.json @@ -0,0 +1,40 @@ +{ + "shape": [ + 27 + ], + "data_type": "int64", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 27 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": 0, + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "zstd", + "configuration": { + "level": 0, + "checksum": false + } + } + ], + "attributes": {}, + "zarr_format": 3, + "node_type": "array", + "storage_transformers": [] +} diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/parent_cell/c/0 b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/parent_cell/c/0 new file mode 100644 index 00000000..4a006dc8 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/parent_cell/c/0 differ diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/parent_cell/zarr.json b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/parent_cell/zarr.json new file mode 100644 index 00000000..10303bd3 --- /dev/null +++ b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/parent_cell/zarr.json @@ -0,0 +1,40 @@ +{ + "shape": [ + 27 + ], + "data_type": "int32", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 27 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": 0, + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "zstd", + "configuration": { + "level": 0, + "checksum": false + } + } + ], + "attributes": {}, + "zarr_format": 3, + "node_type": "array", + "storage_transformers": [] +} diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/slope/c/0 b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/slope/c/0 new file mode 100644 index 00000000..75bdb7d3 Binary files /dev/null and b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/slope/c/0 differ diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/slope/zarr.json b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/slope/zarr.json new file mode 100644 index 00000000..0455e05a --- /dev/null +++ b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/slope/zarr.json @@ -0,0 +1,40 @@ +{ + "shape": [ + 27 + ], + "data_type": "float32", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 27 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": 0.0, + "codecs": [ + { + "name": "bytes", + "configuration": { + "endian": "little" + } + }, + { + "name": "zstd", + "configuration": { + "level": 0, + "checksum": false + } + } + ], + "attributes": {}, + "zarr_format": 3, + "node_type": "array", + "storage_transformers": [] +} diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/values/c/0 b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/values/c/0 new file mode 100644 index 00000000..18b133ed Binary files /dev/null and b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/values/c/0 differ diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/values/zarr.json b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/values/zarr.json new file mode 100644 index 00000000..4074697f --- /dev/null +++ b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/values/zarr.json @@ -0,0 +1,37 @@ +{ + "shape": [ + 26 + ], + "data_type": "uint8", + "chunk_grid": { + "name": "regular", + "configuration": { + "chunk_shape": [ + 26 + ] + } + }, + "chunk_key_encoding": { + "name": "default", + "configuration": { + "separator": "/" + } + }, + "fill_value": 0, + "codecs": [ + { + "name": "bytes" + }, + { + "name": "zstd", + "configuration": { + "level": 0, + "checksum": false + } + } + ], + "attributes": {}, + "zarr_format": 3, + "node_type": "array", + "storage_transformers": [] +} diff --git a/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/zarr.json b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/zarr.json new file mode 100644 index 00000000..a7acf16b --- /dev/null +++ b/examples/juniata_gridded/data/juniata_subreach_adjacency.zarr/zarr.json @@ -0,0 +1,19 @@ +{ + "attributes": { + "format": "COO", + "shape": [ + 27, + 27 + ], + "geodataset": "ddm30", + "subdivided": true, + "cell_id_scheme": "node id = cell id * 1000 + sub index (upstream-most first)", + "data_types": { + "indices_0": "int32", + "indices_1": "int32", + "values": "uint8" + } + }, + "zarr_format": 3, + "node_type": "group" +} diff --git a/examples/juniata_gridded/data/statistics/merit_attribute_statistics_ddm30_conus_attributes.nc.json b/examples/juniata_gridded/data/statistics/merit_attribute_statistics_ddm30_conus_attributes.nc.json new file mode 100644 index 00000000..a97e3c71 --- /dev/null +++ b/examples/juniata_gridded/data/statistics/merit_attribute_statistics_ddm30_conus_attributes.nc.json @@ -0,0 +1,202 @@ +{ + "meanelevation": { + "min": 0.004861295323337654, + "max": 3369.19688318091, + "mean": 703.9322811238138, + "std": 660.8384435790672, + "p10": 69.40527726837776, + "p90": 1759.9149382594505 + }, + "NDVI": { + "min": -0.07139881334372336, + "max": 0.8012113621976084, + "mean": 0.39212428904603336, + "std": 0.16051253610065513, + "p10": 0.18850096432721256, + "p90": 0.6226797447717318 + }, + "snow_fraction": { + "min": 7.450250524098686e-05, + "max": 0.6698027268301546, + "mean": 0.09175650361244392, + "std": 0.09728262891609547, + "p10": 0.0013670597900336228, + "p90": 0.23224604059529871 + }, + "FW": { + "min": 0.0, + "max": 1.0, + "mean": 0.058816597108576625, + "std": 0.14434725529012868, + "p10": 0.000209455598065324, + "p90": 0.13832576788124348 + }, + "Ksat": { + "min": 1.630860676223944, + "max": 360.32950299120023, + "mean": 29.965475362872624, + "std": 36.42079478433762, + "p10": 6.477697036109817, + "p90": 71.25291135094982 + }, + "ALPHA": { + "min": 0.0016555487306490662, + "max": 0.10123577553330375, + "mean": 0.020262717081631222, + "std": 0.012544429442543381, + "p10": 0.007331282362348112, + "p90": 0.03696019671172372 + }, + "N": { + "min": 1.1986454776628077, + "max": 1.6693176918692725, + "mean": 1.3249985062017362, + "std": 0.06297205116060489, + "p10": 1.247119541409867, + "p90": 1.4068147898142649 + }, + "ORMC": { + "min": 0.4428205578327663, + "max": 52.20994461155595, + "mean": 8.741051363357073, + "std": 6.635191044133652, + "p10": 2.4944028294625813, + "p90": 16.503845225470094 + }, + "WCpF2": { + "min": 0.16847541756909717, + "max": 0.6116269056806527, + "mean": 0.37807038847617563, + "std": 0.06387835873786013, + "p10": 0.2910144102103489, + "p90": 0.45110854385255456 + }, + "WCsat": { + "min": 0.3936480963612873, + "max": 0.7172914560492192, + "mean": 0.49274685483597847, + "std": 0.05074944293515737, + "p10": 0.43361513124774137, + "p90": 0.5675219682404052 + }, + "SoilGrids1km_clay": { + "min": 2.0618805485446328, + "max": 53.02753344037019, + "mean": 22.060060387691486, + "std": 7.029637572972272, + "p10": 13.18350418518392, + "p90": 31.155411883878447 + }, + "SoilGrids1km_sand": { + "min": 4.115237064847427, + "max": 94.81344250065135, + "mean": 40.129230147520076, + "std": 14.779979045583923, + "p10": 20.56785571370896, + "p90": 57.42126506346657 + }, + "SoilGrids1km_silt": { + "min": 2.8477707313651184, + "max": 70.83410231479239, + "mean": 37.810677628739946, + "std": 10.427255627718067, + "p10": 24.797150150008086, + "p90": 51.42990197184438 + }, + "meanslope": { + "min": 0.0, + "max": 27.263860282239285, + "mean": 3.397809137240753, + "std": 4.139947075799843, + "p10": 0.4088235262148592, + "p90": 8.639461607208487 + }, + "meanP": { + "min": 57.27565084232705, + "max": 3139.0101362051923, + "mean": 744.2556414051538, + "std": 396.6552667548469, + "p10": 293.9332777773968, + "p90": 1259.2973333326408 + }, + "meanTa": { + "min": -4.214311336403693, + "max": 25.418786005990317, + "mean": 9.397976360355777, + "std": 6.757331808995632, + "p10": 0.6839523293738895, + "p90": 18.991226844959154 + }, + "ETPOT_Hargr": { + "min": 506.3604392672467, + "max": 1979.1649480526967, + "mean": 1082.352619034022, + "std": 328.609970222106, + "p10": 678.405916454727, + "p90": 1554.4528788758657 + }, + "seasonality_P": { + "min": 0.045816066321010616, + "max": 1.018944280140316, + "mean": 0.36544183968232397, + "std": 0.2032374264509107, + "p10": 0.13108854272027207, + "p90": 0.6138264486008584 + }, + "seasonality_PET": { + "min": 0.17664974526005303, + "max": 0.9025561776057203, + "mean": 0.5753634587230027, + "std": 0.16810704414639588, + "p10": 0.33392438075236064, + "p90": 0.8013925131806772 + }, + "snowfall_fraction": { + "min": 0.0, + "max": 0.7009380043020851, + "mean": 0.1738869544314959, + "std": 0.15691078415008627, + "p10": 0.0, + "p90": 0.38890798297778867 + }, + "aridity": { + "min": 0.24188227366112164, + "max": 28.526506740792904, + "mean": 2.2754431321907855, + "std": 2.543380125305846, + "p10": 0.7385919473552908, + "p90": 4.607038839728651 + }, + "Porosity": { + "min": 0.00999999977648258, + "max": 0.2793158438904144, + "mean": 0.13425318477349144, + "std": 0.06511777897495456, + "p10": 0.027224847045667213, + "p90": 0.21367892881640052 + }, + "permeability": { + "min": -16.500000000000004, + "max": -10.916314141407062, + "mean": -13.783503510211885, + "std": 1.0981717690721766, + "p10": -15.275190199691304, + "p90": -12.515690630734188 + }, + "catchsize": { + "min": 1871.0049494176442, + "max": 2818.3185855515253, + "mean": 2311.363148928616, + "std": 260.79175193002754, + "p10": 1955.7352702129808, + "p90": 2670.174390705445 + }, + "log10_uparea": { + "min": 3.2720749363519226, + "max": 6.509329060049167, + "mean": 3.8717224718970327, + "std": 0.6399331225807444, + "p10": 3.3218508865933143, + "p90": 4.795751524059255 + } +} \ No newline at end of file diff --git a/examples/juniata_gridded/ddrs.yaml b/examples/juniata_gridded/ddrs.yaml new file mode 100644 index 00000000..e16b88fa --- /dev/null +++ b/examples/juniata_gridded/ddrs.yaml @@ -0,0 +1,119 @@ +# Gridded (ISIMIP DDM30) Juniata sample — the ddrs mirror of DDR's +# examples/juniata_gridded (DeepGroundwater/ddr PR #194). +# +# Same gauge as examples/juniata (USGS 01567000, 8,657 km²) but the river +# network is the 0.5° DDM30 grid: four cells split into 27 Courant≈1 +# sub-reaches of ~7 km. The 1.6 MB bundle under examples/juniata_gridded/data/ +# is byte-identical to DDR's; nothing external is needed. +# +# Run FROM THE REPO ROOT (data paths below are repo-root-relative): +# +# ddrs --config examples/juniata_gridded/ddrs.yaml plan +# ddrs --config examples/juniata_gridded/ddrs.yaml run --workflow train-and-test --backend cpu +# +# The workspace (.ddrs/) is created beside this file. `plan` relabels the +# sub-reach store into ddrs's subdivided layout and cuts the gauge subgraph +# into .ddrs/adjacency// (parent = grid cell, pieces = sub-reaches). +# +# Hyperparameters mirror DDR's examples/juniata_gridded/train_gridded.py — +# keep them in sync when comparing against the Python implementation. + +mode: training +workflow: train-and-test +geodataset: ddm30 +device: 0 +seed: 42 +np_seed: 42 + +data_sources: + # DDR's sub-reach adjacency zarr (`order`, `parent_cell`, `indices_0/1`, + # `length_m`, `slope`, `lat`, `lon`), 27 rows over 4 cells. This is the + # gridded counterpart of `geospatial_fabric`: ddrs builds the managed + # stores from it on first `plan`. + gridded_network: examples/juniata_gridded/data/juniata_subreach_adjacency.zarr + # 25 attributes for all 5,198 CONUS cells (`COMID` = cell id), so the + # normalisation statistics match a CONUS run. Sidecar: + # data/statistics/merit_attribute_statistics_ddm30_conus_attributes.nc.json + # (generated by scripts/build_attribute_statistics.py; ddrs never recomputes). + attributes: examples/juniata_gridded/data/ddm30_conus_attributes.nc + # Daily lateral inflow per CELL (`Qr(divide_id, time)`, divide_id = cell + # id): dHBV2 UH retrospective area-weighted onto the four cells. Each cell's + # Q' is split evenly across its sub-reaches at routing time. + streamflow: examples/juniata_gridded/data/juniata_qprime.ic + # USGS daily observed discharge at 01567000 (icechunk store). + observations: examples/juniata_gridded/data/juniata_obs.ic + # One-row gauge table. The `cell` column (DDR's name) is the snapped DDM30 + # cell; ddrs reads it as the outlet id and routes to that cell's most + # downstream sub-reach. + gages: examples/juniata_gridded/data/juniata_gage.csv + +experiment: + # Single gauge → one mini-batch (one random rho-day window) per epoch, so + # `epochs` is the optimizer-step count. 30 steps demonstrates that the + # gridded network trains and routes; it is not a converged model. + batch_size: 1 + start_time: 1981/10/01 + end_time: 1995/09/30 + epochs: 30 + rho: 90 + shuffle: true + warmup: 5 + # train_gridded.py uses a constant Adam lr of 1e-3 (no halving). + learning_rate: + 1: 0.001 + grad_clip_max_norm: 1.0 + +kan_head: + hidden_size: 21 + num_hidden_layers: 2 + grid: 50 + k: 2 + input_var_names: + - SoilGrids1km_clay + - aridity + - meanelevation + - meanP + - NDVI + - meanslope + - log10_uparea + - SoilGrids1km_sand + - ETPOT_Hargr + - Porosity + learnable_parameters: + - n + - q_spatial + - p_spatial + +params: + # train_gridded.py::make_config bounds (n is narrower than the MERIT + # Juniata example's [0.015, 0.25]). + parameter_ranges: + n: [0.02, 0.2] + q_spatial: [0.0, 1.0] + p_spatial: [1.0, 200.0] + attribute_minimums: + discharge: 1.0e-4 + slope: 1.0e-3 + velocity: 0.01 + depth: 0.01 + bottom_width: 0.01 + defaults: + p_spatial: 21.0 + log_space_parameters: + - p_spatial + tau: 9 + # The sub-reach store is already split by DDR; ddrs's MERIT reach plan + # must stay off (config load rejects `enabled: true` with gridded_network). + subdivision: + enabled: false + # CPU-friendly: no CUDA requirement anywhere in the sample. + sparse_solver: cpu + use_cuda_graphs: false + +# Evaluation overlay (water years 1996–2010, as train_gridded.py). +# CAUTION — batch_size changes meaning: GAUGES in training, DAYS here. +testing: + start_time: 1995/10/01 + end_time: 2010/09/30 + batch_size: 15 + rho: null diff --git a/scripts/build_attribute_statistics.py b/scripts/build_attribute_statistics.py new file mode 100644 index 00000000..81fe0e23 --- /dev/null +++ b/scripts/build_attribute_statistics.py @@ -0,0 +1,56 @@ +"""Build the attribute-statistics sidecar JSON that ddrs reads at data-load time. + +ddrs never recomputes normalization statistics itself (see +`src/data/dataset.rs::stats_path_from_attrs`): it expects a +`/statistics/merit_attribute_statistics_.json` +file to already exist next to any attributes NetCDF it opens, MERIT or +DDM30-gridded alike (the `merit_` prefix is a fixed template, not a claim +about the network type). This script produces that sidecar, mirroring DDR's +`src/ddr/io/statistics.py::set_statistics` but without needing a DDR `Config` +object, so it can run standalone under DDR's venv against any attributes file. + +Usage: + python build_attribute_statistics.py [--out-dir ] + +Default `--out-dir` is `/statistics`. +""" + +import argparse +import json +from pathlib import Path + +import numpy as np +import xarray as xr + + +def main() -> None: + parser = argparse.ArgumentParser() + parser.add_argument("attributes", type=Path) + parser.add_argument("--out-dir", type=Path, default=None) + args = parser.parse_args() + + out_dir = args.out_dir or (args.attributes.parent / "statistics") + out_dir.mkdir(exist_ok=True, parents=True) + out_file = out_dir / f"merit_attribute_statistics_{args.attributes.name}.json" + + ds = xr.open_dataset(args.attributes) + stats = {} + for attr in list(ds.data_vars.keys()): + data = ds[attr].values + stats[attr] = { + "min": float(np.nanmin(data, axis=0)), + "max": float(np.nanmax(data, axis=0)), + "mean": float(np.nanmean(data, axis=0)), + "std": float(np.nanstd(data, axis=0)), + "p10": float(np.nanpercentile(data, 10, axis=0)), + "p90": float(np.nanpercentile(data, 90, axis=0)), + } + + with open(out_file, "w") as f: + json.dump(stats, f, indent=2) + + print(f"Wrote {out_file}") + + +if __name__ == "__main__": + main() diff --git a/scripts/make_streamflow_fixtures.py b/scripts/make_streamflow_fixtures.py index 4687a808..01a26271 100644 --- a/scripts/make_streamflow_fixtures.py +++ b/scripts/make_streamflow_fixtures.py @@ -11,13 +11,18 @@ qr_daily.ic : 4 divides x 10 days, Qr[j, t] = (j+1)*100 + t qr_hourly.ic : 4 divides x 240 hours, Qr[j, h] = (j+1)*1000 + h qr_minutes.ic : sniff-rejection fixture (units "minutes since ...") + qr_daily_time_major.ic : qr_daily's values stored TIME-MAJOR, Qr(time, divide_id) + -- the layout DDR's gridded (DDM30) Q' stores use; the reader + must return identical values from both. + +Stores that already exist are left alone (icechunk snapshot ids are not +deterministic, so rewriting would churn tracked bytes); delete one to regenerate. Note: xarray normalises "hours since 1981-01-01 00:00:00" to drop the time component on write. We patch the zarr attr back to the full string after to_zarr() so the on-disk CF units string is exactly as documented above. """ from pathlib import Path -import shutil import icechunk import numpy as np @@ -35,14 +40,19 @@ def write_store( time_units: str, *, time_units_on_disk: str | None = None, + time_major: bool = False, ) -> None: - shutil.rmtree(path, ignore_errors=True) + if path.exists(): + print(f"exists, skipping: {path} (delete it to regenerate)") + return + dims = ["time", "divide_id"] if time_major else ["divide_id", "time"] + data = qr.T if time_major else qr storage = icechunk.local_filesystem_storage(str(path)) repo = icechunk.Repository.create(storage) session = repo.writable_session("main") ds = xr.Dataset( data_vars={ - "Qr": (["divide_id", "time"], qr.astype(np.float32), {"units": "m^3/s"}), + "Qr": (dims, data.astype(np.float32), {"units": "m^3/s"}), }, coords={ "divide_id": ("divide_id", DIVIDES), @@ -70,6 +80,10 @@ def main() -> None: daily_times = np.datetime64("1981-01-01") + np.arange(n_days).astype("timedelta64[D]") daily = (np.arange(4)[:, None] + 1) * 100 + np.arange(n_days)[None, :] write_store(FIXTURES / "qr_daily.ic", daily_times, daily, "days since 1981-01-01") + write_store( + FIXTURES / "qr_daily_time_major.ic", daily_times, daily, "days since 1981-01-01", + time_major=True, + ) n_hours = n_days * 24 hourly_times = np.datetime64("1981-01-01T00") + np.arange(n_hours).astype("timedelta64[h]") diff --git a/scripts/snap_gridded_gauges.py b/scripts/snap_gridded_gauges.py new file mode 100644 index 00000000..10409a84 --- /dev/null +++ b/scripts/snap_gridded_gauges.py @@ -0,0 +1,103 @@ +"""Build the CONUS gauge CSV for ddrs's gridded (ISIMIP DDM30) routing path. + +Gridded counterpart of the gages_3000 gauge selection used by the MERIT-fabric +path: takes `gages_3000.csv`, keeps gauges with DRAIN_SQKM above a minimum +(default 2,000 km2 -- a 0.5 degree DDM30 cell is ~2,350 km2, so 2,416 of the +3,211 gages_3000 gauges are sub-cell and cannot be represented by a single +routing element), snaps each surviving gauge to a DDM30 cell via +`ddr_benchmarks.gridded.snap_gauges` (drainage-area-matched 3x3 neighborhood +search against per-cell upstream area from the sub-reach adjacency zarr), then +keeps gauges present in the USGS observation store with daily coverage above a +minimum (default 0.9) over 1981-10-01..1997-09-30 (train_conus.py measures coverage +from `--train-start` to its `--test-end` default), which yields DDR's 620 gauges. +Mirrors `~/projects/ddr/examples/juniata_gridded/train_conus.py::build_network` +(lines 51-91) exactly, including its default thresholds. + +Output columns: STAID,STANAME,DRAIN_SQKM,LAT_GAGE,LNG_GAGE,COMID,da_ratio +`COMID` is the snapped DDM30 cell id (ddrs's gauge reader joins on `COMID`); +`da_ratio` is snap_gauges's cell-upstream-area / gauge-DA ratio, informational. + +Run under DDR's venv (needs `ddr_benchmarks` and `ddr` importable): + cd ~/projects/ddr && \ + ~/projects/ddr/.venv/bin/python ~/projects/ddrs/scripts/snap_gridded_gauges.py +""" + +import argparse +import logging +from pathlib import Path + +import numpy as np +import pandas as pd +import zarr + +SUBREACH = Path("/home/tbindas/projects/ddr/data/ddm30/ddm30_subreach_adjacency.zarr") +OBS = Path("/mnt/ssd1/data/icechunk/usgs_daily_observations") +GAGES = Path("/home/tbindas/projects/ddr/references/gage_info/gages_3000.csv") +OUT_DEFAULT = "/home/tbindas/projects/ddr/references/gage_info/ddm30_conus_gauges.csv" + +MIN_DA_KM2 = 2000.0 +MIN_COVERAGE = 0.9 +START = "1981-10-01" +END = "1997-09-30" # train_conus.py measures coverage over train_start..test_end (default 1997-09-30) + +log = logging.getLogger("snap_gridded_gauges") + + +def main() -> None: + parser = argparse.ArgumentParser(description=__doc__) + parser.add_argument("--out", default=OUT_DEFAULT) + parser.add_argument("--min-da-km2", type=float, default=MIN_DA_KM2) + parser.add_argument("--min-coverage", type=float, default=MIN_COVERAGE) + args = parser.parse_args() + + logging.basicConfig(level=logging.INFO, format="%(message)s") + + from ddr_benchmarks.gridded import cell_areas_km2, snap_gauges, topo_accumulate + + from ddr.io.readers import read_ic + + g = zarr.open_group(str(SUBREACH), mode="r") + parent = g["parent_cell"][:] + node_ids = g["order"][:] + rows, cols = g["indices_0"][:], g["indices_1"][:] + n_all = len(node_ids) + dn = np.full(n_all, -1, dtype=np.int64) + dn[cols] = rows + + _, inv, cnt = np.unique(parent, return_inverse=True, return_counts=True) + k_per_node = cnt[inv] + uparea = topo_accumulate(cell_areas_km2(g["lat"][:]) / k_per_node, dn) + + # the cell's outlet is its most downstream sub-reach + sub_idx = node_ids % 1000 + outlet_of_cell: dict[int, int] = {} + for i in range(n_all): + c = int(parent[i]) + if c not in outlet_of_cell or sub_idx[i] > sub_idx[outlet_of_cell[c]]: + outlet_of_cell[c] = i + ua_cell = {c: float(uparea[i]) for c, i in outlet_of_cell.items()} + + gauges = pd.read_csv(GAGES, dtype={"STAID": str}) + gauges["STAID"] = gauges.STAID.str.zfill(8) + gauges = snap_gauges(gauges[gauges.DRAIN_SQKM > args.min_da_km2], ua_cell) + obs_ds = read_ic(str(OBS)) + gauges = gauges[gauges.STAID.isin(set(obs_ds.gage_id.values.astype(str)))] + cov = ( + obs_ds["streamflow"] + .sel(gage_id=gauges.STAID.tolist(), time=slice(START, END)) + .notnull() + .mean("time") + .values + ) + gauges = gauges[cov > args.min_coverage].reset_index(drop=True) + log.info("%d gauges after DA match and observation coverage", len(gauges)) + + out = gauges[["STAID", "STANAME", "DRAIN_SQKM", "LAT_GAGE", "LNG_GAGE", "cell", "da_ratio"]].rename( + columns={"cell": "COMID"} + ) + out.to_csv(args.out, index=False, quoting=1) # quoting=1 == csv.QUOTE_ALL, keeps STAID's leading zeros + log.info("wrote %d gauges to %s", len(out), args.out) + + +if __name__ == "__main__": + main() diff --git a/src/adjacency/cache.rs b/src/adjacency/cache.rs index 487a8bf8..17bdd3ca 100644 --- a/src/adjacency/cache.rs +++ b/src/adjacency/cache.rs @@ -37,8 +37,12 @@ use serde::{Deserialize, Serialize}; use crate::adjacency::build::{build_conus_adjacency, BuildError, ConusAdjacency}; use crate::adjacency::fabric::{read_fabric_records, resolve_fabric}; use crate::adjacency::gauges::build_gauge_subgraphs; +use crate::adjacency::gridded::GriddedNetwork; use crate::adjacency::subdivide::{plan_reaches, subdivide, ReachPlan, SubdividedAdjacency}; -use crate::adjacency::zarr_write::{write_conus_store_subdivided, write_gauges_store}; +use crate::adjacency::zarr_write::{ + write_conus_store_subdivided, write_conus_store_subdivided_as, write_gauges_store, + write_gauges_store_as, +}; use crate::adjacency::BUILDER_VERSION; use crate::config::Subdivision; use crate::data::error::DataError; @@ -117,6 +121,14 @@ struct CacheManifest { build_duration_secs: f64, /// `git rev-parse HEAD` at build time, empty string if not available. ddrs_git_sha: String, + /// `"merit"` (fabric builder) or `"ddm30"` (`gridded_network`). Defaulted + /// for manifests written before the gridded path existed. + #[serde(default = "default_network_kind")] + network_kind: String, +} + +fn default_network_kind() -> String { + "merit".to_string() } // ── public API ─────────────────────────────────────────────────────────────── @@ -268,21 +280,12 @@ pub fn resolve_or_build( dropped_comids, build_duration_secs, ddrs_git_sha: capture_git_sha(), + network_kind: "merit".to_string(), }; let json = serde_json::to_string_pretty(&manifest)?; fs::write(tmp_dir.join("manifest.json"), json)?; - // Atomic rename into place. If the target already exists (lost race - // between two concurrent builds) treat it as a hit and discard the tmp. - match fs::rename(&tmp_dir, &cache_dir) { - Ok(()) => {} - Err(_) if is_cache_hit(&cache_dir) - && store_paths(&cache_dir, &fabric_stem).conus.is_dir() => { - // Another process won the race; discard our tmp. - let _ = fs::remove_dir_all(&tmp_dir); - } - Err(e) => return Err(AdjacencyCacheError::Io(e)), - } + commit_build(&tmp_dir, &cache_dir, &fabric_stem)?; println!( " adjacency cache written to {} ({:.1}s)", @@ -297,6 +300,148 @@ pub fn resolve_or_build( }) } +/// Resolve the adjacency cache for a **gridded** (DDM30 sub-reach) network, +/// building if necessary — the `data_sources.gridded_network` counterpart of +/// [`resolve_or_build`]. +/// +/// The store is DDR's sub-reach adjacency zarr (`gridded::GriddedNetwork`). +/// It is already a subdivided graph, so there is no fabric read and no reach +/// plan: the build relabels it into ddrs's parent/piece layout, writes it as +/// `_adjacency.zarr`, and cuts the per-gauge subgraphs from the sub-reach +/// rows — `position_lookup` maps a cell id to its **last** row, which is the +/// cell's outlet piece, so a gauge snapped to a cell is read where DDR's +/// trainers read it (`outlet_of_cell`). `` is the store's file stem with +/// a trailing `_adjacency` dropped (`juniata_subreach_adjacency.zarr` → +/// `juniata_subreach_{,gages_}adjacency.zarr`). +/// +/// Key: blake3(store fingerprint ∥ gages fingerprint ∥ "gridded" ∥ +/// BUILDER_VERSION). The store fingerprint hashes every file under the zarr +/// directory (sorted relative path ∥ NUL ∥ bytes), so any array edit rebuilds. +pub fn resolve_or_build_gridded( + workspace_root: &Path, + gridded: &Path, + gages_csv: &Path, +) -> Result { + let stem = gridded + .file_stem() + .map(|s| s.to_string_lossy().into_owned()) + .unwrap_or_else(|| "ddm30".into()); + let stem = stem.strip_suffix("_adjacency").unwrap_or(&stem).to_string(); + + let store_fp = dir_fingerprint(gridded)?; + let gages_fp = gages_fingerprint(gages_csv)?; + let key = content_key_gridded(&store_fp, &gages_fp); + + let cache_dir = adjacency_cache_dir(workspace_root, &key); + if is_cache_hit(&cache_dir) { + let paths = store_paths(&cache_dir, &stem); + if paths.conus.is_dir() && paths.gages.is_dir() { + return Ok(AdjacencyCacheOutcome { paths, key, cache_hit: true }); + } + eprintln!( + " adjacency cache {} is incomplete (manifest without stores) — rebuilding", + cache_dir.display() + ); + fs::remove_dir_all(&cache_dir)?; + } + + println!( + " building DDM30 sub-reach adjacency from {}", + gridded.display() + ); + let tmp_dir = workspace_root + .join("adjacency") + .join(format!(".tmp-{}", key)); + let _ = fs::remove_dir_all(&tmp_dir); + fs::create_dir_all(&tmp_dir)?; + let t0 = Instant::now(); + + let network = GriddedNetwork::open(gridded).map_err(AdjacencyCacheError::Data)?; + let expanded = network.to_subdivided(); + let n = expanded.order.len(); + let n_parent = expanded.parent_order.len(); + let nnz = expanded.rows.len(); + println!( + " gridded network: {n_parent} cells → {n} sub-reaches ({:.2}x), {nnz} edges", + n as f64 / n_parent as f64 + ); + + let conus_dest = tmp_dir.join(format!("{stem}_adjacency.zarr")); + write_conus_store_subdivided_as(&expanded, &conus_dest, GRIDDED_GEODATASET) + .map_err(AdjacencyCacheError::Data)?; + + // Same sub-reach-space view as the fabric path: a gauge's cell resolves to + // that cell's outlet piece (last row). + let expanded_view = ConusAdjacency { + order: expanded.order.clone(), + rows: expanded.rows.clone(), + cols: expanded.cols.clone(), + length_m: expanded.length_m.clone(), + slope: expanded.slope.clone(), + dropped_comids: Vec::new(), + }; + let subgraphs = build_gauge_subgraphs(&expanded_view, gages_csv) + .map_err(AdjacencyCacheError::Data)?; + let n_gauges = subgraphs.len(); + let gages_dest = tmp_dir.join(format!("{stem}_gages_adjacency.zarr")); + write_gauges_store_as(&subgraphs, n, &gages_dest, GRIDDED_GEODATASET) + .map_err(AdjacencyCacheError::Data)?; + + let build_duration_secs = t0.elapsed().as_secs_f64(); + let manifest = CacheManifest { + key: key.clone(), + builder_version: BUILDER_VERSION, + fabric_path: gridded.to_path_buf(), + fabric_resolved_path: gridded.to_path_buf(), + fabric_layer: None, + fabric_fingerprint: store_fp, + gages_path: gages_csv.to_path_buf(), + gages_fingerprint: gages_fp, + n, + n_parent, + nnz, + n_gauges, + dropped_comids: Vec::new(), + build_duration_secs, + ddrs_git_sha: capture_git_sha(), + network_kind: GRIDDED_GEODATASET.to_string(), + }; + fs::write( + tmp_dir.join("manifest.json"), + serde_json::to_string_pretty(&manifest)?, + )?; + + commit_build(&tmp_dir, &cache_dir, &stem)?; + + println!( + " adjacency cache written to {} ({:.1}s)", + cache_dir.display(), + build_duration_secs + ); + Ok(AdjacencyCacheOutcome { + paths: store_paths(&cache_dir, &stem), + key, + cache_hit: false, + }) +} + +/// `geodataset` attr / manifest `network_kind` for the gridded path. +const GRIDDED_GEODATASET: &str = "ddm30"; + +/// Atomically move a finished build into place. If the target already exists +/// (lost race between two concurrent builds) treat it as a hit and discard the +/// temp dir. +fn commit_build(tmp_dir: &Path, cache_dir: &Path, stem: &str) -> Result<(), AdjacencyCacheError> { + match fs::rename(tmp_dir, cache_dir) { + Ok(()) => Ok(()), + Err(_) if is_cache_hit(cache_dir) && store_paths(cache_dir, stem).conus.is_dir() => { + let _ = fs::remove_dir_all(tmp_dir); + Ok(()) + } + Err(e) => Err(AdjacencyCacheError::Io(e)), + } +} + /// Resolve `/adjacency//` for a given key. pub fn adjacency_cache_dir(workspace_root: &Path, key: &str) -> PathBuf { workspace_root.join("adjacency").join(key) @@ -495,6 +640,69 @@ fn content_key( hex.as_str()[..16].to_string() } +/// Content key for [`resolve_or_build_gridded`]: blake3(store_fp ∥ gages_fp ∥ +/// "gridded" ∥ BUILDER_VERSION)[..16]. The literal separates the namespace +/// from fabric keys; `params.subdivision` is absent because config validation +/// rejects it with `gridded_network`. +fn content_key_gridded(store_fp: &str, gages_fp: &str) -> String { + let mut h = blake3::Hasher::new(); + h.update(store_fp.as_bytes()); + h.update(b"\n"); + h.update(gages_fp.as_bytes()); + h.update(b"\ngridded\n"); + h.update(BUILDER_VERSION.to_le_bytes().as_ref()); + let hex = h.finalize().to_hex(); + hex.as_str()[..16].to_string() +} + +#[doc(hidden)] +pub fn content_key_gridded_for_test(store_fp: &str, gages_fp: &str) -> String { + content_key_gridded(store_fp, gages_fp) +} + +/// blake3 over every regular file under `dir`, recursively, in sorted +/// relative-path order (path ∥ NUL ∥ bytes). The fingerprint of a zarr store: +/// metadata and chunks alike, so editing any array — or a store attribute — +/// changes the key. +pub fn dir_fingerprint(dir: &Path) -> Result { + fn walk(root: &Path, dir: &Path, out: &mut Vec) -> std::io::Result<()> { + for entry in fs::read_dir(dir)? { + let p = entry?.path(); + if p.is_dir() { + walk(root, &p, out)?; + } else if p.is_file() { + out.push(p); + } + } + Ok(()) + } + if !dir.is_dir() { + return Err(AdjacencyCacheError::Io(std::io::Error::new( + std::io::ErrorKind::NotFound, + format!("gridded_network store is not a directory: {}", dir.display()), + ))); + } + let mut files = Vec::new(); + walk(dir, dir, &mut files).map_err(AdjacencyCacheError::Io)?; + files.sort(); + let mut hasher = blake3::Hasher::new(); + let mut buf = vec![0u8; 65536]; + for f in &files { + let rel = f.strip_prefix(dir).unwrap_or(f); + hasher.update(rel.to_string_lossy().as_bytes()); + hasher.update(&[0]); + let mut reader = BufReader::new(fs::File::open(f).map_err(AdjacencyCacheError::Io)?); + loop { + let n = reader.read(&mut buf).map_err(AdjacencyCacheError::Io)?; + if n == 0 { + break; + } + hasher.update(&buf[..n]); + } + } + Ok(hasher.finalize().to_hex().to_string()) +} + /// Test hook for [`content_key`]. Integration tests live in their own crate and /// cannot reach a private fn; exposing the real function (rather than /// reimplementing the hash in the test) is what makes diff --git a/src/adjacency/gridded.rs b/src/adjacency/gridded.rs new file mode 100644 index 00000000..4c3f9a22 --- /dev/null +++ b/src/adjacency/gridded.rs @@ -0,0 +1,358 @@ +//! Read a DDR gridded (ISIMIP DDM30) **sub-reach** adjacency zarr and express +//! it as the [`SubdividedAdjacency`] the managed builder already knows how to +//! write and the dataset already knows how to route. +//! +//! The store is the output of DDR's `scripts/build_subdivided_adjacency.py` +//! (DeepGroundwater/ddr #194). Per node (sub-reach) it carries: +//! +//! ```text +//! order int64 node id = cell_id * 1000 + sub_index (upstream-most first) +//! parent_cell int32 DDM30 cell id = row * 720 + col (row 0 = 55.75°S) +//! indices_0/1 int32 COO (downstream row, upstream col), strictly lower-tri +//! length_m float32 piece length (the cell's flow length / pieces) +//! slope float32 inherited from the cell +//! lat, lon float32 cell centre (unused here) +//! ``` +//! +//! and the pieces of a cell form a chain that is **contiguous in `order`**, +//! upstream first. That is exactly ddrs's subdivided layout (`zarr.rs:36-46`): +//! parent `p` owns rows `parent_offset[p]..parent_offset[p+1]` and its outlet — +//! the row a gauge is read at — is the last one. So the translation is a +//! relabelling, not a rebuild: +//! +//! ```text +//! DDM30 store SubdividedAdjacency +//! order [c₀·1000+0, c₀·1000+1, c₁·1000+0, …] +//! parent_cell [c₀, c₀, c₁, …] ──► order = [c₀, c₀, c₁, …] +//! parent_order = [c₀, c₁, …] +//! parent_offset= [0, 2, 3, …] +//! indices_0/1, length_m, slope ──► rows/cols, length_m, slope (copied) +//! ``` +//! +//! Attributes and Q′ are keyed on the cell id (DDR writes both stores with +//! `COMID`/`divide_id` = cell id), so `AttributesStore::open(…, parent_order)` +//! and the streamflow reader need no change, and `pieces_per_row_divisor` +//! splits each cell's Q′ across its pieces the way DDR's trainers divide by +//! `k_per_node`. +//! +//! **Not supported:** DDR's unsplit cell adjacency (`ddm30_adjacency.zarr`, +//! no `parent_cell`). DDR's trainers use the sub-reach store; `open` refuses a +//! store without `parent_cell` with a message saying which script produces one. + +use std::path::{Path, PathBuf}; +use std::sync::Arc; + +use zarrs::array::Array as ZarrArray; +use zarrs::filesystem::FilesystemStore; +use zarrs::group::Group; +use zarrs::storage::ReadableStorage; + +use crate::adjacency::subdivide::SubdividedAdjacency; +use crate::data::error::{DataError, Result}; + +/// A DDM30 sub-reach network as read from disk, validated. +#[derive(Debug, Clone)] +pub struct GriddedNetwork { + pub path: PathBuf, + /// `order`: DDR node ids (`cell * 1000 + sub`). Kept for diagnostics only. + pub node_ids: Vec, + /// DDM30 cell id per sub-reach row. + pub parent_cell: Vec, + /// COO `indices_0` — downstream row. + pub rows: Vec, + /// COO `indices_1` — upstream row. + pub cols: Vec, + pub length_m: Vec, + pub slope: Vec, + /// Row where each cell's contiguous run starts, plus a final `n`. + /// `parent_offset[p]..parent_offset[p+1]` are cell `p`'s rows. + parent_offset: Vec, +} + +impl GriddedNetwork { + /// Open and validate a DDM30 sub-reach store. + pub fn open(path: impl Into) -> Result { + let path = path.into(); + let storage: ReadableStorage = + Arc::new(FilesystemStore::new(&path).map_err(|e| zarr_err(&path, e))?); + let _root = Group::open(storage.clone(), "/").map_err(|e| zarr_err(&path, e))?; + + let node_ids = read_i64_or_i32(&storage, &path, "/order")?; + let n = node_ids.len(); + let parent_cell = match read_i32(&storage, &path, "/parent_cell") { + Ok(v) => v, + Err(_) => { + return Err(DataError::Malformed { + path, + message: "gridded_network store has no `parent_cell` array — this is \ + DDR's unsplit cell adjacency, not a sub-reach store. ddrs \ + routes the sub-reach network only; build one with DDR's \ + scripts/build_subdivided_adjacency.py" + .to_string(), + }); + } + }; + let rows = read_i32(&storage, &path, "/indices_0")?; + let cols = read_i32(&storage, &path, "/indices_1")?; + let length_m = read_f32(&storage, &path, "/length_m")?; + let slope = read_f32(&storage, &path, "/slope")?; + + let malformed = |message: String| DataError::Malformed { + path: path.clone(), + message, + }; + + if parent_cell.len() != n || length_m.len() != n || slope.len() != n { + return Err(malformed(format!( + "order/parent_cell/length_m/slope lengths disagree: {n} / {} / {} / {}", + parent_cell.len(), + length_m.len(), + slope.len() + ))); + } + if rows.len() != cols.len() { + return Err(malformed(format!( + "indices_0/indices_1 lengths disagree: {} / {}", + rows.len(), + cols.len() + ))); + } + // Routing invariant 3 (forward substitution): every edge must run from + // a lower row (upstream) to a strictly higher row (downstream). A + // self-edge would also break `subdivide`'s outlet→inlet argument. + for (k, (&r, &c)) in rows.iter().zip(cols.iter()).enumerate() { + if r < 0 || c < 0 || r as usize >= n || c as usize >= n { + return Err(malformed(format!( + "edge {k} ({r} <- {c}) is outside 0..{n}" + ))); + } + if r <= c { + return Err(malformed(format!( + "edge {k} ({r} <- {c}) is not strictly lower-triangular; the \ + sub-reach store must be topologically ordered (downstream row > \ + upstream row)" + ))); + } + } + for (i, (&l, &s)) in length_m.iter().zip(slope.iter()).enumerate() { + if !(l.is_finite() && l > 0.0) || !s.is_finite() { + return Err(malformed(format!( + "row {i} (node {}) has length_m {l} / slope {s}; both must be finite \ + and length positive", + node_ids[i] + ))); + } + } + + // Each cell's pieces must be one contiguous run: ddrs reads a gauge at + // `parent_offset[p+1] - 1`, and `pieces_per_row_divisor` counts a + // parent's rows from the offsets. A cell that appears in two runs + // would be split into two "parents" with the same id, and the second + // would shadow the first in `IdIndex`. + let mut parent_offset = vec![0_i32]; + let mut seen = std::collections::HashSet::with_capacity(n); + for i in 1..=n { + if i == n || parent_cell[i] != parent_cell[i - 1] { + let cell = parent_cell[i - 1]; + if !seen.insert(cell) { + return Err(malformed(format!( + "cell {cell} occupies more than one contiguous run of rows; the \ + pieces of a cell must be adjacent in `order` (upstream first)" + ))); + } + parent_offset.push(i as i32); + } + } + + Ok(Self { + path, + node_ids, + parent_cell, + rows, + cols, + length_m, + slope, + parent_offset, + }) + } + + /// Number of sub-reach rows. + pub fn n(&self) -> usize { + self.node_ids.len() + } + + /// Number of distinct DDM30 cells (parents). + pub fn n_cells(&self) -> usize { + self.parent_offset.len() - 1 + } + + /// Relabel into ddrs's subdivided layout (see the module diagram). Pure + /// bookkeeping: no edge or length is recomputed. + pub fn to_subdivided(&self) -> SubdividedAdjacency { + let parent_order: Vec = self.parent_offset[..self.n_cells()] + .iter() + .map(|&start| self.parent_cell[start as usize]) + .collect(); + SubdividedAdjacency { + order: self.parent_cell.clone(), + parent_order, + parent_offset: self.parent_offset.clone(), + rows: self.rows.clone(), + cols: self.cols.clone(), + length_m: self.length_m.clone(), + slope: self.slope.clone(), + } + } +} + +// ---------- zarr helpers (mirror src/data/store/zarr.rs's private ones) ---------- + +fn read_i32(storage: &ReadableStorage, store_path: &Path, array_path: &str) -> Result> { + let arr = ZarrArray::open(storage.clone(), array_path).map_err(|e| zarr_err(store_path, e))?; + arr.retrieve_array_subset::>(&arr.subset_all()) + .map_err(|e| zarr_err(store_path, e)) +} + +fn read_f32(storage: &ReadableStorage, store_path: &Path, array_path: &str) -> Result> { + let arr = ZarrArray::open(storage.clone(), array_path).map_err(|e| zarr_err(store_path, e))?; + arr.retrieve_array_subset::>(&arr.subset_all()) + .map_err(|e| zarr_err(store_path, e)) +} + +/// DDR writes `order` as int64 (`build_subdivided_adjacency.py`); accept int32 +/// too so a hand-rolled store still opens. +fn read_i64_or_i32(storage: &ReadableStorage, store_path: &Path, array_path: &str) -> Result> { + let arr = ZarrArray::open(storage.clone(), array_path).map_err(|e| zarr_err(store_path, e))?; + match arr.retrieve_array_subset::>(&arr.subset_all()) { + Ok(v) => Ok(v), + Err(_) => arr + .retrieve_array_subset::>(&arr.subset_all()) + .map(|v| v.into_iter().map(i64::from).collect()) + .map_err(|e| zarr_err(store_path, e)), + } +} + +fn zarr_err(path: &Path, source: E) -> DataError { + DataError::Zarr { + path: path.to_path_buf(), + source: Box::new(source), + } +} + +#[cfg(test)] +mod tests { + use super::*; + use zarrs::array::{data_type, ArrayBuilder}; + use zarrs::group::GroupBuilder; + + /// Two cells, 100 (3 pieces) → 200 (2 pieces): a 5-row chain. + /// + /// ```text + /// 100·0 → 100·1 → 100·2 → 200·0 → 200·1 + /// row 0 1 2 3 4 + /// ``` + struct Synthetic { + order: Vec, + parent: Vec, + rows: Vec, + cols: Vec, + length: Vec, + slope: Vec, + } + + fn chain() -> Synthetic { + Synthetic { + order: vec![100_000, 100_001, 100_002, 200_000, 200_001], + parent: vec![100, 100, 100, 200, 200], + rows: vec![1, 2, 3, 4], + cols: vec![0, 1, 2, 3], + length: vec![3000.0; 5], + slope: vec![0.001, 0.001, 0.001, 0.002, 0.002], + } + } + + fn write_store(dir: &Path, s: &Synthetic, with_parent: bool) { + let storage = Arc::new(FilesystemStore::new(dir).unwrap()); + GroupBuilder::new().build(storage.clone(), "/").unwrap().store_metadata().unwrap(); + macro_rules! put { + ($name:expr, $dt:expr, $fill:expr, $data:expr) => {{ + let data = $data; + let a = ArrayBuilder::new(vec![data.len() as u64], vec![data.len().max(1) as u64], $dt, $fill) + .build(storage.clone(), $name) + .unwrap(); + a.store_metadata().unwrap(); + if !data.is_empty() { + a.store_chunk(&[0], &data).unwrap(); + } + }}; + } + put!("/order", data_type::int64(), 0_i64, s.order.clone()); + if with_parent { + put!("/parent_cell", data_type::int32(), 0_i32, s.parent.clone()); + } + put!("/indices_0", data_type::int32(), 0_i32, s.rows.clone()); + put!("/indices_1", data_type::int32(), 0_i32, s.cols.clone()); + put!("/length_m", data_type::float32(), 0.0_f32, s.length.clone()); + put!("/slope", data_type::float32(), 0.0_f32, s.slope.clone()); + } + + fn open_synthetic(tag: &str, s: &Synthetic, with_parent: bool) -> Result { + let dir = tempfile::tempdir().unwrap(); + let store = dir.path().join(format!("{tag}_subreach_adjacency.zarr")); + write_store(&store, s, with_parent); + GriddedNetwork::open(&store) + } + + #[test] + fn chain_relabels_into_two_parents_with_contiguous_pieces() { + let g = open_synthetic("chain", &chain(), true).expect("open"); + assert_eq!(g.n(), 5); + assert_eq!(g.n_cells(), 2); + let sub = g.to_subdivided(); + assert_eq!(sub.order, vec![100, 100, 100, 200, 200]); + assert_eq!(sub.parent_order, vec![100, 200]); + assert_eq!(sub.parent_offset, vec![0, 3, 5]); + assert_eq!(sub.pieces(0), 3); + assert_eq!(sub.pieces(1), 2); + // The gauge row for a cell is its LAST piece. + assert_eq!(sub.outlet(0), 2); + assert_eq!(sub.outlet(1), 4); + assert_eq!(sub.rows, vec![1, 2, 3, 4]); + assert_eq!(sub.cols, vec![0, 1, 2, 3]); + assert_eq!(sub.length_m, vec![3000.0; 5]); + assert_eq!(sub.slope[3], 0.002); + } + + #[test] + fn store_without_parent_cell_is_refused_with_the_fix() { + let err = open_synthetic("noparent", &chain(), false).unwrap_err(); + let msg = err.to_string(); + assert!(msg.contains("parent_cell"), "{msg}"); + assert!(msg.contains("build_subdivided_adjacency.py"), "{msg}"); + } + + #[test] + fn upper_triangular_edge_is_refused() { + let mut s = chain(); + s.rows[0] = 0; + s.cols[0] = 1; // 0 <- 1: downstream row below upstream row + let msg = open_synthetic("upper", &s, true).unwrap_err().to_string(); + assert!(msg.contains("lower-triangular"), "{msg}"); + } + + #[test] + fn non_contiguous_cell_is_refused() { + let mut s = chain(); + s.parent = vec![100, 200, 100, 200, 200]; // cell 100 split in two runs + let msg = open_synthetic("split", &s, true).unwrap_err().to_string(); + assert!(msg.contains("contiguous"), "{msg}"); + } + + #[test] + fn non_positive_length_is_refused() { + let mut s = chain(); + s.length[2] = 0.0; + let msg = open_synthetic("zerolen", &s, true).unwrap_err().to_string(); + assert!(msg.contains("length_m 0"), "{msg}"); + } +} diff --git a/src/adjacency/mod.rs b/src/adjacency/mod.rs index 1dc92273..0f73d0b4 100644 --- a/src/adjacency/mod.rs +++ b/src/adjacency/mod.rs @@ -35,6 +35,7 @@ pub mod cache; pub mod dbf; pub mod fabric; pub mod gauges; +pub mod gridded; pub mod gpkg; pub mod subdivide; pub mod validate; diff --git a/src/adjacency/zarr_write.rs b/src/adjacency/zarr_write.rs index b7c51a28..d70d7919 100644 --- a/src/adjacency/zarr_write.rs +++ b/src/adjacency/zarr_write.rs @@ -76,7 +76,7 @@ pub fn write_conus_store(adj: &ConusAdjacency, dest: &Path) -> Result<()> { // Root group attrs — mirrors coo_to_zarr (zarr_io.py:131-138) and the real // store's root zarr.json. - let root_attrs = coo_root_attrs(n); + let root_attrs = coo_root_attrs(n, "merit"); let root = GroupBuilder::new() .attributes(root_attrs) .build(storage.clone(), "/") @@ -107,13 +107,24 @@ pub fn write_conus_store(adj: &ConusAdjacency, dest: &Path) -> Result<()> { /// [`write_conus_store`] plus the two identity arrays, which the reader would /// have synthesized anyway — so the disabled path stays a true no-op. pub fn write_conus_store_subdivided(sub: &SubdividedAdjacency, dest: &Path) -> Result<()> { + write_conus_store_subdivided_as(sub, dest, "merit") +} + +/// [`write_conus_store_subdivided`] with an explicit `geodataset` root attr +/// (`"merit"` for the fabric builder, `"ddm30"` for `gridded_network`). The +/// attr is provenance only — no reader dispatches on it. +pub fn write_conus_store_subdivided_as( + sub: &SubdividedAdjacency, + dest: &Path, + geodataset: &str, +) -> Result<()> { let storage = Arc::new(FilesystemStore::new(dest).map_err(|e| zarr_err(dest, e))?); let n_sub = sub.order.len(); let nnz = sub.rows.len(); let root = GroupBuilder::new() - .attributes(coo_root_attrs(n_sub)) + .attributes(coo_root_attrs(n_sub, geodataset)) .build(storage.clone(), "/") .map_err(|e| zarr_err(dest, e))?; root.store_metadata().map_err(|e| zarr_err(dest, e))?; @@ -141,6 +152,16 @@ pub fn write_conus_store_subdivided(sub: &SubdividedAdjacency, dest: &Path) -> R /// /// Round-trips through [`crate::data::store::zarr::GagesAdjacencyStore::open`]. pub fn write_gauges_store(subgraphs: &[GaugeSubgraph], conus_n: usize, dest: &Path) -> Result<()> { + write_gauges_store_as(subgraphs, conus_n, dest, "merit") +} + +/// [`write_gauges_store`] with an explicit `geodataset` attr on every subgroup. +pub fn write_gauges_store_as( + subgraphs: &[GaugeSubgraph], + conus_n: usize, + dest: &Path, + geodataset: &str, +) -> Result<()> { let storage = Arc::new(FilesystemStore::new(dest).map_err(|e| zarr_err(dest, e))?); // Empty root group (matches the real merit_gages store's root zarr.json). @@ -152,7 +173,7 @@ pub fn write_gauges_store(subgraphs: &[GaugeSubgraph], conus_n: usize, dest: &Pa for sg in subgraphs { let group_path = format!("/{}", sg.staid.as_str()); - let mut attrs = coo_root_attrs(conus_n); + let mut attrs = coo_root_attrs(conus_n, geodataset); attrs.insert( "gage_catchment".into(), serde_json::Value::Number(sg.gage_catchment.into()), @@ -181,7 +202,7 @@ pub fn write_gauges_store(subgraphs: &[GaugeSubgraph], conus_n: usize, dest: &Pa /// The COO `format`/`shape`/`geodataset`/`data_types` attr block shared by the /// CONUS root and every gauge subgroup (zarr_io.py:131-138, 382-392). -fn coo_root_attrs(dim: usize) -> serde_json::Map { +fn coo_root_attrs(dim: usize, geodataset: &str) -> serde_json::Map { let mut m = serde_json::Map::new(); m.insert("format".into(), serde_json::Value::String("COO".into())); m.insert( @@ -193,7 +214,7 @@ fn coo_root_attrs(dim: usize) -> serde_json::Map { ); m.insert( "geodataset".into(), - serde_json::Value::String("merit".into()), + serde_json::Value::String(geodataset.into()), ); let mut data_types = serde_json::Map::new(); data_types.insert("indices_0".into(), serde_json::Value::String("int32".into())); diff --git a/src/baseline/cache.rs b/src/baseline/cache.rs index a75400cf..83eb20f5 100644 --- a/src/baseline/cache.rs +++ b/src/baseline/cache.rs @@ -369,6 +369,7 @@ mod tests { gages_adjacency: Some(PathBuf::from("/dev/null/gages_adj.zarr")), geospatial_fabric: None, geospatial_fabric_layer: None, + gridded_network: None, streamflow: PathBuf::from("/dev/null/sf.ic"), observations: PathBuf::from("/dev/null/obs.ic"), gages: PathBuf::from("/dev/null/gages.csv"), @@ -496,6 +497,7 @@ mod tests { gages_adjacency: None, geospatial_fabric: Some(PathBuf::from("/dev/null/fabric.shp")), geospatial_fabric_layer: None, + gridded_network: None, streamflow: PathBuf::from("/dev/null/sf.ic"), observations: PathBuf::from("/dev/null/obs.ic"), gages: PathBuf::from("/dev/null/gages.csv"), diff --git a/src/cli/plan.rs b/src/cli/plan.rs index 6f5754a6..9fefccbf 100644 --- a/src/cli/plan.rs +++ b/src/cli/plan.rs @@ -205,6 +205,9 @@ pub fn plan(input: PlanInput, workspace: &Workspace) -> Result { + let gridded = ds.gridded_network.as_ref().expect("checked by guard"); + let outcome = crate::adjacency::cache::resolve_or_build_gridded( + workspace.root(), + gridded, + &ds.gages, + ) + .map_err(|e| CliError::ConfigInvalid { + path: config_path.into(), + source: Box::new(e), + })?; + Ok(ResolvedAdjacency { + conus: outcome.paths.conus, + gages: outcome.paths.gages, + cache_key: Some(outcome.key), + cache_hit: Some(outcome.cache_hit), + }) + } // Fabric-only (managed build). `validate_data_sources` (config.rs) has // already rejected the partial-adjacency and neither-source cases at // load time, so a missing fabric here is an internal invariant break. _ => { let fabric = ds.geospatial_fabric.as_ref().ok_or_else(|| CliError::ConfigInvalid { path: config_path.into(), - source: "data_sources: no adjacency paths and no geospatial_fabric".into(), + source: "data_sources: no adjacency paths and no geospatial_fabric \ + or gridded_network" + .into(), })?; let outcome = crate::adjacency::cache::resolve_or_build( workspace.root(), diff --git a/src/cli/show.rs b/src/cli/show.rs index 40bf1e79..ed97bcd1 100644 --- a/src/cli/show.rs +++ b/src/cli/show.rs @@ -1,8 +1,11 @@ +use std::path::Path; + use crate::cli::{manifest::Manifest, workspace::Workspace}; use crate::error::CliError; pub fn run_show(ws: &Workspace, run_id: &str, as_json: bool) -> Result<(), CliError> { - let path = ws.runs_dir().join(run_id).join("manifest.json"); + let run_dir = ws.runs_dir().join(run_id); + let path = run_dir.join("manifest.json"); let m = Manifest::read(&path)?; if as_json { println!("{}", serde_json::to_string_pretty(&m) @@ -28,6 +31,123 @@ pub fn run_show(ws: &Workspace, run_id: &str, as_json: bool) -> Result<(), CliEr } } if let Some(p) = &m.outputs.plot { println!("plot {}", p.display()); } + print_metrics(&m.metrics, &run_dir); } Ok(()) } + +/// Render the manifest's `metrics` map, then the summed-Q' baseline comparison +/// when the run carries one. +/// +/// Text mode used to print nothing from `metrics`, so `ddrs show | grep nse` +/// came back empty on a perfectly good `train-and-test` run and read like a +/// failure. The numbers were only ever in `--json`. +fn print_metrics(metrics: &serde_json::Value, run_dir: &Path) { + let Some(map) = metrics.as_object() else { return }; + if map.is_empty() { + return; + } + println!("metrics"); + // Sorted for a stable, greppable order (serde_json preserves insertion + // order only with the preserve_order feature). + let mut keys: Vec<&String> = map.keys().collect(); + keys.sort(); + for k in keys { + println!(" {:<22} {}", k, render_scalar(&map[k])); + } + print_baseline_comparison(map, run_dir); +} + +fn render_scalar(v: &serde_json::Value) -> String { + match v { + serde_json::Value::Number(n) => match n.as_f64() { + // Integers (epoch counts, gauge counts) shouldn't gain decimals. + Some(f) if f.fract() == 0.0 && f.abs() < 1e15 => format!("{}", f as i64), + Some(f) => format!("{f:.4}"), + None => n.to_string(), + }, + serde_json::Value::String(s) => s.clone(), + other => other.to_string(), + } +} + +/// `train-and-test` copies the summed-Q' baseline to +/// `/baseline/manifest.json`, but its metrics are **per-gauge arrays**, +/// not reduced medians, so a reader comparing them to the run's +/// `median_nse_finite` was comparing a scalar against a list. Reduce here. +fn print_baseline_comparison( + run_metrics: &serde_json::Map, + run_dir: &Path, +) { + let text = match std::fs::read_to_string(run_dir.join("baseline/manifest.json")) { + Ok(t) => t, + Err(_) => return, // train-only run, or baseline copy unavailable + }; + let Ok(b): std::result::Result = serde_json::from_str(&text) else { + return; + }; + let bm = &b["metrics"]; + println!("baseline (summed Q', median over gauges)"); + for (label, key, run_key) in [ + ("NSE", "nse", "median_nse_finite"), + ("KGE", "kge", "median_kge_finite"), + ] { + let Some(base) = median_finite(&bm[key]) else { continue }; + match run_metrics.get(run_key).and_then(|v| v.as_f64()) { + Some(routed) => println!( + " {:<22} {:.4} routed {:.4} Δ {:+.4}{}", + label, + base, + routed, + routed - base, + if routed > base { "" } else { " ← routed does NOT beat baseline" }, + ), + None => println!(" {:<22} {:.4}", label, base), + } + } +} + +/// Median of the finite entries of a JSON array, or `None` if there are none. +fn median_finite(v: &serde_json::Value) -> Option { + let mut xs: Vec = v + .as_array()? + .iter() + .filter_map(|x| x.as_f64()) + .filter(|x| x.is_finite()) + .collect(); + if xs.is_empty() { + return None; + } + xs.sort_by(|a, b| a.partial_cmp(b).expect("finite values compare")); + let n = xs.len(); + Some(if n % 2 == 1 { + xs[n / 2] + } else { + (xs[n / 2 - 1] + xs[n / 2]) / 2.0 + }) +} + +#[cfg(test)] +mod tests { + use super::*; + use serde_json::json; + + #[test] + fn median_finite_handles_odd_even_and_nonfinite() { + assert_eq!(median_finite(&json!([3.0, 1.0, 2.0])), Some(2.0)); + assert_eq!(median_finite(&json!([1.0, 2.0, 3.0, 4.0])), Some(2.5)); + // nulls (serde's NaN encoding here) and missing entries are skipped. + assert_eq!(median_finite(&json!([1.0, null, 3.0])), Some(2.0)); + assert_eq!(median_finite(&json!([])), None); + assert_eq!(median_finite(&json!([null])), None); + assert_eq!(median_finite(&json!("not an array")), None); + } + + #[test] + fn render_scalar_keeps_counts_integral() { + assert_eq!(render_scalar(&json!(30)), "30"); + assert_eq!(render_scalar(&json!(620.0)), "620"); + assert_eq!(render_scalar(&json!(0.7510778)), "0.7511"); + assert_eq!(render_scalar(&json!("ok")), "ok"); + } +} diff --git a/src/cli/sources.rs b/src/cli/sources.rs index 48a4ae9e..90938edd 100644 --- a/src/cli/sources.rs +++ b/src/cli/sources.rs @@ -239,6 +239,9 @@ fn lock_sources_from_config(cfg_path: &Path, ws: &Workspace) -> Result<(), CliEr if let Some(p) = &ds.geospatial_fabric { pairs.push(("geospatial_fabric".into(), p.clone())); } + if let Some(p) = &ds.gridded_network { + pairs.push(("gridded_network".into(), p.clone())); + } // Parallel reachability + fingerprint. std::thread::scope is fine — these // are I/O-bound and the count is small. diff --git a/src/config.rs b/src/config.rs index f9225e2c..815028a9 100644 --- a/src/config.rs +++ b/src/config.rs @@ -94,6 +94,7 @@ where } #[derive(Debug, Clone, Deserialize)] +#[serde(deny_unknown_fields)] pub struct DataSources { #[serde(deserialize_with = "deserialize_one_or_many_paths")] pub attributes: Vec, @@ -127,9 +128,20 @@ pub struct DataSources { /// `.shp`/`.dbf` fabrics. #[serde(default)] pub geospatial_fabric_layer: Option, + /// Path to a DDR gridded (ISIMIP DDM30) **sub-reach** adjacency zarr — + /// the output of DDR's `scripts/build_subdivided_adjacency.py` (`order`, + /// `parent_cell`, `indices_0/1`, `length_m`, `slope`). A third, + /// mutually-exclusive adjacency source next to `geospatial_fabric` and + /// the explicit zarr pair: `ddrs plan` turns it into ddrs's subdivided + /// store layout (parent = grid cell, pieces = sub-reaches) plus per-gauge + /// subgraphs under `.ddrs/adjacency//`. See + /// `docs/superpowers/specs/2026-09-08-ddrs-gridded-routing-design.md`. + #[serde(default)] + pub gridded_network: Option, } #[derive(Debug, Clone, Deserialize)] +#[serde(deny_unknown_fields)] pub struct Experiment { pub batch_size: usize, pub start_time: String, @@ -218,6 +230,7 @@ pub enum OptimizerKind { /// component weights. See `src/training/loss.rs`. #[derive(Debug, Clone, Deserialize)] #[serde(default, rename_all = "kebab-case")] +#[serde(deny_unknown_fields)] pub struct LossConfig { /// Which objective to optimize. pub kind: LossKind, @@ -289,6 +302,7 @@ pub enum LossKind { } #[derive(Debug, Clone, Deserialize)] +#[serde(deny_unknown_fields)] pub struct KanHeadConfigSection { pub hidden_size: usize, pub num_hidden_layers: usize, @@ -677,6 +691,7 @@ impl Default for Params { #[derive(Debug, Default, Deserialize)] #[serde(default)] +#[serde(deny_unknown_fields)] struct ParamsRaw { parameter_ranges: HashMap, attribute_minimums: HashMap, @@ -813,6 +828,7 @@ pub struct Config { /// Fields are all optional so absent keys inherit from `experiment:`. #[derive(Debug, Default, Deserialize, Clone)] #[serde(default)] +#[serde(deny_unknown_fields)] struct TestingOverridesRaw { pub start_time: Option, pub end_time: Option, @@ -839,6 +855,7 @@ where /// YAML-shaped intermediate; the public `Config` has nicer types. #[derive(Debug, Default, Deserialize)] #[serde(default)] +#[serde(deny_unknown_fields)] struct ConfigRaw { mode: Option, geodataset: Option, @@ -858,13 +875,25 @@ struct ConfigRaw { impl From for Config { fn from(r: ConfigRaw) -> Self { + // Resolved before the struct literal moves `data_sources` out of `r`. + let geodataset = r.geodataset.clone().unwrap_or_else(|| { + match r.data_sources.as_ref() { + Some(ds) if ds.gridded_network.is_some() => "ddm30".to_string(), + _ => "merit".to_string(), + } + }); Self { params: r.params.into(), data_sources: r.data_sources, experiment: r.experiment, kan_head: r.kan_head, mode: r.mode.unwrap_or_else(|| "training".to_string()), - geodataset: r.geodataset.unwrap_or_else(|| "merit".to_string()), + // Absent `geodataset:` is INFERRED from the adjacency source rather + // than defaulting blindly to "merit", so a gridded config that + // omits the key is labelled `ddm30` instead of being mislabelled. + // An explicitly wrong value is then a real disagreement, which + // `validate_geodataset` rejects. + geodataset, seed: r.seed.unwrap_or(42), np_seed: r.np_seed.unwrap_or(42), workflow: r.workflow, @@ -915,6 +944,10 @@ impl Config { path: path.to_path_buf(), source: serde_yaml::Error::custom(msg), })?; + validate_geodataset(&cfg).map_err(|msg| DataError::Yaml { + path: path.to_path_buf(), + source: serde_yaml::Error::custom(msg), + })?; validate_leakance(&cfg).map_err(|msg| DataError::Yaml { path: path.to_path_buf(), source: serde_yaml::Error::custom(msg), @@ -963,6 +996,26 @@ fn validate_data_sources(cfg: &Config) -> std::result::Result<(), String> { let has_conus = ds.conus_adjacency.is_some(); let has_gages = ds.gages_adjacency.is_some(); let has_fabric = ds.geospatial_fabric.is_some(); + let has_gridded = ds.gridded_network.is_some(); + + // The gridded network is its own managed-build input; combining it with + // either other adjacency source would leave one of them silently unused. + if has_gridded && has_fabric { + return Err( + "data_sources: `gridded_network` and `geospatial_fabric` are both set; \ + they are alternative adjacency sources — keep exactly one." + .to_string(), + ); + } + if has_gridded && (has_conus || has_gages) { + return Err( + "data_sources: `gridded_network` is set alongside explicit \ + `conus_adjacency`/`gages_adjacency` paths; the gridded store is built \ + into the managed cache by `ddrs plan` — remove the explicit paths \ + (or drop `gridded_network` to use them as-is)." + .to_string(), + ); + } // Layer selection is a gpkg-only concept; reject it for dBASE fabrics // (or with no fabric at all) at load time rather than silently ignoring. @@ -982,10 +1035,10 @@ fn validate_data_sources(cfg: &Config) -> std::result::Result<(), String> { } } - match (has_conus, has_gages, has_fabric) { + match (has_conus, has_gages, has_fabric || has_gridded) { // Both explicit zarr paths: valid. (true, true, _) => Ok(()), - // Managed build: neither zarr, fabric present: valid. + // Managed build: neither zarr, fabric or gridded network present: valid. (false, false, true) => Ok(()), // Partial adjacency: exactly one of the two zarr keys provided. (true, false, _) => Err( @@ -998,16 +1051,59 @@ fn validate_data_sources(cfg: &Config) -> std::result::Result<(), String> { provide both adjacency paths or remove both and set `geospatial_fabric`." .to_string(), ), - // Neither zarr and no fabric. + // Neither zarr and no fabric / gridded network. (false, false, false) => Err( "data_sources: adjacency sources are missing — either set both \ `conus_adjacency` and `gages_adjacency`, or set `geospatial_fabric` \ - for a managed adjacency build." + (MERIT) or `gridded_network` (DDM30) for a managed adjacency build." .to_string(), ), } } +/// `geodataset:` must agree with the adjacency source it is labelling. +/// +/// Nothing in the routing path reads `geodataset` — it is provenance, written +/// into the managed-build zarr attrs and copied into the run's `config.yaml` +/// snapshot. That is exactly why a wrong value is dangerous: the run succeeds +/// and its audit trail names the wrong network. `ddrs sources use` replaces the +/// `data_sources:` block only, so switching MERIT ↔ gridded leaves a stale label +/// behind unless the user remembers to edit it. This turns that silent +/// mislabelling into a load-time error. +/// +/// Only the two managed-build sources are decidable. Explicit +/// `conus_adjacency`/`gages_adjacency` paths can point at either kind of +/// pre-built store, so any label is accepted there. +fn validate_geodataset(cfg: &Config) -> std::result::Result<(), String> { + let Some(ds) = cfg.data_sources.as_ref() else { + return Ok(()); + }; + let label = cfg.geodataset.trim(); + let expected = if ds.gridded_network.is_some() { + "ddm30" + } else if ds.geospatial_fabric.is_some() { + "merit" + } else { + return Ok(()); // explicit adjacency paths: undecidable, allow any label + }; + if label.eq_ignore_ascii_case(expected) { + return Ok(()); + } + let (source_key, other) = if expected == "ddm30" { + ("gridded_network", "merit") + } else { + ("geospatial_fabric", "ddm30") + }; + Err(format!( + "geodataset: `{label}` contradicts `data_sources.{source_key}` — that source \ + builds a `{expected}` network. Nothing reads `geodataset` at runtime, so this \ + would not fail the run; it would only mislabel the run's config snapshot and \ + the adjacency store attrs. Fix: set `geodataset: {expected}` (or omit the key \ + and it is inferred). If you meant to route the {other} network, change \ + `data_sources` instead — `ddrs sources use ` swaps the block for you." + )) +} + fn validate_leakance(cfg: &Config) -> std::result::Result<(), String> { if cfg.params.use_leakance && cfg.params.use_cuda_graphs { return Err( @@ -1068,6 +1164,22 @@ fn validate_subdivision(cfg: &Config) -> std::result::Result<(), String> { ); } if s.enabled { + if cfg + .data_sources + .as_ref() + .is_some_and(|ds| ds.gridded_network.is_some()) + { + // The DDM30 sub-reach store arrives already split (DDR sizes its + // pieces from celerity at mean flow); ddrs's MERIT reach plan has + // no `catchsize`/fabric to work from and would be silently inert. + return Err( + "params.subdivision: `enabled: true` does not apply to \ + `data_sources.gridded_network` — the DDM30 sub-reach store is \ + already subdivided by DDR's build_subdivided_adjacency.py. \ + Set `enabled: false` (or resize the sub-reaches in DDR)." + .to_string(), + ); + } validate_subdivision_reaches_the_builder(cfg)?; } Ok(()) @@ -1521,9 +1633,12 @@ workflow: train // ── data_sources validation matrix ────────────────────────────────────── + /// `geodataset:` is deliberately omitted so it is inferred from the block + /// under test: stamping `merit` here would contradict any gridded + /// `data_sources` and trip `validate_geodataset`. fn write_yaml_with_data_sources(name: &str, ds_block: &str) -> std::path::PathBuf { let yaml = format!( - "mode: training\ngeodataset: merit\nseed: 1\nnp_seed: 1\n{ds_block}" + "mode: training\nseed: 1\nnp_seed: 1\n{ds_block}" ); let path = std::env::temp_dir().join(name); std::fs::write(&path, yaml).unwrap(); @@ -1588,6 +1703,224 @@ data_sources: ); } + #[test] + fn gridded_network_alone_is_a_valid_managed_build() { + let ds_block = r#" +data_sources: + attributes: /dev/null/attrs.nc + gridded_network: /dev/null/ddm30_subreach_adjacency.zarr + streamflow: /dev/null/sf.ic + observations: /dev/null/obs.ic + gages: /dev/null/gages.csv +"#; + let path = write_yaml_with_data_sources("ddrs_ds_gridded_only.yaml", ds_block); + let cfg = Config::from_yaml_file(&path).expect("gridded_network alone must be valid"); + let ds = cfg.data_sources.as_ref().unwrap(); + assert!(ds.gridded_network.is_some()); + assert!(ds.conus_adjacency.is_none() && ds.gages_adjacency.is_none()); + } + + /// A misspelled key used to take its default in silence, which the docs + /// call the single most common cause of "my config change did nothing". + /// Every config section now sets `deny_unknown_fields`, so it is a load + /// error naming the key. + #[test] + fn typod_keys_are_rejected_in_every_section() { + let base = |extra: &str| { + format!( + "mode: training\ngeodataset: merit\nseed: 1\nnp_seed: 1\n\ + data_sources:\n attributes: /dev/null/a.nc\n \ + conus_adjacency: /dev/null/c.zarr\n gages_adjacency: /dev/null/g.zarr\n \ + streamflow: /dev/null/s.ic\n observations: /dev/null/o.ic\n \ + gages: /dev/null/g.csv\n{extra}" + ) + }; + // (label, yaml fragment, the misspelled key that must be named) + let cases = [ + ("top level", "epocs: 5\n", "epocs"), + ("experiment", "experiment:\n epocs: 5\n", "epocs"), + ("kan_head", "kan_head:\n hidden_sizes: 21\n", "hidden_sizes"), + ("params", "params:\n taus: 9\n", "taus"), + ("testing", "testing:\n batchsize: 15\n", "batchsize"), + ]; + for (label, extra, key) in cases { + let path = std::env::temp_dir() + .join(format!("ddrs_typo_{}.yaml", label.replace(' ', "_"))); + std::fs::write(&path, base(extra)).unwrap(); + let err = Config::from_yaml_file(&path) + .expect_err(&format!("{label}: typo `{key}` must be rejected")); + let msg = format!("{err}"); + assert!( + msg.contains(key), + "{label}: error should name the offending key `{key}`, got: {msg}" + ); + } + } + + /// The guard must not reject keys that really exist, including the + /// `mlp:` backward-compat alias for `kan_head:`. + #[test] + fn known_keys_and_the_mlp_alias_still_parse() { + let yaml = "mode: training\ngeodataset: merit\nseed: 1\nnp_seed: 1\n\ + data_sources:\n attributes: /dev/null/a.nc\n \ + conus_adjacency: /dev/null/c.zarr\n gages_adjacency: /dev/null/g.zarr\n \ + streamflow: /dev/null/s.ic\n observations: /dev/null/o.ic\n \ + gages: /dev/null/g.csv\n\ + mlp:\n hidden_size: 21\n num_hidden_layers: 2\n grid: 50\n k: 2\n \ + input_var_names: [meanP]\n learnable_parameters: [n]\n"; + let path = std::env::temp_dir().join("ddrs_mlp_alias.yaml"); + std::fs::write(&path, yaml).unwrap(); + let cfg = Config::from_yaml_file(&path).expect("mlp alias must still parse"); + assert_eq!(cfg.kan_head.as_ref().unwrap().hidden_size, 21); + } + + #[test] + fn absent_geodataset_is_inferred_from_the_adjacency_source() { + // Gridded source, no `geodataset:` key → inferred `ddm30`, not the + // blind "merit" default that used to mislabel every gridded run. + let ds_block = r#" +data_sources: + attributes: /dev/null/attrs.nc + gridded_network: /dev/null/ddm30_subreach_adjacency.zarr + streamflow: /dev/null/sf.ic + observations: /dev/null/obs.ic + gages: /dev/null/gages.csv +"#; + let yaml = format!("mode: training\nseed: 1\nnp_seed: 1\n{ds_block}"); + let path = std::env::temp_dir().join("ddrs_geodataset_infer.yaml"); + std::fs::write(&path, yaml).unwrap(); + let cfg = Config::from_yaml_file(&path).expect("infers ddm30"); + assert_eq!(cfg.geodataset, "ddm30"); + + // Fabric source, no key → merit. + let ds_block = r#" +data_sources: + attributes: /dev/null/attrs.nc + geospatial_fabric: /dev/null/rivers.shp + streamflow: /dev/null/sf.ic + observations: /dev/null/obs.ic + gages: /dev/null/gages.csv +"#; + let yaml = format!("mode: training\nseed: 1\nnp_seed: 1\n{ds_block}"); + let path = std::env::temp_dir().join("ddrs_geodataset_infer_merit.yaml"); + std::fs::write(&path, yaml).unwrap(); + assert_eq!(Config::from_yaml_file(&path).unwrap().geodataset, "merit"); + } + + #[test] + fn geodataset_contradicting_the_source_is_rejected() { + // The `ddrs sources use conus` footgun: block swapped to MERIT, stale + // `geodataset: ddm30` left behind above it. + let yaml = "mode: training\ngeodataset: ddm30\nseed: 1\nnp_seed: 1\n +data_sources: + attributes: /dev/null/attrs.nc + geospatial_fabric: /dev/null/rivers.shp + streamflow: /dev/null/sf.ic + observations: /dev/null/obs.ic + gages: /dev/null/gages.csv +"; + let path = std::env::temp_dir().join("ddrs_geodataset_stale.yaml"); + std::fs::write(&path, yaml).unwrap(); + let msg = format!("{}", Config::from_yaml_file(&path).unwrap_err()); + assert!( + msg.contains("geodataset") && msg.contains("geospatial_fabric"), + "expected the mislabel error, got: {msg}" + ); + + // And the mirror: gridded source labelled merit. + let yaml = "mode: training\ngeodataset: merit\nseed: 1\nnp_seed: 1\n +data_sources: + attributes: /dev/null/attrs.nc + gridded_network: /dev/null/ddm30.zarr + streamflow: /dev/null/sf.ic + observations: /dev/null/obs.ic + gages: /dev/null/gages.csv +"; + let path = std::env::temp_dir().join("ddrs_geodataset_stale2.yaml"); + std::fs::write(&path, yaml).unwrap(); + let msg = format!("{}", Config::from_yaml_file(&path).unwrap_err()); + assert!(msg.contains("ddm30"), "expected the mislabel error, got: {msg}"); + } + + #[test] + fn geodataset_is_unconstrained_with_explicit_adjacency_paths() { + // A pre-built store can be either network, so the label is undecidable + // and must not be second-guessed. + let ds_block = r#" +data_sources: + attributes: /dev/null/attrs.nc + conus_adjacency: /dev/null/conus.zarr + gages_adjacency: /dev/null/gages.zarr + streamflow: /dev/null/sf.ic + observations: /dev/null/obs.ic + gages: /dev/null/gages.csv +"#; + let path = write_yaml_with_data_sources("ddrs_geodataset_explicit.yaml", ds_block); + // write_yaml_with_data_sources stamps `geodataset: merit`. + assert!(Config::from_yaml_file(&path).is_ok()); + } + + #[test] + fn gridded_network_with_fabric_rejected() { + let ds_block = r#" +data_sources: + attributes: /dev/null/attrs.nc + gridded_network: /dev/null/ddm30_subreach_adjacency.zarr + geospatial_fabric: /dev/null/rivers.shp + streamflow: /dev/null/sf.ic + observations: /dev/null/obs.ic + gages: /dev/null/gages.csv +"#; + let path = write_yaml_with_data_sources("ddrs_ds_gridded_fabric.yaml", ds_block); + let msg = format!("{}", Config::from_yaml_file(&path).unwrap_err()); + assert!( + msg.contains("gridded_network") && msg.contains("geospatial_fabric"), + "expected the both-sources error, got: {msg}" + ); + } + + #[test] + fn gridded_network_with_explicit_adjacency_rejected() { + let ds_block = r#" +data_sources: + attributes: /dev/null/attrs.nc + gridded_network: /dev/null/ddm30_subreach_adjacency.zarr + conus_adjacency: /dev/null/conus.zarr + gages_adjacency: /dev/null/gages.zarr + streamflow: /dev/null/sf.ic + observations: /dev/null/obs.ic + gages: /dev/null/gages.csv +"#; + let path = write_yaml_with_data_sources("ddrs_ds_gridded_explicit.yaml", ds_block); + let msg = format!("{}", Config::from_yaml_file(&path).unwrap_err()); + assert!( + msg.contains("gridded_network") && msg.contains("conus_adjacency"), + "expected the explicit-paths error, got: {msg}" + ); + } + + #[test] + fn gridded_network_with_subdivision_enabled_rejected() { + let block = r#" +data_sources: + attributes: /dev/null/attrs.nc + gridded_network: /dev/null/ddm30_subreach_adjacency.zarr + streamflow: /dev/null/sf.ic + observations: /dev/null/obs.ic + gages: /dev/null/gages.csv +params: + use_cuda_graphs: false + subdivision: + enabled: true +"#; + let path = write_yaml_with_data_sources("ddrs_ds_gridded_subdiv.yaml", block); + let msg = format!("{}", Config::from_yaml_file(&path).unwrap_err()); + assert!( + msg.contains("params.subdivision") && msg.contains("gridded_network"), + "expected the subdivision-vs-gridded error, got: {msg}" + ); + } + #[test] fn partial_adjacency_conus_only_rejected() { // Only conus_adjacency, no gages_adjacency → config error. diff --git a/src/data/store/gage_csv.rs b/src/data/store/gage_csv.rs index b2134755..b5167041 100644 --- a/src/data/store/gage_csv.rs +++ b/src/data/store/gage_csv.rs @@ -51,7 +51,10 @@ struct RawRow { lat_gage: f64, #[serde(rename = "LNG_GAGE")] lng_gage: f64, - #[serde(rename = "COMID", default)] + /// Outlet reach id. DDR's gridded (DDM30) gauge CSVs write the snapped + /// grid cell under `cell` instead of `COMID`; the value plays the same role + /// (the id the adjacency `order` is keyed on), so it is accepted as an alias. + #[serde(rename = "COMID", alias = "cell", default)] comid: Option, #[serde(rename = "COMID_DRAIN_SQKM", default)] comid_drain_sqkm: Option, @@ -230,6 +233,21 @@ STAID,STANAME,DRAIN_SQKM,LAT_GAGE,LNG_GAGE,COMID,COMID_DRAIN_SQKM,COMID_UNITAREA assert_eq!(r1.da_valid, Some(false)); } + /// DDR's gridded gauge CSV (`examples/juniata_gridded/data/juniata_gage.csv` + /// layout): the snapped DDM30 cell is under `cell`, and there is no COMID, + /// DA_VALID or FLOW_SCALE column. + #[test] + fn gridded_csv_cell_column_is_read_as_comid() { + let csv = "STAID,STANAME,DRAIN_SQKM,LAT_GAGE,LNG_GAGE,cell,da_ratio\n\ + 01567000,\"JUNIATA RIVER AT NEWPORT, PA\",8657.0,40.4785,-77.1294,138445,1.088\n"; + let m = GageMetadata::from_reader(csv.as_bytes(), std::path::PathBuf::from("")) + .expect("parse gridded csv"); + let r = &m.rows[0]; + assert_eq!(r.comid, Some(138445)); + assert_eq!(r.da_valid, None); + assert_eq!(r.flow_scale, None); + } + #[test] fn lookup_by_staid_uses_padded_form() { let m = GageMetadata::from_reader( diff --git a/src/data/store/icechunk.rs b/src/data/store/icechunk.rs index 58dbccc0..5806ad45 100644 --- a/src/data/store/icechunk.rs +++ b/src/data/store/icechunk.rs @@ -190,6 +190,53 @@ pub struct StreamflowStore { #[allow(dead_code)] storage: Arc, qr: ZarrArray, + /// `true` when `Qr` is stored `(time, divide_id)` — DDR's gridded (DDM30) + /// stores; `false` for the contract's `(divide_id, time)`. Decided at open + /// from the array's dimension names, falling back to its shape. + pub time_major: bool, +} + +/// Which axis of `Qr` is the divide axis. The contract +/// (`docs/nh-qprime-store-contract.md`) is `(divide_id, time)`, but DDR's +/// `build_gridded_qprime.py` writes `(time, divide_id)`, and zarrs answers an +/// out-of-range subset with fill values rather than an error — so a layout +/// guess that is wrong reads NaN silently. Dimension names decide when present; +/// otherwise the shape does, and a shape matching neither is refused. +fn detect_time_major( + qr: &ZarrArray, + n_divide: usize, + n_time: usize, + path: &Path, +) -> Result { + let shape = qr.shape(); + if shape.len() != 2 { + return Err(DataError::Malformed { + path: path.to_path_buf(), + message: format!("Qr must be 2-D, got shape {shape:?}"), + }); + } + if let Some(names) = qr.dimension_names() { + let named: Vec<&str> = names.iter().map(|n| n.as_deref().unwrap_or("")).collect(); + match named.as_slice() { + ["divide_id", "time"] => return Ok(false), + ["time", "divide_id"] => return Ok(true), + _ => {} // unnamed or unexpected names: fall through to the shape + } + } + let (a, b) = (shape[0] as usize, shape[1] as usize); + if a == n_divide && b == n_time { + Ok(false) + } else if a == n_time && b == n_divide { + Ok(true) + } else { + Err(DataError::Malformed { + path: path.to_path_buf(), + message: format!( + "Qr shape {shape:?} matches neither (divide_id, time) = ({n_divide}, {n_time}) \ + nor (time, divide_id)" + ), + }) + } } impl StreamflowStore { @@ -239,8 +286,9 @@ impl StreamflowStore { // 3. Open the `/Qr` data var (Task 3 will actually read from it). let qr = ZarrArray::open(readable.clone(), "/Qr") .map_err(|e| ic_err(&path, e))?; + let time_major = detect_time_major(&qr, index.len(), n_time, &path)?; - Ok(Self { path, index, time_start, n_time, resolution, storage, qr }) + Ok(Self { path, index, time_start, n_time, resolution, storage, qr, time_major }) } } @@ -508,17 +556,28 @@ impl StreamflowStore { let div_range_end = max_pos + 1; let div_count = div_range_end - min_pos; - // Qr is stored as (divide_id, time). Subset: axis 0 = divide, axis 1 = time. - let subset = zarrs::array::ArraySubset::new_with_ranges(&[ - (min_pos as u64)..(div_range_end as u64), - (start_step as u64)..(end_step as u64), - ]); + // Contract layout is (divide_id, time): axis 0 = divide, axis 1 = time. + // Time-major stores (DDR's gridded Q') swap the axes; the raw buffer is + // then (n_steps, div_count) row-major and is addressed accordingly. + let div_axis = (min_pos as u64)..(div_range_end as u64); + let time_axis = (start_step as u64)..(end_step as u64); + let subset = if self.time_major { + zarrs::array::ArraySubset::new_with_ranges(&[time_axis, div_axis]) + } else { + zarrs::array::ArraySubset::new_with_ranges(&[div_axis, time_axis]) + }; let raw_f32: Vec = self .qr .retrieve_array_subset(&subset) .map_err(|e| ic_err(&self.path, e))?; - // raw_f32 is row-major: shape (div_count, n_steps). debug_assert_eq!(raw_f32.len(), div_count * n_steps); + let raw_idx = |local_div: usize, t: usize| -> usize { + if self.time_major { + t * div_count + local_div + } else { + local_div * n_steps + t + } + }; // Scatter into the output. Walk `comids` in order; for each // non-missing entry consume the next element of `positions`. @@ -531,8 +590,7 @@ impl StreamflowStore { next_present += 1; let local_div = div_pos - min_pos; for t in 0..n_steps { - let raw_idx = local_div * n_steps + t; - out[(t, out_col)] = raw_f32[raw_idx]; + out[(t, out_col)] = raw_f32[raw_idx(local_div, t)]; } } diff --git a/src/data/store/zarr.rs b/src/data/store/zarr.rs index 51f3adf8..d9250e07 100644 --- a/src/data/store/zarr.rs +++ b/src/data/store/zarr.rs @@ -216,13 +216,21 @@ impl GageSubgraph { /// node appears as an edge endpoint. Single-divide catchments have NO /// edges, so this returns empty — callers must filter with /// [`GageSubgraph::is_headwater`] first, as training does. + /// + /// On a **subdivided** store `order` repeats a parent's id once per piece, + /// so the position set is deduplicated by id: each divide (MERIT reach or + /// DDM30 cell) appears once, in the position order of its first piece. + /// The baseline sums Q′ per divide, and Q′ is stored per divide, so a + /// per-piece list would count every divide `m` times. pub fn upstream_comids(&self, conus: &ConusAdjacencyStore) -> Vec { let mut positions: std::collections::BTreeSet = std::collections::BTreeSet::new(); positions.extend(self.indices_0.iter().copied()); positions.extend(self.indices_1.iter().copied()); + let mut seen: std::collections::HashSet = std::collections::HashSet::new(); positions .into_iter() .map(|pos| conus.order[pos as usize]) + .filter(|c| seen.insert(*c)) .collect() } } @@ -376,6 +384,26 @@ mod tests { assert_eq!(comids, vec![Comid(100), Comid(200), Comid(300), Comid(400)]); } + /// Subdivided store: parent 100 owns rows 0..3, parent 200 owns rows 3..5 + /// (`order` repeats the parent id per piece). The upstream set of a gauge + /// on the outlet piece names each parent ONCE — the baseline reads Q′ per + /// parent, and a per-piece list would triple-count parent 100. + #[test] + fn upstream_comids_names_each_subdivided_parent_once() { + let mut conus = fake_conus(vec![100, 100, 100, 200, 200]); + conus.parent_order = vec![Comid(100), Comid(200)]; + conus.parent_offset = vec![0, 3, 5]; + conus.index = IdIndex::new(conus.parent_order.clone()); + let sg = GageSubgraph { + staid: Staid::from("00000004"), + gage_idx: 4, + gage_catchment: "200".into(), + indices_0: vec![1, 2, 3, 4], + indices_1: vec![0, 1, 2, 3], + }; + assert_eq!(sg.upstream_comids(&conus), vec![Comid(100), Comid(200)]); + } + #[test] fn single_divide_subgraph_is_headwater_and_has_empty_upstream() { // Single-divide catchments have NO edges in the gages store — only a diff --git a/tests/fixtures/qr_daily_time_major.ic/manifests/1ZQAR89FZ5DSH6AJMKEG b/tests/fixtures/qr_daily_time_major.ic/manifests/1ZQAR89FZ5DSH6AJMKEG new file mode 100644 index 00000000..ac540958 Binary files /dev/null and b/tests/fixtures/qr_daily_time_major.ic/manifests/1ZQAR89FZ5DSH6AJMKEG differ diff --git a/tests/fixtures/qr_daily_time_major.ic/manifests/FHNTWK6FAXKF41VN8NK0 b/tests/fixtures/qr_daily_time_major.ic/manifests/FHNTWK6FAXKF41VN8NK0 new file mode 100644 index 00000000..ea039567 Binary files /dev/null and b/tests/fixtures/qr_daily_time_major.ic/manifests/FHNTWK6FAXKF41VN8NK0 differ diff --git a/tests/fixtures/qr_daily_time_major.ic/manifests/QFRFNDESAR9XZ41HEKWG b/tests/fixtures/qr_daily_time_major.ic/manifests/QFRFNDESAR9XZ41HEKWG new file mode 100644 index 00000000..a65bb947 Binary files /dev/null and b/tests/fixtures/qr_daily_time_major.ic/manifests/QFRFNDESAR9XZ41HEKWG differ diff --git a/tests/fixtures/qr_daily_time_major.ic/overwritten/repo.30714756847413.HJEX80MSFK7C39YMWW30 b/tests/fixtures/qr_daily_time_major.ic/overwritten/repo.30714756847413.HJEX80MSFK7C39YMWW30 new file mode 100644 index 00000000..e76dfbae Binary files /dev/null and b/tests/fixtures/qr_daily_time_major.ic/overwritten/repo.30714756847413.HJEX80MSFK7C39YMWW30 differ diff --git a/tests/fixtures/qr_daily_time_major.ic/repo b/tests/fixtures/qr_daily_time_major.ic/repo new file mode 100644 index 00000000..fde40110 Binary files /dev/null and b/tests/fixtures/qr_daily_time_major.ic/repo differ diff --git a/tests/fixtures/qr_daily_time_major.ic/snapshots/1CECHNKREP0F1RSTCMT0 b/tests/fixtures/qr_daily_time_major.ic/snapshots/1CECHNKREP0F1RSTCMT0 new file mode 100644 index 00000000..a96b1b85 Binary files /dev/null and b/tests/fixtures/qr_daily_time_major.ic/snapshots/1CECHNKREP0F1RSTCMT0 differ diff --git a/tests/fixtures/qr_daily_time_major.ic/snapshots/HGT03GVPCD6B4ZVTDQ9G b/tests/fixtures/qr_daily_time_major.ic/snapshots/HGT03GVPCD6B4ZVTDQ9G new file mode 100644 index 00000000..f4e8fc51 Binary files /dev/null and b/tests/fixtures/qr_daily_time_major.ic/snapshots/HGT03GVPCD6B4ZVTDQ9G differ diff --git a/tests/fixtures/qr_daily_time_major.ic/transactions/1CECHNKREP0F1RSTCMT0 b/tests/fixtures/qr_daily_time_major.ic/transactions/1CECHNKREP0F1RSTCMT0 new file mode 100644 index 00000000..2a1bf888 Binary files /dev/null and b/tests/fixtures/qr_daily_time_major.ic/transactions/1CECHNKREP0F1RSTCMT0 differ diff --git a/tests/fixtures/qr_daily_time_major.ic/transactions/HGT03GVPCD6B4ZVTDQ9G b/tests/fixtures/qr_daily_time_major.ic/transactions/HGT03GVPCD6B4ZVTDQ9G new file mode 100644 index 00000000..ed3ab6be Binary files /dev/null and b/tests/fixtures/qr_daily_time_major.ic/transactions/HGT03GVPCD6B4ZVTDQ9G differ diff --git a/tests/gridded_acceptance.rs b/tests/gridded_acceptance.rs new file mode 100644 index 00000000..d890ec0b --- /dev/null +++ b/tests/gridded_acceptance.rs @@ -0,0 +1,105 @@ +//! End-to-end routing acceptance on the committed gridded (ISIMIP DDM30) +//! Juniata bundle: full `train-and-test` on CPU through the gridded managed +//! build, then metric floors on the manifest. +//! +//! Reference result (2026-09-09, seed 42, 30 epochs): routed NSE 0.751 / +//! KGE 0.730 vs summed-Q' baseline NSE 0.594 / KGE 0.672. Seeds {42, 7, 123} +//! spanned NSE 0.744–0.758 / KGE 0.728–0.732. DDR's Python trainer reports +//! 0.795 / 0.737 on the same bundle: the gap is the different RNG window +//! stream plus ddrs's parent-space KAN head (one parameter set per cell; +//! DDR varies `log10_uparea` per sub-reach) — see +//! `docs/superpowers/specs/2026-09-08-ddrs-gridded-routing-design.md` §7. +//! The baseline has no RNG and matches DDR's 0.594 / 0.671 to rounding, so it +//! gets a tight band and doubles as the cross-implementation reader check — +//! it is what caught the time-major `Qr(time, divide_id)` layout of DDR's +//! gridded Q' store. +//! +//! Runs only at opt-level 3 (~30 s; a debug build would take minutes): +//! +//! cargo test --release --test gridded_acceptance -- --nocapture +//! +//! Must run from the repo root (bundle data paths in +//! `examples/juniata_gridded/ddrs.yaml` are repo-root-relative). The +//! workspace goes to a tempdir, so `examples/juniata_gridded/.ddrs/` from +//! manual runs is untouched. + +use std::path::PathBuf; + +use ddrs::cli::run::{run, RunInput}; +use ddrs::cli::workspace::Workspace; + +/// Routed-metric floors: below the three-seed range (NSE 0.744–0.758 / +/// KGE 0.728–0.732) by a margin, above the baseline (NSE 0.594 / KGE 0.672). +const NSE_FLOOR: f64 = 0.70; +const KGE_FLOOR: f64 = 0.69; +/// Deterministic summed-Q' NSE (0.594; DDR 0.594). +const BASELINE_NSE_RANGE: (f64, f64) = (0.57, 0.62); + +fn json_f64(v: &serde_json::Value, key: &str) -> f64 { + v.get(key) + .and_then(|x| x.as_f64()) + .unwrap_or_else(|| panic!("manifest metrics missing numeric `{key}`: {v}")) +} + +#[test] +fn gridded_juniata_train_and_test_meets_metric_floors_and_beats_baseline() { + if cfg!(debug_assertions) { + eprintln!( + "skipping: gridded acceptance needs opt-level 3 — \ + run `cargo test --release --test gridded_acceptance -- --nocapture`" + ); + return; + } + + let tmp = tempfile::tempdir().unwrap(); + let run_dir = run(RunInput { + workspace: Workspace::with_root(tmp.path().join(".ddrs")), + config_path: PathBuf::from("examples/juniata_gridded/ddrs.yaml"), + workflow: None, // config says train-and-test + plot: false, + strict: false, + max_mini_batches: None, + batch_order_from: None, + backend: "cpu".into(), + }) + .expect("gridded juniata train-and-test run failed"); + + let manifest: serde_json::Value = serde_json::from_str( + &std::fs::read_to_string(run_dir.join("manifest.json")).expect("read manifest.json"), + ) + .expect("parse manifest.json"); + let metrics = &manifest["metrics"]; + assert_eq!( + metrics["n_gauges_total"].as_u64(), + Some(1), + "expected the single Juniata gauge, got: {metrics}" + ); + let nse = json_f64(metrics, "median_nse_finite"); + let kge = json_f64(metrics, "median_kge_finite"); + + let baseline: serde_json::Value = serde_json::from_str( + &std::fs::read_to_string(run_dir.join("baseline/manifest.json")) + .expect("read baseline/manifest.json"), + ) + .expect("parse baseline/manifest.json"); + let baseline_nse = baseline["metrics"]["nse"][0] + .as_f64() + .expect("baseline metrics.nse[0] missing or null"); + + eprintln!( + "gridded juniata acceptance: routed NSE {nse:.4} / KGE {kge:.4}, baseline NSE {baseline_nse:.4}" + ); + + assert!( + (BASELINE_NSE_RANGE.0..=BASELINE_NSE_RANGE.1).contains(&baseline_nse), + "summed-Q' baseline NSE {baseline_nse:.4} outside {BASELINE_NSE_RANGE:?} — \ + no RNG here, so this is a reader/baseline regression (check the Qr axis \ + order: DDR's gridded store is (time, divide_id)), not noise" + ); + assert!(nse >= NSE_FLOOR, "routed NSE {nse:.4} < floor {NSE_FLOOR} (seed-42 reference 0.751)"); + assert!(kge >= KGE_FLOOR, "routed KGE {kge:.4} < floor {KGE_FLOOR} (seed-42 reference 0.730)"); + assert!( + nse > baseline_nse, + "routed NSE {nse:.4} does not beat the summed-Q' baseline {baseline_nse:.4}" + ); +} diff --git a/tests/gridded_bundle.rs b/tests/gridded_bundle.rs new file mode 100644 index 00000000..6996b80d --- /dev/null +++ b/tests/gridded_bundle.rs @@ -0,0 +1,157 @@ +//! Contract test for the committed gridded (ISIMIP DDM30) Juniata bundle +//! (`examples/juniata_gridded/data/`, byte-identical to DeepGroundwater/ddr +//! PR #194's `examples/juniata_gridded/data/`). +//! +//! Exercises the whole gridded ingestion path through the REAL code: +//! `GriddedNetwork` reads DDR's sub-reach zarr, the managed build relabels +//! it into ddrs's subdivided layout and cuts the gauge subgraph, and +//! `MeritGagesDataset::open` accepts the result. Never skips — the bundle +//! is committed in-repo. + +use ddrs::adjacency::cache::resolve_or_build_gridded; +use ddrs::adjacency::gridded::GriddedNetwork; +use ddrs::config::Config; +use ddrs::data::{ + AttrStats, ConusAdjacencyStore, GageMetadata, GagesAdjacencyStore, MeritGagesDataset, Staid, +}; + +const BUNDLE: &str = "examples/juniata_gridded/data"; +const CONFIG: &str = "examples/juniata_gridded/ddrs.yaml"; +const GAGE: &str = "01567000"; +/// The four DDM30 cells of the Juniata network, upstream first, and the +/// number of sub-reaches DDR cut each into (`length_m` ≈ 7.1–7.8 km). +const CELLS: [(i32, usize); 4] = [(139163, 9), (138443, 6), (138444, 6), (138445, 6)]; +/// The gauge's snapped cell (`juniata_gage.csv`'s `cell` column). +const GAGE_CELL: i64 = 138445; + +#[test] +fn subreach_store_reads_as_27_rows_over_4_cells() { + let g = GriddedNetwork::open(format!("{BUNDLE}/juniata_subreach_adjacency.zarr")) + .expect("open bundle sub-reach store"); + assert_eq!(g.n(), 27); + assert_eq!(g.n_cells(), 4); + assert_eq!(g.rows.len(), 26, "a 27-node tree has 26 edges"); + + let sub = g.to_subdivided(); + let cells: Vec = CELLS.iter().map(|c| c.0).collect(); + assert_eq!(sub.parent_order, cells); + for (p, &(cell, pieces)) in CELLS.iter().enumerate() { + assert_eq!(sub.pieces(p), pieces, "cell {cell} piece count"); + // Every piece of a cell carries the cell id in `order`. + for r in sub.inlet(p)..=sub.outlet(p) { + assert_eq!(sub.order[r], cell); + } + } + // DDR's `length_m` is the cell flow length / pieces: all pieces of a cell + // are equal, and the cell is ~7 km per piece. + for p in 0..4 { + let l0 = sub.length_m[sub.inlet(p)]; + assert!((7000.0..8000.0).contains(&l0), "piece length {l0} m"); + for r in sub.inlet(p)..=sub.outlet(p) { + assert_eq!(sub.length_m[r], l0); + } + } +} + +#[test] +fn gauge_csv_cell_column_names_the_outlet_cell() { + let meta = GageMetadata::open(format!("{BUNDLE}/juniata_gage.csv")).expect("open gauge csv"); + assert_eq!(meta.rows.len(), 1); + assert_eq!(meta.rows[0].staid, Staid::from(GAGE)); + assert_eq!(meta.rows[0].comid, Some(GAGE_CELL)); + // No FLOW_SCALE / COMID_DRAIN_SQKM: the gauge reach's inflow scale stays 1, + // matching DDR's train_gridded.py (flow_scale=None). + assert!(meta.rows[0].flow_scale.is_none()); +} + +/// The managed build into a scratch workspace: the CONUS store round-trips +/// through `ConusAdjacencyStore` as a subdivided store, and the gauge is read +/// at cell 138445's most downstream sub-reach (the last row). +#[test] +fn managed_build_cuts_the_gauge_at_the_cells_outlet_piece() { + let tmp = tempfile::tempdir().unwrap(); + let out = resolve_or_build_gridded( + tmp.path(), + std::path::Path::new(&format!("{BUNDLE}/juniata_subreach_adjacency.zarr")), + std::path::Path::new(&format!("{BUNDLE}/juniata_gage.csv")), + ) + .expect("gridded managed build"); + assert!(!out.cache_hit); + assert!(out.paths.conus.ends_with("juniata_subreach_adjacency.zarr")); + assert!(out.paths.gages.ends_with("juniata_subreach_gages_adjacency.zarr")); + + let conus = ConusAdjacencyStore::open(&out.paths.conus).expect("open built conus store"); + assert_eq!(conus.n, 27); + assert_eq!(conus.nnz, 26); + assert_eq!(conus.n_parent(), 4); + assert_eq!(conus.parent_offset, vec![0, 9, 15, 21, 27]); + let violations = conus + .indices_0 + .iter() + .zip(&conus.indices_1) + .filter(|(r, c)| r <= c) + .count(); + assert_eq!(violations, 0, "{violations} edges are not strictly lower-triangular"); + + let gages = GagesAdjacencyStore::open(&out.paths.gages, &[Staid::from(GAGE)]) + .expect("open built gages store"); + let sg = gages.subgraphs.get(&Staid::from(GAGE)).expect("01567000 subgraph present"); + assert!(!sg.is_headwater()); + assert_eq!(sg.gage_catchment, GAGE_CELL.to_string()); + // Outlet piece of the last cell = row 26 (DDR's `outlet_of_cell`). + assert_eq!(sg.gage_idx, 26); + // Single-gauge bundle: the upstream set is every cell, each once (the + // baseline sums a cell's Q' once, not once per piece). + let upstream = sg.upstream_comids(&conus); + assert_eq!(upstream.len(), 4); + + // Second resolve is a cache hit with the same key. + let again = resolve_or_build_gridded( + tmp.path(), + std::path::Path::new(&format!("{BUNDLE}/juniata_subreach_adjacency.zarr")), + std::path::Path::new(&format!("{BUNDLE}/juniata_gage.csv")), + ) + .expect("second resolve"); + assert!(again.cache_hit); + assert_eq!(again.key, out.key); +} + +#[test] +fn statistics_cover_all_ten_kan_inputs() { + let stats = AttrStats::open(format!( + "{BUNDLE}/statistics/merit_attribute_statistics_ddm30_conus_attributes.nc.json" + )) + .expect("open bundle attribute statistics"); + let cfg = Config::from_yaml_file(CONFIG).expect("load example config"); + let vars = cfg.kan_head.as_ref().unwrap().input_var_names.clone(); + assert_eq!(vars.len(), 10); + let means = stats.means_f32(&vars); + let stds = stats.stds_f32(&vars); + for (i, v) in vars.iter().enumerate() { + assert!(means[i].is_finite(), "mean for {v} not finite"); + assert!(stds[i].is_finite() && stds[i] > 0.0, "std for {v} not usable"); + } +} + +/// End-to-end: the shipped config resolves its adjacency through the gridded +/// managed build and opens through the full dataset pipeline, exactly as +/// `ddrs run` would. +#[test] +fn example_config_opens_the_dataset() { + let mut cfg = Config::from_yaml_file(CONFIG).expect("load example config"); + let ds = cfg.data_sources.as_mut().unwrap(); + assert!(ds.gridded_network.is_some()); + let tmp = tempfile::tempdir().unwrap(); + let out = resolve_or_build_gridded( + tmp.path(), + ds.gridded_network.as_ref().unwrap(), + &ds.gages, + ) + .expect("gridded managed build"); + ds.conus_adjacency = Some(out.paths.conus); + ds.gages_adjacency = Some(out.paths.gages); + + let dataset = MeritGagesDataset::open(&cfg).expect("open gridded Juniata dataset"); + assert_eq!(dataset.len(), 1); + assert_eq!(dataset.staids(), &[Staid::from(GAGE)]); +} diff --git a/tests/hourly_streamflow.rs b/tests/hourly_streamflow.rs index b4f47fd2..2e50e4a6 100644 --- a/tests/hourly_streamflow.rs +++ b/tests/hourly_streamflow.rs @@ -143,6 +143,40 @@ fn streamflow_source_reports_resolution() { assert_eq!(hourly.resolution(), Frequency::Hourly); } +/// DDR's gridded (DDM30) Q' stores are written `Qr(time, divide_id)`; the +/// reader must detect the layout and return the same values as the +/// `(divide_id, time)` contract layout. Before this, the wrong-axis subset +/// came back as fill values (zarrs does not error on an out-of-range read), +/// so the summed-Q' baseline on the gridded Juniata bundle summed nothing. +#[test] +fn time_major_store_reads_identically_to_divide_major() { + let a = StreamflowStore::open(fixture("qr_daily.ic")).expect("open divide-major"); + let b = StreamflowStore::open(fixture("qr_daily_time_major.ic")).expect("open time-major"); + assert!(!a.time_major); + assert!(b.time_major); + assert_eq!(b.resolution, Frequency::Daily); + assert_eq!(b.n_time, 10); + + let w = RhoWindow { + start_day_idx: 2, + rho_days: 4, + window_start: d(1981, 1, 3), + }; + // Reversed and partial COMID list, so the divide-axis span is not the + // whole array and the scatter order differs from storage order. + let comids = [Comid(104), Comid(102), Comid(101)]; + let qa = a.read_window(&w, &comids).expect("read_window divide-major"); + let qb = b.read_window(&w, &comids).expect("read_window time-major"); + assert_eq!(qa, qb); + assert_eq!(qb[(0, 0)], 402.0); // divide 104 (j=3), day 2 + assert_eq!(qb[(71, 2)], 104.0); // divide 101 (j=0), day 4 + + let da = a.read_window_daily(d(1981, 1, 2), 5, &COMIDS).expect("daily divide-major"); + let db = b.read_window_daily(d(1981, 1, 2), 5, &COMIDS).expect("daily time-major"); + assert_eq!(da, db); + assert_eq!(db[(4, 1)], 205.0); // divide 102, day 5 +} + #[test] fn daily_fixture_read_window_keeps_repeat24_semantics() { // Pins the daily path: values repeat 24x per day with the trailing-day trim.