diff --git a/AGENTS.md b/AGENTS.md index 55706b53..8ca7996f 100644 --- a/AGENTS.md +++ b/AGENTS.md @@ -42,9 +42,55 @@ Paths inside `request.json` are relative to the request file's directory. and sw_stage_m (m) instead /io/input/static.json read-only: catchment attributes /io/output/result.csv write: one row per forcing row, spinup included +/io/output/routing.csv write when requested: one row per forcing row and reach /io/output/run.json write: {"status": "ok"} ``` +When `request.request.routing` is non-empty, the adapter must write +`output/routing.csv` as a long table with one row for every forcing time and +every reach declared by `static.json` under `routing_network.reaches`. The +required identifier columns are `time` and `reach_id`; the remaining required +columns are the names listed in `request.request.routing`. + +A probe that declares `requires.routing` must also declare +`requires.static: [routing_network]`. `routing_network` is an object whose +`reaches` member is a non-empty list of unique, non-empty string IDs: + +```json +{ + "routing_network": { + "reaches": ["A", "B", "C"] + } +} +``` + +Keep topology metadata such as junctions, headwaters, outlets, lengths or +areas in separate fields under `routing_network`; do not place reach objects +inside `reaches`. Reach IDs are strings and leading zeros are significant. + +For example, the routing table is: + +```csv +time,reach_id,reach_q_in,q_out,channel_storage +2001-01-01,A,12.3,8.4,336960 +2001-01-01,B,7.1,4.9,190080 +2001-01-01,C,13.3,9.2,354240 +``` + +Each `(time, reach_id)` pair must be unique. Every declared reach must have +exactly one row for every forcing row, including spinup, and each reach's time +axis must match `forcing.csv`. `reach_q_in` and `q_out` are interval-mean discharge +rates in m3/s. `channel_storage` is the absolute volume stored in that reach at +the end of the interval, in m3. It is not a tendency and must not be converted +to catchment-average depth. A missing requested routing column, reach or time +row is a contract error. Models that do not declare `emits.routing` do not +write `routing.csv`. + +A routing table has `n_steps * n_reaches` rows. Probe authors must set +`case.max_output_mb` high enough for that table; the normal output size +allow-list and limit still apply. + + The model seed is deterministic for reproducible stochastic inference, but it is not the generator seed recorded in the host-side report. Exit 0 on success. There is no network. Do not attempt to download weights or data at run time; @@ -171,6 +217,9 @@ numbers in both runs; keep it that way and do not reseed from the clock. | `evspsbl` | evapotranspiration | mm/day | | `mrro` | total runoff | mm/day | | `dis` | river discharge | m3/s | +| `reach_q_in` | interval-mean total inflow into the identified reach control volume, including upstream, lateral, and prescribed boundary inflows; positive into the reach; reach-indexed in `routing.csv`, never a catchment-average flux | m3/s | +| `q_out` | interval-mean flow across the identified reach's downstream boundary; positive downstream and negative for reverse or backwater flow; reach-indexed in `routing.csv`, never a catchment-average flux | m3/s | +| `channel_storage` | absolute water volume stored in the routed river channel of the identified reach at the end of the row's time interval; excludes land-surface, soil, groundwater, floodplain, and overland storage; reach-indexed in `routing.csv`, not a tendency and not catchment-average `channel` storage | m3 | | `gwex` | a declared exchange with the outside: regional groundwater, inter-basin transfer; positive into the catchment | mm/day | | `gw_sw_exchange` | net river-aquifer exchange, positive into the aquifer; unlike `gwex`, this moves water between two stores inside the control volume (`gw` and `channel`), so it is never added to `gwex` or counted as a `closure` source | mm/day | | `gw_to_sw` | groundwater-to-river exchange component: the aquifer losing to the river, so it is never positive; a component of `gw_sw_exchange`, not of `gwex`, and not an addition to it | mm/day | @@ -255,6 +304,28 @@ emits: states: [] ``` +A model that also exposes native reach-indexed routing output declares it +separately from the catchment-level result columns: + +```yaml +emits: + fluxes: [mrro, dis] + states: [channel] + routing: [reach_q_in, q_out, channel_storage] +needs_static: [routing_network] +``` + +Variables declared under `emits.routing` belong in `routing.csv`, not in +`result.csv`, `emits.fluxes` or `emits.states`. Declaring them means the +adapter must write the requested long table; do not construct `reach_q_in` or +`channel_storage` from a budget residual merely to satisfy the contract. + + +A model that declares `emits.routing` must also list `routing_network` in +`needs_static` or `uses_static`; otherwise the manifest is invalid because the +routing output has no declared network topology. + + Every probe that needs more than that will be `N/A` with reason `INCOMPLETE`, which is the honest outcome: not a fail, but a probe that cannot be put to the model, so it counts neither way. Fabricating an @@ -297,7 +368,10 @@ ht run --model # the actual evaluation ``` `verify-adapter` runs a single seed, on the same window the evaluation will -use, and checks the shape of what came back. Get that green before looking +use, and checks the shape of what came back. For a model that declares +`emits.routing`, this includes the presence of `routing.csv`, its requested +columns, reach IDs, complete time coverage, unique `(time, reach_id)` pairs +and finite numeric values. Get that green before looking at any residual. Without `--probe` it checks the closure probe, or the first probe the model can consume when it cannot consume that one: a step it does not declare, a forcing or static input the probe does not generate, or a window that drops a diff --git a/docs/adapting-a-model.md b/docs/adapting-a-model.md index 64735a6d..5d63167e 100644 --- a/docs/adapting-a-model.md +++ b/docs/adapting-a-model.md @@ -64,6 +64,33 @@ diagnostic but omitting it from the output is a protocol error. The separate `ts` diagnostic is instantaneous surface temperature for radiative checks; it cannot replace the interval-end mean-layer `tsoil_layer`. +Reach-indexed river-network outputs use the optional `routing` group and are +written separately from the catchment-level result table: + +```yaml +emits: + fluxes: [mrro, dis] + states: [channel] + routing: [reach_q_in, q_out, channel_storage] +needs_static: [routing_network] +``` + +A model that declares `emits.routing` must also list `routing_network` in +`needs_static` or `uses_static`, so its adapter has the network topology that +defines the declared reaches. + +Variables declared under `emits.routing` belong in `output/routing.csv`, not +`output/result.csv`. The harness requests them only on a probe that declares +`requires.routing` and supplies `routing_network.reaches` in `static.json`. +The routing table contains one row for every forcing time and declared reach. + +`reach_q_in` is the interval-mean rate of total inflow into the reach control +volume, including upstream, lateral, and prescribed boundary inflows, and is +positive into the reach. `q_out` is the interval-mean flow across the reach's +downstream boundary, positive downstream and negative for reverse or backwater +flow. `channel_storage` is the absolute end-of-step channel water volume and +excludes land-surface, soil, groundwater, floodplain, and overland storage. + The `stage` diagnostic is a water-surface elevation in metres on the fixed vertical datum declared by the case. `bed_elevation_m` is the rigid reach bed on that same datum, so water depth is `stage - bed_elevation_m` at every step. @@ -91,8 +118,10 @@ river/reach-inflow path whenever the column is present. ## 2. Write the adapter -Read `/io/request.json`, read the forcing, call your model, write the table. -The plumbing is about thirty lines and does not depend on your model. +Read `/io/request.json`, read the forcing and static inputs, call your model, +and write `output/result.csv`. When `request.request.routing` is non-empty, +also write `output/routing.csv` using the paths declared under +`request.output`. The plumbing does not otherwise depend on your model. ```python request = json.loads(Path(sys.argv[-1]).read_text()) diff --git a/docs/writing-a-probe.md b/docs/writing-a-probe.md index aa5fbbd9..9277a8bf 100644 --- a/docs/writing-a-probe.md +++ b/docs/writing-a-probe.md @@ -85,6 +85,26 @@ Requirements: Determinism is not a nicety. A benchmark whose failures cannot be reproduced is unusable the first time a result is disputed. +## Required variables + +`requires` names the inputs and outputs that the verdict depends on. Ask for +only what the probe actually uses: a model that cannot supply a requested +output is N/A with reason INCOMPLETE, not a physics failure. + +A reach-indexed routing probe declares its routing outputs separately from the +catchment-level `fluxes` and `states`: + +```yaml +requires: + routing: [reach_q_in, q_out, channel_storage] + static: [routing_network] +``` + +Its generator must supply `static["routing_network"]["reaches"]` as a +non-empty list of unique, non-empty string IDs. The adapter writes the +requested variables to `output/routing.csv`, with one row for every forcing +time and declared reach. + ## Criteria Each entry in `criteria` names a registered criterion and its parameters. diff --git a/schemas/model.schema.json b/schemas/model.schema.json index be1eb500..ea9a1d50 100644 --- a/schemas/model.schema.json +++ b/schemas/model.schema.json @@ -40,6 +40,10 @@ "type": "array", "items": {"type": "string"}, "uniqueItems": true, "description": "Absolute storage values, never tendencies. The harness differences them." }, + "routing": { + "type": "array", "items": {"type": "string"}, "uniqueItems": true, + "description": "Reach-indexed variables written to output/routing.csv. reach_q_in is total inflow into the reach control volume (upstream, lateral, and prescribed), positive into the reach; q_out crosses the downstream boundary, positive downstream and negative for reverse or backwater flow. Both are interval-mean discharges in m3 s-1. channel_storage is absolute end-of-step channel water volume in m3 and excludes land-surface, soil, groundwater, floodplain, and overland storage." + }, "diagnostics": { "type": "array", "items": {"type": "string"}, "uniqueItems": true, "description": "Diagnostics such as a surface or layer temperature, or the liquid share of a snowpack. Never added to any water-storage sum: `lwsnl` is water, but it is water already counted inside `snw`, so summing it beside `snw` would count the same kilogram twice. Each criterion defines their time handling." diff --git a/schemas/probe.schema.json b/schemas/probe.schema.json index a6c270d5..c08463c4 100644 --- a/schemas/probe.schema.json +++ b/schemas/probe.schema.json @@ -103,6 +103,14 @@ "uniqueItems": true, "description": "Diagnostics read by criteria, such as a surface or layer temperature, or the liquid share of a snowpack. Never added to any water-storage sum: `lwsnl` is water, but it is water already counted inside `snw`, so summing it beside `snw` would count the same kilogram twice. Each criterion defines their time handling." }, + "routing": { + "type": "array", + "items": { + "type": "string" + }, + "uniqueItems": true, + "description": "Reach-indexed variables read from output/routing.csv. reach_q_in is total inflow into the reach control volume (upstream, lateral, and prescribed), positive into the reach; q_out crosses the downstream boundary, positive downstream and negative for reverse or backwater flow. Both are interval-mean discharges in m3 s-1. channel_storage is absolute end-of-step channel water volume in m3 and excludes land-surface, soil, groundwater, floodplain, and overland storage." + }, "forcing": { "type": "array", "items": { diff --git a/src/hydroturing/cli.py b/src/hydroturing/cli.py index 91302798..cf90d69a 100644 --- a/src/hydroturing/cli.py +++ b/src/hydroturing/cli.py @@ -170,7 +170,7 @@ def cmd_list(args) -> int: models = registry.all_models() print(f"probes ({len(probes)}):") for p in probes: - req = ", ".join(p.required_vars) + req = ", ".join((*p.required_vars, *p.requires_routing)) steps = "/".join(p.timesteps) print(f" {p.id:<34} {p.law:<9} {p.track:<10} {steps:<10} seeds={p.n_seeds} needs: {req}") print(f"\nmodels ({len(models)}):") @@ -179,9 +179,10 @@ def cmd_list(args) -> int: days = resolve_window_days(m, matching) if matching else m.window_days window = "full" if days is None else f"{days}d" steps = "/".join(m.timesteps) + emits = ", ".join((*m.emitted, *m.emits_routing)) print( f" {m.name:<24} v{m.version:<8} runner={m.runner:<11} {steps:<22} " - f"window={window:<5} emits: {', '.join(m.emitted)}" + f"window={window:<5} emits: {emits}" ) return 0 @@ -233,12 +234,29 @@ def cmd_verify_adapter(args) -> int: if args.probe: candidates = [registry.find_probe(args.probe, roots)] else: - # The closure probe is the reference implementation and the natural - # smoke test. A model that cannot consume it, at its step, with its - # forcing or on its window, is checked on the first probe it can - # consume instead: one it cannot is N/A for it and asks the adapter - # nothing. - candidates = sorted(probes, key=lambda p: (p.id != "mass/catchment-closure", p.id)) + if model.emits_routing: + # Reach-indexed output can only be verified on a probe that + # supplies a routing network. Prefer one when the model declares + # routing output, then fall back to the ordinary probe order. + candidates = sorted( + probes, + key=lambda p: ( + not bool(p.requires_routing), + p.id != "mass/catchment-closure", + p.id, + ), + ) + else: + # The closure probe is the reference implementation and the + # natural smoke test. A model that cannot consume it is checked + # on the first probe it can consume instead. + candidates = sorted( + probes, + key=lambda p: ( + p.id != "mass/catchment-closure", + p.id, + ), + ) workdir = Path(args.workdir) if args.workdir else None result = None diff --git a/src/hydroturing/harness.py b/src/hydroturing/harness.py index 2b3715a0..37e796bd 100644 --- a/src/hydroturing/harness.py +++ b/src/hydroturing/harness.py @@ -512,6 +512,7 @@ def verify_adapter_contract( requires_fluxes=model.emits_fluxes, requires_states=model.emits_states, requires_diagnostics=model.emits_diagnostics, + requires_routing=(model.emits_routing if probe.requires_routing else ()), variants=(), criteria=(), ) diff --git a/src/hydroturing/protocol.py b/src/hydroturing/protocol.py index 16ede4e4..3eaed4be 100644 --- a/src/hydroturing/protocol.py +++ b/src/hydroturing/protocol.py @@ -9,6 +9,7 @@ input/forcing.csv <- harness writes, read only input/static.json <- harness writes, read only output/result.csv -> adapter writes + output/routing.csv -> adapter writes, when reach-indexed output is requested output/run.json -> adapter writes CSV is the required format. NetCDF is accepted when present but never @@ -35,10 +36,11 @@ STATIC_FILE = "input/static.json" RESULT_CSV = "output/result.csv" RESULT_NC = "output/result.nc" +ROUTING_CSV = "output/routing.csv" RUN_FILE = "output/run.json" TIME_COL = "time" - +REACH_COL = "reach_id" class ProtocolError(RuntimeError): """The adapter did not honour the contract.""" @@ -93,6 +95,33 @@ def after_spinup(self, frame: pd.DataFrame) -> pd.DataFrame: return frame.iloc[self.spinup_steps :].reset_index(drop=True) +def _routing_reaches(case: Case) -> tuple[str, ...]: + """Return and validate the reach IDs declared by a routing case.""" + network = case.static.get("routing_network") + if not isinstance(network, dict) or "reaches" not in network: + raise ValueError( + "routing probes require static.json routing_network.reaches" + ) + + reaches = network["reaches"] + if ( + not isinstance(reaches, list) + or not reaches + or any(not isinstance(reach, str) or not reach for reach in reaches) + ): + raise ValueError( + "static.json routing_network.reaches must be a non-empty " + "list of non-empty string IDs" + ) + + declared = tuple(reaches) + if len(declared) != len(set(declared)): + raise ValueError( + "static.json routing_network.reaches contains duplicates" + ) + return declared + + @dataclass class RunResult: """What a model returned for one case, after contract validation.""" @@ -101,6 +130,7 @@ class RunResult: table: pd.DataFrame meta: dict[str, Any] wall_seconds: float + routing: pd.DataFrame | None = None # The manifest of the model that produced the table. A criterion whose # verdict depends on what the model declared it consumes — a prescribed # driver it may or may not have read — needs this; every other criterion @@ -111,6 +141,10 @@ class RunResult: def stage(io_dir: Path, case: Case, probe: ProbeSpec, model: ModelManifest) -> Path: """Write request.json and the read-only inputs. Returns the request path.""" + if probe.requires_routing: + _routing_reaches(case) + + (io_dir / "input").mkdir(parents=True, exist_ok=True) (io_dir / "output").mkdir(parents=True, exist_ok=True) @@ -125,21 +159,33 @@ def stage(io_dir: Path, case: Case, probe: ProbeSpec, model: ModelManifest) -> P json.dump(case.static, fh, indent=2) case_id, model_seed = _opaque_case_metadata(case) + requested_routing = ( + list(model.emits_routing) if probe.requires_routing else [] + ) request = { "case_id": case_id, "seed": model_seed, "timestep": case.timestep, "n_steps": case.n_steps, "request": { - # Asking for every declared output keeps this part of the request - # invariant across probes and prevents it identifying the criterion. + # Catchment-level groups remain invariant across probes. + # Routing is requested only when the probe supplies a network. "fluxes": list(model.emits_fluxes), "states": list(model.emits_states), "diagnostics": list(model.emits_diagnostics), + "routing": requested_routing, }, "input": {"forcing": FORCING_FILE, "static": STATIC_FILE}, - "output": {"table": RESULT_CSV, "run": RUN_FILE}, - "units": {v: UNITS[v] for v in model.emitted if v in UNITS}, + "output": { + "table": RESULT_CSV, + "routing": ROUTING_CSV, + "run": RUN_FILE, + }, + "units": { + v: UNITS[v] + for v in (*model.emitted, *requested_routing) + if v in UNITS + }, "notes": ( "States are absolute storages, not tendencies. The harness " "differences them itself. Row i's ts and rlus are instantaneous " @@ -147,6 +193,13 @@ def stage(io_dir: Path, case: Case, probe: ProbeSpec, model: ModelManifest) -> P "when supplied. Row i's tsoil_layer is the mean temperature of " "the specified layer at the interval end; hfg and hfg_bottom " "are interval-mean boundary fluxes for the soil-storage check." + " Routing reach_q_in is total inflow into the reach control " + "volume (upstream, lateral, and prescribed), positive into the " + "reach; q_out is flow across the downstream boundary, positive " + "downstream and negative for reverse or backwater flow. Both are " + "interval-mean m3 s-1. channel_storage is absolute end-of-step " + "channel water volume in m3 and excludes land-surface, soil, " + "groundwater, floodplain, and overland storage." ), } request_path = io_dir / REQUEST_FILE @@ -161,7 +214,7 @@ def _validate_output_files(io_dir: Path, probe: ProbeSpec) -> None: if not output_dir.exists(): return - allowed = {RESULT_CSV, RESULT_NC, RUN_FILE} + allowed = {RESULT_CSV, RESULT_NC, ROUTING_CSV, RUN_FILE} total = 0 for path in output_dir.iterdir(): relative = f"output/{path.name}" @@ -180,7 +233,12 @@ def _validate_output_files(io_dir: Path, probe: ProbeSpec) -> None: ) -def _validate_time_axis(table: pd.DataFrame, case: Case) -> None: +def _validate_time_axis( + table: pd.DataFrame, + case: Case, + *, + label: str = "result", +) -> None: expected = case.forcing[TIME_COL].reset_index(drop=True) actual = table[TIME_COL].reset_index(drop=True) @@ -200,11 +258,90 @@ def _validate_time_axis(table: pd.DataFrame, case: Case) -> None: if not equal: first = int(mismatch.to_numpy().nonzero()[0][0]) raise ProtocolError( - f"result time axis differs from the forcing at row {first}: " + f"{label} time axis differs from the forcing at row {first}: " f"got {actual.iloc[first]!r}, expected {expected.iloc[first]!r}" ) +def _read_routing( + io_dir: Path, + case: Case, + probe: ProbeSpec, +) -> pd.DataFrame | None: + """Read and validate the optional reach-indexed routing result.""" + if not probe.requires_routing: + return None + + path = io_dir / ROUTING_CSV + if not path.exists(): + raise ProtocolError(f"adapter did not write required {ROUTING_CSV}") + + table = pd.read_csv(path, dtype={REACH_COL: str}, keep_default_na=False) + required_columns = [TIME_COL, REACH_COL, *probe.requires_routing] + missing = [ + column for column in required_columns if column not in table.columns + ] + if missing: + raise ProtocolError(f"routing result is missing columns: {missing}") + + if table[REACH_COL].isna().any(): + raise ProtocolError(f"routing result has missing '{REACH_COL}' values") + table[REACH_COL] = table[REACH_COL].astype(str) + + duplicate = table.duplicated([TIME_COL, REACH_COL]) + if duplicate.any(): + raise ProtocolError( + "routing result has duplicate (time, reach_id) rows: " + f"{int(duplicate.sum())} duplicates" + ) + + expected_reaches = _routing_reaches(case) + + actual_reaches = set(table[REACH_COL]) + expected_set = set(expected_reaches) + if actual_reaches != expected_set: + missing_reaches = sorted(expected_set - actual_reaches) + extra_reaches = sorted(actual_reaches - expected_set) + raise ProtocolError( + "routing reach IDs differ from static.json: " + f"missing={missing_reaches}, extra={extra_reaches}" + ) + + expected_rows = case.n_steps * len(expected_reaches) + if len(table) != expected_rows: + raise ProtocolError( + f"routing result has {len(table)} rows, expected {expected_rows} " + "(one row per forcing step and reach, spinup included)" + ) + + for reach in expected_reaches: + reach_table = table.loc[ + table[REACH_COL] == reach + ].reset_index(drop=True) + if len(reach_table) != case.n_steps: + raise ProtocolError( + f"routing reach '{reach}' has {len(reach_table)} rows, " + f"expected {case.n_steps}" + ) + _validate_time_axis( + reach_table, + case, + label=f"routing.csv reach {reach!r}", + ) + + for var in probe.requires_routing: + column = pd.to_numeric(table[var], errors="coerce") + invalid = ~np.isfinite(column) + if invalid.any(): + raise ProtocolError( + f"routing column '{var}' has " + f"{int(invalid.sum())} non-finite values" + ) + table[var] = column + + return table + + def read_result( io_dir: Path, case: Case, probe: ProbeSpec, wall_seconds: float, model: ModelManifest | None = None, @@ -238,6 +375,7 @@ def read_result( ) _validate_time_axis(table, case) + routing = _read_routing(io_dir, case, probe) missing = [v for v in probe.required_vars if v not in table.columns] if missing: @@ -272,4 +410,11 @@ def read_result( f"run.json reports n_steps={meta['n_steps']!r}, expected {case.n_steps}" ) - return RunResult(case=case, table=table, meta=meta, wall_seconds=wall_seconds, model=model) + return RunResult( + case=case, + table=table, + meta=meta, + wall_seconds=wall_seconds, + routing=routing, + model=model, + ) diff --git a/src/hydroturing/spec.py b/src/hydroturing/spec.py index 618e33e4..28456e52 100644 --- a/src/hydroturing/spec.py +++ b/src/hydroturing/spec.py @@ -39,11 +39,19 @@ # pack rather than evidence of a phase change. DIAG_VARS = ("ts", "tsoil_layer", "stage", "lwsnl", "csnow") +# Reach-indexed routing variables live in output/routing.csv rather than the +# one-row-per-timestep result table. +ROUTING_VARS = ("reach_q_in", "q_out", "channel_storage") + UNITS = { "pr": "mm day-1", "evspsbl": "mm day-1", "mrro": "mm day-1", "dis": "m3 s-1", + "q_in": "m3 s-1", + "reach_q_in": "m3 s-1", + "q_out": "m3 s-1", + "channel_storage": "m3", "gwex": "mm day-1", # Net river-aquifer exchange, positive into the aquifer. Unlike gwex # (a source or sink crossing the catchment boundary), this moves water @@ -270,6 +278,8 @@ class ProbeSpec: min_scored_fraction: float = 1.0 # Missing diagnostic outputs cause INCOMPLETE, as for missing fluxes. requires_diagnostics: tuple[str, ...] = () + # Reach-indexed variables requested from output/routing.csv. + requires_routing: tuple[str, ...] = () # Case-supplied inputs the verdict rests on: forcing columns and # static.json keys a model must declare it consumes, or it is judged # against values it never read and is INCOMPATIBLE instead. @@ -379,6 +389,8 @@ class ModelManifest: window_days: int | str | None = None # Diagnostics the model reports in addition to its fluxes and states. emits_diagnostics: tuple[str, ...] = () + # Reach-indexed variables written to output/routing.csv. + emits_routing: tuple[str, ...] = () # Static inputs the adapter cannot run without. needs_static: tuple[str, ...] = () # Optional inputs the adapter consumes whenever the case supplies them. @@ -404,7 +416,11 @@ def missing_for(self, probe: ProbeSpec) -> list[str]: INCOMPLETE: the model cannot demonstrate conservation because it never says enough to be checked, and it has not violated it either. """ - return [v for v in probe.required_vars if v not in self.emitted] + missing = [v for v in probe.required_vars if v not in self.emitted] + missing.extend( + v for v in probe.requires_routing if v not in self.emits_routing + ) + return missing def load_probe(path: str | Path) -> ProbeSpec: @@ -497,12 +513,24 @@ def load_probe(path: str | Path) -> ProbeSpec: # asked for as a state would be met by the flux and never noticed. unknown = [ name - for key, known in (("fluxes", FLUX_VARS), ("states", STATE_VARS), ("diagnostics", DIAG_VARS)) + for key, known in ( + ("fluxes", FLUX_VARS), + ("states", STATE_VARS), + ("diagnostics", DIAG_VARS), + ("routing", ROUTING_VARS), + ) for name in requires.get(key, []) if name not in known ] if unknown: raise SpecError(f"{spec_file}: unknown variables in requires: {unknown}") + if requires.get("routing") and "routing_network" not in requires.get( + "static", [] + ): + raise SpecError( + f"{spec_file}: requires.routing also needs requires.static " + "to include routing_network" + ) return ProbeSpec( id=raw["id"], title=raw["title"], @@ -514,6 +542,7 @@ def load_probe(path: str | Path) -> ProbeSpec: requires_fluxes=tuple(requires.get("fluxes", [])), requires_states=tuple(requires.get("states", [])), requires_diagnostics=tuple(requires.get("diagnostics", [])), + requires_routing=tuple(requires.get("routing", [])), requires_forcing=tuple(requires.get("forcing", [])), requires_static=tuple(requires.get("static", [])), generator=case["generator"], @@ -591,9 +620,19 @@ def load_model(path: str | Path) -> ModelManifest: unknown = [v for v in raw["emits"]["fluxes"] if v not in FLUX_VARS] unknown += [v for v in raw["emits"]["states"] if v not in STATE_VARS] unknown += [v for v in raw["emits"].get("diagnostics", []) if v not in DIAG_VARS] + unknown += [v for v in raw["emits"].get("routing", []) if v not in ROUTING_VARS] if unknown: raise SpecError(f"{spec_file}: unknown variables in emits: {unknown}") + if raw["emits"].get("routing") and "routing_network" not in { + *raw.get("needs_static", []), + *raw.get("uses_static", []), + }: + raise SpecError( + f"{spec_file}: emits.routing also needs needs_static or " + "uses_static to include routing_network" + ) + declared = raw["timestep"] timesteps = tuple(declared) if isinstance(declared, list) else (declared,) if len(set(timesteps)) != len(timesteps): @@ -607,6 +646,7 @@ def load_model(path: str | Path) -> ModelManifest: emits_fluxes=tuple(raw["emits"]["fluxes"]), emits_states=tuple(raw["emits"]["states"]), emits_diagnostics=tuple(raw["emits"].get("diagnostics", [])), + emits_routing=tuple(raw["emits"].get("routing", [])), runner=runner, needs_forcing=tuple(raw.get("needs_forcing", [])), needs_static=tuple(raw.get("needs_static", [])), diff --git a/templates/probe.template.yaml b/templates/probe.template.yaml index ccda50b5..4b404525 100644 --- a/templates/probe.template.yaml +++ b/templates/probe.template.yaml @@ -53,12 +53,14 @@ description: > #! fluxes: pr, evspsbl, mrro, dis (rates) #! states: mrso, snw, canopy (absolute storages, never tendencies) #! diagnostics: ts (read by a criterion, never integrated) +#! routing: reach_q_in, q_out, channel_storage (reach-indexed outputs) #! forcing, static: case inputs the verdict rests on; a model must declare #! it consumes them (needs_* or uses_*) or it is #! N/A with reason INCOMPATIBLE requires: fluxes: [pr, evspsbl, mrro] states: [mrso, snw, canopy] + routing: [] case: #! Must define generate(seed) -> (DataFrame, dict) and be deterministic. diff --git a/tests/test_docs_in_sync.py b/tests/test_docs_in_sync.py index a0ba694c..19c66480 100644 --- a/tests/test_docs_in_sync.py +++ b/tests/test_docs_in_sync.py @@ -15,7 +15,7 @@ import pytest import yaml -from hydroturing.spec import DIAG_VARS, FLUX_VARS, REPO_ROOT, STATE_VARS, load_model, load_probe +from hydroturing.spec import DIAG_VARS, FLUX_VARS, REPO_ROOT, ROUTING_VARS, STATE_VARS, load_model, load_probe PROBES = {p.id: p for p in (load_probe(x) for x in REPO_ROOT.glob("probes/*/*/probe.yaml"))} PROBE_IDS = sorted(PROBES) @@ -229,7 +229,7 @@ def test_site_charts_are_drawn_from_the_archive(): assert (passed, scored) == expected, f"{model}: chart says {passed} of {scored}, archive {expected}" -@pytest.mark.parametrize("variable", FLUX_VARS + STATE_VARS + DIAG_VARS) +@pytest.mark.parametrize("variable", FLUX_VARS + STATE_VARS + DIAG_VARS + ROUTING_VARS) def test_agents_doc_defines_every_variable(variable): assert f"| `{variable}` |" in read("AGENTS.md"), ( f"{variable} is accepted by spec.py but the table in AGENTS.md does not define it" diff --git a/tests/test_routing_protocol.py b/tests/test_routing_protocol.py new file mode 100644 index 00000000..3f311dfc --- /dev/null +++ b/tests/test_routing_protocol.py @@ -0,0 +1,441 @@ +"""Reach-indexed routing declarations and /io tables obey the contract.""" + +from __future__ import annotations + +import json +from dataclasses import replace +from types import SimpleNamespace + +import numpy as np +import pandas as pd +import pytest +import yaml + +from hydroturing import registry +from hydroturing.harness import verify_adapter_contract +from hydroturing.protocol import Case, ProtocolError, read_result, stage +from hydroturing.spec import SpecError, load_model, load_probe + + +@pytest.fixture +def base_probe(): + return registry.find_probe("mass/catchment-closure") + + +@pytest.fixture +def routing_probe(base_probe): + return replace( + base_probe, + requires_fluxes=(), + requires_states=(), + requires_diagnostics=(), + requires_routing=("reach_q_in", "q_out", "channel_storage"), + requires_static=("routing_network",), + variants=(), + ) + + +@pytest.fixture +def case(): + return Case( + probe_id="mass/routing-contract-test", + seed=11, + forcing=pd.DataFrame( + { + "time": pd.date_range("2000-01-01", periods=3, freq="D"), + "pr": [0.0, 1.0, 0.0], + } + ), + static={"routing_network": {"reaches": ["A", "B"]}}, + spinup_steps=1, + ) + + +def _write_probe(path, base_probe, routing): + path.mkdir(parents=True) + raw = yaml.safe_load((base_probe.path / "probe.yaml").read_text()) + raw["id"] = "mass/routing-contract-test" + raw["requires"] = {"routing": routing, "static": ["routing_network"]} + (path / "generate.py").write_text("# Only manifest loading is tested.\n") + (path / "probe.yaml").write_text(yaml.safe_dump(raw)) + + +def _write_model(path, routing): + path.mkdir() + manifest = { + "name": path.name, + "version": "1", + "entrypoint": ["python", "adapter.py"], + "timestep": "PT1D", + "emits": {"fluxes": [], "states": [], "routing": routing}, + "needs_static": ["routing_network"], + } + (path / "model.yaml").write_text(yaml.safe_dump(manifest)) + + +def _routing_table(case): + rows = [] + for time in case.forcing["time"]: + rows.extend( + [ + { + "time": time.strftime("%Y-%m-%d"), + "reach_id": "A", + "reach_q_in": 2.0, + "q_out": 1.0, + "channel_storage": 86400.0, + }, + { + "time": time.strftime("%Y-%m-%d"), + "reach_id": "B", + "reach_q_in": 3.0, + "q_out": 2.0, + "channel_storage": 86400.0, + }, + ] + ) + return pd.DataFrame(rows) + + +def _write_outputs(path, case, routing=None): + output = path / "output" + output.mkdir(parents=True, exist_ok=True) + pd.DataFrame({"time": case.forcing["time"]}).to_csv( + output / "result.csv", index=False + ) + if routing is not None: + routing.to_csv(output / "routing.csv", index=False) + (output / "run.json").write_text(json.dumps({"status": "ok"})) + + +def test_routing_declarations_load_and_old_defaults_remain_empty( + base_probe, tmp_path +): + probe_dir = tmp_path / "mass" / "routing-contract-test" + variables = ["reach_q_in", "q_out", "channel_storage"] + _write_probe(probe_dir, base_probe, variables) + probe = load_probe(probe_dir) + + model_dir = tmp_path / "routing_model" + _write_model(model_dir, variables) + model = load_model(model_dir) + + assert probe.requires_routing == model.emits_routing == tuple(variables) + assert model.missing_for(probe) == [] + assert base_probe.requires_routing == () + assert registry.find_model("reference_bucket").emits_routing == () + + +def test_routing_probe_must_require_network(base_probe, tmp_path): + path = tmp_path / "mass" / "routing-contract-test" + _write_probe(path, base_probe, ["q_out"]) + raw = yaml.safe_load((path / "probe.yaml").read_text()) + raw["requires"].pop("static") + (path / "probe.yaml").write_text(yaml.safe_dump(raw)) + + with pytest.raises(SpecError, match="routing_network"): + load_probe(path) + + +def test_routing_model_must_declare_network(tmp_path): + path = tmp_path / "routing_model" + _write_model(path, ["q_out"]) + raw = yaml.safe_load((path / "model.yaml").read_text()) + raw.pop("needs_static") + (path / "model.yaml").write_text(yaml.safe_dump(raw)) + + with pytest.raises(SpecError, match="routing_network"): + load_model(path) + + +@pytest.mark.parametrize("kind", ["probe", "model"]) +def test_unknown_routing_variable_is_rejected(base_probe, tmp_path, kind): + if kind == "probe": + path = tmp_path / "mass" / "routing-contract-test" + _write_probe(path, base_probe, ["not_a_routing_variable"]) + loader = load_probe + else: + path = tmp_path / "routing_model" + _write_model(path, ["not_a_routing_variable"]) + loader = load_model + + with pytest.raises(SpecError, match="unknown variables"): + loader(path) + + +def test_stage_requests_routing_path_and_units(case, routing_probe, tmp_path): + model = replace( + registry.find_model("reference_bucket"), + emits_routing=("reach_q_in", "q_out", "channel_storage"), + ) + request = json.loads(stage(tmp_path, case, routing_probe, model).read_text()) + + assert request["request"]["routing"] == [ + "reach_q_in", + "q_out", + "channel_storage", + ] + assert request["output"]["routing"] == "output/routing.csv" + assert request["units"]["reach_q_in"] == "m3 s-1" + assert request["units"]["q_out"] == "m3 s-1" + assert request["units"]["channel_storage"] == "m3" + + +def test_stage_omits_routing_for_nonrouting_probe( + case, base_probe, tmp_path +): + model = replace( + registry.find_model("reference_bucket"), + emits_routing=("reach_q_in",), + ) + request = json.loads( + stage(tmp_path, case, base_probe, model).read_text() + ) + + assert request["request"]["routing"] == [] + assert "reach_q_in" not in request["units"] + + +def test_unrequested_routing_output_is_ignored( + case, routing_probe, tmp_path +): + probe = replace( + routing_probe, + requires_routing=(), + requires_static=(), + ) + case = replace(case, static={}) + _write_outputs(tmp_path, case, _routing_table(case)) + + result = read_result(tmp_path, case, probe, 0.0) + + assert result.routing is None + + +def test_valid_routing_table_roundtrips(case, routing_probe, tmp_path): + expected = _routing_table(case) + _write_outputs(tmp_path, case, expected) + + result = read_result(tmp_path, case, routing_probe, 0.0) + + assert result.routing is not None + pd.testing.assert_frame_equal(result.routing, expected) + + +def test_routing_output_requires_declared_reaches( + case, routing_probe, tmp_path +): + case = replace(case, static={}) + _write_outputs(tmp_path, case, _routing_table(case)) + + with pytest.raises(ValueError, match="routing_network.reaches"): + read_result(tmp_path, case, routing_probe, 0.0) + + +@pytest.mark.parametrize( + ("reaches", "message"), + [ + ([], "non-empty"), + (["A", "A"], "duplicates"), + ([{"id": "A"}, {"id": "B"}], "string IDs"), + ], +) +def test_invalid_declared_reach_ids_are_rejected( + case, routing_probe, tmp_path, reaches, message +): + case = replace( + case, + static={"routing_network": {"reaches": reaches}}, + ) + _write_outputs(tmp_path, case, _routing_table(case)) + + with pytest.raises(ValueError, match=message): + read_result(tmp_path, case, routing_probe, 0.0) + + +def test_stage_rejects_malformed_network(case, routing_probe, tmp_path): + case = replace( + case, + static={"routing_network": {"reaches": [1, 2]}}, + ) + model = registry.find_model("reference_bucket") + + with pytest.raises(ValueError, match="string IDs"): + stage(tmp_path, case, routing_probe, model) + + assert not (tmp_path / "request.json").exists() + + +def test_leading_zero_reach_ids_are_preserved( + case, routing_probe, tmp_path +): + case = replace( + case, + static={"routing_network": {"reaches": ["007", "012"]}}, + ) + table = _routing_table(case) + table["reach_id"] = table["reach_id"].replace( + {"A": "007", "B": "012"} + ) + _write_outputs(tmp_path, case, table) + + result = read_result(tmp_path, case, routing_probe, 0.0) + + assert result.routing is not None + assert set(result.routing["reach_id"]) == {"007", "012"} + + +@pytest.mark.parametrize( + "reach_id", + ["NA", "N/A", "n/a", "nan", "None", "null", "NULL"], +) +def test_na_like_reach_ids_are_preserved( + case, routing_probe, tmp_path, reach_id +): + case = replace( + case, + static={"routing_network": {"reaches": [reach_id, "B"]}}, + ) + table = _routing_table(case) + table["reach_id"] = table["reach_id"].replace({"A": reach_id}) + _write_outputs(tmp_path, case, table) + + result = read_result(tmp_path, case, routing_probe, 0.0) + + assert result.routing is not None + assert set(result.routing["reach_id"]) == {reach_id, "B"} + + +def test_missing_required_routing_file_is_rejected(case, routing_probe, tmp_path): + _write_outputs(tmp_path, case) + with pytest.raises(ProtocolError, match="required output/routing.csv"): + read_result(tmp_path, case, routing_probe, 0.0) + + +@pytest.mark.parametrize( + ("mutation", "message"), + [ + (lambda table: table.drop(columns="q_out"), "missing columns"), + ( + lambda table: pd.concat([table, table.iloc[[0]]], ignore_index=True), + "duplicate", + ), + (lambda table: table.loc[table["reach_id"] != "B"], "reach IDs"), + ], +) +def test_malformed_routing_shape_is_rejected( + case, routing_probe, tmp_path, mutation, message +): + _write_outputs(tmp_path, case, mutation(_routing_table(case))) + with pytest.raises(ProtocolError, match=message): + read_result(tmp_path, case, routing_probe, 0.0) + + +def test_routing_time_axis_must_match_forcing(case, routing_probe, tmp_path): + table = _routing_table(case) + table.loc[table["reach_id"] == "A", "time"] = pd.date_range( + "2001-01-01", periods=case.n_steps, freq="D" + ).strftime("%Y-%m-%d") + _write_outputs(tmp_path, case, table) + + with pytest.raises(ProtocolError, match=r"routing\.csv reach 'A' time axis"): + read_result(tmp_path, case, routing_probe, 0.0) + + +@pytest.mark.parametrize("bad_value", [np.nan, np.inf, "not-a-number"]) +def test_routing_values_must_be_finite( + case, routing_probe, tmp_path, bad_value +): + table = _routing_table(case) + table["reach_q_in"] = table["reach_q_in"].astype(object) + table.loc[0, "reach_q_in"] = bad_value + _write_outputs(tmp_path, case, table) + + with pytest.raises(ProtocolError, match="non-finite"): + read_result(tmp_path, case, routing_probe, 0.0) + + +def test_verify_adapter_skips_routing_on_nonrouting_probe( + base_probe, tmp_path +): + model = replace( + registry.find_model("reference_bucket"), + emits_routing=("reach_q_in",), + ) + + result = verify_adapter_contract( + model, base_probe, 11, workdir=tmp_path + ) + + assert result.routing is None + + +def test_verify_adapter_checks_routing_on_network_probe( + routing_probe, tmp_path, monkeypatch +): + import hydroturing.harness as harness + + original_build_case = harness.build_case + + def build_case_with_network(*args, **kwargs): + built = original_build_case(*args, **kwargs) + static = dict(built.static) + static["routing_network"] = {"reaches": ["A"]} + return replace(built, static=static) + + monkeypatch.setattr( + harness, + "build_case", + build_case_with_network, + ) + model = replace( + registry.find_model("reference_bucket"), + emits_routing=("reach_q_in",), + ) + + with pytest.raises(ProtocolError, match="routing.csv"): + verify_adapter_contract( + model, routing_probe, 11, workdir=tmp_path + ) + + +def test_verify_adapter_prefers_routing_probe( + base_probe, routing_probe, monkeypatch +): + from hydroturing import cli + + model = replace( + registry.find_model("reference_bucket"), + emits_routing=("reach_q_in",), + ) + attempted = [] + + def stop_after_selection(model, probe, seed, **kwargs): + attempted.append(probe) + raise RuntimeError("stop after selecting the probe") + + monkeypatch.setattr( + cli.registry, + "find_model", + lambda name: model, + ) + monkeypatch.setattr( + cli, + "_all_probes", + lambda args: [base_probe, routing_probe], + ) + monkeypatch.setattr( + cli, + "verify_adapter_contract", + stop_after_selection, + ) + args = SimpleNamespace( + model="reference_bucket", + probe=None, + workdir=None, + window=None, + csv=None, + ) + + assert cli.cmd_verify_adapter(args) == 2 + assert attempted == [routing_probe]