feat(leakance): temperature-annealed gate replaces the continuous factor

.claude/skills/ddrs-dev/references/config.md

@@ -229,6 +229,7 @@ and bounds every topology from above.
 | `use_leakance` | false | |
 | `leakance_losing_only` | **true** | Clamps `head = max(0, depth − d_gw)`, so gaining reaches produce `zeta ≡ 0` |
 | `leakance_impervious_threshold` | 0.7 | Masks reaches whose `corridor_impervious` is **`>`** this value (not `≥`) |
+| `leakance_gate` | absent | Temperature-annealed 0/1 gate on the head's `leakance_factor` output; see §Leakance gate. Absent ⇒ byte-identical to the ungated readers. Requires `use_leakance: true` |
 | `tau` | 9 | **Not** a routing sub-step count. Since 2026-08-08: hours the routed output is ADVANCED before daily scoring (translation-only inverse routing, dMC-Juniata's sign). Slice `[tau : -(24-tau)]`, pooled day i ↔ obs day i, valid range [0, 24) (`src/training/loss.rs`). Nothing in `src/routing/` reads it. Default 9 = the measured CONUS optimum (findings §5g). **Legacy scale (pre-2026-08-08): old = new + 11**, slice `[13+tau : -11+tau]`, day i ↔ obs day i+1; DDR-Python still uses it, all older configs/checkpoints carry it (old shipped 3 ≡ new −8, wrong direction). Never copy a `tau:` value across the convention boundary. |
 | `log_space_parameters` | `["p_spatial"]` | |
 | `defaults` | `{p_spatial: 21.0}` | Value used when a parameter is not in `learnable_parameters` |
@@ -265,6 +266,9 @@ when subdivision is enabled), plus one at dataset open.
 | | nested `validate_subdivision_reaches_the_builder`: `enabled: true` + explicit `conus_adjacency`/`gages_adjacency` pointing at a non-subdivided store | `"params.subdivision"` + `"conflicts with the explicit"` |
 | `validate_geodataset` | `geodataset:` contradicting the adjacency source (`ddm30` with `geospatial_fabric`, `merit` with `gridded_network`) | `"geodataset"` + the source key. Absent ⇒ inferred; explicit adjacency paths ⇒ any label allowed |
 | `validate_leakance` | `use_leakance` + `use_cuda_graphs` | both key names |
+| `validate_leakance_gate` | `leakance_gate` block without `use_leakance: true` | `"leakance_gate"` + `"use_leakance"` |
+| | empty `leakance_gate.temperature` | `"temperature"` |
+| | a temperature that is not positive and finite | `"epoch N"` + `"positive"` |
 | `validate_ddr_match` | `use_cuda_graphs: true` without the deprecated `ddr_match: true` | `"use_cuda_graphs: true` requires the DEPRECATED `ddr_match: true"` |
 | `validate_enforce_positivity` | `enforce_positivity: true` + `ddr_match: true` | `"requires \`ddr_match: false\`"` |
 | `validate_disagg_pretrained` | `freeze: true` without `pretrained_checkpoint` | `"freeze: true requires pretrained_checkpoint"` |
@@ -316,3 +320,58 @@ target/release/eval --config config/experiments/leakance_hourly_on.yaml \
 ```
 
 Takes about 10 minutes, no retrain.
+
+## Leakance gate (`params.leakance_gate`, added 2026-09-16)
+
+`zeta = leakance_factor · area_z · K_D · (depth − d_gw)` multiplies
+`leakance_factor` and `K_D`, so only their product reaches the physics: an
+exact scaling degeneracy (measured `rho(leakance_factor, K_D) = +0.9986` over
+346,321 CONUS reaches on a trained arm, with 97 % of the four outputs' variance
+in one direction). Physically a reach either sits above a losing aquifer or it
+does not, so `leakance_factor` is meant to be a SELECTOR of where leakance
+acts, not a continuous multiplier. The gate makes it one without a
+straight-through estimator: the head's normalized output `u ∈ (0, 1)` becomes
+
+```text
+g = sigmoid( logit(clamp(u, 1e-6, 1 − 1e-6)) / tau ),   logit(x) = ln(x / (1 − x))
+```
+
+`tau = 1` is the identity (`g = u`; the code returns `u` untouched, so it is
+BIT-exact, not round-tripped through `ln`/`exp`). Smaller `tau` sharpens toward
+a step at `u = 0.5`; `0.5` is a fixed point at every temperature. The clamp
+margin `1e-6` sits ~17 f32 ulps below 1.0 so `1 − u` keeps ~3 % relative
+accuracy at the edge, and it only alters values the head has already saturated
+past `|pre-activation| > 13.8`. Autograd through `clamp → log → div → sigmoid`
+gives the exact gradient; `TimestepLeakanceOp` and `src/routing/` are untouched
+because the gate is applied to the head OUTPUT before `setup_inputs`
+denormalizes it (`src/training/gate.rs::leakance_gate`).
+
+```yaml
+params:
+  use_leakance: true
+  leakance_gate:
+    temperature: {1: 1.0, 6: 0.5, 11: 0.2, 16: 0.05}   # keyed by 1-indexed epoch
+```
+
+| Key | Default | Notes |
+|---|---|---|
+| `temperature` | required | Epoch-keyed schedule resolved like `experiment.learning_rate` (`LeakanceGate::resolve`: largest key `<= epoch`, first value before it). Every value must be positive and finite; `1.0` = identity, `> 1` softens |
+
+Which temperature each reader applies:
+
+| Reader | Temperature | Log line |
+|---|---|---|
+| `src/training/forward.rs::forward` (training) | `resolve(epoch)`, threaded by the driver as `gate_tau` | `epoch E lr=… leakance_gate_tau=…` |
+| `src/training/forward.rs::forward_eval_core` (eval) | `final_temperature()` | `leakance gate: final temperature tau=…` |
+| `src/training/probe.rs::probe_forward` | `final_temperature()` (after lifting, so leaf grads carry the gate's Jacobian) | same, from `probe_zeta_gradient` |
+| `src/dump_parameters.rs::dump` / `dump_init` | `final_temperature()` (the dumped factor is the gated value) | same |
+
+`forward` asserts `gate_tau.is_some() == params.leakance_gate.is_some()`, so a
+caller cannot train an ungated model against a gated config or vice versa. In
+`forward_eval_core` the gate is applied BEFORE `LeakanceOverride`, so an
+override still replaces the value that reaches denormalization. Tests:
+`tests/leakance_gate.rs` (tau = 1 bit-exact identity through all three readers,
+central-difference gradcheck at five temperatures, monotonicity, limits,
+saturation without NaN/Inf) and the gated rows of `tests/gamma_eval_parity.rs`.
+This does not re-open the leakance verdict in `research-status.md`; it removes
+one degeneracy so a future arm can be judged on the gate, not on the product.

.claude/skills/ddrs-dev/references/testing.md

@@ -74,7 +74,10 @@ cargo test --lib && cargo test && \
   mkdir -p output && cargo run --release --example compare_ddr_sandbox
 ```
 If you touched `src/training/forward.rs` (disagg / leakance threading), also run
-`cargo test --test leakance_off_parity`.
+`cargo test --test leakance_off_parity`. If you touched any of the three head
+readers (`forward`, `forward_eval_core`, `probe_forward`) or
+`src/training/gate.rs`, also run
+`cargo test --features fixtures --test gamma_eval_parity --test leakance_gate`.
 
 If you touched `src/routing/leakance.rs`, `src/geometry.rs`'s depth inversion, or
 the S6/S25 leakance call site in `src/routing/mmc_op.rs::forward_chain_inner`, run
@@ -148,7 +151,7 @@ line citation outright rather than trusting it to stay pinned.
 | Gridded (DDM30) ingestion | `gridded_bundle` (sub-reach store → subdivided layout, gauge at the cell's last piece, cache hit, dataset opens), `cargo test --lib adjacency::gridded` (synthetic-store validation), `hourly_streamflow::time_major_store_reads_identically_to_divide_major`, `cargo test --lib zarr::tests::upstream_comids_names_each_subdivided_parent_once` |
 | Sparse / autograd | `sparse_gradcheck`, `sp8_gradcheck` |
 | KAN head | the 4 `kan_head_*` fixture tests (need `--features fixtures`) |
-| Leakance | `leakance_reference_match` (cross-implementation, vs DDR `_compute_zeta` @ `c2bd0f9`), `leakance_gradcheck`, `leakance_off_parity`, `zeta_accum` (incl. multi-timestep volume accounting) |
+| Leakance | `leakance_reference_match` (cross-implementation, vs DDR `_compute_zeta` @ `c2bd0f9`), `leakance_gradcheck`, `leakance_off_parity`, `zeta_accum` (incl. multi-timestep volume accounting), `leakance_gate` (tau = 1 bit-exact identity through all three readers, gate gradcheck, saturation) |
 | Subdivision | `subdivide`, `subdivision_integration`, `gauge_mass_conservation`; `compare_ddr_sandbox` must still report ABSOLUTE MATCH |
 | Adjacency | `adjacency_parity` (managed builder byte-identical to the petgraph engine on `order`/`indices_0`/`indices_1`), `adjacency_build`, `data_zarr_store::conus_adjacency_loads_real_merit_zarr` (invariant 3 on real CONUS data) |
 | CLI / data | `data_dataset`, `data_static`, `cli_manifest`, `cli_lockfile`, `cli_json_contract` |

experiments/leakance/bias_decomposition.py

@@ -0,0 +1,60 @@
+"""Decompose where each learned term acts. KGE splits skill into three independent
+error modes: r (timing and shape), alpha = sd_sim/sd_obs (amplitude), beta =
+mean_sim/mean_obs (volume). Roughness conserves mass, so it can move r and alpha
+but not beta. Leakance removes mass, so beta is the only mode it can reach that
+roughness cannot. Comparing the same three numbers for the unrouted inflow, the
+leakance-off arm and the leakance-on arm says which bias each term actually absorbed."""
+import numpy as np, zarr, json
+B="/home/tbindas/projects/ddrs/.ddrs/runs"
+ON=f"{B}/2026-09-16T04-44-15Z-train-and-test"; OFF=f"{B}/2026-09-16T06-11-35Z-train-and-test"
+
+def comps(sim, obs):
+    out=[]
+    for i in range(sim.shape[0]):
+        s,o = sim[i], obs[i]
+        m = np.isfinite(s)&np.isfinite(o)&(o>=0)
+        if m.sum() < 365: out.append((np.nan,)*4); continue
+        s,o = s[m], o[m]
+        so, oo = s.std(), o.std()
+        if oo <= 0 or so <= 0: out.append((np.nan,)*4); continue
+        r = np.corrcoef(s,o)[0,1]; a = so/oo; b = s.mean()/o.mean()
+        kge = 1 - np.sqrt((r-1)**2 + (a-1)**2 + (b-1)**2)
+        out.append((r,a,b,kge))
+    return np.array(out)
+
+def load_eval(run):
+    g = zarr.open_group(f"{run}/eval/predictions.zarr", mode="r")
+    return np.asarray(g["predictions"][:],float), np.asarray(g["observations"][:],float), np.asarray(g["gage_ids"][:])
+
+res={}
+for lab, run in [("leakance ON", ON), ("leakance OFF", OFF)]:
+    sim, obs, ids = load_eval(run)
+    if sim.shape[0] < sim.shape[1] and sim.shape[0] < 5000: pass
+    else: sim, obs = sim.T, obs.T
+    res[lab] = (comps(sim, obs), ids)
+    print(f"{lab}: {sim.shape} gauges x days")
+
+# unrouted summed inflow, from the shared baseline
+bm = json.load(open(f"{ON}/baseline/manifest.json"))
+ng = len(bm["gage_ids"]); nd = len(bm["time_range_daily"])
+bp = np.fromfile(f"{ON}/baseline/predictions.f32", dtype=np.float32).reshape(ng, nd).astype(float)
+bo = np.fromfile(f"{ON}/baseline/observations.f32", dtype=np.float32).reshape(ng, nd).astype(float)
+res["summed inflow (no routing)"] = (comps(bp, bo), np.asarray(bm["gage_ids"]))
+print(f"summed inflow: {bp.shape} gauges x days")
+
+print(f"\n{'arm':28s} {'r':>8s} {'alpha':>8s} {'beta':>8s} {'KGE':>8s}   n")
+order=["summed inflow (no routing)","leakance OFF","leakance ON"]
+med={}
+for k in order:
+    c,_ = res[k]; ok=np.isfinite(c[:,3])
+    med[k]=np.nanmedian(c[ok],axis=0)
+    print(f"{k:28s} " + " ".join(f"{v:8.4f}" for v in med[k]) + f"   {ok.sum()}")
+
+print("\nwhat routing moved (leakance OFF minus unrouted):")
+d=med["leakance OFF"]-med["summed inflow (no routing)"]
+print("   " + " ".join(f"{n}={v:+.4f}" for n,v in zip(["r","alpha","beta","KGE"],d)))
+print("what leakance added on top (ON minus OFF):")
+d2=med["leakance ON"]-med["leakance OFF"]
+print("   " + " ".join(f"{n}={v:+.4f}" for n,v in zip(["r","alpha","beta","KGE"],d2)))
+print("\nvolume bias |beta-1| at the median gauge:")
+for k in order: print(f"   {k:28s} {abs(med[k][2]-1)*100:6.2f} %")

experiments/leakance/clamp_saturation.py

@@ -0,0 +1,66 @@
+#!/usr/bin/env python
+"""How often do the geometry clamps saturate? Each saturated clamp is a gradient
+sink: the parameter behind it receives zero gradient at that reach, so no amount
+of training can identify it there. The equations reference flags this as the
+measurement the repo has never taken.
+
+Reconstructs the geometry chain from a parameter dump plus the eval diagnostic's
+mean depth, and reports the saturated fraction for each clamp in
+params.attribute_minimums and for the side-slope bounds.
+
+Caveat stated up front: this uses each reach's MEAN eval depth, so it is a
+lower bound on saturation. A clamp that fires only at low flow is invisible here.
+
+usage: clamp_saturation.py <params.nc> <kan_parameters.nc> [--config C]
+"""
+import argparse, json
+import numpy as np, netCDF4 as nc, pandas as pd
+
+def main():
+    ap = argparse.ArgumentParser()
+    ap.add_argument("params"); ap.add_argument("zeta")
+    ap.add_argument("--depth-lb", type=float, default=0.01)
+    ap.add_argument("--velocity-lb", type=float, default=0.01)
+    ap.add_argument("--bottom-width-lb", type=float, default=0.01)
+    ap.add_argument("--side-slope-bounds", type=float, nargs=2, default=(0.5, 50.0))
+    a = ap.parse_args()
+
+    P = nc.Dataset(a.params); Z = nc.Dataset(a.zeta)
+    pc = np.asarray(P["COMID"][:], dtype=np.int64)
+    zc = np.asarray(Z["COMID_eval"][:], dtype=np.int64)
+    pos = pd.Index(pc).get_indexer(zc); ok = pos >= 0
+    g = lambda d, v: np.asarray(d[v][:], dtype=float)
+    depth = g(Z, "depth_mean")[ok]
+    q = g(P, "q_spatial")[pos[ok]] if "q_spatial" in P.variables else np.full(ok.sum(), 0.65)
+    p = g(P, "p_spatial")[pos[ok]] if "p_spatial" in P.variables else np.full(ok.sum(), 21.0)
+    n = g(P, "n")[pos[ok]]
+
+    # geometry chain, mirroring src/geometry.rs
+    top_width   = p * depth**q
+    side_slope  = top_width * q / (2.0 * depth)
+    side_clamped = np.clip(side_slope, *a.side_slope_bounds)
+    bottom_width = top_width - 2.0 * side_clamped * depth
+
+    rep = {
+        "n_eval_reaches": int(ok.sum()),
+        "caveat": "computed from each reach's MEAN eval depth, so every figure is a LOWER bound on saturation",
+        "depth at its floor": float((depth <= a.depth_lb * 1.001).mean()),
+        "side slope at its lower bound 0.5": float((side_slope <= a.side_slope_bounds[0]).mean()),
+        "side slope at its upper bound 50": float((side_slope >= a.side_slope_bounds[1]).mean()),
+        "side slope clamped either way": float(((side_slope <= a.side_slope_bounds[0]) | (side_slope >= a.side_slope_bounds[1])).mean()),
+        "bottom width at or below its floor": float((bottom_width <= a.bottom_width_lb).mean()),
+        "roughness at its range floor 0.015": float((n <= 0.015 * 1.001).mean()),
+        "roughness at its range ceiling 0.25": float((n >= 0.25 * 0.999).mean()),
+    }
+    # the side-slope clamp is the interesting one: it is a pure function of q
+    rep["note_side_slope"] = ("side_slope = p*d^q * q / (2d) = (q/2)*p*d^(q-1); with q fixed at 0.65 "
+                              "it saturates the 0.5 lower bound whenever d exceeds (q*p/(2*0.5))^(1/(1-q))")
+    if np.isfinite(q).all() and np.allclose(q, q[0]):
+        qq, pp = q[0], p[0]
+        d_crit = (qq * pp / (2.0 * a.side_slope_bounds[0]))**(1.0 / (1.0 - qq))
+        rep["depth above which the side slope pins at 0.5 (m)"] = float(d_crit)
+        rep["fraction of reaches above that depth"] = float((depth > d_crit).mean())
+    print(json.dumps(rep, indent=2))
+
+if __name__ == "__main__":
+    main()