From dd37883d3ce6235a20eca1bee8f741e59478f81a Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 30 Oct 2026 12:07:45 +0000 Subject: [PATCH] =?UTF-8?q?research:=20D=C3=BCtsch=20analyzer=20=E2=80=94?= =?UTF-8?q?=20horizon-limited=20coverage/events,=20Parseval=20HF=20weights?= MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit 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 --- research/duetsch_re100_kc5/analyze.py | 19 ++++++++++-- tests/unit/test_duetsch_protocol.py | 42 +++++++++++++++++++++++++++ 2 files changed, 58 insertions(+), 3 deletions(-) diff --git a/research/duetsch_re100_kc5/analyze.py b/research/duetsch_re100_kc5/analyze.py index f089769..0f21cfd 100644 --- a/research/duetsch_re100_kc5/analyze.py +++ b/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) diff --git a/tests/unit/test_duetsch_protocol.py b/tests/unit/test_duetsch_protocol.py index 945633f..16c0259 100644 --- a/tests/unit/test_duetsch_protocol.py +++ b/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