diff --git a/experiments/levers/blend_mechanism.py b/experiments/levers/blend_mechanism.py new file mode 100644 index 00000000..7c6af2ba --- /dev/null +++ b/experiments/levers/blend_mechanism.py @@ -0,0 +1,96 @@ +"""Why the multi-product blend scores higher: exact per-gauge ambiguity decomposition, test vs training years. + +For weights w_i summing to 1 and blend pbar = sum w_i p_i, at every time step + sum_i w_i (p_i - o)^2 = (pbar - o)^2 + sum_i w_i (p_i - pbar)^2, +so, dividing the time sums by sum (o - obar)^2, + NSE_blend = sum_i w_i NSE_i + SPREAD, SPREAD = sum_t sum_i w_i (p_i - pbar)^2 / sum_t (o - obar)^2 >= 0. +SPREAD is how much the members disagree, in NSE units; it is the whole gain over the weighted-mean member. +This is the ambiguity decomposition of Krogh & Vedelsby (1995, NIPS 7, "Neural network ensembles, cross validation, +and active learning") divided by the observed variance: a standard multi-model ensemble gain, not a routing result. +SPREAD is split by time scale on each member's deviation d_i = p_i - pbar: water-year mean (year-to-year volume), +30-day centred moving average of the rest (seasonal / monthly), and the remainder (event scale, < ~30 d). +Writes blend_mechanism.txt. +""" +import numpy as np +import pandas as pd + +import common as C + +TEST = {"UH": C.S42, "DIST": "2026-09-13T17-22-30Z-train-and-test", "LSTM": "2026-09-13T13-55-03Z-train-and-test"} +TRAIN = {"UH": "2026-09-30T13-13-59Z-train-and-test", "DIST": "2026-09-30T14-56-59Z-train-and-test", + "LSTM": "2026-09-30T14-30-08Z-train-and-test"} +WEIGHTS = {"equal": {"UH": 1 / 3, "DIST": 1 / 3, "LSTM": 1 / 3}, "training-year": {"UH": 0.5, "DIST": 0.25, "LSTM": 0.25}} + + +def load(runs, t0, t1, ref_ids=None): + """Members aligned to ref_ids on the days every run has inside [t0, t1]; observations from the UH run.""" + raw = {k: C.load_preds(r) for k, r in runs.items()} + ref = ref_ids or raw["UH"][0] + days = pd.DatetimeIndex(sorted(set.intersection(*[set(v[1]) for v in raw.values()]))) + days = days[(days >= t0) & (days <= t1)] + P = {} + for k, (i, t, Pr, Or) in raw.items(): + pos = t.get_indexer(days) + P[k] = C.align(i, Pr, ref)[:, pos] + if k == "UH": + O = C.align(i, Or, ref)[:, pos] + return ref, days, P, O + + +def bands(d, tt): + """Split d[g, t] into (water-year mean, 30-day MA of the rest, remainder).""" + wy = C.water_year(tt) + low = np.empty_like(d) + for y in np.unique(wy): + s = wy == y + low[:, s] = d[:, s].mean(1, keepdims=True) + rest = d - low + mid = pd.DataFrame(rest.T).rolling(30, center=True, min_periods=1).mean().to_numpy().T + return low, mid, rest - mid + + +def decompose(P, O, tt, w): + m = np.isfinite(O) & np.all([np.isfinite(P[k]) for k in P], axis=0) + nv = m.sum(1) + om = np.where(m, O, 0.0) + mo = om.sum(1) / np.maximum(nv, 1) + sst = (np.where(m, O - mo[:, None], 0.0) ** 2).sum(1) + pbar = sum(w[k] * P[k] for k in P) + nse = {k: 1 - (np.where(m, P[k] - O, 0.0) ** 2).sum(1) / sst for k in P} + nse_blend = 1 - (np.where(m, pbar - O, 0.0) ** 2).sum(1) / sst + spread = sum(w[k] * (np.where(m, P[k] - pbar, 0.0) ** 2).sum(1) for k in P) / sst + parts = np.zeros((3, len(sst))) + for k in P: + for j, b in enumerate(bands(P[k] - pbar, tt)): + parts[j] += w[k] * (np.where(m, b, 0.0) ** 2).sum(1) / sst + wmean = sum(w[k] * nse[k] for k in P) + ok = (nv >= 365) & np.isfinite(nse_blend) + df = pd.DataFrame(dict(blend=nse_blend, wmean=wmean, spread=spread, low=parts[0], mid=parts[1], high=parts[2], + **{f"nse_{k}": nse[k] for k in P}))[ok] + assert np.allclose(df.blend, df.wmean + df.spread, atol=1e-8), "identity failed" + return df + + +out = [] +ids, tt_te, P_te, O_te = load(TEST, pd.Timestamp("1995-10-01"), pd.Timestamp("2010-09-30")) +_, tt_tr, P_tr, O_tr = load(TRAIN, pd.Timestamp("1985-11-01"), pd.Timestamp("1995-09-30"), ids) +for period, (P, O, tt) in {"test WY1996-2010": (P_te, O_te, tt_te), "training WY1986-1995": (P_tr, O_tr, tt_tr)}.items(): + for wname, w in WEIGHTS.items(): + D = decompose(P, O, tt, w) + med = D.median() + tot = D[["low", "mid", "high"]].sum(1) + share = (D[["low", "mid", "high"]].div(tot, axis=0)).median() + out.append(f"{period}, {wname} weights, {len(D)} gauges (medians over gauges)") + out.append(f" member NSE: UH {med.nse_UH:.4f} DIST {med.nse_DIST:.4f} LSTM {med.nse_LSTM:.4f}") + out.append(f" weighted mean of member NSEs {med.wmean:.4f} + spread {med.spread:.4f} = blend; blend median " + f"{med.blend:.4f}; blend - UH: median {D.blend.median() - D.nse_UH.median():+.4f}, " + f"paired {np.median(D.blend - D.nse_UH):+.4f}") + out.append(f" paired per gauge: blend - UH = (wmean - UH) + spread: median {np.median(D.wmean - D.nse_UH):+.4f} " + f"+ {med.spread:+.4f}") + out.append(f" spread by time scale, median share: year-to-year volume {share.low:.2f}, seasonal/monthly " + f"{share.mid:.2f}, event (< 30 d) {share.high:.2f}; band sum / spread median " + f"{np.median(tot / D.spread):.3f}") + out.append("") +txt = "\n".join(out) +(C.HERE / "blend_mechanism.txt").write_text(txt + "\n") +print(txt) diff --git a/experiments/levers/blend_mechanism.txt b/experiments/levers/blend_mechanism.txt new file mode 100644 index 00000000..29994d38 --- /dev/null +++ b/experiments/levers/blend_mechanism.txt @@ -0,0 +1,24 @@ +test WY1996-2010, equal weights, 2365 gauges (medians over gauges) + member NSE: UH 0.7391 DIST 0.7300 LSTM 0.6359 + weighted mean of member NSEs 0.6957 + spread 0.0590 = blend; blend median 0.7637; blend - UH: median +0.0246, paired +0.0182 + paired per gauge: blend - UH = (wmean - UH) + spread: median -0.0352 + +0.0590 + spread by time scale, median share: year-to-year volume 0.06, seasonal/monthly 0.16, event (< 30 d) 0.76; band sum / spread median 0.975 + +test WY1996-2010, training-year weights, 2365 gauges (medians over gauges) + member NSE: UH 0.7391 DIST 0.7300 LSTM 0.6359 + weighted mean of member NSEs 0.7064 + spread 0.0511 = blend; blend median 0.7650; blend - UH: median +0.0259, paired +0.0194 + paired per gauge: blend - UH = (wmean - UH) + spread: median -0.0264 + +0.0511 + spread by time scale, median share: year-to-year volume 0.06, seasonal/monthly 0.16, event (< 30 d) 0.76; band sum / spread median 0.975 + +training WY1986-1995, equal weights, 2365 gauges (medians over gauges) + member NSE: UH 0.7861 DIST 0.7801 LSTM 0.7079 + weighted mean of member NSEs 0.7528 + spread 0.0439 = blend; blend median 0.8038; blend - UH: median +0.0177, paired +0.0148 + paired per gauge: blend - UH = (wmean - UH) + spread: median -0.0241 + +0.0439 + spread by time scale, median share: year-to-year volume 0.06, seasonal/monthly 0.18, event (< 30 d) 0.74; band sum / spread median 0.972 + +training WY1986-1995, training-year weights, 2365 gauges (medians over gauges) + member NSE: UH 0.7861 DIST 0.7801 LSTM 0.7079 + weighted mean of member NSEs 0.7605 + spread 0.0381 = blend; blend median 0.8043; blend - UH: median +0.0182, paired +0.0153 + paired per gauge: blend - UH = (wmean - UH) + spread: median -0.0180 + +0.0381 + spread by time scale, median share: year-to-year volume 0.06, seasonal/monthly 0.18, event (< 30 d) 0.74; band sum / spread median 0.972 + diff --git a/experiments/levers/page/runoff_blend_explainer.html b/experiments/levers/page/runoff_blend_explainer.html new file mode 100644 index 00000000..e6a0342a --- /dev/null +++ b/experiments/levers/page/runoff_blend_explainer.html @@ -0,0 +1,226 @@ +Multi-Product Runoff Forcing + + + + + +
+
+
ddrs · 30 September 2026 · what is running and why
+

Multi-Product Runoff Forcing

+

Two experiments are running. One builds a new way to feed runoff into the routing: a weighted blend of three runoff products instead of one, the only lever we found that can lift the population median NSE by about +0.03. The other measures how well the gauges can pin down each dam parameter. This page explains both.

+
+ +
+
Running now
+

Two experiments

+
+
+ building, then full-population runs +

Multi-product runoff forcing, with γ fixed at 0

+
+ QuestionDoes routing a blend of three runoff products, with the roughness exponent γ fixed at 0, lift the median NSE of all 2,365 evaluation gauges by about +0.03? + StepsTeach the engine to read several runoff stores and average them per catchment; pass the regression gates; then train on the full population: blend with γ = 0 at seeds 42 and 43, and the blend alone at seed 42. + Compared toThe existing no-dam runs with the same training setup (2026-09-27T07-29-47Z seed 42, 10-31-30Z seed 43), paired gauge by gauge. + ExpectAbout +0.026 at the population median from the blend with γ = 0 (offline estimate); KGE about flat. + WhenGates in about 1 to 2 hours, then about 2.5 hours of training and testing. +
+
+
+ offline landscapes, then engine checks +

How learnable is each dam parameter?

+
+ QuestionFor each dam parameter (residence time, seasonal rule curve, flood pool size, release target, capture share), can the downstream gauge pin it down, and why could training not grow the flood pools? + StepsMap the training loss along each parameter at every on-reach dam gauge; measure its curvature at the optimum, how wide the good region is, and the slope at the value training starts from; confirm the flood-pool plateau with a few engine runs. + OutputA verdict per parameter: stiff (learnable), sloppy (flat), or learnable only from a good start, on the same scale as the routing parameters in the paper. + WhenAbout 3 hours. +
+
+
+

Also open: pull request #52, the records from after #51 merged (no code changes). Its checks are running; it merges when you say so.

+
+ +
+
Background
+

What runoff forcing is

+

ddrs does not make runoff. For each of the 346,321 MERIT unit catchments it takes the runoff that drains into that reach, called Q′ (lateral inflow, m³/s), from a hydrologic model run beforehand, and routes it down the river network with Muskingum-Cunge. Routing can delay and smooth that water, and the dam law can store and release it, but whatever is wrong in Q′ reaches the gauge. Today every run uses one Q′ product: dHBV2 with a unit hydrograph (UH).

+ +
+ +
+
The idea
+

Three products, averaged

+
+ + + + + + + +
ProductWhat it isWeightError correlation with UHTest-period volume, predicted / observed
UHdHBV2 with a unit hydrograph; today's forcing0.501.001.045
DISTdHBV2, distributed version, AORC forcing0.250.931.065
LSTMA daily NeuralHydrology CudaLSTM, trained on October 1980 to September 19950.250.620.965
Equal-weight blendFor reference1.024
+
Q′blend(catchment, day) = Σi wi · Q′i(catchment, day) / Σi covers the catchment wi
+

Why it helps. This is the standard multi-model ensemble gain (the ambiguity decomposition of Krogh and Vedelsby, 1995), not a routing result. At every gauge, exactly: blend NSE = weighted mean of the members' NSEs + how much the members disagree (their spread around the blend, divided by the variance of observed flow). The weaker members cost 0.026 against UH; their disagreement removes 0.051 of squared error; net +0.019 per gauge, +0.026 at the population median. 76 % of that disagreement is at the event scale (under 30 days), 16 % seasonal, and only 6 % year-to-year volume. So the products mostly miss different storms, and averaging cancels those misses. The LSTM helps because its errors correlate only 0.62 with UH's, not because it is good (alone it scores 0.636). Weighted, 75 % of gauges improve.

+

How the members and weights were chosen, without looking at the test years. Starting from UH, products were added one at a time while they improved the training-year (1986 to 1995) summed-runoff score: DIST, then LSTM, then no further product helped (dHBV2 lumped and HydroDL were rejected). The weights 0.50 / 0.25 / 0.25 were fitted on routed training-year runs (1986 to 1995). The test years played no part in either choice.

+

What the engine build changes. The offline estimate averaged three separately routed runs. The engine version averages the runoff first and trains one routing head on the blend, so there is still one set of channel parameters. Routing is mildly nonlinear, so these two can differ; the new runs measure that.

+
+ +
+
Where it stands against the target
+

Every lever against +0.03

+
+
Hypothesis

If the population median is limited by runoff errors that no routing or dam law can fix, only a lever that changes the runoff itself can move the median by about +0.03, while the dam and routing levers stay in the thousandths.

+
Plain English

Each row is one change and the median NSE change it brings over all 2,365 gauges. Rows marked measured come from full training runs; the others are offline estimates, fitted on training years and scored on test years.

+
Graph

The dot is the change in the population median, the whisker its 95 % interval, or for the dam recipe the two seeds. The dashed line is the +0.03 target. Orange rows are measured, grey rows are offline estimates.

+
+
+
Measured, full training runsOffline estimateTarget +0.03
+
+
+

Result. The dam recipe moves the population median by about +0.004 (two seeds, +0.009 and −0.002). A runoff multiplier fitted honestly gives +0.0005, and fixing γ at 0 about +0.002. The blend gives +0.026 on its own. Inside the blend, γ = 0 adds about +0.001 with p and q fixed (as in the running experiment), and about +0.005 if p and q are learned too. Conclusion. The blend is the only change that gets near +0.03. With p and q fixed it needs the dam recipe's few thousandths on top to get there. The running experiment will show whether the engine version, one head trained on the blended runoff, matches these offline estimates.

+
+ +
+
Before relying on it
+

What to watch for

+
+

KGE barely moves. The blend smooths the hydrographs (variability ratio 0.914 to 0.869), which NSE rewards and KGE penalises. The weighted blend lifts the median KGE by +0.004 (interval includes 0) and leaves it unchanged gauge by gauge; the equal-weight blend loses 0.005 gauge by gauge.

+

A little of the gain may be period-specific. The blend gains +0.018 on the training years and +0.026 on the test years. In 1996 to 2010, DIST runs high and the LSTM low on volume, but year-to-year volume differences make up only 6 % of the members' disagreement, so most of the gain is not drift cancelling.

+

All three stores must exist for any new domain or period, which rules out the global runs today.

+

The model changes meaning. One routing head trained on a blend of three products is not the same as a head trained on one product. For the paper's inflow-product comparisons, the single-product runs remain the reference.

+

One seed each so far for the offline estimates; the two no-dam seeds differ by 0.005 at the median, which the new two-seed runs will account for.

+
+

Sources: experiments/levers/report.txt, combine_ens.txt, blend_mechanism.txt, store_selection.txt; the full-population dam runs in experiments/reservoir/full_run/paired_full_run_full_s7_s{42,43}.json; branch multi-product-forcing.

+
+
+ +