research: Dütsch analysis — Morison fit per force candidate, explicit --reference

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

research/duetsch_re100_kc5/analyze.py

@@ -62,11 +62,18 @@ def common_points(dicts):
     return sorted(keys)
 
 
-def cycle_morison(case, ks):
+def candidates(case):
+    """Per-step Fx of each reported force candidate (none is selected as an RHS)."""
+    d = case["diag"]
+    return {"cme": d[:, 0], "cme_wen": d[:, 0] + d[:, 2], "cme_node": d[:, 0] + d[:, 6] - d[:, 9]}
+
+
+def cycle_morison(case, ks, fx=None):
     T = case["T"]
-    idx = np.concatenate([np.arange(len(case["force"]))[cycle_slices(case)[k]] for k in ks])
+    fx = case["force"][:, 0] if fx is None else fx
+    idx = np.concatenate([np.arange(len(fx))[cycle_slices(case)[k]] for k in ks])
     b = np.array([p.interval_bases(int(m), T) for m in idx])
-    return p.morison_fit(case["force"][idx, 0] / case["F_ref"], b[:, 0], b[:, 1])
+    return p.morison_fit(fx[idx] / case["F_ref"], b[:, 0], b[:, 1])
 
 
 def pair_metrics(case, k):
@@ -91,7 +98,14 @@ def ledger(case):
     return pairs, p.settled_start({k: v["pass"] for k, v in pairs.items()})
 
 
-def main(paths):
+def main(argv):
+    """argv: CASE_DIR ... [--reference CASE_DIR]; the reference defaults to the
+    largest D (resolution ladder); domain runs name the 60 D case explicitly."""
+    ref_path = None
+    if "--reference" in argv:
+        i = argv.index("--reference")
+        ref_path, argv = argv[i + 1], argv[:i] + argv[i + 2:]
+    paths = argv if ref_path is None or ref_path in argv else argv + [ref_path]
     cases = [load(path) for path in paths]
     ledgers = [ledger(c) for c in cases]
     window = p.common_window([k for _, k in ledgers],
@@ -104,15 +118,19 @@ def main(paths):
                      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()))
         if window:
-            entry["morison"] = cycle_morison(case, window)
+            entry["morison_obs"] = cycle_morison(case, window)
+            entry["morison"] = {name: cycle_morison(case, window, fx)
+                                for name, fx in candidates(case).items()}
             entry["fy_mean"] = float(np.mean(np.concatenate(
                 [case["force"][cycle_slices(case)[k], 1] for k in window])) / case["F_ref"])
         report["cases"].append(entry)
     if window and len(cases) > 1:  # cross-resolution waveform against the finest case
-        fine = max(cases, key=lambda c: c["D"])
+        fine = (next(c for c in cases if c["path"] == ref_path) if ref_path
+                else max(cases, key=lambda c: c["D"]))
+        report["reference"] = fine["path"]
         ref = np.concatenate([coarse_wave(fine, k) for k in window])
         for case, entry in zip(cases, report["cases"]):
-            r = fine["T"] // case["T"]
+            r = fine["T"] // case["T"]  # 1 for equal-D domain runs
             ref_c = ref.reshape(-1, r).mean(axis=1)
             own = np.concatenate([coarse_wave(case, k) for k in window])
             entry["wave_rms_vs_finest"] = float(np.sqrt(np.mean((own - ref_c) ** 2)))