Muskingum-Cunge physics corrections behind ddr_match (+0.0296 NSE over baseline) - #33
Merged
Conversation
Builds gages_2000_area_balanced.csv (1,841 gauges) from GAGES-II per the 2026-08-02 small-basin-domination handoff: DA_VALID recomputed as ABS_DIFF/DRAIN_SQKM <= 10%, >=80% obs coverage in both train and eval windows, non-headwater subgraph required, all >=5,000 km2 basins kept. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Documents the gages_2000_area_balanced.csv method: relative DA_VALID, coverage-in-both-windows filter, headwater exclusion, area-balanced subsample, and the baseline-invalidation consequences of changing the gauge population. Verified via fresh-context retrieval test. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Adds `probe_n_slope`, which answers "does the routing objective actually
want a stronger n-vs-basin-scale relationship?" with forward passes only.
Training an answer costs ~13 h per attempt; this costs ~0.6 h per field.
Motivation: across the 2026-07-30..08-02 CONUS series every trained head
landed on a nearly scale-independent n field with a faint positive tilt
against log10_uparea (Spearman +0.205 best run, +0.076 gentlest, -0.230 at
init). Two explanations fit equally well — under-trained, or the flat field
is already this objective's optimum — and they imply opposite next steps.
Reuses the existing `EvalParams::Frozen` path, so no training code is
touched. Reads the trained field from a checkpoint, rewrites n, re-scores.
p_spatial/q_spatial stay at trained values so the slope is the only
manipulated variable.
Two modes, because the first is provably insufficient alone:
--amplify K scales the trained field's log10_uparea-linear component.
Useless when the field has collapsed onto a bound: on the
lr=1e-2 run (47.3% of reaches at the 0.015 floor) K=5 moved
median n by 0.9% and dNSE by 1.6e-5, because clamping eats
the downward half and the residual linear term is tiny.
--impose L:H builds n from scratch as a rank-monotone ramp in
log10_uparea, independent of where training landed. This is
the mode that actually tests the hypothesis.
Plus `flat` (trained median everywhere) and `floor` (range minimum = identity
routing) controls, both on by default.
Two caveats documented in the module header:
- `EvalParams::Frozen` has no cross-chunk state injection (eval.rs), so it
cold-restarts each chunk. `slope x1` therefore will NOT reproduce a run's
headline NSE — every field goes through the identical path, so use it as
the internal reference, never the manifest number.
- Full 15-year eval is ~4.5 h/field (measured: phase2 4.23 h and 5.00 h in
the two most recent runs), so the window defaults to two water years.
Validated end-to-end on CPU against the lr=1e-2 run: all four fields scored,
table and CSV emitted. First `--impose` result is still pending.
Co-Authored-By: Claude <noreply@anthropic.com>
Four corrections gated behind params.ddr_match (default true = DDR parity): trapezoidal celerity beta, Cunge-derived Muskingum X, negative-discharge instrumentation, and Courant sub-stepping. Each with hand-derived backward branches and finite-difference gradchecks. Co-Authored-By: Claude <noreply@anthropic.com>
Adds params.ddr_match (default true). No behavioural change: the flag is threaded to forward_chain_inner and SavedState but not yet branched on. Later tasks use it to select corrected trapezoidal celerity and Cunge X. Co-Authored-By: Claude <noreply@anthropic.com>
…param Review of 461761d found two hazards to foreclose before physics branches land: - cuda_graph/geometry_kernel.rs hardcodes DDR's 5/3 celerity, so ddr_match=false + use_cuda_graphs=true would give a corrected backward against a DDR forward -- a silent gradient mismatch the planned CPU gradchecks would not catch. Now rejected at load, mirroring validate_leakance. - forward_chain_inner already takes cfg; the added positional bool was a round trip that risked desyncing from the cfg used to build TimestepState. Also records that forward_chain_inner_pinned is intentionally DDR-only until revived, so a later task does not silently diverge. Co-Authored-By: Claude <noreply@anthropic.com>
Completes the flag added in c09bb0e: the two Experiment literals in baseline/cache.rs need the new field, and the load-path test was untracked. Without these, cargo test --lib fails to compile at HEAD. Co-Authored-By: Claude <noreply@anthropic.com>
The S28 clamp_min(1e-4) silently rewrites negative solve output to +1e-4, creating mass and hiding Courant instability. Muskingum coefficients are only non-negative for 2X <= Cr <= 2(1-X); with X = 0.3 that window is [0.6, 1.4], and 69.8% of CONUS reaches fall outside it at mean flow. This counts the negatives before the clamp and reports them per forward. No numerics change, identical in both ddr_match modes. Co-Authored-By: Claude <noreply@anthropic.com>
Three fixes to the negative-Muskingum-solve diagnostic introduced in e5a513b: 1. CUDA-graph path (the shipped default: use_cuda_graphs=true in config/merit_training.yaml) never enters forward_chain_inner, so both counters stayed at zero — indistinguishable from "measured, found zero". MuskingumCunge::forward now prints an explicit UNAVAILABLE notice when graphs were requested, so silence can never be misread as a zero-negative measurement. 2. The per-timestep host sync (primitive_to_vec → Tensor::to_data) is now opt-in, mirroring enable_zeta_accumulation. A new enable_negative_discharge_tracking setter (default OFF) gates the device→host transfer; when off the forward pays zero added cost. The training driver (forward_with_kan, src/training/forward.rs) enables it so the per-micro-batch count still appears. Forward-chain helpers (__spike_*) and graph fallback pass false. 3. The test previously reset and asserted zero — it would have passed if the entire counting block were deleted. The rewritten test adds: - Case B (positive control): very short/steep/low-n reaches with large q_t and tiny q_prime_t drive c3 < 0 and produce actual negative x_sol values before the S28 clamp. Asserts neg_b > 0. - Case C (negative control): long reaches / gentle slope / typical flow keep c3 >= 0. Asserts neg_c == 0. Both cases share one #[test] function to avoid process-global counter races (no --test-threads=1 required). Co-Authored-By: Claude <noreply@anthropic.com>
S17 computed `celerity = velocity_cl * 5/3`. The 5/3 Kleitz-Seddon factor
is dQ/dA for a *wide rectangular* channel under Manning, but S7-S13 build a
trapezoid (top_width, side_slope, bottom_width, area, wetted perimeter).
For the sections this solver produces, kappa = b/y ~ 0.7-1.8, so the true
ratio is ~1.30-1.36 and the hardcoded 5/3 is 22-27% too high — the routed
wave travels that much too fast.
Under `ddr_match: false`, S17 now uses the exact trapezoidal celerity
c = dQ/dA = (dQ/dy)/T = v * beta
beta = 5/3 - (4/3)*A*sqrt(1+z^2)/(T*P)
which reduces to 5/3 as b/y -> inf (wide rectangular) and to exactly 4/3 as
b -> 0 (triangular). `ddr_match: true` (the default) is untouched, so
examples/compare_ddr_sandbox still reports ABSOLUTE MATCH
(max abs diff 1.53e-5 m3/s).
The hand-written backward (invariant 4 — no autograd-tape unrolling) gained
five terms. B17's dcelerity/dvelocity_cl becomes `beta` instead of the
constant 5/3, and with G = 5/3 - beta = (4/3)*A*u/(T*P), u = sqrt(1+z^2):
dbeta/dA = -G/A dbeta/dT = +G/T
dbeta/dP = +G/P dbeta/dz = -G*z/(1+z^2)
Each is ADDED to the accumulator that already carries that quantity's
gradient through the hyd_radius chain: gA before the S12 decomposition, gP
before S13, gT into gtw_total, gz into gss_combined. The wp/area
recomputation moved above B18 so B17 can read it; B14 reuses the same
tensors.
tests/celerity_beta.rs covers both halves: f64 helpers assert the wide
rectangular and triangular limits, agreement with a finite-difference dQ/dA
to 1e-6 on four realistic (b, z, y) triples, and a regression guard on the
22-27% overshoot; then four finite-difference gradchecks (n, q_spatial,
p_spatial, q_t) with ddr_match = false, worst relative error 1.4e-3 against
a 5e-3 tolerance.
Note the physics comment in the task description said beta is bounded below
by the triangular 4/3. It is not: beta is non-monotone in b/y and decays
toward 1 for narrow near-rectangular sections (b/y -> 0, z -> 0). The test
asserts the true bound, beta in (1, 5/3).
forward_chain_inner_pinned (the CUDA-graph capture path) is deliberately
left DDR-only; config load already rejects ddr_match: false with
use_cuda_graphs: true.
Co-Authored-By: Claude <noreply@anthropic.com>
The Muskingum storage weight X arrived as a caller-supplied constant
(forward.rs sets 0.3). A constant X severs the link between the scheme's
NUMERICAL diffusion and the channel's PHYSICAL hydraulic diffusivity --
the entire point of Muskingum-*Cunge*. It is Cunge-optimal only on the
measure-zero locus Q/(B*S*c*L) == 0.4; at a CONUS-median reach it
over-diffuses ~30x (documented median 28x).
S19 now branches:
ddr_match=true -> the caller's constant (unchanged; invariant 1)
ddr_match=false -> X = clamp(0.5*(1 - Q/(B*S*c*L)), 0, 0.5)
which makes D_num = c*L*(0.5-X) equal D_phys = Q/(2*B*S).
Backward (B19). X feeds only two_kx = 2k*X and two_k_1mx = 2k*(1-X), so
gX = 2k*(g_2kx_total - g_2k1mx_total). With W = Q/(B*S*c*L) and
X_raw = 0.5*(1-W) this adds THREE hand-derived terms, each masked to zero
where the [0, 0.5] clamp saturates:
dX/dQ = -0.5*W/Q dX/dB = +0.5*W/B dX/dc = +0.5*W/c
S and L are constants, not parents. q_t therefore gains a FOURTH gradient
path (S25 RHS, S24 SpMV, S2 depth chain, and now S19). dX/dQ is written
in its cancelled form -0.5/(B*S*c*L) because W/Q is 0/0 when the raw
per-step q_t is zero.
Ordering is the hazard: dX/dc must join gcelerity before B17 consumes it,
and dX/dB must join gtw_total before B7 consumes it. Both are done at the
existing fold-in points alongside the Task 3 beta terms.
x_effective is saved on TimestepState via a new forward_chain_inner
out-param rather than the forward_saved_idx array, which is index-locked
to cuda_graph::PersistentScratch. The CUDA-graph path is ddr_match=true
only by config validation, so x_effective == x_storage there.
Verification:
compare_ddr_sandbox ABSOLUTE MATCH (max abs 1.53e-5 m3/s)
gradcheck worst rel n 2.7e-4, q_spatial 1.0e-3,
p_spatial 2.3e-3, q_t 2.8e-4
falsification all four gradchecks FAIL with the B19
terms disabled (the 1000 m fixture from
celerity_beta.rs saturates X and would
have made them vacuous; cunge_x.rs uses
5000 m and guards this explicitly)
Negative-solve rate on config/experiments/gradaccum_smoke.yaml (cpu),
per mini-batch, before the S28 clamp:
ddr_match: true 0.004% 0.006% 0.007% 0.012%
ddr_match: false 0.027% 0.082% 0.025% 0.048%
The ~7-14x rise is EXPECTED, not a regression: Cunge X ~ 0.49 narrows the
non-negative-coefficient window 2X <= Cr <= 2(1-X) from [0.6, 1.4] to
roughly [0.98, 1.02]. This is the documented reason Task 5 (Courant
sub-stepping) exists; Tasks 4 and 5 must be evaluated together.
Co-Authored-By: Claude <noreply@anthropic.com>
…bled Replaces experiment.use_frozen_kan_head (ce26b77) — misleadingly named (it toggled the disaggregation head, not a frozen main KAN head, and "frozen" collided with disaggregation.freeze) and living in the wrong section. The switch is now `enabled: true|false` inside the block it controls; default true preserves the presence-enables contract for every existing config. `enabled: false` strips the block at load, giving flat repeat-24 (nearest) upsampling — the one-line ablation switch. No run ever consumed the old key, so no migration is needed. Also adds #[serde(deny_unknown_fields)] to DisaggregationSection: phantom keys (e.g. use_precip, removed in 334f0fe) previously silently built a different head than intended; they now fail load with "unknown field". All tracked configs carrying a disagg block verified clean. Tests: tests/disagg_enabled.rs (default-on, off-strips, on-keeps, unknown-key-rejected). Full suite green; compare_ddr_sandbox ABSOLUTE MATCH. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…cs result First run with ddr_match: false beats the summed-Q' baseline by +0.0296 median NSE on 1,841 area-balanced gauges (0.6736 vs 0.6440) -- 3x the previous best margin -- with the most physically defensible parameter field yet: median n 0.0467, 76.3% inside the NLCD natural-channel band, and the strongest scale dependence recorded (rho(n, log10_uparea) = +0.323). Records the two corrections (trapezoidal celerity beta; Cunge-derived X), the negative-discharge instrumentation and its first measurement, and Courant sub-stepping as attempted-and-abandoned with the evidence -- a 44x spread in K cannot be forced into a 4%-wide window by a global timestep, and variable dx (reach subdivision) is the real fix. States plainly that four variables changed at once, so the gain is NOT attributable to the physics without the matched ddr_match: true control. Co-Authored-By: Claude <noreply@anthropic.com>
…: false)
`collate::compress` built `outflow_idx` — the reaches whose routed
discharge is summed to form a gauge's prediction — by filtering the COO
for edges whose ROW is the gauge reach and taking the COLS. In this
adjacency `indices_0` = rows = DOWNSTREAM segment and `indices_1` = cols
= UPSTREAM segment (`src/data/store/zarr.rs:39-42`), so the prediction
was the sum of the gauge's UPSTREAM NEIGHBOURS, silently excluding the
gauge's own reach and therefore that reach's entire local drainage.
A USGS gauge measures ALL drainage above it, and we do not know where
along its reach the gauge physically sits, so the gauge reach's flow must
be included. Mass conservation means the MC solve at that reach already
accumulates every upstream contribution PLUS its own lateral inflow —
the correct answer is the single index `[gauge_compressed[g]]`.
Evidence, gauge 01457000:
- DRAIN_SQKM 366.8; the gauge reach's own COMID_UNITAREA_SQKM 250.1
= 68% of the basin
- subgraph order [73006562, 73006585, 73005764], gauge COMID 73005764
(the outlet); edges indices_0 = [229475, 229475] (both the gauge
reach), indices_1 = [228316, 228315] (two headwaters)
- observed mean 7.60 m3/s; summed-Q' baseline (which sums ALL subgraph
divides) 7.38; ddrs routed 1.58
- a constant 0.166-0.230 suppression across all 15 eval years, on peaks
as well as means
26 of 1841 gauges fell below 0.5x baseline, all small (139-453 km2, 3-9
reaches) where the gauge reach is a large area fraction; median ratio
across all gauges 0.952. Because the omitted mass is always positive,
this biased EVERY ddrs-vs-baseline comparison against ddrs — the
benchmark understated the routing's skill.
This is a faithful port of a DDR defect
(`~/projects/ddr/src/ddr/geodatazoo/merit.py:226-234`), so the correction
is gated behind the existing `params.ddr_match` flag:
- `ddr_match: true` (default) — DDR's upstream-cols convention, with
its empty-cols fallback to the gauge's own index. Preserved exactly.
- `ddr_match: false` — the gauge's own reach.
DDR's Lynker path validates `outflow_idx` against the flowpath `toid`
column (`lynker_hydrofabric.py:239-250`) while the MERIT path does not;
that is where this would have been caught upstream.
Plumbed as `compress(unioned, conus_order, ddr_match)`, with
`MeritGagesDataset` capturing `cfg.params.ddr_match` at open() and
passing it at both production call sites — mirroring the existing
`want_gauge_std` precedent.
Gates:
- tests/gauge_mass_conservation.rs (new): drives a constant q_prime to
steady state on the 01457000 topology and asserts the gauge sees the
whole network's lateral inflow. ddr_match=false gives 30.0000 m3/s
against an expected 30.0000 (ratio 1.0000); ddr_match=true gives
20.0000 (ratio 0.6667 — exactly the gauge reach's own 10 m3/s
dropped). The gauge reach itself carries 30.0 in both, so the defect
is in extraction, not routing.
- collate.rs::outflow_idx_includes_the_gauge_reach_when_not_ddr_match
asserts both conventions; the pre-existing outflow_idx tests are kept
and re-commented to say they pin the ddr_match=true path.
`outflow_idx` is downstream of the solver and is never built by
`examples/compare_ddr_sandbox`, which still reports ABSOLUTE MATCH
(max abs diff 1.526e-5 m3/s). cargo test --lib: 261 passed.
Co-Authored-By: Claude <noreply@anthropic.com>
…heck Adds a ddrs-eval-plots reference that turns learned n/p/q + post-clamp slope into physical width and depth at baseflow for all 346,321 CONUS reaches, maps them, and tests the geometry parameterisation against downstream hydraulic geometry. The attributes NetCDF has no width or depth variable, so this exponent check is the only internal validation available for p_spatial and q_spatial. Measured on run 2026-08-03T13-11-00Z: width exponent b = 0.226 against Leopold & Maddock's ~0.50, depth f = 0.600 against ~0.40 -- channels modelled too narrow and too deep, increasingly so downstream. Records that this is STRUCTURAL, not a training failure: with w = p*d^q and d proportional to Q^(3/(5+3q)), the width exponent is 3q/(5+3q), which maxes at 0.375 when q = 1 -- below L&M's 0.50 for ANY admissible q. Reaching 0.5 needs q > 1 (outside parameter_ranges) or a p that grows fast enough with Q. Also notes b+f+m = 1 is an identity, not a validation, and that the exponents are invariant to the assumed baseflow specific discharge while absolute widths are not. Co-Authored-By: Claude <noreply@anthropic.com>
Adds `params.enforce_positivity: bool` (default false), the gate for the Muskingum positivity clamp: floor K at dt(1+d)/2 and cap X at min(0.5*Cr, 1-0.5*Cr)*(1-d) so 2X <= Cr <= 2(1-X) holds on every reach-timestep, hence c1,c3 >= 0 and the S27 forward substitution can never go negative. The clamp targets the INPUTS, not the coefficients: c1+c2+c3=1 holds for any (K, X), so mass is preserved exactly, while clamping c3 would break the partition identity. Rejected at load together with `ddr_match: true`, because raising K and lowering X changes forward output and would break compare_ddr_sandbox's ABSOLUTE MATCH (invariant 1). `ddr_match: false` already forbids `use_cuda_graphs: true`, so the CUDA-graph path is excluded transitively. Plumbing mirrors `ddr_match` exactly: Params field + Default, ParamsRaw Option<bool> + From, and a `validate_enforce_positivity` sibling called from the same site as `validate_ddr_match`. Co-Authored-By: Claude <noreply@anthropic.com>
Muskingum's coefficients are non-negative exactly on 2X <= Cr <= 2(1-X)
(Cr = dt/K). c2 and c4 are positive unconditionally; c1 >= 0 <=> Cr >= 2X and
c3 >= 0 <=> Cr <= 2(1-X). Since S27 is forward substitution in topological
order, x[i] = b[i] + c1[i]*sum_up x[j] with b = c2*(N q_t) + c3*q_t + c4*q',
and q_t, q' > 0, c1,c3 >= 0 makes every x[i] >= 0 by induction. So enforcing
the window drives "negative solves before clamp" to exactly zero and S28's
clamp_min(1e-4) stops creating mass.
forward_chain_inner gains two flag-gated steps (enforce_pos = !ddr_match &&
cfg.params.enforce_positivity, so ddr_match:true is untouched):
S18' k_musk = max(L/c, dt(1+d)/2) => Cr in (0, 2/(1+d)]
S19' x_eff = min(x_cunge, hi_a, hi_b)
hi_a = 0.5*Cr*(1-d) (c1 >= 0)
hi_b = (1 - 0.5*Cr)*(1-d) (c3 >= 0)
The clamp targets the INPUTS, never the coefficients: c1+c2+c3 = 1 holds
identically for any (K, X), so mass is preserved for free, whereas clamping c3
would break the partition. The K floor is what guarantees hi_b > 0, so the
three-way min needs no clamp_min. d = POSITIVITY_DELTA = 1e-2 is mandatory:
at d = 0 the cap lands exactly on c1 = 0 / c3 = 0 and f32 roundoff crosses it.
A reach with K < dt/2 is sub-grid -- the timestep cannot resolve its transit,
and the unclamped scheme expressed that as oscillation S28 clamped to 1e-4
anyway. The floor makes that coarse-graining explicit.
tests/positivity_clamp.rs: 10-reach chain spanning Cr_raw 0.051 .. 7.899.
branch mix cunge 20% hi_a 30% hi_b 50%
K floor 5/10 floored (straddles)
negative solves 3/10 clamp OFF -> 0/10 clamp ON
min c1 / min c3 +2.54e-4 / +5.00e-5 (matches the d=1e-2 prediction)
off-parity x_sol bit-identical to the pre-change path (golden
captured from mmc_op.rs at fa5bcb4)
falsifiability forcing enforce_pos=false fails 4 of the 7 tests
compare_ddr_sandbox ABSOLUTE MATCH, max abs 1.53e-5 m3/s (unchanged)
fixture_is_not_vacuous is a standing guard modelled on the cunge_x.rs lesson
(a 1000 m fixture saturated X's clamp on every reach and made that gradcheck
pass with the backward deleted). It asserts all three min branches win, the
K floor straddles, Cr spans both sides of the window, and at least one Cunge
win is interior to X's own [0, 0.5] clamp.
NOTE: the backward is NOT yet updated (plan Task 3). Gradients under
enforce_positivity: true are wrong until B18'/B19' land; the flag defaults to
false and no shipped config sets it.
Co-Authored-By: Claude <noreply@anthropic.com>
B18'/B19' close the gradient hole Task 2 opened. The forward's two new steps
are both non-smooth, so each contributes a mask, and one contributes a path
that did not exist before:
S18' k_musk = max(k_raw, k_floor) -> gradient masked where floored
S19' x_eff = min(x_cunge, hi_a, hi_b) -> gradient to ONE branch, and
hi_a/hi_b depend on cr = dt/k_musk, so x_eff now reaches celerity
through K as well as through the Cunge W = Q/(B*S*c*L).
Accumulation order (mirrors the existing B18/B19 ordering constraint --
both new paths reach celerity through the SAME clamp_min, so the cap term
must land BEFORE the floor mask, and the floor mask before B18 divides by
celerity^2):
gk_musk = g_2k_total*2 (c1..c4)
+ gx_eff*[mask_a]*(+0.5(1-d)) * (-dt/k_musk^2) NEW
+ gx_eff*[mask_b]*(-0.5(1-d)) * (-dt/k_musk^2) NEW
gk_raw = gk_musk * mask(k_raw > k_floor) NEW
gcelerity = -gk_raw*length/celerity^2 + XGrads.g_celerity
XGrads (dX/dQ, dX/dB, dX/dc) is the CUNGE branch's, so all three components
are additionally masked by mask_cunge. The [0, 0.5] clamp mask likewise
belongs to x_cunge, not to x_eff: when hi_a/hi_b win, x_eff can sit in the
interior while x_cunge is saturated, and vice versa.
Tie-break priority for the three-way min, deterministic and total:
Cunge > hi_a > hi_b
mask_cunge = (x_cunge <= hi_a) && (x_cunge <= hi_b)
mask_a = !mask_cunge && (hi_a <= hi_b)
mask_b = !mask_cunge && !mask_a
so the three masks partition exactly. x_cunge/hi_a/hi_b are not saved; the
backward recomputes them from the saved length, slope, top_width, celerity,
q_t and k_muskingum with the same ops the forward used.
recomputed_min_branches_reproduce_the_forward_x pins that recomputation
against the X the forward's coefficients were actually built from (recovered
from the saved denom), on BOTH fixtures.
TimestepState gains `enforce_pos`, storing the already-combined
`!ddr_match && params.enforce_positivity` so the backward cannot drift from
the forward's gate. The CUDA-graph replay path hard-codes false (its captured
graph is ddr_match: true).
FALSIFIABILITY (the cunge_x.rs lesson: a 1000 m fixture saturated X's clamp
on every reach and made that gradcheck pass with the backward deleted).
Deleting each new piece in turn and re-running:
intact worst rel 3.3e-4 (n) / 3.2e-3 (q_spatial)
2.3e-4 (p_spatial) / 7.2e-4 (q_t)
two new gk_musk terms cut worst rel 3.9e-1 -- all 4 FAIL
mask_cunge cut worst rel 1.0e0 -- all 4 FAIL
K-floor mask cut worst rel 1.0e0 -- all 4 FAIL
grad_fixture(): same 10 reaches as fixture() with a GRADED q_t. Flat q_t made
reaches 6..9 gradient-dead -- not by masking but because c1+c2+c3 = 1 makes
q_next ~ q_t when x_up ~ i_t ~ q_t, which cancelled every S18'/S19' effect
into f32 noise and would have left the hi_a branch untested. fixture() itself
is untouched (GOLDEN_X_SOL_BITS is keyed to it).
branch mix cunge 20% hi_a 30% hi_b 50%
K floor 5/10 floored (straddles); Cr_raw 0.035 .. 7.899
resolved 6/10 reaches (10/10 for q_t) above the FD noise floor
Two f32 hazards had to be handled explicitly rather than papered over:
* the sibling gradchecks' absolute step max(1e-3*x, 1e-3) is a 2.9% swing in
n = 0.035, enough to walk reach 0 across the x_cunge/hi_b boundary
(0.0044 vs 0.0099); central differences then average two slopes and
disagree by 30%. REL_STEP is a pure relative step instead.
* an O(600) unweighted sum() puts l_plus - l_minus at the round-off floor,
so conditioning_weights rescales each reach to O(1) and compare_grads
states the surviving quantum (NOISE_ULPS*ulp(loss)/(2*eps)) explicitly
instead of hiding it in a magic ABS_TOL. It also asserts >= 4 reaches stay
resolved, so the noise gate cannot swallow the test.
Gate set: positivity_clamp 12, cunge_x 11, celerity_beta 9, sparse_gradcheck 1,
mmc 13, leakance_gradcheck 16, ddr_match_flag 5, --lib 261 -- all pass.
compare_ddr_sandbox ABSOLUTE MATCH, max abs 1.53e-5 m3/s (unchanged).
Co-Authored-By: Claude <noreply@anthropic.com>
…es on real data
Measurement harness for `params.enforce_positivity` (S18'/S19'). Drives
`forward_chain_inner` in the same q_next-fed-back loop as
`MuskingumCunge::forward`, so the numbers come from the production kernel
sequence, not a re-derivation:
* `negative solves before clamp` — exact over every routed timestep
(the same atomic counters the training log prints),
* `Cr = dt/K` before and after the S18' floor,
* `x_cunge` (pre-cap) vs `x_eff` (used), and the fraction where the cap binds,
* `k_musk / k_raw` over the floored reaches,
* `min c1` / `min c3`.
Result on the full 1841-gauge CONUS network (92,488 reaches, 2,135 hourly
steps, trained head from run 2026-08-03T13-11-00Z epoch_10_mb_27):
off = 55,181/197,461,880 (0.0279%), on = 0/197,461,880 (0.0000%).
`MuskingumCunge::probe_inputs` is the only lib change — a diagnostic
accessor bundle; no numerics touched.
Co-Authored-By: Claude <noreply@anthropic.com>
Correct the X_max cost table: evaluating the non-monotone X_max(Cr) at each Cr percentile does not give the percentiles of X. Measured on 1,841 gauges, the cap binds on 95.3% of reach-timesteps and median X falls 0.4976 -> 0.0794 (6.3x), not to 0.45 as first claimed. Also record: median Cr = 0.226 means the typical MERIT reach is ~4.4x too LONG for an hourly step, so the variable-dx fix is mostly subdivision, not merging. Co-Authored-By: Claude <noreply@anthropic.com>
…C Q' store Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Every physics gradcheck for the three `ddr_match: false` corrections (trapezoidal beta celerity, Cunge X, positivity clamp) declares `type I = NdArray<f32>` — CPU only. But `ddr_match: false` + `use_cuda_graphs: false` + a CUDA backend is a legal, actively-used configuration, so those hand-written `Backward<I, N>` impls shipped un-exercised on the hardware that runs them. `tests/cuda_backward_parity.rs` closes that gap: * native central-difference gradcheck with `Cuda<f32, i32>` as the inner backend, over both configs (beta+CungeX, and +positivity clamp) and all four learnable parents — correctness, not merely CPU agreement; * CUDA-vs-CPU analytic gradient parity at a derived 1e-5 tolerance (measured worst 2.784e-7 ~ 2 f32 ulp), with a zero-pattern assertion that catches a mask disagreeing across backends; * the same gradcheck through `sparse_solver: cuda`, the production configuration, which swaps in the cuSPARSE triangular solve; * per-op CUDA checks for `min_pair` / `clamp` / `greater_elem` / `lower_elem` / `lower_equal` / `bool_and` / `bool_not` / `mask_fill` / `powf` / `powf_scalar` / `recip` / `sqrt`, at the exact edges the physics hits; * the transitive `enforce_positivity => !use_cuda_graphs` guard, which nothing asserted before — it falls out of `validate_enforce_positivity` and `validate_ddr_match` separately and which one fires depends on `ddr_match`. Falsifiability is proven, not assumed. Non-vacuity guards run as a PRECONDITION of every gradcheck (beta must vary and differ from 5/3, side_slope must be off its clamp, Cunge X must be interior, all three min branches must win, the fixture must straddle the K floor), and two permanent negative controls pin the tolerances' sharpness. Five temporary mutations of `mmc_op.rs` were each confirmed to fail the suite and then reverted — numbers in `.claude/PHYSICS-CORRECTIONS.md`. Co-Authored-By: Claude <noreply@anthropic.com>
…ort ones Two-sided rule targeting dx = c_ref*dt so Cr ~ 1 network-wide. Splitting alone leaves 17.3% of reaches over-Courant (they are too SHORT); merging them is rejected because it destroys junction structure. Clamping their length instead drives frac(Cr>2) to 0.00% and frac(Cr<0.5) to 0.05%. Critically the clamp is a BUILD-TIME constant, unlike enforce_positivity's runtime K floor: no gradient path, so it cannot recreate the X ~ Cr ~ 1/n coupling that pinned n at its floor on 98% of reaches. Records an erratum: subdivision does NOT restore the Cunge X dynamic range (median moves only 0.4973 -> 0.4815). Uncapped subdivision is infeasible on solver critical path (9.2x), not memory. Co-Authored-By: Claude <noreply@anthropic.com>
…x limit K, X (with qo as unit-width discharge), and C0/C1/C2 all match ddrs exactly. Ponce also confirms the failure mode: too-large dx -> negative outflows, and negative C2 (our c3) -> dips on the rising limb. But the Ponce-Theurer limit dx <= (c*dt + qo/(So*c))/2 is REJECTED here. D ~ 0.012 on MERIT, so it collapses to ~c*dt/2 (C ~ 2), which violates the non-negativity ceiling C <= 2(1-X) and puts 57.3%% of reaches above Cr=2 -- worse, not better. C*D >= xi is unsatisfiable at any dx that also keeps the coefficients non-negative; it is a diffusion-routing criterion and MERIT is advection-dominated. Not a slope-floor artifact (1e-3 -> 1e-6 moves median D only 0.0112 -> 0.0120). Co-Authored-By: Claude <noreply@anthropic.com>
Adds the `params.subdivision:` block that will gate static reach
subdivision (variable dx) so Cr = c*dt/dx lands near 1 network-wide.
Config-only: nothing consumes it yet.
Fields and their reasons:
* `enabled` (false) — default off, byte-identical to today.
* `max_pieces` (8) — uncapped subdivision is infeasible (13.2x reaches,
9.2x solver critical path) AND unestimable: Sum(m) swings 2.3M-10.5M
across defensible reference-flow choices, while Sum(min(m,8)) is stable
to +/-12%. The cap bounds the cost and makes it measurable.
* `reference_n` (0.05), `reference_discharge_coefficient` (0.01),
`reference_discharge_exponent` (0.9) — the reference celerity that sets
the piece count is config-specified, deliberately NOT checkpoint-derived:
the graph must not depend on training state.
* `min_length_fraction` (1.0) — short reaches get their length clamped UP
rather than merged (merging would destroy junction structure). This is a
BUILD-TIME constant, unlike `enforce_positivity`'s runtime K floor, so it
has no gradient path and cannot recreate the X ~ Cr ~ 1/n coupling that
drove n to its bound.
`validate_subdivision` rejects `max_pieces: 0` and `enabled: true` together
with `use_cuda_graphs: true` — a captured graph is sized to a fixed reach
count and subdivision changes that count.
Plumbing mirrors `ddr_match`/`enforce_positivity`: Params field + Default,
ParamsRaw field + From, and a validator called from the same site.
ParamsRaw carries `Subdivision` directly rather than `Option<Subdivision>`;
its struct-level `#[serde(default)]` already yields `Subdivision::default()`
when the block is absent, which is exactly what `Params::default()` holds.
Co-Authored-By: Claude <noreply@anthropic.com>
Co-Authored-By: Claude <noreply@anthropic.com>
Task 2 review surfaced two holes. 1. min_length_fraction: 0.0 on a zero-length reach gave l_eff = 0, hence K = L/c = 0 and c1 = 1. MERIT has 11 sub-10 m reaches so this is reachable. Added ABS_MIN_LENGTH_M = 1.0 applied last, after the clamp. 2. The clamp was unbounded. reference_celerity uses r = Q_ref^0.4 with no slope dependence while v ~ sqrt(S), so steep small catchments got big-river depth AND steep-slope velocity — 8.9 m/s at slope 1e-2, i.e. a 32 km dx_target that would stretch every short steep headwater to 32 km. That is the source of the measured p99 = 36x / max = 48,597x clamp tail. Added max_clamp_factor (default 4.0) plus a tighter [0.05, 5.0] flood-wave celerity band. Bounding the clamp means the worst reaches stay over-Courant instead of being silently rewritten into 30 km channels. Task 8 must report that residual. Also fixes the Task 4 cache key in the plan: it hashed only enabled + max_pieces, but all six fields change dx_target and hence the built graph. Co-Authored-By: Claude <noreply@anthropic.com>
Task 3 of the reach-subdivision plan. `subdivide(&ConusAdjacency, &ReachPlan) -> SubdividedAdjacency` expands each parent reach into its planned pieces, chaining them upstream->downstream over a contiguous row block and rewiring every external edge to outlet(upstream) -> inlet(downstream). Lower-triangularity (routing invariant 3) survives: internal links run base+k-1 -> base+k, and for an external edge r > c we get inlet(r) = off[r] >= off[c+1] > off[c+1]-1 = outlet(c). The r != c assertion is backed by the real CONUS COO (346,321 reaches / 338,814 edges), which carries zero self-edges and zero upper-triangular entries. Nothing is wired up yet; this is pure preprocessing over plain Vecs. Co-Authored-By: Claude <noreply@anthropic.com>
Task 4 of the reach-subdivision plan. - `write_conus_store_subdivided` writes `/parent_order` (int32 [n_parent]) and `/parent_offset` (int32 [n_parent+1]) alongside the existing arrays. - `ConusAdjacencyStore` gains `parent_order`/`parent_offset` and builds its `IdIndex` from `parent_order` — `order` has duplicate COMIDs once subdivided, so a lookup on it would be ambiguous. Missing arrays synthesize the identity, so every pre-existing store (engine exports, caches built at BUILDER_VERSION 1) keeps loading unchanged. - `content_key` hashes all seven `Subdivision` fields after BUILDER_VERSION; each one feeds `reference_celerity` -> `dx_target` -> the built graph, so hashing a subset would silently reuse a stale adjacency. - BUILDER_VERSION 1 -> 2. - `resolve_or_build` runs `plan_reaches` + `subdivide` between `build_conus_adjacency` and the gauge subgraphs, so subgraphs are cut from the expanded graph and each gauge resolves to its parent's outlet piece. Disabled (the default) is an exact no-op and never opens the attributes NetCDF. Drainage area: `catchsize` is the LOCAL divide area, not upstream area, so it is accumulated downstream over the topological order. Verified against the fabric's own `log10_uparea` on real CONUS — ratio p5/p50/p95 all 1.000 over 346,321 reaches. (`log10_uparea` cannot be the source: NaN on 88% of the global attributes file.) Gates: subdivide 28 passed, ddr_match_flag 5, lib 262, adjacency_parity --ignored (release) passed unchanged, cargo check --all-targets clean. Co-Authored-By: Claude <noreply@anthropic.com>
A reach split into m pieces of length L/m must receive q'/m per piece: HEC-HMS's lateral term is C4*(q_L*dx) with q_L an inflow per unit length, and the pieces chain in series so the parent's outlet piece still carries the whole reach's runoff. Total q' per parent reach is conserved exactly. `SparseAdjacency` gains `parent_offset: Option<Vec<i32>>` (the sub-reach parent map written by `adjacency::subdivide`). `setup_inputs` turns it into a per-row divisor tensor once; `forward` divides q' by it immediately AFTER the `discharge_lb` clamp. Order matters: clamping first floors the parent's inflow once, whereas dividing first would floor each piece independently and inject m*discharge_lb on a dry reach. `pieces_per_row_divisor` returns None when no parent owns more than one row, so an un-subdivided network (or an identity map) skips the op entirely and stays bit-identical -- asserted by `identity_parent_map_is_bit_identical_to_none` and by compare_ddr_sandbox (max abs diff 1.53e-5 m3/s, ABSOLUTE MATCH). Measured: 1 piece = 10.000000 m3/s, 4 pieces = 10.000001 m3/s at steady state with q' = 10; 4 pieces without the parent map = 40.000000 m3/s. Co-Authored-By: Claude <noreply@anthropic.com>
Thread the reach-subdivision parent map from ConusAdjacencyStore through collate::compress and into SparseAdjacency, so Task 5's lateral-inflow split actually fires in production instead of seeing parent_offset: None. compress gains a `conus_parent_offset` argument and re-expresses it in COMPRESSED space: gauge_compressed holds compressed sub-reach positions, not parent indices, so indexing the CONUS parent map with one is a category error. A parent's pieces are a contiguous ascending run of CONUS rows and a gauge subgraph enters a parent at its outlet then walks the internal chain upstream, so every present parent must be present in full; that is asserted rather than assumed (a partial chain would give the engine m smaller than the one the piece lengths were cut for). With that map in hand the `ddr_match: false` branch resolves each gauge's parent and reads parent_offset[p + 1] - 1 — the outlet piece. Any earlier piece drops the downstream fraction of the reach's own lateral inflow: the same class of bug as 2fe6bee, in a new form. The `ddr_match: true` branch is untouched. Measured on the 3-reach fixture (two headwaters into the gauge reach, q' = 10 m3/s each): 1 piece = 30.000000 m3/s, 4 pieces = 30.000006. The gauge reach's four pieces form the ramp 22.5 / 25.0 / 27.5 / 30.0, so pointing outflow_idx at the inlet fails the test at 22.5. Co-Authored-By: Claude <noreply@anthropic.com>
Task 6 of docs/superpowers/plans/2026-08-05-reach-subdivision.md, plus the
length/slope indexing defect Task 7 flagged.
Fix: per-row geometry was read through the PARENT index
-------------------------------------------------------
`dataset.rs` sliced `conus.length_m` / `conus.slope` with
`conus.index.position(comid)`. Since Task 4 that index is built from
`parent_order` and returns a PARENT position, while the geometry arrays are
SUB-REACH arrays — so under subdivision every row got the geometry of an
unrelated parent-numbered row (all 1000 m / 0.001 on the chain3 fixture
instead of 1000/1000/1000/3000/3000/900 and three distinct slopes). Without
subdivision the two spaces coincide, which is what made it silent.
`compress` now returns `conus_positions` (the compression's inverse map) and
the new `slice_reach_geometry` indexes with it.
Gather: run the KAN at parent resolution
----------------------------------------
Sub-reaches inherit their parent's hydraulics — MERIT carries no within-reach
variation, so there is nothing better available. Rather than duplicate
attribute rows, the head runs once per MERIT reach and
`gather_params_to_subreaches` expands every output onto the routing rows via
`select`, whose backward is a scatter-add: each parent receives the SUMMED
gradient from all its pieces, the correct semantics for a shared parameter.
~5x less head compute at `max_pieces: 8`, identical numbers.
`CompressedAdj::parent_comids` supplies the parent COMID list that
`finalize_attrs` now slices, so `spatial_attributes` is `(N_parent, F)`.
Consistent index space elsewhere: `AttributesStore::open{,_multi}` and
`dump_parameters` now take `conus.parent_order` (the dump's netCDF COMID
dimension is a coordinate — duplicates would be malformed), and
`probe_forward` / `probe_n_slope` gather before consuming head outputs.
Disabled path is untouched: with one row per parent `parent_row_index`
returns `None` and the map is returned as-is — no `select`, no tape entry.
Verified
--------
- `compare_ddr_sandbox`: ABSOLUTE MATCH, max abs 1.53e-5 m3/s.
- Break-and-confirm, geometry: reverting to `index.position` makes
`subdivided_rows_get_their_own_length_and_slope` fail at row 3
(length 1000 != 3000); all six rows read 1000 m / 0.001.
- Break-and-confirm, gather: a no-op gather fails
`every_piece_inherits_its_parents_parameters` ([0.02, 0.05, 0.1] vs
[0.02, 0.02, 0.02, 0.05, 0.05, 0.1]) and `gradient_sums_back_to_the_parent`.
Severing the tape inside the gather leaves the forward test passing and
fails only the gradient test — so the two discriminate independently.
- `cargo test --tests --no-fail-fast`: 80 binaries, exit 0.
- `cargo test --lib` 266; subdivision_integration 7 (was 1); subdivide 28,
gauge_mass_conservation 4, mmc 13, sparse_gradcheck 1, kan_head 7 unchanged.
Co-Authored-By: Claude <noreply@anthropic.com>
…the builder Subdivision runs inside `adjacency::cache::resolve_or_build`, which `cli::plan::resolve_adjacency` only reaches when `data_sources` supplies no `conus_adjacency`/`gages_adjacency`. With those keys set — as the stock `ddrs.yaml` has them — `params.subdivision.enabled: true` was SILENTLY INERT: no subdivision, no warning, and a run whose manifest claimed subdivision while routing the un-split network. Same failure class as the 2026-07-01 stale-binary disaggregation 2x2. Config load now rejects the combination, naming the conflict and pointing at `geospatial_fabric`. The one legitimate case — an explicit path to a store that WAS built subdivided — is admitted: `adjacency::validate::store_is_subdivided` compares the declared shapes of `/order` and `/parent_order` (`n_parent < n`) from two small zarr.json files, so no array data is read. An identity parent map (managed build with subdivision off) is correctly treated as un-subdivided. tests/subdivide.rs: 28 -> 33. Co-Authored-By: Claude <noreply@anthropic.com>
probe_courant gains `--fabric` (switch to the managed adjacency build,
the only path that subdivides), `--max-pieces`, `--reference-n`,
`--min-length-fraction`, `--clamp-report`, `--divide-hotstart` and
`--trace-steps`, plus `frac Cr < 0.5`, `frac c1>=0 AND c3>=0`, the
non-negativity window width `2(1-2X)`, and per-arm timing.
Two correctness fixes the measurement exposed:
* `ProbeInputs` now carries `pieces_per_row`, so the probe routes `q'/m`
exactly as `forward` does (clamp first, then divide). It was feeding a
subdivided network m-times too much water.
* probe_courant gathers KAN outputs parent -> sub-reach via
`gather_params_to_subreaches`, as `training::forward` does.
HOT START (plan Task 8 Step 4b): `setup_inputs` cold-started with the
UNDIVIDED `q'_0`, so under subdivision every parent outlet began at m x
its true steady state. Measured on 1,841 CONUS gauges / 184,676 rows at
`max_pieces: 8`: 2.94x total network discharge, decaying below 10% only
after 221 hourly steps and below 5% after 282 — against `warmup` = 5 d =
120 steps, where it was still 41.7% off. So `divide_hotstart_by_pieces`
defaults to true. Exact no-op without subdivision (no divisor exists),
so `compare_ddr_sandbox` stays an ABSOLUTE MATCH at 1.53e-5 m3/s and
`adjacency_parity` still passes element-for-element.
`subdivide::plan_stats` + `cache::{reach_plan, parent_adjacency_from_fabric}`
expose the reach-plan cost (piece histogram, clamped fraction,
clamp-factor percentiles, length inflation) without writing a store.
tests/subdivision_integration.rs: 7 -> 12. The plan's sketched
`subdivision_makes_coefficients_non_negative_without_the_clamp`
(frac c1<0 < 1%) is deliberately NOT added: the measurement refutes it —
frac c1<0 goes 93.0% -> 98.8% at cap 8. The added
`both_coefficients_are_non_negative_only_inside_a_window_that_collapses_at_x_half`
records why: both coefficients are non-negative only for
`2X <= Cr <= 2(1-X)`, a window of width `2(1-2X)` measuring 0.013 at the
CONUS-median X of 0.4966. A build-time piece count cannot land a
flow-varying Cr inside a 1-3% window.
Co-Authored-By: Claude <noreply@anthropic.com>
Task 9 of the reach-subdivision plan, with the conclusion the Task-8 measurement actually produced. `.claude/REACH-SUBDIVISION.md` (new): STATUS: NO-GO for "non-negative by construction". Both Muskingum coefficients are non-negative only inside `[2X, 2(1-X)]`, a window of width `2(1-2X)` — 80% at X=0.3 but 1.4% at the measured CONUS median X=0.4966. A build-time piece count fixes dx from a reference flow while Cr tracks the routed celerity, which varies severalfold within a storm, so no cap lands Cr inside that window. Measured (1,841 gauges, enforce_positivity OFF): frac c1<0 gets WORSE, 93.0% -> 98.79% at cap 8; negative solves fall only 35% for a 2.05x network, 1.5x step time and +23.9% total channel length. What it does buy — Cr>2 2.10% -> 0.16%, c3<0 3.93% -> 0.31% — comes from the short-reach length clamp, not the splitting. Also records: the plan's cap-sweep cost table is materially wrong for the shipped `reference_n: 0.05` (Sigma m 709,974 = 2.05x, not 4.77x; 51.0% clamped, not 34.4%) because the trained median n is 0.130; the root cause X -> 0.5 (cell Reynolds D ~ 0.012, advection-dominated), which retroactively vindicates the constant X = 0.3; the hot-start decision (undivided cold start = 2.94x discharge, 41.7% off at the configured warmup of t=120, so divide_hotstart_by_pieces defaults true); and the design that stays in-tree — two-sided rule, two index spaces, per-file plumbing, and the catchsize/cache-key/silent-inertness gotchas. PHYSICS-CORRECTIONS.md: the "sub-stepping does not substitute" section claimed subdivision would fix this. Replaced with the measured refutation in the file's existing erratum style, plus a second erratum — subdivision does not restore the Cunge X's dynamic range either (raw X median moves only 0.4973 -> 0.4815 even uncapped). Added the `params.subdivision` row to the blast-radius table (no gradient path). CLAUDE.md: new `params.subdivision` section — default off, what it does, the NO-GO status, and the three enable-time requirements (`geospatial_fabric`, `use_cuda_graphs: false`, retraining), including why enabled + explicit adjacency paths is a config error rather than a warning. Co-Authored-By: Claude <noreply@anthropic.com>
… key Corrections found during execution: catchsize is LOCAL divide area (median 36.7 km2), not upstream drainage — must be accumulated downstream; the cache key must hash all seven Subdivision fields; and the hotstart division question carried from Task 5 is now a Task 8 step. Co-Authored-By: Claude <noreply@anthropic.com>
DDRS_HOURLY_DUMP=<path> writes the (n_gauges, n_hours) f32 hourly series before tau_trim_and_downsample so params.tau can be swept exactly offline. Spec: docs/superpowers/specs/2026-08-05-per-gauge-tau-sweep-design.md Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Reads the DDRS_HOURLY_DUMP raw f32 + eval predictions.zarr, sweeps tau 0..23 (block-mean pooling, obs day i+1 convention), reports per-gauge NSE(tau) with area-bin medians, baseline comparison, and longitude correlation. Pilot restricted to WY1996 per the spec. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
All three method gates PASS. NSE(tau) plateaus at tau=14-19 (local-midnight pooling) in every bin < 30,000 km2; global tau=19 gains +0.114 median NSE (1,841 gauges, WY1996) and tau=16 ties the summed-Q' baseline in the <1,000 km2 bin. Skill references updated per the maintenance rule; sweep script gains an all-NaN-curve guard. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…-midnight phase Fable subagent review (sound-with-corrections): disagg was OFF in the swept run, hourly signal 97.6% UTC-day-constant, so the sweep resolves day-pairing + blend weight, not sub-daily phase. Robust covariate is drainage area (accumulated-lag signature); longitude fingerprint sign-flips uncensored. 596/661 tau=23 pins are real optima beyond the sweep edge. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…pling Day-center interpolation with +/-1-day context reads so 15-day eval chunks tile without edge clamps; nearest path byte-identical; global zarr-v2 reader rejects the env var loudly. Findings 5b records the AORC-UTC vs USGS-LST mechanistic prior (predicts tau in [16,19] a priori) and that src/data has no timezone logic anywhere. Adds the 3-arm gages_3000 driver + config. Gates: cargo test --lib (270), cargo test (0 fail), compare_ddr_sandbox ABSOLUTE MATCH (1.5e-5 m3/s). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…step-function artifact 3 arms (nearest/linear/quadratic) x 2,365 gauges: optimum tau=18-20 in every arm, curves neither sharpen nor shift, nearest >= linear >= quadratic at each optimum. Small-basin global tau=19 beats baseline 0.674 vs 0.645 (WY1996). Area-correlated late tail points at convention offset + double-routing lag. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
… constants by ~0.002 c=0.28 @ b=0.3 fitted on nearest-arm WY1996 curves. Area term absorbs the area signal; tz term does not track per-gauge optima (residual-vs-tz -0.45). Constant tau ~19-20 captures nearly all recoverable skill at daily-info resolution. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
… 2001 Reframe from selective-equifinality/bias-absorber thesis to an empirical answer to Beven's uncertainty and equifinality challenges, with a normative judgment-criteria payload. Full two-axis matrix (lumped/distributed dHBV, daily/hourly LSTM); open non-directional questions; keep Intro+Methods, rewrite Results/Discussion/Conclusion. Gates Results on timezone alignment, spatial-axis parameter analysis, a replicate-seed noise floor, and the double-routing confound. Co-Authored-By: Claude Opus 4.8 <noreply@anthropic.com>
User preference: deterministic NdArray for diagnostics, GPU stays free. Comment notes the no-mixed-backends-within-a-set rule; today's completed sets ran entirely on cuda. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…_lumped left-censored outlier Same epoch-30 checkpoint and 2,365-gauge network, streamflow store swapped. UH retro / daily LSTM / hourly-native LSTM peak within 2h of the aorc2f distributed reference (per-gauge best-tau median 18-19 on all four), so the lag lives in the shared forcing/obs day convention, not any runoff model. aorc2f_lumped peaks at tau~0-3 with 47% of optima pinned at tau=0: aligned ~1 day differently, CF day convention needs inspection. Hourly-native arm does not sharpen the curve but shows the strongest correct-sign longitude correlation (-0.220, censored subset). Adds tau_src_* experiment configs. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
Corrections: attribution softened (replication cannot separate day convention from common-mode MC routing); lumped CF-convention candidate refuted (metadata byte-identical to distributed — shift is in the written data, obs-free waveform lead ~23h vs NSE-optimum shift ~17h); hourly-arm longitude claim cut as a censoring artifact (interior sign flips to +0.075). New measurements: obs-free cross-arm hydrograph cross-correlation, extended tau sweep -13..47 (constant-tau optimum interior at 20; 33% of per-gauge optima beyond 23), censoring asymmetry, flatness-not-floor check. Proposed discriminator: tau sweep on routing-free summed q'. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…d, MC over-delays ~2x Summed daily q' repeat-24'd in the arms' phase, identical sweep machinery (tau=11 algebraically reproduces the cached day-aligned baseline NSE, median |diff| 0.0003). Optimum tau=6 vs routed 20. Per-bin: summed q' leads the gauge by an area-growing travel time (2h to 19h), routed lags by the mirror image, so MC adds ~2x the needed delay (the double routing of DDR's tau docstring, now measured). Zero-area intercept ~tau 10-11 puts the UTC-vs-LST convention offset at 0-2h, refuting §5b's prior as dominant and explaining every failed longitude fingerprint. Routed still beats summed at both optima in every bin. Extended curves + comparison plot under output/tau_sweep/. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…g structural in the store (§5h) dump_gamma_uh_params.py reads the CONUS2717_AORC2F_v3_gradaccum ep100 Ann head (the checkpoint behind the 2026-07-29 UH-routed re-export of the distributed aorc2f store) and computes a_eff*theta_eff per divide across all 197,088 exported divides. Median UH kernel mean 1.50 days (IQR 0.87-4.22), 71.6% of divides above 1 day, spearman +0.383 with log uparea: the per-divide lateral inflows carry area-dependent network-scale travel time before MC routes anything, confirming §5g's measured 2x over-delay at the parameter level. Clean fix target is a sub-grid-only UH on lateral inflows (col-7 unrouted is not it: 0.29 routed). Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…↔ obs day i, default 9 BREAKING: params.tau is now the number of hours the routed output is advanced before daily scoring (dMC-Juniata's sign; 0 = day-aligned; valid [0,24)). Legacy scale = new + 11: old shipped 3 ≡ new -8 (wrong-direction shift), old measured optimum 20 ≡ new 9 (the new default). Pooled day i is scored against obs day i (was i+1): driver obs pairing updated, eval drops the first pooled day instead of the last (zarr day axis unchanged), probe binary aligned. tau_sweep.py moved to the new axis (start 24+tau vs zarr days, TAUS -12..23); pre-change dumps need the old script from history. DDR-parity fixture tests pin the legacy window explicitly (legacy tau=3 cuts the same hours as new tau=16) — fixtures unchanged, all green. Equivalence unit test: new tau=9 ≡ legacy tau=20 offset one pooled day. Adds config/experiments/tau9_train_*.yaml (epoch-30 snapshot with only streamflow + sparse_solver:cpu changed) and run_tau9_source_trains.sh: five CPU train-and-test retrains, one per streamflow store. Findings §5i. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
…ats the baseline Three arms complete (30 epochs, 1,841 gauges, full test window): aorc2f distributed 0.620 -> 0.706 median NSE (old tau vs tau=9), UH retro 0.707, daily LSTM 0.578; summed-q' baseline 0.642. Adds the remaining-arms driver (aorc2f_lumped + hourly_lstm legacy-eval chain after the original driver was killed mid-eval), the failed zero-batch eval-resume config as a record, and the aorc_dhbv_distributed source group used by the campaign. Co-Authored-By: Claude Fable 5 <noreply@anthropic.com>
This file contains hidden or bidirectional Unicode text that may be interpreted or compiled differently than what appears below. To review, open the file in an editor that reveals hidden Unicode characters.
Learn more about bidirectional Unicode characters
Sign up for free
to join this conversation on GitHub.
Already have an account?
Sign in to comment
1 participant
Add this suggestion to a batch that can be applied as a single commit.This suggestion is invalid because no changes were made to the code.Suggestions cannot be applied while the pull request is closed.Suggestions cannot be applied while viewing a subset of changes.Only one suggestion per line can be applied in a batch.Add this suggestion to a batch that can be applied as a single commit.Applying suggestions on deleted lines is not supported.You must change the existing code in this line in order to create a valid suggestion.Outdated suggestions cannot be applied.This suggestion has been applied or marked resolved.Suggestions cannot be applied from pending reviews.Suggestions cannot be applied on multi-line comments.Suggestions cannot be applied while the pull request is queued to merge.Suggestion cannot be applied right now. Please check back later.
Corrects two physics defects in the Muskingum-Cunge core and instruments a third, all behind
params.ddr_match: bool(defaulttrue, so DDR parity andcompare_ddr_sandbox's ABSOLUTE MATCH are untouched).Full writeup:
docs/2026-08-03-ddr-match-findings.md. Dataflow diagram:.claude/PHYSICS-CORRECTIONS.md.What was wrong
Celerity used the wide-rectangular limit.
c = v · 5/3isdQ/dAfor a wide rectangular channel, but the solver builds a trapezoid. The exact kinematic celerity isc = v·βwithβ = 5/3 − (4/3)·A·√(1+z²)/(T·P). These channels haveκ = b/y ≈ 0.7–1.8, givingβ ≈ 1.30–1.36— so5/3was 22–27% too high.Muskingum X was a constant 0.3, not Cunge-derived. That severs the link between numerical and physical diffusion, which is the defining feature of Muskingum-Cunge. Corrected to
X = clamp(0.5(1 − Q/(B·S·c·L)), 0, 0.5), which makesD_num = c·L·(0.5−X)equalD_phys = Q/(2·B·S). The constant was injecting a median 28× excess diffusion.Negative discharge was clamped without ever being measured.
q_next = x_sol.clamp_min(1e-4)at S28 silently rewrote negatives to+1e-4, creating mass and zeroing gradients. Now counted before the clamp — first measurement: 0.004–0.012% of reach-timesteps on the DDR path, 0.027–0.082% with Cunge X. This corrects the audit's framing: ~70% of reaches carry a negative Muskingum coefficient, but only ~0.01–0.05% actually produce negative discharge.Result
First run with
ddr_match: false— 1,841 area-balanced gauges, Adam without gradient accumulation (280 optimizer steps in 10 epochs), disagg head off:Largest margin over the summed-Q' baseline this project has produced — 3× the previous best (+0.0099). Observations byte-identical between the two series (max diff 0.000e+00 over 10,026,921 cells).
Area-stratified, where the gain actually lives:
Learned parameter field over 346,321 CONUS reaches: median
n0.0467, 76.3% inside the NLCD natural-channel band 0.025–0.15, and ρ(n, log10_uparea) = +0.323 — the strongest scale dependence recorded, on all three learnable parameters.Caveats — please read before quoting the headline
Four variables changed at once versus every prior run: corrected celerity, Cunge X, disagg head off, and the rebalanced gauge set. The +0.0296 is not attributable to the physics on this evidence. The matched
ddr_match: truecontrol on the same gauge set has not been run and is the first follow-up.The baseline differs (0.6440 on this 1,841-gauge population vs 0.6754 on
gages_3000), so none of these absolutes compare to earlier runs.FHV moved +5.4 → −6.0 — peaks under-predicted ~6%. Mechanistically consistent: corrected celerity is ~20% lower, so
K = L/cis longer and the router attenuates more.Floor fraction rose to 6.56% from 0.73% in the 50-epoch run. Far below the 47.3% collapse, but worth checking whether the pinned reaches concentrate in small headwaters.
Courant sub-stepping: attempted, abandoned, documented
Cunge
X ≈ 0.49narrows the non-negative window2X ≤ Cr ≤ 2(1−X)to[0.98, 1.02]. Sub-stepping was planned to fix it and cannot work:Kspans 425 s to 18,551 s (44×), and even ideal per-reach integer sub-stepping only lands 6.3% of reaches in the window. The correct fix is variable Δx — reach subdivision, as HEC-HMS does — which changes adjacency topology. Recorded in the findings doc so it is not re-attempted.This also reframes
X = 0.3: window width is2 − 4X, so the constant was the most numerically forgiving choice available, plausibly a deliberate stability trade.Verification
compare_ddr_sandboxABSOLUTE MATCH (1.53e-5 m³/s) ·cunge_x11/11 ·celerity_beta9/9 ·sp8_gradcheck5/5 ·sparse_gradcheck1/1 ·mmc13/13 ·leakance_gradcheck16/16 ·zeta_accum8/8 ·cargo test --lib260 passed.Both new backward branches are falsifiable, not vacuous — they FAIL with their terms disabled (celerity rel 1.3e-1 vs 1.44e-3 passing; X rel 8.9e-2…5.7e-1 vs 2.84e-3). The Cunge-X gradcheck fixture was raised from 1000 m to 5000 m because at 1000 m
W ≈ 1.6saturates the clamp on every reach, maskinggXto zero and letting all four tests pass with the terms deleted.Config guard:
ddr_match: false+use_cuda_graphs: trueis rejected at load —cuda_graph/geometry_kernel.rs:296hardcodes DDR's 5/3, so a captured graph would replay a DDR forward against a corrected backward. That mismatch would be invisible to the CPU gradchecks.Not changed (found, deferred)
attribute_minimums.slope: 1e-3clamps 33.2% of reaches; it is an exact invariance absorbable by scalingnby√f. Undoing it puts 94.2% of implied physicalnin the NLCD band. Needs a retrain.leakance.rs:35-36uses(p·d)^qinstead ofp·d^q— dimensionally incoherent, inherited from DDR. Leakance is closed/NO-GO.The tau fix: routing was scored half a day out of phase
This branch also lands the diagnosis and fix for
params.tau, the hourly-to-daily pooling offset intau_trim_and_downsample. Full evidence chain indocs/2026-08-06-tau-sweep-pilot-findings.md(sections 1 through 5i).Diagnosis. The shipped
tau: 3scored the routed hydrograph with a pooling window roughly half a local day out of phase. ADDRS_HOURLY_DUMPinstrument plus an offline sweep (scripts/tau_sweep.py) showed the NSE-optimal offset is 17 to 20 on the legacy scale in every drainage-area bin, on every q' store tested (aorc2f distributed, UH retrospective, daily LSTM, hourly-native LSTM), while nothing insrc/data/or the stores explains it as a data or timezone convention:scripts/dump_gamma_uh_params.py). MC then routes the network on top of inflows that already carry network-scale delay.Fix.
tauis redefined on a signed-at-zero convention (commit54cd386): the slice is[tau : -(24-tau)], pooled dayiis scored against observation dayi, andtauis literally the number of hours the routed output is advanced before daily scoring, the same shift (same sign, same magnitude) as the dMC-Juniata paper's inverse-routing tau.tau = 0is exactly day-aligned; the legacy scale wasold = new + 11, so the old shipped 3 was a shift in the wrong direction (new −8) and the measured optimum 20 becomes the new defaulttau: 9. Old configs and checkpoints carry legacy-scale values and must not be reused verbatim. DDR-parity fixture tests pin the legacy window explicitly (legacy tau=3 cuts the same hours as new tau=16), so all parity gates stay green; an equivalence unit test proves new tau=9 reproduces legacy tau=20 exactly.Result. Retraining the identical config at
tau: 9(30 epochs, 1,841 gauges, full 1995-2010 test window):The timing correction alone is worth +0.086 median NSE and flips the routing from losing to the summed-q' baseline to beating it by +0.064. Two arms (hourly-native LSTM, aorc2f lumped) are still evaluating; the lumped store is a known outlier (its data is aligned about a day differently, see findings section 5e) and runs as the consistency control. tau remains a translation-only patch: the root cause (double-carried travel time in the exported q', plus the slow trained celerity) is documented for a follow-up store re-export.
🤖 Generated with Claude Code