research: Dütsch analyzer — horizon-limited coverage/events, Parseval HF weights

Coverage and no_donor/eq_only totals now stop at the horizon; HF energy
uses one-sided Parseval weights (interior bins x2, DC/Nyquist x1). Focused
regressions for both. run.py, src and thresholds unchanged.

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>

research/duetsch_re100_kc5/analyze.py

@@ -107,13 +107,23 @@ def hf_component(case, ks):
     spec = np.fft.rfft(resid)
     high = spec.copy()
     high[:21] = 0.0
-    energy = np.abs(spec) ** 2
+    energy = one_sided_energy(spec, T)
     return dict(label="unqualified phase-locked HF force component, not checked",
                 cycles=list(ks), hf_rms_over_F_ref=float(np.sqrt(np.mean(np.fft.irfft(high, T) ** 2))),
                 hf_energy_fraction=float(energy[21:].sum() / max(energy.sum(), 1e-300)),
                 top_harmonics=[int(j) for j in np.argsort(energy[21:])[::-1][:5] + 21])
 
 
+def one_sided_energy(spec, n):
+    """Parseval-consistent energy per rfft bin: interior bins count twice, DC and
+    (for even n) Nyquist once."""
+    w = np.full(len(spec), 2.0)
+    w[0] = 1.0
+    if n % 2 == 0:
+        w[-1] = 1.0
+    return w * np.abs(spec) ** 2
+
+
 def ledger(case, horizon=p.LAST_SCORED):
     ks = [k for k in sorted(cycle_slices(case)) if k <= horizon]
     pairs = {k: pair_metrics(case, k) for k in ks if k - 1 > p.DISCARD}
@@ -144,11 +154,14 @@ def main(argv):
                   horizon_mode="default-cap" if horizon == p.LAST_SCORED else
                   "extended duration diagnostic (declared amendment, D40/60D only)")
     for case, (pairs, k_case) in zip(cases, ledgers):
+        limit = round((0.25 + horizon) * case["T"])  # events inside the horizon only
         entry = {key: case[key] for key in ("path", "D", "T", "box", "n", "u_max", "nu", "tau", "A",
                                             "steps", "cycles", "F_ref")}
         entry.update(duration_T=case["steps"] / case["T"], K_case=k_case, pairs=pairs,
-                     coverage={k: coverage(case, k) for k in sorted(cycle_slices(case)) if k > p.DISCARD},
-                     no_donor=float(case["diag"][:, 13].sum()), eq_only=float(case["diag"][:, 12].sum()),
+                     coverage={k: coverage(case, k) for k in sorted(cycle_slices(case))
+                               if p.DISCARD < k <= horizon},
+                     no_donor=float(case["diag"][:limit, 13].sum()),
+                     eq_only=float(case["diag"][:limit, 12].sum()),
                      hf_component=hf_component(case, window or tuple(
                          [k for k in sorted(cycle_slices(case)) if p.DISCARD < k <= horizon][-3:])))
         if window:  # residual drift inside the scored window (pairs K..K+2)

tests/unit/test_duetsch_protocol.py

@@ -108,3 +108,45 @@ def test_generated_config_passes_the_real_loader(tmp_path):
     assert cfg.collision == "MRT" and cfg.streaming == "double_buffer" and not cfg.ddf_shifting
     assert cfg.velocity == 0.0 and cfg.viscosity == pytest.approx(u_max * 20 / 100)
     assert cfg.y_wall_bc == "bounce_back"
+
+
+def _analyze():
+    import sys
+    sys.path.insert(0, str(Path(__file__).resolve().parents[2] / "research/duetsch_re100_kc5"))
+    try:
+        import analyze
+    finally:
+        sys.path.pop(0)
+    return analyze
+
+
+def test_hf_energy_is_parseval_weighted():
+    a = _analyze()
+    n = 360
+    t = np.arange(n)
+    x = 1.0 + np.cos(2 * np.pi * 5 * t / n) + 2.0 * np.cos(2 * np.pi * 30 * t / n) + 0.5 * np.cos(np.pi * t)
+    e = a.one_sided_energy(np.fft.rfft(x), n) / n
+    assert e.sum() == pytest.approx(np.sum(x ** 2))               # Parseval
+    expected = (2.0 ** 2 / 2 + 0.5 ** 2) / np.mean(x ** 2)       # bins >= 21 (30 and Nyquist)
+    assert e[21:].sum() / e.sum() == pytest.approx(expected)
+
+
+def test_horizon_excludes_future_coverage_and_events(monkeypatch):
+    a = _analyze()
+    T, horizon = 360, 4
+    steps = round((0.25 + 6) * T)                                 # cycles 1..6 recorded
+    diag = np.zeros((steps, 15))
+    diag[round((0.25 + horizon) * T):, 13] = 1.0                  # synthetic future no_donor
+    rows = [(s, p.cycle(s, T), 180.0, 0.0, 0.0, 0.1, 0.0, 1.0) for s in range(steps)
+            if p.bin_phase(s, T) is not None and p.cycle(s, T) > p.DISCARD]
+    case = dict(path="c", D=20, T=T, box=10, n=200, u_max=.04, nu=.008, tau=.524, A=1.0,
+                steps=steps, cycles=6, F_ref=1.0, force=np.zeros((steps, 2)), diag=diag,
+                samples=np.asarray(rows, float))
+    monkeypatch.setattr(a, "load", lambda path: case)
+    out = []
+    monkeypatch.setattr("builtins.print", lambda text: out.append(text))
+    a.main(["c", "--horizon", str(horizon)])
+    import json
+    entry = json.loads(out[0])["cases"][0]
+    assert sorted(int(k) for k in entry["coverage"]) == [3, 4]   # nothing from cycles 5, 6
+    assert entry["no_donor"] == 0.0