feat(theory): establish steady-cloak analytical baseline

Freeze the source-derived action contract and direct q-in objective so the new MFS strip model and serial CFD runner provide a compact, auditable basis for the steady-cloak study.

Co-authored-by: Cursor <cursoragent@cursor.com>

src/SR_analysis/data/v5/ill_075L/config.json

@@ -1,5 +1,5 @@
 {
-  "case_name": "ill_075L_sc",
+  "case_name": "ill_075L",
   "scene_type": "illusion",
   "seed": 43,
   "SI": 400,

src/SR_analysis/data/v5/ill_15L/config.json

@@ -1,5 +1,5 @@
 {
-  "case_name": "ill_15L_sc",
+  "case_name": "ill_15L",
   "scene_type": "illusion",
   "seed": 43,
   "SI": 800,

src/SR_analysis/data/v5/ill_1L/config.json

@@ -1,5 +1,5 @@
 {
-  "case_name": "ill_1L_sc",
+  "case_name": "ill_1L",
   "scene_type": "illusion",
   "seed": 43,
   "SI": 600,

src/SR_analysis/data/v5/ill_2L/config.json

@@ -1,5 +1,5 @@
 {
-  "case_name": "ill_2L_sc",
+  "case_name": "ill_2L",
   "scene_type": "illusion",
   "seed": 43,
   "SI": 800,

src/SR_analysis/data/v5/kar_d075/config.json

@@ -1,5 +1,5 @@
 {
-  "case_name": "kar_d075_sc",
+  "case_name": "kar_d075",
   "scene_type": "karman",
   "seed": 44,
   "SI": 800,

src/SR_analysis/data/v5/kar_d15/config.json

@@ -1,5 +1,5 @@
 {
-  "case_name": "kar_d15_sc",
+  "case_name": "kar_d15",
   "scene_type": "karman",
   "seed": 45,
   "SI": 800,

src/SR_analysis/data/v5/kar_d2/config.json

@@ -1,5 +1,5 @@
 {
-  "case_name": "kar_d2_sc",
+  "case_name": "kar_d2",
   "scene_type": "karman",
   "seed": 45,
   "SI": 800,

src/SR_analysis/data/v5/kar_re100/config.json

@@ -1,5 +1,5 @@
 {
-  "case_name": "kar_re100_sc",
+  "case_name": "kar_re100",
   "scene_type": "karman",
   "seed": 45,
   "SI": 800,

src/SR_analysis/experiments/v5/bridge_v5.py

@@ -9,19 +9,19 @@ For each effective V5 case:
   5. Creates SR_analysis/data/v5/{case}/ with controlled.npz, target.npz, etc.
 
 Effective cases (with best seed and DTW):
-  - kar_re100_sc (seed45, dtw=0.918) — baseline
-  - kar_d075_sc (seed44, dtw=0.911) — vardist
-  - kar_d15_sc (seed45, dtw=0.892) — vardist
-  - kar_d2_sc (seed45, dtw=0.816) — vardist, borderline
-  - ill_075L_sc (seed43, dtw=0.807)
-  - ill_1L_sc (seed43, dtw=0.900)
-  - ill_15L_sc (seed43, dtw=0.906)
-  - ill_2L_sc (seed43, dtw=0.800)
+  - kar_re100 (seed45, dtw=0.918) — baseline
+  - kar_d075 (seed44, dtw=0.911) — vardist
+  - kar_d15 (seed45, dtw=0.892) — vardist
+  - kar_d2 (seed45, dtw=0.816) — vardist, borderline
+  - ill_075L (seed43, dtw=0.807)
+  - ill_1L (seed43, dtw=0.900)
+  - ill_15L (seed43, dtw=0.906)
+  - ill_2L (seed43, dtw=0.800)
 
 Usage:
     PYTHONPATH=src conda run -n pycuda_3_10 python -m SR_analysis.experiments.v5.bridge_v5
     PYTHONPATH=src conda run -n pycuda_3_10 python -m SR_analysis.experiments.v5.bridge_v5 --include-failed
-    PYTHONPATH=src conda run -n pycuda_3_10 python -m SR_analysis.experiments.v5.bridge_v5 --cases kar_re100_sc,ill_1L_sc
+    PYTHONPATH=src conda run -n pycuda_3_10 python -m SR_analysis.experiments.v5.bridge_v5 --cases kar_re100,ill_1L
 """
 from __future__ import annotations
 
@@ -32,6 +32,8 @@ from pathlib import Path
 
 import numpy as np
 
+from drl_pinball.case_registry import CASE_REGISTRY
+
 _SR_ROOT = Path(__file__).resolve().parents[2]
 _REPO = _SR_ROOT.parents[1]
 _EVAL_OUT = _REPO / "src" / "drl_pinball" / "eval" / "output" / "train"
@@ -42,71 +44,44 @@ _SR_DATA = _SR_ROOT / "data" / "v5"
 L0 = 20.0
 U0 = 0.01
 
-# Best seeds from eval metrics.json
-EFFECTIVE_CASES = {
-    "kar_re100_sc": {"seed": 45, "dtw": 0.918, "scene_type": "karman",
-                     "cal_name": "kar_re100",
-                     "si": 800, "action_scale": 12.0, "action_bias": (0.0, 0.0, 0.0),
-                     "s_dim": 12, "conv_len": 30,
-                     "obs_slice": (2, 14), "n_objects": 7},
-    "kar_d075_sc": {"seed": 44, "dtw": 0.911, "scene_type": "karman",
-                    "cal_name": "kar_d075_sc",
-                    "si": 800, "action_scale": 12.0, "action_bias": (0.0, 0.0, 0.0),
-                    "s_dim": 12, "conv_len": 30,
-                    "obs_slice": (2, 14), "n_objects": 7},
-    "kar_d15_sc": {"seed": 45, "dtw": 0.892, "scene_type": "karman",
-                   "cal_name": "kar_d15_sc",
-                   "si": 800, "action_scale": 12.0, "action_bias": (0.0, 0.0, 0.0),
-                   "s_dim": 12, "conv_len": 30,
-                   "obs_slice": (2, 14), "n_objects": 7},
-    "kar_d2_sc": {"seed": 45, "dtw": 0.816, "scene_type": "karman",
-                  "cal_name": "kar_d2_sc",
-                  "si": 800, "action_scale": 12.0, "action_bias": (0.0, 0.0, 0.0),
-                  "s_dim": 12, "conv_len": 30,
-                  "obs_slice": (2, 14), "n_objects": 7},
-    "ill_075L_sc": {"seed": 43, "dtw": 0.807, "scene_type": "illusion",
-                    "cal_name": "ill_075L",
-                    "si": 400, "action_scale": 12.0, "action_bias": (0.0, 0.0, 0.0),
-                    "s_dim": 14, "conv_len": 36,
-                    "target_diam": 0.75,
-                    "obs_slice": (0, 12), "n_objects": 6},
-    "ill_1L_sc": {"seed": 43, "dtw": 0.900, "scene_type": "illusion",
-                  "cal_name": "ill_1L",
-                  "si": 600, "action_scale": 12.0, "action_bias": (0.0, 0.0, 0.0),
-                  "s_dim": 14, "conv_len": 36,
-                  "target_diam": 1.0,
-                  "obs_slice": (0, 12), "n_objects": 6},
-    "ill_15L_sc": {"seed": 43, "dtw": 0.906, "scene_type": "illusion",
-                   "cal_name": "ill_15L",
-                   "si": 800, "action_scale": 12.0, "action_bias": (0.0, 0.0, 0.0),
-                   "s_dim": 14, "conv_len": 36,
-                   "target_diam": 1.5,
-                   "obs_slice": (0, 12), "n_objects": 6},
-    "ill_2L_sc": {"seed": 43, "dtw": 0.800, "scene_type": "illusion",
-                  "cal_name": "ill_2L",
-                  "si": 800, "action_scale": 12.0, "action_bias": (0.0, 0.0, 0.0),
-                  "s_dim": 14, "conv_len": 36,
-                  "target_diam": 2.0,
-                  "obs_slice": (0, 12), "n_objects": 6},
+# Evaluation attributes not supplied by the canonical registry.
+_CASE_METRICS = {
+    "kar_re100": (45, 0.918, 30, (2, 14), 7),
+    "kar_d075": (44, 0.911, 30, (2, 14), 7),
+    "kar_d15": (45, 0.892, 30, (2, 14), 7),
+    "kar_d2": (45, 0.816, 30, (2, 14), 7),
+    "ill_075L": (43, 0.807, 36, (0, 12), 6),
+    "ill_1L": (43, 0.900, 36, (0, 12), 6),
+    "ill_15L": (43, 0.906, 36, (0, 12), 6),
+    "ill_2L": (43, 0.800, 36, (0, 12), 6),
+    "kar_re60": (43, 0.364, 30, (2, 14), 7),
+    "kar_re200": (43, 0.712, 30, (2, 14), 7),
+    "kar_re400": (43, 0.565, 30, (2, 14), 7),
 }
 
-FAILED_CASES = {
-    "kar_re60_sc": {"seed": 43, "dtw": 0.364, "scene_type": "karman",
-                    "cal_name": "kar_re60_sc",
-                    "si": 800, "action_scale": 12.0, "action_bias": (0.0, 0.0, 0.0),
-                    "s_dim": 12, "conv_len": 30,
-                    "obs_slice": (2, 14), "n_objects": 7},
-    "kar_re200_sc": {"seed": 43, "dtw": 0.712, "scene_type": "karman",
-                     "cal_name": "kar_re200_sc",
-                     "si": 500, "action_scale": 12.0, "action_bias": (0.0, 0.0, 0.0),
-                     "s_dim": 12, "conv_len": 30,
-                     "obs_slice": (2, 14), "n_objects": 7},
-    "kar_re400_sc": {"seed": 43, "dtw": 0.565, "scene_type": "karman",
-                     "cal_name": "kar_re400_sc",
-                     "si": 400, "action_scale": 12.0, "action_bias": (0.0, 0.0, 0.0),
-                     "s_dim": 12, "conv_len": 30,
-                     "obs_slice": (2, 14), "n_objects": 7},
-}
+
+def _bridge_config(case_id: str) -> dict:
+    case = CASE_REGISTRY[case_id]
+    seed, dtw, conv_len, obs_slice, n_objects = _CASE_METRICS[case_id]
+    cfg = {
+        "seed": seed, "dtw": dtw, "scene_type": case.scene_type,
+        "cal_name": case.calibration, "si": case.si,
+        "action_scale": 12.0, "action_bias": (0.0, 0.0, 0.0),
+        "s_dim": 14 if case.scene_type == "illusion" else 12,
+        "conv_len": conv_len, "obs_slice": obs_slice, "n_objects": n_objects,
+    }
+    if case.target_diam is not None:
+        cfg["target_diam"] = case.target_diam
+    return cfg
+
+
+EFFECTIVE_CASES = {case_id: _bridge_config(case_id) for case_id in (
+    "kar_re100", "kar_d075", "kar_d15", "kar_d2",
+    "ill_075L", "ill_1L", "ill_15L", "ill_2L",
+)}
+FAILED_CASES = {case_id: _bridge_config(case_id) for case_id in (
+    "kar_re60", "kar_re200", "kar_re400",
+)}
 
 
 def bridge_case(case_name: str, cfg: dict, out_dir: Path) -> bool:

src/SR_analysis/experiments/v5/stage_2_fit_v5.py

@@ -9,10 +9,10 @@ Supports per-scene individual and cross-scene joint fitting.
 
 Usage:
     # Per-scene (fast, niter=40)
-    PYTHONPATH=src conda run -n sr_env python -m SR_analysis.experiments.v5.stage_2_fit_v5 --scene kar_re100_sc
+    PYTHONPATH=src conda run -n sr_env python -m SR_analysis.experiments.v5.stage_2_fit_v5 --scene kar_re100
 
     # Per-scene deep search (niter=120)
-    PYTHONPATH=src conda run -n sr_env python -m SR_analysis.experiments.v5.stage_2_fit_v5 --scene kar_re100_sc --deep
+    PYTHONPATH=src conda run -n sr_env python -m SR_analysis.experiments.v5.stage_2_fit_v5 --scene kar_re100 --deep
 
     # Joint: all effective Karman cases
     PYTHONPATH=src conda run -n sr_env python -m SR_analysis.experiments.v5.stage_2_fit_v5 --group karman_v5 --mode joint
@@ -21,7 +21,7 @@ Usage:
     PYTHONPATH=src conda run -n sr_env python -m SR_analysis.experiments.v5.stage_2_fit_v5 --group illusion_v5 --mode joint
 
     # Joint: mixed V5 + legacy (cross-solver)
-    PYTHONPATH=src conda run -n sr_env python -m SR_analysis.experiments.v5.stage_2_fit_v5 --scenes kar_re100_sc,kar_d075_sc,kar_d15_sc,karman_re100,karman_re200 --mode joint --output-label mixed_karman
+    PYTHONPATH=src conda run -n sr_env python -m SR_analysis.experiments.v5.stage_2_fit_v5 --scenes kar_re100,kar_d075,kar_d15,karman_re100,karman_re200 --mode joint --output-label mixed_karman
 """
 from __future__ import annotations
 
@@ -46,6 +46,7 @@ from SR_analysis.utils.feature_builder import (
     PHASE_STATE_KEYS, ILLUSION_PHASE_KEYS, CORE_FEAT_KEYS_V2, MU_FEAT_KEYS,
 )
 from SR_analysis.configs import get_scene  # for legacy scenes
+from drl_pinball.case_registry import CASE_REGISTRY
 
 ALL_MU = list(MU_FEAT_KEYS) + ["mu_Cl_tot"]
 PHYS = [k for k in CORE_FEAT_KEYS_V2 if not k.startswith(("aF_","aB_","aT_","daF","daB","daT"))]
@@ -68,15 +69,10 @@ V5_ILLUSION_CFG = {
     "action_bias": V5_ACTION_BIAS,
 }
 
-# SI overrides per case
-V5_SI_OVERRIDE = {
-    "ill_075L_sc": 400, "ill_1L_sc": 600, "ill_15L_sc": 800, "ill_2L_sc": 800,
-    "kar_re60_sc": 800, "kar_re200_sc": 500, "kar_re400_sc": 400,
-}
-
+# V5 case attributes come from drl_pinball.case_registry.
 
 def is_v5_scene(name: str) -> bool:
-    return "_sc" in name or "_tr" in name
+    return name in CASE_REGISTRY
 
 
 def load_controlled_v5(case_name: str):
@@ -97,8 +93,8 @@ def load_controlled_v5(case_name: str):
         cfg["si_actual"] = v5_cfg.get("SI", 800)
         cfg["scene_type"] = v5_cfg.get("scene_type", "")
 
-    # Apply SI override
-    si = V5_SI_OVERRIDE.get(case_name, cfg.get("si_actual", 800))
+    # Apply canonical training/evaluation SI.
+    si = CASE_REGISTRY[case_name].si
     cfg["sample_interval"] = si
     cfg["si_actual"] = si  # for feature builder
 
@@ -148,8 +144,8 @@ FORMULA_DIR = _SR_ROOT / "results" / "formulas_v5"
 os.makedirs(FORMULA_DIR, exist_ok=True)
 
 FIT_GROUPS_V5 = {
-    "karman_v5": ["kar_re100_sc", "kar_d075_sc", "kar_d15_sc", "kar_d2_sc"],
-    "illusion_v5": ["ill_075L_sc", "ill_1L_sc", "ill_15L_sc", "ill_2L_sc"],
+    "karman_v5": ["kar_re100", "kar_d075", "kar_d15", "kar_d2"],
+    "illusion_v5": ["ill_075L", "ill_1L", "ill_15L", "ill_2L"],
 }
 
 

src/SR_analysis/results/README.md

@@ -59,8 +59,11 @@ It adds, without refitting formulas:
 - `term_contributions_sr_closed_loop_400.csv`: additive terms on 400-step SR trajectories, with exact next-action reconstruction;
 - `ablation_summary.csv` and `scaling_summary.csv`;
 - seven constant-rotation steady time-series exports;
-- publication-starting-point PNG/PDF diagnostics under `figures/`;
+- four publication PNG/PDF figures under `figures/` (performance/duration, pointwise generalization, steady calibration, and a phase-aligned Re100 example);
+- publication tables under `tables/` for offline next-action RMSE and 40-step term deletion, replacing the former figure versions;
+- package-level `README.md` and `phase_alignment.json` documenting scope and the exact example-window alignment;
 - a SHA-256 manifest for the complete derived package.
+- a GPU-generated `07_flow_field_comparison_karman_re100` exporter/output contract: target, PPO, and canonical SR vorticity are selected by one-run-per-condition exact same-sample center-sensor phase matching and share one physical crop/color scale. The retained field comparison is one matched snapshot per controller, not an ensemble or robustness result; no field artifact exists until the exporter is run.
 
 The 400-step PPO similarities are `0.954856`, `0.946973`, `0.900280`, `0.845378` for Kármán Re50/100/200/400 and `0.976328`, `0.975728`, `0.926701` for Illusion 0.75L/1L/1.5L. All completed with finite telemetry. PPO inference runs on CPU while physical GPU 2 remains dedicated to PyCUDA CFD, avoiding PyTorch/PyCUDA context conflicts.
 

src/SR_analysis/results/formulas_v5/ill_075L_ill_1L_ill_15L_ill_2L_front.json

@@ -1,5 +1,5 @@
 {
-  "scene": "ill_075L_sc_ill_1L_sc_ill_15L_sc_ill_2L_sc",
+  "scene": "ill_075L_ill_1L_ill_15L_ill_2L",
   "channel": "front",
   "output": "alpha",
   "feature_keys": [

src/SR_analysis/results/formulas_v5/ill_1L_front.json

@@ -1,5 +1,5 @@
 {
-  "scene": "ill_1L_sc",
+  "scene": "ill_1L",
   "channel": "front",
   "output": "alpha",
   "feature_keys": [

src/SR_analysis/results/formulas_v5/ill_1L_top.json

@@ -1,5 +1,5 @@
 {
-  "scene": "ill_1L_sc",
+  "scene": "ill_1L",
   "channel": "top",
   "output": "alpha",
   "feature_keys": [

src/SR_analysis/results/formulas_v5/kar_d075_front.json

@@ -1,5 +1,5 @@
 {
-  "scene": "kar_d075_sc",
+  "scene": "kar_d075",
   "channel": "front",
   "output": "alpha",
   "feature_keys": [

src/SR_analysis/results/formulas_v5/kar_d075_top.json

@@ -1,5 +1,5 @@
 {
-  "scene": "kar_d075_sc",
+  "scene": "kar_d075",
   "channel": "top",
   "output": "alpha",
   "feature_keys": [

src/SR_analysis/results/formulas_v5/kar_d15_front.json

@@ -1,5 +1,5 @@
 {
-  "scene": "kar_d15_sc",
+  "scene": "kar_d15",
   "channel": "front",
   "output": "alpha",
   "feature_keys": [

src/SR_analysis/results/formulas_v5/kar_d15_top.json

@@ -1,5 +1,5 @@
 {
-  "scene": "kar_d15_sc",
+  "scene": "kar_d15",
   "channel": "top",
   "output": "alpha",
   "feature_keys": [

src/SR_analysis/results/formulas_v5/kar_d2_front.json

@@ -1,5 +1,5 @@
 {
-  "scene": "kar_d2_sc",
+  "scene": "kar_d2",
   "channel": "front",
   "output": "alpha",
   "feature_keys": [

src/SR_analysis/results/formulas_v5/kar_d2_top.json

@@ -1,5 +1,5 @@
 {
-  "scene": "kar_d2_sc",
+  "scene": "kar_d2",
   "channel": "top",
   "output": "alpha",
   "feature_keys": [

src/SR_analysis/results/formulas_v5/kar_re100_front.json

@@ -1,5 +1,5 @@
 {
-  "scene": "kar_re100_sc",
+  "scene": "kar_re100",
   "channel": "front",
   "output": "alpha",
   "feature_keys": [

src/SR_analysis/results/formulas_v5/kar_re100_kar_d075_kar_d15_kar_d2_front.json

@@ -1,5 +1,5 @@
 {
-  "scene": "kar_re100_sc_kar_d075_sc_kar_d15_sc_kar_d2_sc",
+  "scene": "kar_re100_kar_d075_kar_d15_kar_d2",
   "channel": "front",
   "output": "alpha",
   "feature_keys": [

src/SR_analysis/results/runs/article2-plotting-package-20260721/README.md

@@ -0,0 +1,13 @@
+# SR publication plotting package
+
+This package contains five generated publication figures, two replacement tables, and two PPT-ready 16:9 summary pages. All values are derived from frozen canonical artifacts without formula refitting.
+
+- `01_training_case_performance`: panel (a) compares PPO and SR absolute closed-loop legacy DTW similarity at **400 control steps** for all seven training cases; panel (b) reports each controller's 400-minus-200 change. Here 200 and 400 are run durations in control steps, not Reynolds numbers.
+- `02_pointwise_generalization`: 200-step SR results at training conditions and sampled unseen interpolation/extrapolation conditions. Each unseen point is one realization. Categories are unconnected because the evidence does not establish a continuous parameter law or statistical robustness.
+- `03_steady_rotation_calibration`: the retained disturbance-free constant rear-rotation sweep diagnostic.
+- `06_example_timeseries_karman_re100`: late-window phase portraits and PPO/SR actions over approximately three target cycles. Target indices 96–145 (`t_D=38.4`–`58.0`) are used. PPO and SR segments at the same indices are circularly aligned to target `sensors_center_uy` by mean-centered normalized cross-correlation; deterministic roll lags are PPO=0 and SR=-4 samples. No target actions exist or are plotted.
+- `07_flow_field_comparison_karman_re100` (generated by `tools/export_flow_comparison.py`): target wake, PPO-controlled pinball, and canonical SR-controlled pinball vorticity at one deterministic stable-cycle phase. Over indices [96,146), each trajectory standardizes downstream center-sensor (ux,uy) separately and defines theta=atan2(v_z,u_z). An exhaustive joint search minimizes the two wrapped angular errors plus 0.001 rad/sample times total temporal separation, with a same-direction branch check. Each condition runs once through index 145 while all 50 stable-window fields are held as float32 host-memory candidates (375 MiB total for three 512x1280 buffers); the selected field and phase diagnostic are therefore the exact same CFD sample from the same run. This matches center-sensor limit-cycle phase; it does not assert six-sensor state equality or exact full-field identity. It stores one compressed NPZ, JSON provenance, and PNG/PDF. This is a single snapshot, not an ensemble or uncertainty estimate.
+- `tables/offline_action_rmse.*`: offline SR-vs-PPO next-action RMSE on causally aligned PPO-visited states. This is not closed-loop performance.
+- `tables/term_deletion.*`: actual 40-step closed-loop deletion results, canonical formula identities, same-window parent comparisons, and aggregate mean/min summaries.
+
+The `presentation/` directory contains two white-background 16:9 PNG/PDF pages and its source-artifact README. PNG and PDF files share each retained figure stem. `phase_alignment.json` records the exact Figure 06 selection/alignment contract. `manifest.json` hashes every package CSV, Markdown, JSON/NPZ metadata or field artifact, and publication figure.

src/SR_analysis/results/runs/article2-plotting-package-20260721/phase_alignment.json

@@ -0,0 +1,19 @@
+{
+  "scene": "karman_re100",
+  "source_package": "article2-timeseries-csv-20260720",
+  "target_index_start_inclusive": 96,
+  "target_index_stop_exclusive": 146,
+  "target_t_D_start": 38.4,
+  "target_t_D_end": 58.0,
+  "samples": 50,
+  "displayed_target_cycles": 3,
+  "alignment_signal": "sensors_center_uy",
+  "alignment_method": "mean-centered normalized circular cross-correlation; smallest maximizing roll selected",
+  "roll_lags_samples": {
+    "Target": 0,
+    "PPO": 0,
+    "SR": -4
+  },
+  "control_dt_D_over_U0": 0.4,
+  "target_actions_plotted": false
+}

src/SR_analysis/results/runs/article2-plotting-package-20260721/presentation/README.md

@@ -0,0 +1,10 @@
+# SR presentation graphics
+
+- `01_sr_method_and_formulas`: a 16:9 methodology-to-formula summary. Formula coefficients and expressions are read dynamically from the canonical Kármán and Illusion JSON artifacts. The page states the front odd projection, lower-cylinder symmetry map, and evidence-bounded physical interpretation.
+- `02_sr_quantitative_evidence`: a 16:9 summary of frozen closed-loop evidence: 400-control-step SR/PPO training-case ranges and means, 200-step pointwise generalization ranges, and 40-step parent-relative term-deletion effects. It deliberately excludes offline imitation RMSE from article-performance claims.
+
+Reproduce both PNG and PDF pages with:
+
+`conda run -n base python src/SR_analysis/tools/plot_sr_presentation.py`
+
+Sources: canonical formula JSONs in `article-refit-{karman,illusion}-topology-a-20260718`; 400-step SR and PPO DTW convergence CSVs; `article2-generalization-summary-20260720/generalization.csv`; and `tables/term_deletion.csv`. All similarities use `legacy_dtw_v1_abs_n_unclipped`. Windows and limitations are printed on the pages.

src/SR_analysis/results/runs/article2-plotting-package-20260721/tables/offline_action_rmse.csv

@@ -0,0 +1,8 @@
+objective,training_case,n_causally_aligned_states,rmse_alpha_front,rmse_alpha_upper,rmse_alpha_lower,aggregate_rmse_alpha
+karman,karman_re50,197,1.73860463,1.39458231,1.13278257,1.44346830
+karman,karman_re100,197,0.51356311,1.24024221,1.32638969,1.08953467
+karman,karman_re200,197,1.47993917,2.29449643,2.86751651,2.28601434
+karman,karman_re400,197,0.40632824,2.55717318,2.72429918,2.16994956
+illusion,illusion_0.75L,197,0.82228198,1.77041216,1.43962577,1.40036058
+illusion,illusion_1L,197,0.85035894,1.20967232,0.58682215,0.91847295
+illusion,illusion_1.5L,197,3.20534050,2.22584755,1.74754340,2.46863748

src/SR_analysis/results/runs/article2-plotting-package-20260721/tables/offline_action_rmse.md

@@ -0,0 +1,15 @@
+# Offline SR-vs-PPO next-action RMSE
+
+RMSE is evaluated on causally aligned PPO-visited states: each post-action state predicts the next PPO action. This is an offline imitation diagnostic, not closed-loop performance.
+
+| Objective | Training case | States | Front RMSE | Upper RMSE | Lower RMSE | Aggregate RMSE |
+|---|---|---:|---:|---:|---:|---:|
+| karman | karman_re50 | 197 | 1.73860463 | 1.39458231 | 1.13278257 | 1.44346830 |
+| karman | karman_re100 | 197 | 0.51356311 | 1.24024221 | 1.32638969 | 1.08953467 |
+| karman | karman_re200 | 197 | 1.47993917 | 2.29449643 | 2.86751651 | 2.28601434 |
+| karman | karman_re400 | 197 | 0.40632824 | 2.55717318 | 2.72429918 | 2.16994956 |
+| illusion | illusion_0.75L | 197 | 0.82228198 | 1.77041216 | 1.43962577 | 1.40036058 |
+| illusion | illusion_1L | 197 | 0.85035894 | 1.20967232 | 0.58682215 | 0.91847295 |
+| illusion | illusion_1.5L | 197 | 3.20534050 | 2.22584755 | 1.74754340 | 2.46863748 |
+
+All action errors are in dimensionless surface-speed α. Aggregate RMSE pools the three action channels and all aligned states within a case.

src/SR_analysis/results/runs/article2-plotting-package-20260721/tables/term_deletion.csv

@@ -0,0 +1,25 @@
+objective,head,variant,parent_formula,deleted_term,retained_formula,case,steps,legacy_dtw,parent_legacy_dtw_40_step,delta_vs_parent_40_step,variant_mean,variant_min
+karman,front odd projection,k_front0,-0.3813914338074549 * Cd_rear_a,-0.3813914338074549*Cd_rear_a,0,karman_re50,40,0.93936105,0.94588134,-0.00652029,0.88656666,0.81328922
+karman,front odd projection,k_front0,-0.3813914338074549 * Cd_rear_a,-0.3813914338074549*Cd_rear_a,0,karman_re100,40,0.91218229,0.91160360,+0.00057869,0.88656666,0.81328922
+karman,front odd projection,k_front0,-0.3813914338074549 * Cd_rear_a,-0.3813914338074549*Cd_rear_a,0,karman_re200,40,0.88143408,0.88661335,-0.00517927,0.88656666,0.81328922
+karman,front odd projection,k_front0,-0.3813914338074549 * Cd_rear_a,-0.3813914338074549*Cd_rear_a,0,karman_re400,40,0.81328922,0.83848301,-0.02519379,0.88656666,0.81328922
+karman,upper/lower shared-symmetry rear,k_rear0,1.3077817865589976 * Cl_rear_s - 3.431208680510616,1.3077817865589976*Cl_rear_s,-3.431208680510616,karman_re50,40,0.91712633,0.94588134,-0.02875501,0.85046363,0.78953465
+karman,upper/lower shared-symmetry rear,k_rear0,1.3077817865589976 * Cl_rear_s - 3.431208680510616,1.3077817865589976*Cl_rear_s,-3.431208680510616,karman_re100,40,0.86994335,0.91160360,-0.04166025,0.85046363,0.78953465
+karman,upper/lower shared-symmetry rear,k_rear0,1.3077817865589976 * Cl_rear_s - 3.431208680510616,1.3077817865589976*Cl_rear_s,-3.431208680510616,karman_re200,40,0.82525019,0.88661335,-0.06136315,0.85046363,0.78953465
+karman,upper/lower shared-symmetry rear,k_rear0,1.3077817865589976 * Cl_rear_s - 3.431208680510616,1.3077817865589976*Cl_rear_s,-3.431208680510616,karman_re400,40,0.78953465,0.83848301,-0.04894836,0.85046363,0.78953465
+karman,upper/lower shared-symmetry rear,k_rear1,1.3077817865589976 * Cl_rear_s - 3.431208680510616,-3.431208680510616,1.3077817865589976*Cl_rear_s,karman_re50,40,0.83934413,0.94588134,-0.10653720,0.75314220,0.64740091
+karman,upper/lower shared-symmetry rear,k_rear1,1.3077817865589976 * Cl_rear_s - 3.431208680510616,-3.431208680510616,1.3077817865589976*Cl_rear_s,karman_re100,40,0.81773500,0.91160360,-0.09386860,0.75314220,0.64740091
+karman,upper/lower shared-symmetry rear,k_rear1,1.3077817865589976 * Cl_rear_s - 3.431208680510616,-3.431208680510616,1.3077817865589976*Cl_rear_s,karman_re200,40,0.70808876,0.88661335,-0.17852458,0.75314220,0.64740091
+karman,upper/lower shared-symmetry rear,k_rear1,1.3077817865589976 * Cl_rear_s - 3.431208680510616,-3.431208680510616,1.3077817865589976*Cl_rear_s,karman_re400,40,0.64740091,0.83848301,-0.19108210,0.75314220,0.64740091
+illusion,front odd projection,i_front0,-1.8266041890847688 * Cd_rear_a + 2.064492712013321 * Cl_F,-1.8266041890847688*Cd_rear_a,2.064492712013321*Cl_F,illusion_0.75L,40,0.93937493,0.95545614,-0.01608120,0.91931368,0.88508469
+illusion,front odd projection,i_front0,-1.8266041890847688 * Cd_rear_a + 2.064492712013321 * Cl_F,-1.8266041890847688*Cd_rear_a,2.064492712013321*Cl_F,illusion_1L,40,0.93348141,0.94330800,-0.00982659,0.91931368,0.88508469
+illusion,front odd projection,i_front0,-1.8266041890847688 * Cd_rear_a + 2.064492712013321 * Cl_F,-1.8266041890847688*Cd_rear_a,2.064492712013321*Cl_F,illusion_1.5L,40,0.88508469,0.87939766,+0.00568703,0.91931368,0.88508469
+illusion,front odd projection,i_front1,-1.8266041890847688 * Cd_rear_a + 2.064492712013321 * Cl_F,2.064492712013321*Cl_F,-1.8266041890847688*Cd_rear_a,illusion_0.75L,40,0.96006748,0.95545614,+0.00461134,0.90788838,0.83656314
+illusion,front odd projection,i_front1,-1.8266041890847688 * Cd_rear_a + 2.064492712013321 * Cl_F,2.064492712013321*Cl_F,-1.8266041890847688*Cd_rear_a,illusion_1L,40,0.92703454,0.94330800,-0.01627347,0.90788838,0.83656314
+illusion,front odd projection,i_front1,-1.8266041890847688 * Cd_rear_a + 2.064492712013321 * Cl_F,2.064492712013321*Cl_F,-1.8266041890847688*Cd_rear_a,illusion_1.5L,40,0.83656314,0.87939766,-0.04283452,0.90788838,0.83656314
+illusion,upper/lower shared-symmetry rear,i_rear0,1.2544396437730243 * Cd_rear_a - 1.5280742226055937 * Cl_F,1.2544396437730243*Cd_rear_a,-1.5280742226055937*Cl_F,illusion_0.75L,40,0.95448757,0.95545614,-0.00096857,0.92631380,0.88760304
+illusion,upper/lower shared-symmetry rear,i_rear0,1.2544396437730243 * Cd_rear_a - 1.5280742226055937 * Cl_F,1.2544396437730243*Cd_rear_a,-1.5280742226055937*Cl_F,illusion_1L,40,0.93685081,0.94330800,-0.00645719,0.92631380,0.88760304
+illusion,upper/lower shared-symmetry rear,i_rear0,1.2544396437730243 * Cd_rear_a - 1.5280742226055937 * Cl_F,1.2544396437730243*Cd_rear_a,-1.5280742226055937*Cl_F,illusion_1.5L,40,0.88760304,0.87939766,+0.00820537,0.92631380,0.88760304
+illusion,upper/lower shared-symmetry rear,i_rear1,1.2544396437730243 * Cd_rear_a - 1.5280742226055937 * Cl_F,-1.5280742226055937*Cl_F,1.2544396437730243*Cd_rear_a,illusion_0.75L,40,0.95464065,0.95545614,-0.00081549,0.91197030,0.85081581
+illusion,upper/lower shared-symmetry rear,i_rear1,1.2544396437730243 * Cd_rear_a - 1.5280742226055937 * Cl_F,-1.5280742226055937*Cl_F,1.2544396437730243*Cd_rear_a,illusion_1L,40,0.93045444,0.94330800,-0.01285357,0.91197030,0.85081581
+illusion,upper/lower shared-symmetry rear,i_rear1,1.2544396437730243 * Cd_rear_a - 1.5280742226055937 * Cl_F,-1.5280742226055937*Cl_F,1.2544396437730243*Cd_rear_a,illusion_1.5L,40,0.85081581,0.87939766,-0.02858185,0.91197030,0.85081581

src/SR_analysis/results/runs/article2-plotting-package-20260721/tables/term_deletion.md

@@ -0,0 +1,34 @@
+# Forty-step closed-loop term deletion
+
+Every comparison uses the same 40-control-step window and the legacy DTW metric. Parent values come from the scientifically comparable parent L2 runs; deltas are deletion minus parent. Formula identities are read from the canonical parent and deletion JSON artifacts.
+
+Front formulas are scalar generators deployed through the odd projection `α_F(x) = ½[f(x) − f(Gx)]`. Rear formulas define the upper action; the lower action is mapped by `α_L(x) = −α_U(Gx)`.
+
+| Objective | Head | Variant | Parent formula | Deleted term | Retained formula | Case | DTW | Parent DTW | Δ | Mean | Min |
+|---|---|---|---|---|---|---|---:|---:|---:|---:|---:|
+| karman | front odd projection | k_front0 | `-0.3813914338074549 · Cd_rear_a` | `-0.3813914338074549·Cd_rear_a` | `0` | karman_re50 | 0.93936105 | 0.94588134 | -0.00652029 | 0.88656666 | 0.81328922 |
+| karman | front odd projection | k_front0 | `-0.3813914338074549 · Cd_rear_a` | `-0.3813914338074549·Cd_rear_a` | `0` | karman_re100 | 0.91218229 | 0.91160360 | +0.00057869 | 0.88656666 | 0.81328922 |
+| karman | front odd projection | k_front0 | `-0.3813914338074549 · Cd_rear_a` | `-0.3813914338074549·Cd_rear_a` | `0` | karman_re200 | 0.88143408 | 0.88661335 | -0.00517927 | 0.88656666 | 0.81328922 |
+| karman | front odd projection | k_front0 | `-0.3813914338074549 · Cd_rear_a` | `-0.3813914338074549·Cd_rear_a` | `0` | karman_re400 | 0.81328922 | 0.83848301 | -0.02519379 | 0.88656666 | 0.81328922 |
+| karman | upper/lower shared-symmetry rear | k_rear0 | `1.3077817865589976 · Cl_rear_s - 3.431208680510616` | `1.3077817865589976·Cl_rear_s` | `-3.431208680510616` | karman_re50 | 0.91712633 | 0.94588134 | -0.02875501 | 0.85046363 | 0.78953465 |
+| karman | upper/lower shared-symmetry rear | k_rear0 | `1.3077817865589976 · Cl_rear_s - 3.431208680510616` | `1.3077817865589976·Cl_rear_s` | `-3.431208680510616` | karman_re100 | 0.86994335 | 0.91160360 | -0.04166025 | 0.85046363 | 0.78953465 |
+| karman | upper/lower shared-symmetry rear | k_rear0 | `1.3077817865589976 · Cl_rear_s - 3.431208680510616` | `1.3077817865589976·Cl_rear_s` | `-3.431208680510616` | karman_re200 | 0.82525019 | 0.88661335 | -0.06136315 | 0.85046363 | 0.78953465 |
+| karman | upper/lower shared-symmetry rear | k_rear0 | `1.3077817865589976 · Cl_rear_s - 3.431208680510616` | `1.3077817865589976·Cl_rear_s` | `-3.431208680510616` | karman_re400 | 0.78953465 | 0.83848301 | -0.04894836 | 0.85046363 | 0.78953465 |
+| karman | upper/lower shared-symmetry rear | k_rear1 | `1.3077817865589976 · Cl_rear_s - 3.431208680510616` | `-3.431208680510616` | `1.3077817865589976·Cl_rear_s` | karman_re50 | 0.83934413 | 0.94588134 | -0.10653720 | 0.75314220 | 0.64740091 |
+| karman | upper/lower shared-symmetry rear | k_rear1 | `1.3077817865589976 · Cl_rear_s - 3.431208680510616` | `-3.431208680510616` | `1.3077817865589976·Cl_rear_s` | karman_re100 | 0.81773500 | 0.91160360 | -0.09386860 | 0.75314220 | 0.64740091 |
+| karman | upper/lower shared-symmetry rear | k_rear1 | `1.3077817865589976 · Cl_rear_s - 3.431208680510616` | `-3.431208680510616` | `1.3077817865589976·Cl_rear_s` | karman_re200 | 0.70808876 | 0.88661335 | -0.17852458 | 0.75314220 | 0.64740091 |
+| karman | upper/lower shared-symmetry rear | k_rear1 | `1.3077817865589976 · Cl_rear_s - 3.431208680510616` | `-3.431208680510616` | `1.3077817865589976·Cl_rear_s` | karman_re400 | 0.64740091 | 0.83848301 | -0.19108210 | 0.75314220 | 0.64740091 |
+| illusion | front odd projection | i_front0 | `-1.8266041890847688 · Cd_rear_a + 2.064492712013321 · Cl_F` | `-1.8266041890847688·Cd_rear_a` | `2.064492712013321·Cl_F` | illusion_0.75L | 0.93937493 | 0.95545614 | -0.01608120 | 0.91931368 | 0.88508469 |
+| illusion | front odd projection | i_front0 | `-1.8266041890847688 · Cd_rear_a + 2.064492712013321 · Cl_F` | `-1.8266041890847688·Cd_rear_a` | `2.064492712013321·Cl_F` | illusion_1L | 0.93348141 | 0.94330800 | -0.00982659 | 0.91931368 | 0.88508469 |
+| illusion | front odd projection | i_front0 | `-1.8266041890847688 · Cd_rear_a + 2.064492712013321 · Cl_F` | `-1.8266041890847688·Cd_rear_a` | `2.064492712013321·Cl_F` | illusion_1.5L | 0.88508469 | 0.87939766 | +0.00568703 | 0.91931368 | 0.88508469 |
+| illusion | front odd projection | i_front1 | `-1.8266041890847688 · Cd_rear_a + 2.064492712013321 · Cl_F` | `2.064492712013321·Cl_F` | `-1.8266041890847688·Cd_rear_a` | illusion_0.75L | 0.96006748 | 0.95545614 | +0.00461134 | 0.90788838 | 0.83656314 |
+| illusion | front odd projection | i_front1 | `-1.8266041890847688 · Cd_rear_a + 2.064492712013321 · Cl_F` | `2.064492712013321·Cl_F` | `-1.8266041890847688·Cd_rear_a` | illusion_1L | 0.92703454 | 0.94330800 | -0.01627347 | 0.90788838 | 0.83656314 |
+| illusion | front odd projection | i_front1 | `-1.8266041890847688 · Cd_rear_a + 2.064492712013321 · Cl_F` | `2.064492712013321·Cl_F` | `-1.8266041890847688·Cd_rear_a` | illusion_1.5L | 0.83656314 | 0.87939766 | -0.04283452 | 0.90788838 | 0.83656314 |
+| illusion | upper/lower shared-symmetry rear | i_rear0 | `1.2544396437730243 · Cd_rear_a - 1.5280742226055937 · Cl_F` | `1.2544396437730243·Cd_rear_a` | `-1.5280742226055937·Cl_F` | illusion_0.75L | 0.95448757 | 0.95545614 | -0.00096857 | 0.92631380 | 0.88760304 |
+| illusion | upper/lower shared-symmetry rear | i_rear0 | `1.2544396437730243 · Cd_rear_a - 1.5280742226055937 · Cl_F` | `1.2544396437730243·Cd_rear_a` | `-1.5280742226055937·Cl_F` | illusion_1L | 0.93685081 | 0.94330800 | -0.00645719 | 0.92631380 | 0.88760304 |
+| illusion | upper/lower shared-symmetry rear | i_rear0 | `1.2544396437730243 · Cd_rear_a - 1.5280742226055937 · Cl_F` | `1.2544396437730243·Cd_rear_a` | `-1.5280742226055937·Cl_F` | illusion_1.5L | 0.88760304 | 0.87939766 | +0.00820537 | 0.92631380 | 0.88760304 |
+| illusion | upper/lower shared-symmetry rear | i_rear1 | `1.2544396437730243 · Cd_rear_a - 1.5280742226055937 · Cl_F` | `-1.5280742226055937·Cl_F` | `1.2544396437730243·Cd_rear_a` | illusion_0.75L | 0.95464065 | 0.95545614 | -0.00081549 | 0.91197030 | 0.85081581 |
+| illusion | upper/lower shared-symmetry rear | i_rear1 | `1.2544396437730243 · Cd_rear_a - 1.5280742226055937 · Cl_F` | `-1.5280742226055937·Cl_F` | `1.2544396437730243·Cd_rear_a` | illusion_1L | 0.93045444 | 0.94330800 | -0.01285357 | 0.91197030 | 0.85081581 |
+| illusion | upper/lower shared-symmetry rear | i_rear1 | `1.2544396437730243 · Cd_rear_a - 1.5280742226055937 · Cl_F` | `-1.5280742226055937·Cl_F` | `1.2544396437730243·Cd_rear_a` | illusion_1.5L | 0.85081581 | 0.87939766 | -0.02858185 | 0.91197030 | 0.85081581 |
+
+Interpretation: Kármán deletion effects are term-dependent, with deleting the rear constant producing the largest degradation; the tested front term is weak over this short window. Illusion deletions remain stable and comparatively close to their parent, so these runs do not establish uniqueness of every term. Means and minima summarize the listed training cases only.

src/SR_analysis/tests/test_export_flow_comparison.py

@@ -0,0 +1,96 @@
+from __future__ import annotations
+
+import importlib.util
+import json
+from pathlib import Path
+
+import numpy as np
+
+MODULE_PATH = Path(__file__).resolve().parents[1] / "tools/export_flow_comparison.py"
+spec = importlib.util.spec_from_file_location("export_flow_comparison", MODULE_PATH)
+module = importlib.util.module_from_spec(spec)
+assert spec.loader
+spec.loader.exec_module(module)
+
+
+def test_wrapped_phase_difference_crosses_branch_cut():
+    delta = module.wrapped_phase_difference(-np.pi + 0.02, np.pi - 0.03)
+    assert np.isclose(delta, 0.05)
+
+
+def test_joint_phase_selection_is_deterministic_and_nearby():
+    n = 20
+    phase = np.linspace(0, 4 * np.pi, n, endpoint=False)
+    def trace(offset: float) -> np.ndarray:
+        out = np.zeros((n, 6))
+        out[:, 2] = 2.0 + 3.0 * np.cos(phase + offset)
+        out[:, 3] = -1.0 + 0.5 * np.sin(phase + offset)
+        return out
+    first = module.select_joint_phase_match(trace(0), trace(0.05), trace(-0.04), 0, n, temporal_weight=0.001)
+    second = module.select_joint_phase_match(trace(0), trace(0.05), trace(-0.04), 0, n, temporal_weight=0.001)
+    assert first == second
+    assert len(set(first["indices"].values())) == 1
+    assert np.isclose(first["wrapped_angle_errors_rad"]["PPO"], 0.05)
+    assert np.isclose(first["wrapped_angle_errors_rad"]["SR"], 0.04)
+    assert first["objective_terms"]["temporal_distance_samples"] == 0
+
+
+def test_exact_candidate_selection_uses_trace_index_mapping():
+    candidates = np.empty((4, 2, 3), dtype=np.float32)
+    for candidate_index, trace_index in enumerate(range(96, 100)):
+        candidates[candidate_index].fill(trace_index)
+    field, candidate_index = module.select_exact_candidate_field(
+        candidates, np.arange(96, 100), 98, expected_start=96, expected_stop=100
+    )
+    assert candidate_index == 2
+    np.testing.assert_array_equal(field, np.full((2, 3), 98, dtype=np.float32))
+
+
+def test_exact_candidate_selection_rejects_mapping_mismatch():
+    candidates = np.zeros((4, 2, 3), dtype=np.float32)
+    mismatched = np.array([96, 97, 99, 100])
+    with np.testing.assert_raises_regex(ValueError, "does not exactly cover"):
+        module.select_exact_candidate_field(
+            candidates, mismatched, 99, expected_start=96, expected_stop=100
+        )
+
+
+def test_exact_candidate_selection_requires_float32():
+    candidates = np.zeros((4, 2, 3), dtype=np.float64)
+    with np.testing.assert_raises_regex(ValueError, "float32"):
+        module.select_exact_candidate_field(
+            candidates, np.arange(96, 100), 98, expected_start=96, expected_stop=100
+        )
+
+
+def test_orientation_transposes_xy_to_image_yx():
+    omega_xy = np.arange(12).reshape(4, 3)
+    oriented = module.orient_vorticity_xy_to_yx(omega_xy, (4, 3))
+    assert oriented.shape == (3, 4)
+    np.testing.assert_array_equal(oriented[:, 2], omega_xy[2, :])
+
+
+def test_crop_uses_physical_x_and_centered_y_coordinates():
+    field = np.arange(7 * 10).reshape(7, 10)
+    cropped, meta = module.crop_yx(field, (1.0, 3.0), (-1.0, 1.0), d_lattice=2.0, center_y_lattice=3.0)
+    np.testing.assert_array_equal(cropped, field[1:6, 2:7])
+    assert meta["x_slice"] == [2, 7]
+    assert meta["y_slice"] == [1, 6]
+    assert meta["extent_xD_yD"] == [1.0, 3.0, -1.0, 1.0]
+
+
+def test_manifest_includes_npz(tmp_path: Path):
+    repo = tmp_path
+    package = repo / "package"
+    package.mkdir()
+    (package / "data.npz").write_bytes(b"npz")
+    (package / "figure.png").write_bytes(b"png")
+    (package / "ignored.txt").write_text("ignored", encoding="utf-8")
+    (package / "manifest.json").write_text(json.dumps({"source_policy": "test", "summary": {}}), encoding="utf-8")
+    module.update_package_manifest(package, repo)
+    manifest = json.loads((package / "manifest.json").read_text(encoding="utf-8"))
+    paths = {entry["path"] for entry in manifest["artifacts"]}
+    assert "package/data.npz" in paths
+    assert "package/figure.png" in paths
+    assert "package/ignored.txt" not in paths
+    assert manifest["summary"]["flow_field_npz"] == 1

src/SR_analysis/tests/test_formula_provenance.py

@@ -20,10 +20,10 @@ from SR_analysis.utils.provenance import atomic_write_json_below, hash_json, saf
 
 FORMULA_DIR = Path(__file__).parents[1] / "results" / "formulas_v5"
 CANONICAL_FORMULAS = (
-    "ill_1L_sc_front.json",
-    "ill_1L_sc_top.json",
-    "kar_d075_sc_front.json",
-    "kar_d075_sc_top.json",
+    "ill_1L_front.json",
+    "ill_1L_top.json",
+    "kar_d075_front.json",
+    "kar_d075_top.json",
 )
 
 

src/SR_analysis/tests/test_plot_sr_diagnostics.py

@@ -0,0 +1,22 @@
+from __future__ import annotations
+import importlib.util
+from pathlib import Path
+import numpy as np
+
+MODULE_PATH = Path(__file__).resolve().parents[1] / "tools/plot_sr_diagnostics.py"
+spec = importlib.util.spec_from_file_location("plot_sr_diagnostics", MODULE_PATH)
+module = importlib.util.module_from_spec(spec); assert spec.loader; spec.loader.exec_module(module)
+
+
+def test_circular_lag_recovers_roll():
+    reference = np.array([0.0, 1.0, 0.0, -1.0, 0.2])
+    signal = np.roll(reference, 2)
+    assert module.circular_lag(reference, signal) == -2
+
+
+def test_generalization_includes_all_training_conditions():
+    repo = Path(__file__).resolve().parents[3]
+    points = module.generalization_points(repo)
+    training = [p for p in points if p["classification"] == "training"]
+    assert {p["scene"] for p in training} == set(module.TRAINING_SCENES)
+    assert len(training) == 7

src/SR_analysis/tools/export_flow_comparison.py

@@ -0,0 +1,433 @@
+#!/usr/bin/env python3
+"""Export a phase-matched Target/PPO/SR Kármán vorticity comparison.
+
+CUDA-backed modules are imported only after argument validation.  The exporter
+uses the canonical Stage-3 environment and policy constructors for controlled
+runs and mirrors ``build_karman_cloak_env`` exactly for the target-only run.
+"""
+from __future__ import annotations
+
+import argparse
+import hashlib
+import json
+import os
+import platform
+import subprocess
+import sys
+from datetime import datetime, timezone
+from pathlib import Path
+from typing import Any, Mapping, Sequence
+
+import numpy as np
+
+REPO_ROOT = Path(__file__).resolve().parents[3]
+SRC_ROOT = REPO_ROOT / "src"
+for root in (REPO_ROOT, SRC_ROOT):
+    if str(root) not in sys.path:
+        sys.path.insert(0, str(root))
+
+from SR_analysis.configs import FIFO_LEN, LEGACY_CFG_DIR, get_scene  # noqa: E402
+from SR_analysis.stage_3_validate import (  # noqa: E402
+    DATA_TYPE,
+    ValidationPlan,
+    build_karman_environment,
+    build_policy,
+    load_formula_pair,
+    prepare_plan,
+)
+from SR_analysis.utils.provenance import atomic_write_json, hash_file, hash_json  # noqa: E402
+
+STEM = "07_flow_field_comparison_karman_re100"
+SCHEMA = "sr-flow-comparison-v2"
+CANDIDATE_DTYPE = np.dtype(np.float32)
+DEFAULT_PACKAGE = REPO_ROOT / "src/SR_analysis/results/runs/article2-plotting-package-20260721"
+DEFAULT_FORMULAS = REPO_ROOT / "src/SR_analysis/results/runs/article-refit-karman-topology-a-20260718/formulas"
+DEFAULT_ALIGNMENT = DEFAULT_PACKAGE / "phase_alignment.json"
+D_LATTICE = 20.0
+SENSOR_LAYOUT = ("upper_ux", "upper_uy", "center_ux", "center_uy", "lower_ux", "lower_uy")
+
+
+def sha256_bytes(data: bytes) -> str:
+    return hashlib.sha256(data).hexdigest()
+
+
+def git_sha() -> str | None:
+    try:
+        return subprocess.run(["git", "-C", str(REPO_ROOT), "rev-parse", "HEAD"], check=True, capture_output=True, text=True).stdout.strip()
+    except (OSError, subprocess.CalledProcessError):
+        return None
+
+
+def wrapped_phase_difference(a: np.ndarray | float, b: np.ndarray | float) -> np.ndarray:
+    """Signed shortest angular difference ``a-b`` in [-pi, pi]."""
+    delta = np.asarray(a, dtype=np.float64) - np.asarray(b, dtype=np.float64)
+    return np.arctan2(np.sin(delta), np.cos(delta))
+
+
+def standardized_center_phase(trace: np.ndarray, start: int, stop: int) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
+    """Return center-sensor phase from per-trajectory stable-window z scores."""
+    sensors = np.asarray(trace, dtype=np.float64)
+    if sensors.ndim != 2 or sensors.shape[1] < 4 or not 0 <= start < stop <= len(sensors):
+        raise ValueError("trace/window must provide center (ux,uy) channels")
+    center = sensors[start:stop, 2:4]
+    mean = center.mean(axis=0)
+    scale = center.std(axis=0, ddof=0)
+    if np.any(~np.isfinite(scale)) or np.any(scale <= np.finfo(float).eps):
+        raise ValueError("stable-window center-sensor scale is zero or non-finite")
+    z = (center - mean) / scale
+    return np.arctan2(z[:, 1], z[:, 0]), mean, scale
+
+
+def select_joint_phase_match(
+    target_trace: np.ndarray,
+    ppo_trace: np.ndarray,
+    sr_trace: np.ndarray,
+    start: int,
+    stop: int,
+    *,
+    temporal_weight: float = 0.001,
+) -> dict[str, Any]:
+    """Jointly match downstream center-sensor limit-cycle phase.
+
+    Exhaustive stable-window search minimizes the two wrapped phase errors plus
+    ``temporal_weight * (|i_ppo-i_target| + |i_sr-i_target|)``.  Requiring the
+    same local phase direction rejects branch-reversed phase-portrait matches.
+    """
+    if temporal_weight < 0:
+        raise ValueError("temporal_weight must be non-negative")
+    phases, means, scales, directions = [], {}, {}, []
+    for label, trace in (("Target", target_trace), ("PPO", ppo_trace), ("SR", sr_trace)):
+        theta, mean, scale = standardized_center_phase(trace, start, stop)
+        phases.append(theta); means[label] = mean.tolist(); scales[label] = scale.tolist()
+        directions.append(np.sign(np.gradient(np.unwrap(theta))))
+    best = None
+    for ti, target_index in enumerate(range(start, stop)):
+        for pi, ppo_index in enumerate(range(start, stop)):
+            if directions[1][pi] != directions[0][ti]:
+                continue
+            for si, sr_index in enumerate(range(start, stop)):
+                if directions[2][si] != directions[0][ti]:
+                    continue
+                ppo_error = abs(float(wrapped_phase_difference(phases[1][pi], phases[0][ti])))
+                sr_error = abs(float(wrapped_phase_difference(phases[2][si], phases[0][ti])))
+                temporal_distance = abs(ppo_index-target_index) + abs(sr_index-target_index)
+                penalty = temporal_weight * temporal_distance
+                objective = ppo_error + sr_error + penalty
+                key = (objective, temporal_distance, target_index, ppo_index, sr_index)
+                if best is None or key < best[0]:
+                    best = (key, ti, pi, si, ppo_error, sr_error, penalty)
+    if best is None:
+        raise ValueError("no same-direction phase match in stable window")
+    key, ti, pi, si, ppo_error, sr_error, penalty = best
+    indices = {"Target": start+ti, "PPO": start+pi, "SR": start+si}
+    return {
+        "indices": indices,
+        "phase_angles_rad": {"Target": float(phases[0][ti]), "PPO": float(phases[1][pi]), "SR": float(phases[2][si])},
+        "wrapped_angle_errors_rad": {"PPO": ppo_error, "SR": sr_error},
+        "local_phase_direction": {"Target": int(directions[0][ti]), "PPO": int(directions[1][pi]), "SR": int(directions[2][si])},
+        "stable_window_center_mean": means,
+        "stable_window_center_std": scales,
+        "objective_terms": {"angle_error_sum_rad": ppo_error+sr_error, "temporal_distance_samples": int(key[1]), "temporal_weight_rad_per_sample": temporal_weight, "temporal_penalty_rad": penalty, "objective": float(key[0])},
+        "candidate_count_per_trajectory": stop-start,
+    }
+
+
+def orient_vorticity_xy_to_yx(omega_xy: np.ndarray, field_shape: Sequence[int]) -> np.ndarray:
+    """Convert vorticity_from_ddf's (NX, NY) result to image (NY, NX)."""
+    nx, ny = map(int, field_shape[:2])
+    omega = np.asarray(omega_xy)
+    if omega.shape != (nx, ny):
+        raise ValueError(f"expected vorticity shape {(nx, ny)}, got {omega.shape}")
+    return omega.T.copy()
+
+
+def crop_yx(field_yx: np.ndarray, xlim_d: Sequence[float], ylim_d: Sequence[float], *, d_lattice: float = D_LATTICE, center_y_lattice: float | None = None) -> tuple[np.ndarray, dict[str, Any]]:
+    """Crop an image-oriented field using physical x/D and centered y/D."""
+    field = np.asarray(field_yx)
+    if field.ndim != 2:
+        raise ValueError("field must be 2-D in (y, x) order")
+    ny, nx = field.shape
+    cy = (ny - 1) / 2 if center_y_lattice is None else float(center_y_lattice)
+    x0 = max(0, int(np.ceil(float(xlim_d[0]) * d_lattice)))
+    x1 = min(nx, int(np.floor(float(xlim_d[1]) * d_lattice)) + 1)
+    y0 = max(0, int(np.ceil(cy + float(ylim_d[0]) * d_lattice)))
+    y1 = min(ny, int(np.floor(cy + float(ylim_d[1]) * d_lattice)) + 1)
+    if x0 >= x1 or y0 >= y1:
+        raise ValueError("requested crop does not intersect field")
+    extent = ((x0 / d_lattice), ((x1 - 1) / d_lattice), ((y0 - cy) / d_lattice), ((y1 - 1 - cy) / d_lattice))
+    return field[y0:y1, x0:x1].copy(), {"x_slice": [x0, x1], "y_slice": [y0, y1], "extent_xD_yD": list(extent)}
+
+
+def manifest_artifacts(package_dir: Path, repo_root: Path = REPO_ROOT) -> list[dict[str, str]]:
+    suffixes = {".csv", ".png", ".pdf", ".md", ".json", ".npz"}
+    return [
+        {"path": str(path.relative_to(repo_root)), "sha256": hash_file(path)}
+        for path in sorted(package_dir.rglob("*"))
+        if path.is_file() and path.suffix.lower() in suffixes and path.name != "manifest.json"
+    ]
+
+
+def update_package_manifest(package_dir: Path, repo_root: Path = REPO_ROOT) -> None:
+    path = package_dir / "manifest.json"
+    old = json.loads(path.read_text(encoding="utf-8")) if path.is_file() else {}
+    summary = dict(old.get("summary", {}))
+    summary["publication_figures"] = max(5, int(summary.get("publication_figures", 0)))
+    summary["flow_field_npz"] = 1
+    summary["presentation_pages"] = 2
+    atomic_write_json(path, {"schema_version": "sr-plotting-package-v3", "source_policy": old.get("source_policy", "immutable article artifacts; no scientific refit"), "summary": summary, "artifacts": manifest_artifacts(package_dir, repo_root)})
+
+
+def _runtime_cfd() -> tuple[Any, Any, Any]:
+    from LegacyCelerisLab import FlowField
+    from SR_analysis.utils.cfd_interface import load_legacy_configs, vorticity_from_ddf
+    return FlowField, load_legacy_configs, vorticity_from_ddf
+
+
+def _candidate_index_array(start: int, stop: int) -> np.ndarray:
+    """Return the exact trace indices represented by stable-window fields."""
+    if not 0 <= start < stop:
+        raise ValueError("candidate window must satisfy 0 <= start < stop")
+    return np.arange(start, stop, dtype=np.int64)
+
+
+def select_exact_candidate_field(
+    candidates: np.ndarray,
+    candidate_trace_indices: np.ndarray,
+    selected_trace_index: int,
+    *,
+    expected_start: int,
+    expected_stop: int,
+) -> tuple[np.ndarray, int]:
+    """Select a field by its same-run trace index, rejecting mapping drift."""
+    fields = np.asarray(candidates)
+    indices = np.asarray(candidate_trace_indices)
+    expected = _candidate_index_array(expected_start, expected_stop)
+    if fields.ndim != 3 or fields.dtype != CANDIDATE_DTYPE:
+        raise ValueError("candidate fields must be a float32 (sample,y,x) array")
+    if indices.ndim != 1 or not np.issubdtype(indices.dtype, np.integer):
+        raise ValueError("candidate trace indices must be a one-dimensional integer array")
+    if fields.shape[0] != len(indices):
+        raise ValueError("candidate field/index counts differ")
+    if not np.array_equal(indices, expected):
+        raise ValueError("candidate trace-index mapping does not exactly cover the stable window")
+    matches = np.flatnonzero(indices == int(selected_trace_index))
+    if len(matches) != 1:
+        raise IndexError("selected trace index has no unique same-run candidate field")
+    candidate_index = int(matches[0])
+    if candidate_index != int(selected_trace_index) - expected_start:
+        raise AssertionError("candidate offset and trace index disagree")
+    return fields[candidate_index], candidate_index
+
+
+def _allocate_candidates(field_shape: Sequence[int], start: int, stop: int) -> np.ndarray:
+    nx, ny = map(int, field_shape[:2])
+    return np.empty((stop - start, ny, nx), dtype=CANDIDATE_DTYPE)
+
+
+def _store_candidate(
+    candidates: np.ndarray,
+    trace_index: int,
+    start: int,
+    omega_xy: np.ndarray,
+    field_shape: Sequence[int],
+) -> None:
+    candidate_index = trace_index - start
+    if not 0 <= candidate_index < len(candidates):
+        raise IndexError("candidate trace index lies outside allocated stable window")
+    field = D_LATTICE * orient_vorticity_xy_to_yx(omega_xy, field_shape)
+    candidates[candidate_index] = np.asarray(field, dtype=CANDIDATE_DTYPE)
+
+
+def _target_trace_and_candidates(
+    cfg: Mapping[str, Any], device: int, n_samples: int, candidate_start: int
+) -> tuple[np.ndarray, np.ndarray, np.ndarray, tuple[int, int]]:
+    """Run Target once and couple every stable trace sample to its field."""
+    FlowField, load_configs, vorticity_from_ddf = _runtime_cfd()
+    cuda_cfg, field_cfg = load_configs(LEGACY_CFG_DIR)
+    ff = FlowField(field_cfg._replace(viscosity=float(cfg["nu"])), cuda_cfg, device_id=device)
+    try:
+        cy = (ff.FIELD_SHAPE[1] - 1) / 2
+        ff.add_cylinder((10.0 * D_LATTICE, cy, 0.0), D_LATTICE)
+        for y_off in (2.0, 0.0, -2.0):
+            ff.add_sensor((40.0 * D_LATTICE, cy + y_off * D_LATTICE, 0.0), D_LATTICE / 4.0)
+        n_obj = ff.obs.size // 2
+        zero = np.zeros(n_obj, dtype=DATA_TYPE)
+        ff.run(int(4 * ff.FIELD_SHAPE[0] / float(cfg["u0"])), zero)
+        rows: list[np.ndarray] = []
+        candidates = _allocate_candidates(ff.FIELD_SHAPE, candidate_start, n_samples)
+        candidate_indices = _candidate_index_array(candidate_start, n_samples)
+        for index in range(n_samples):
+            ff.run(int(cfg["sample_interval"]), zero)
+            rows.append(ff.obs.copy()[2:8].astype(np.float64))
+            if index >= candidate_start:
+                _store_candidate(candidates, index, candidate_start, vorticity_from_ddf(ff, float(cfg["u0"])), ff.FIELD_SHAPE)
+        if len(rows) != n_samples or len(candidates) != n_samples - candidate_start:
+            raise AssertionError("Target trace/candidate collection is incomplete")
+        return np.asarray(rows), candidates, candidate_indices, tuple(map(int, ff.FIELD_SHAPE[:2]))
+    finally:
+        del ff
+
+
+def _controlled_trace_and_candidates(
+    plan: ValidationPlan, device: int, n_samples: int, candidate_start: int
+) -> tuple[np.ndarray, np.ndarray, np.ndarray, tuple[int, int]]:
+    """Run one controller once and couple stable trace samples to fields."""
+    env = build_karman_environment(plan, device)
+    try:
+        policy = build_policy(plan)
+        raw = np.asarray(env.current_raw, dtype=np.float64)
+        rows: list[np.ndarray] = []
+        candidates = _allocate_candidates(env.ff.FIELD_SHAPE, candidate_start, n_samples)
+        candidate_indices = _candidate_index_array(candidate_start, n_samples)
+        _, _, vorticity_from_ddf = _runtime_cfd()
+        for index in range(n_samples):
+            omega, _, _ = policy.action(raw, index)
+            raw = env.step(omega)
+            policy.observe(omega)
+            rows.append(raw[:6].copy())
+            if index >= candidate_start:
+                _store_candidate(candidates, index, candidate_start, vorticity_from_ddf(env.ff, float(plan.cfg["u0"])), env.ff.FIELD_SHAPE)
+        if len(rows) != n_samples or len(candidates) != n_samples - candidate_start:
+            raise AssertionError("controlled trace/candidate collection is incomplete")
+        return np.asarray(rows), candidates, candidate_indices, tuple(map(int, env.ff.FIELD_SHAPE[:2]))
+    finally:
+        env.close()
+
+
+def _make_plan(scene: str, mode: str, n_steps: int, pair: Any | None, model_device: str) -> ValidationPlan:
+    plan = prepare_plan(scene=scene, mode=mode, n_steps=n_steps, run_id="flow-comparison-in-memory", output_root=REPO_ROOT / ".flow-comparison-unused", formula_pair=pair)
+    return ValidationPlan(**{**plan.__dict__, "cfg": {**plan.cfg, "model_device": model_device}})
+
+
+def _plot(fields: Mapping[str, np.ndarray], crop_meta: Mapping[str, Any], cfg: Mapping[str, Any], output_dir: Path, vmax: float) -> None:
+    import matplotlib
+    matplotlib.use("Agg")
+    import matplotlib.pyplot as plt
+    from matplotlib.patches import Circle
+
+    fig, axes = plt.subplots(1, 3, figsize=(14.5, 4.3), sharex=True, sharey=True, constrained_layout=True)
+    extent = crop_meta["extent_xD_yD"]
+    for ax, label in zip(axes, ("Target", "PPO", "SR")):
+        ax.imshow(fields[label], origin="lower", extent=extent, aspect="equal", cmap="RdBu_r", vmin=-vmax, vmax=vmax, interpolation="nearest")
+        geometry = [(10.0, 0.0, 1.0, "disturbance")] if label == "Target" else [
+            (10.0, 0.0, 1.0, "disturbance"),
+            (float(cfg["pinball_front_x"]), 0.0, 0.5, "front"),
+            (float(cfg["pinball_rear_x"]), 0.75, 0.5, "upper"),
+            (float(cfg["pinball_rear_x"]), -0.75, 0.5, "lower"),
+        ]
+        for x, y, radius, name in geometry:
+            ax.add_patch(Circle((x, y), radius, facecolor="white", edgecolor="black", linewidth=0.8, zorder=4))
+        ax.scatter([40.0] * 3, [2.0, 0.0, -2.0], s=12, marker="x", color="black", linewidths=0.8, zorder=5)
+        ax.set_title(label)
+        ax.set_xlabel(r"$x/D$")
+    axes[0].set_ylabel(r"$y/D$")
+    for suffix in ("png", "pdf"):
+        fig.savefig(output_dir / f"{STEM}.{suffix}", dpi=300 if suffix == "png" else None)
+    plt.close(fig)
+
+
+def build_parser() -> argparse.ArgumentParser:
+    parser = argparse.ArgumentParser(description=__doc__)
+    parser.add_argument("--scene", default="karman_re100", choices=("karman_re100",), help="explicit canonical scene")
+    parser.add_argument("--package-dir", type=Path, default=DEFAULT_PACKAGE)
+    parser.add_argument("--alignment-metadata", type=Path, default=DEFAULT_ALIGNMENT)
+    parser.add_argument("--formula-front", type=Path, default=DEFAULT_FORMULAS / "joint_front.json")
+    parser.add_argument("--formula-rear", type=Path, default=DEFAULT_FORMULAS / "joint_rear_shared_upper.json")
+    parser.add_argument("--temporal-weight", type=float, default=0.001, help="phase objective penalty in radians per sample of temporal separation")
+    parser.add_argument("--device", type=int, default=0, help="logical CFD device after CUDA_VISIBLE_DEVICES masking")
+    parser.add_argument("--model-device", choices=("cpu",), default="cpu")
+    parser.add_argument("--xlim", nargs=2, type=float, default=(7.0, 48.0), metavar=("XMIN_D", "XMAX_D"))
+    parser.add_argument("--ylim", nargs=2, type=float, default=(-6.0, 6.0), metavar=("YMIN_D", "YMAX_D"))
+    parser.add_argument("--replace", action="store_true")
+    return parser
+
+
+def main(argv: Sequence[str] | None = None) -> int:
+    args = build_parser().parse_args(argv)
+    package = args.package_dir.resolve()
+    outputs = [package / f"{STEM}.npz", package / f"{STEM}.json", package / "figures" / f"{STEM}.png", package / "figures" / f"{STEM}.pdf"]
+    existing = [p for p in outputs if p.exists()]
+    if existing and not args.replace:
+        raise FileExistsError("refusing to overwrite: " + ", ".join(map(str, existing)))
+    alignment = json.loads(args.alignment_metadata.read_text(encoding="utf-8"))
+    start, stop = int(alignment["target_index_start_inclusive"]), int(alignment["target_index_stop_exclusive"])
+    cfg = get_scene(args.scene)
+    if cfg["scene_id"] != "karman":
+        raise ValueError("only the Kármán environment is supported")
+    pair = load_formula_pair(args.formula_front, args.formula_rear)
+    n_samples = stop
+    ppo_plan = _make_plan(args.scene, "ppo", n_samples, None, args.model_device)
+    sr_plan = _make_plan(args.scene, "pysr", n_samples, pair, args.model_device)
+
+    # Each condition is initialized and run exactly once. The three float32
+    # stable-window buffers remain in host memory until joint phase selection.
+    target_trace, target_candidates, target_candidate_indices, shape = _target_trace_and_candidates(
+        cfg, args.device, n_samples, start
+    )
+    ppo_trace, ppo_candidates, ppo_candidate_indices, ppo_shape = _controlled_trace_and_candidates(
+        ppo_plan, args.device, n_samples, start
+    )
+    sr_trace, sr_candidates, sr_candidate_indices, sr_shape = _controlled_trace_and_candidates(
+        sr_plan, args.device, n_samples, start
+    )
+    if shape != ppo_shape or shape != sr_shape:
+        raise ValueError("field shapes differ between cases")
+
+    phase_match = select_joint_phase_match(target_trace, ppo_trace, sr_trace, start, stop, temporal_weight=args.temporal_weight)
+    selected_indices = phase_match["indices"]
+    candidate_sets = {
+        "Target": (target_candidates, target_candidate_indices),
+        "PPO": (ppo_candidates, ppo_candidate_indices),
+        "SR": (sr_candidates, sr_candidate_indices),
+    }
+    full: dict[str, np.ndarray] = {}
+    selected_candidate_indices: dict[str, int] = {}
+    for label, (candidates, trace_indices) in candidate_sets.items():
+        field, candidate_index = select_exact_candidate_field(
+            candidates, trace_indices, selected_indices[label], expected_start=start, expected_stop=stop
+        )
+        full[label] = field.copy()
+        selected_candidate_indices[label] = candidate_index
+    candidate_field_bytes = int(sum(candidates.nbytes for candidates, _ in candidate_sets.values()))
+    cropped: dict[str, np.ndarray] = {}
+    crop_meta = None
+    for label, field in full.items():
+        cropped[label], meta = crop_yx(field, args.xlim, args.ylim)
+        crop_meta = crop_meta or meta
+        if meta != crop_meta:
+            raise AssertionError("inconsistent crop metadata")
+    finite = np.concatenate([np.abs(value[np.isfinite(value)]) for value in cropped.values()])
+    if finite.size == 0:
+        raise FloatingPointError("captured fields contain no finite vorticity")
+    vmax = float(np.percentile(finite, 99.5))
+    if not np.isfinite(vmax) or vmax <= 0:
+        raise FloatingPointError("invalid shared color normalization")
+
+    package.mkdir(parents=True, exist_ok=True)
+    (package / "figures").mkdir(parents=True, exist_ok=True)
+    npz_path = outputs[0]
+    np.savez_compressed(npz_path, target_vorticity_yx=full["Target"], ppo_vorticity_yx=full["PPO"], sr_vorticity_yx=full["SR"], target_sensors=target_trace, ppo_sensors=ppo_trace, sr_sensors=sr_trace, selected_indices=np.asarray([selected_indices[x] for x in ("Target", "PPO", "SR")]), phase_angles_rad=np.asarray([phase_match["phase_angles_rad"][x] for x in ("Target", "PPO", "SR")]), wrapped_angle_errors_rad=np.asarray([0.0, phase_match["wrapped_angle_errors_rad"]["PPO"], phase_match["wrapped_angle_errors_rad"]["SR"]]), local_phase_direction=np.asarray([phase_match["local_phase_direction"][x] for x in ("Target", "PPO", "SR")]), selected_candidate_indices=np.asarray([selected_candidate_indices[x] for x in ("Target", "PPO", "SR")]), extent_xD_yD=np.asarray(crop_meta["extent_xD_yD"]), crop_x_slice=np.asarray(crop_meta["x_slice"]), crop_y_slice=np.asarray(crop_meta["y_slice"]), field_shape_xy=np.asarray(shape), sensor_layout=np.asarray(SENSOR_LAYOUT))
+    _plot(cropped, crop_meta, cfg, package / "figures", vmax)
+
+    dt = float(cfg["control_dt"])
+    model_path = Path(ppo_plan.model_path) if ppo_plan.model_path else None
+    metadata = {
+        "schema_version": SCHEMA, "scene": args.scene, "description": "single deterministic downstream center-sensor limit-cycle phase snapshot; not full-field identity or an ensemble",
+        "phase_matching": {"method": "one run per condition with same-sample stable-window field capture; exhaustive joint phase search", "phase_definition": "theta=atan2(z(center_uy), z(center_ux)); each trajectory standardized separately over stable window", "objective_formula": "|wrap(theta_PPO-theta_Target)| + |wrap(theta_SR-theta_Target)| + w*(|i_PPO-i_Target|+|i_SR-i_Target|)", "target_window": [start, stop], "direction_check": "all selected local unwrapped-phase derivatives have the same sign", "figure06_alignment_context_only": alignment, **phase_match},
+        "selected": {label: {"index": int(selected_indices[label]), "candidate_index": int(selected_candidate_indices[label]), "t_D_over_U0": float(selected_indices[label] * dt), "phase_angle_rad": float(phase_match["phase_angles_rad"][label]), "wrapped_phase_error_rad": 0.0 if label == "Target" else float(phase_match["wrapped_angle_errors_rad"][label])} for label in ("Target", "PPO", "SR")},
+        "sample_coupling": {"contract": "for every condition, the selected field and phase-diagnostic sensor values come from the same CFD sample in the same run", "candidate_trace_index_mapping": "candidate_index = trace_index - stable_window_start; exact contiguous mapping asserted before selection", "candidate_window": [start, stop], "candidate_count_per_condition": stop-start, "candidate_dtype": CANDIDATE_DTYPE.name, "candidate_field_shape_yx": [int(shape[1]), int(shape[0])], "candidate_host_memory_bytes": candidate_field_bytes, "candidate_host_memory_mib": candidate_field_bytes / (1024**2), "disk_contract": "only the three selected full fields and complete sensor traces are stored; stable-window candidates are memory-only"},
+        "field": {"quantity": "omega_z D/U0", "source": "vorticity_from_ddf", "source_shape_order": "(NX,NY)", "stored_shape_order": "(NY,NX)", "crop": crop_meta, "shared_symmetric_vmax_percentile": {"percentile": 99.5, "vmax": vmax}, "geometry_D": {"disturbance": [10.0, 0.0, 1.0], "pinball": [[cfg["pinball_front_x"], 0.0, 0.5], [cfg["pinball_rear_x"], 0.75, 0.5], [cfg["pinball_rear_x"], -0.75, 0.5]], "sensors": [[40.0, 2.0], [40.0, 0.0], [40.0, -2.0]]}},
+        "contracts": {"target": "exact build_karman_cloak_env geometry, stabilization, and sample interval; target has no pinball", "controlled": "stage_3_validate.prepare_plan/build_karman_environment/build_policy", "ppo_model_device": args.model_device, "cfd_logical_device": args.device},
+        "hashes": {"config": hash_json(cfg), "formula_front": hash_file(args.formula_front), "formula_rear": hash_file(args.formula_rear), "formula_pair": pair.pair_hash, "ppo_model": hash_file(model_path) if model_path and model_path.is_file() else None, "alignment_metadata": hash_file(args.alignment_metadata), "npz": hash_file(npz_path), "png": hash_file(outputs[2]), "pdf": hash_file(outputs[3]), "exporter": hash_file(Path(__file__))},
+        "paths": {"npz": str(npz_path.relative_to(REPO_ROOT)), "formula_front": str(args.formula_front.resolve()), "formula_rear": str(args.formula_rear.resolve()), "ppo_model": str(model_path) if model_path else None},
+        "provenance": {"git_sha": git_sha(), "command": " ".join(sys.argv), "created_utc": datetime.now(timezone.utc).isoformat(), "python": sys.version, "platform": platform.platform(), "numpy": np.__version__, "CUDA_VISIBLE_DEVICES": os.environ.get("CUDA_VISIBLE_DEVICES")},
+    }
+    metadata["record_hash"] = hash_json(metadata)
+    atomic_write_json(outputs[1], metadata)
+    update_package_manifest(package)
+    print(json.dumps({"outputs": [str(p) for p in outputs], "selected": metadata["selected"], "vmax": vmax}, indent=2))
+    return 0
+
+
+if __name__ == "__main__":
+    raise SystemExit(main())

src/SR_analysis/tools/plot_sr_diagnostics.py

@@ -1,11 +1,13 @@
 #!/usr/bin/env python3
-"""Render compact diagnostic figures from the SR plotting CSV packages."""
+"""Render publication SR figures and tables from frozen article CSV/JSON artifacts."""
 from __future__ import annotations
 
 import argparse
 import csv
+import hashlib
+import json
 from pathlib import Path
-from typing import Sequence
+from typing import Any, Mapping, Sequence
 
 import matplotlib.pyplot as plt
 import numpy as np
@@ -15,6 +17,11 @@ TRAINING_SCENES = (
     "illusion_0.75L", "illusion_1L", "illusion_1.5L",
 )
 LABELS = ("K50", "K100", "K200", "K400", "I0.75L", "I1L", "I1.5L")
+ACTIONS = ("front", "upper", "lower")
+PACKAGE = "article2-plotting-package-20260721"
+STANDARD = "article2-timeseries-csv-20260720"
+LONG_SR = "article2-long-timeseries-csv-20260720"
+COLORS = {"Target": "#222222", "PPO": "#2878B5", "SR": "#D95319"}
 
 
 def read_rows(path: Path) -> list[dict[str, str]]:
@@ -22,187 +29,249 @@ def read_rows(path: Path) -> list[dict[str, str]]:
         return list(csv.DictReader(handle))
 
 
+def write_rows(path: Path, fieldnames: Sequence[str], rows: Sequence[Mapping[str, Any]]) -> None:
+    path.parent.mkdir(parents=True, exist_ok=True)
+    with path.open("w", newline="", encoding="utf-8") as handle:
+        writer = csv.DictWriter(handle, fieldnames=fieldnames, lineterminator="\n")
+        writer.writeheader()
+        writer.writerows(rows)
+
+
 def final_similarity(root: Path, package: str, scene: str, source: str) -> float:
-    rows = read_rows(root / package / scene / f"{source}_dtw_convergence.csv")
-    return float(rows[-1]["similarity"])
+    return float(read_rows(root / package / scene / f"{source}_dtw_convergence.csv")[-1]["similarity"])
 
 
 def save(fig: plt.Figure, output: Path, stem: str) -> None:
     output.mkdir(parents=True, exist_ok=True)
-    fig.savefig(output / f"{stem}.png", dpi=220, bbox_inches="tight")
+    fig.savefig(output / f"{stem}.png", dpi=300, bbox_inches="tight")
     fig.savefig(output / f"{stem}.pdf", bbox_inches="tight")
     plt.close(fig)
 
 
 def plot_performance(root: Path, output: Path) -> None:
-    standard = "article2-timeseries-csv-20260720"
-    long_sr = "article2-long-timeseries-csv-20260720"
-    ppo200 = [final_similarity(root, standard, scene, "ppo") for scene in TRAINING_SCENES]
-    sr200 = [final_similarity(root, standard, scene, "sr") for scene in TRAINING_SCENES]
-    sr400 = [final_similarity(root, long_sr, scene, "sr") for scene in TRAINING_SCENES]
-    long_ppo_root = root / "article2-plotting-package-20260721/long_ppo"
-    ppo400 = []
-    for scene in TRAINING_SCENES:
-        path = long_ppo_root / scene / "ppo_dtw_convergence.csv"
-        ppo400.append(float(read_rows(path)[-1]["similarity"]) if path.is_file() else np.nan)
+    ppo200 = np.array([final_similarity(root, STANDARD, s, "ppo") for s in TRAINING_SCENES])
+    sr200 = np.array([final_similarity(root, STANDARD, s, "sr") for s in TRAINING_SCENES])
+    ppo400 = np.array([final_similarity(root / PACKAGE / "long_ppo", "", s, "ppo") for s in TRAINING_SCENES])
+    sr400 = np.array([final_similarity(root, LONG_SR, s, "sr") for s in TRAINING_SCENES])
+    x = np.arange(len(LABELS)); width = 0.36
+    fig, axes = plt.subplots(1, 2, figsize=(12.2, 4.7), gridspec_kw={"width_ratios": (1.15, 1)})
+    axes[0].bar(x - width / 2, ppo400, width, label="PPO — 400 control steps", color=COLORS["PPO"])
+    axes[0].bar(x + width / 2, sr400, width, label="SR — 400 control steps", color=COLORS["SR"])
+    axes[0].set_ylabel("Legacy DTW similarity"); axes[0].set_ylim(0.75, 1.005)
+    axes[0].set_title("(a) Closed-loop performance at 400 control steps")
+    axes[0].legend(frameon=False, fontsize=9)
+    axes[1].axhline(0, color="0.25", linewidth=1)
+    axes[1].scatter(x - 0.08, ppo400 - ppo200, label="PPO: 400 − 200 control steps", marker="o", s=48, color=COLORS["PPO"])
+    axes[1].scatter(x + 0.08, sr400 - sr200, label="SR: 400 − 200 control steps", marker="s", s=48, color=COLORS["SR"])
+    axes[1].set_ylabel("Change in legacy DTW similarity")
+    axes[1].set_title("(b) Duration stability from 200 to 400 steps")
+    axes[1].legend(frameon=False, fontsize=9)
+    for ax in axes:
+        ax.set_xticks(x, LABELS); ax.set_xlabel("Training case"); ax.grid(axis="y", alpha=0.25)
+    fig.tight_layout(); save(fig, output, "01_training_case_performance")
 
-    x = np.arange(len(LABELS))
-    width = 0.2
-    fig, ax = plt.subplots(figsize=(11, 4.8))
-    ax.bar(x - 1.5 * width, ppo200, width, label="PPO 200")
-    ax.bar(x - 0.5 * width, sr200, width, label="SR 200")
-    ax.bar(x + 0.5 * width, ppo400, width, label="PPO 400")
-    ax.bar(x + 1.5 * width, sr400, width, label="SR 400")
-    ax.set_ylabel("Legacy DTW similarity")
-    ax.set_xlabel("Training case")
-    ax.set_xticks(x, LABELS)
-    ax.set_ylim(0.65, 1.01)
-    ax.grid(axis="y", alpha=0.25)
-    ax.legend(ncol=4, frameon=False)
-    ax.set_title("Closed-loop PPO and symbolic-controller performance")
-    fig.tight_layout()
-    save(fig, output, "01_training_case_performance")
+
+def generalization_points(repo_root: Path) -> list[dict[str, Any]]:
+    path = repo_root / "src/SR_analysis/results/runs/article2-generalization-summary-20260720/generalization.csv"
+    points: list[dict[str, Any]] = [dict(row) for row in read_rows(path)]
+    run_root = repo_root / "src/SR_analysis/results/runs"
+    for scene in TRAINING_SCENES:
+        objective = "karman" if scene.startswith("karman") else "illusion"
+        parameter = float(scene.removeprefix("karman_re").removeprefix("illusion_").removesuffix("L"))
+        points.append({"objective": objective, "scene": scene, "parameter": parameter,
+                       "classification": "training", "legacy_dtw": final_similarity(run_root, STANDARD, scene, "sr")})
+    return points
 
 
 def plot_generalization(repo_root: Path, output: Path) -> None:
-    path = repo_root / "src/SR_analysis/results/runs/article2-generalization-summary-20260720/generalization.csv"
-    rows = read_rows(path)
-    fig, axes = plt.subplots(1, 2, figsize=(10.5, 4.2))
-    for ax, objective, xlabel in zip(axes, ("karman", "illusion"), ("Code Reynolds label", "Target-size label L")):
-        subset = sorted((row for row in rows if row["objective"] == objective), key=lambda row: float(row["parameter"]))
-        interpolation = [row for row in subset if row["classification"] == "interpolation"]
-        extrapolation = [row for row in subset if row["classification"] == "extrapolation"]
-        ax.plot([float(r["parameter"]) for r in subset], [float(r["legacy_dtw"]) for r in subset], color="0.6", linewidth=1)
-        ax.scatter([float(r["parameter"]) for r in interpolation], [float(r["legacy_dtw"]) for r in interpolation], label="Interpolation", s=48)
-        ax.scatter([float(r["parameter"]) for r in extrapolation], [float(r["legacy_dtw"]) for r in extrapolation], label="Extrapolation", marker="s", s=48)
-        ax.set_xlabel(xlabel)
-        ax.set_ylabel("Legacy DTW similarity")
-        ax.set_ylim(0.7, 1.01)
-        ax.grid(alpha=0.25)
-        ax.legend(frameon=False)
-        ax.set_title("Kármán unseen Re" if objective == "karman" else "Illusion unseen target size")
-    fig.tight_layout()
-    save(fig, output, "02_pointwise_generalization")
+    rows = generalization_points(repo_root)
+    fig, axes = plt.subplots(1, 2, figsize=(11, 4.4))
+    styles = {"training": ("o", "#222222", "Training condition"),
+              "interpolation": ("^", "#2E8B57", "Unseen interpolation"),
+              "extrapolation": ("s", "#A23B72", "Unseen extrapolation")}
+    for ax, objective, xlabel in zip(axes, ("karman", "illusion"), ("Reynolds-number label", "Target-size label, L")):
+        subset = [r for r in rows if r["objective"] == objective]
+        for category, (marker, color, label) in styles.items():
+            selected = sorted((r for r in subset if r["classification"] == category), key=lambda r: float(r["parameter"]))
+            ax.scatter([float(r["parameter"]) for r in selected], [float(r["legacy_dtw"]) for r in selected],
+                       marker=marker, color=color, edgecolor="white", linewidth=0.6, s=65, label=label, zorder=3)
+        ax.set_xlabel(xlabel); ax.set_ylabel("Legacy DTW similarity"); ax.set_ylim(0.7, 1.01)
+        ax.grid(alpha=0.25); ax.legend(frameon=False, fontsize=8.5)
+        ax.set_title("Kármán conditions" if objective == "karman" else "Illusion conditions")
+    fig.suptitle("Frozen SR controller: training and pointwise unseen 200-step realizations", y=1.01)
+    fig.tight_layout(); save(fig, output, "02_pointwise_generalization")
 
 
 def plot_steady(repo_root: Path, output: Path) -> None:
-    path = repo_root / "src/SR_analysis/results/runs/article2-steady-analysis-20260720/steady_rotation_sweep.csv"
-    rows = read_rows(path)
-    magnitude = np.asarray([float(row["rear_alpha_magnitude"]) for row in rows])
-    similarity = np.asarray([float(row["legacy_dtw"]) for row in rows])
-    error = np.asarray([float(row["mean_abs_sensor_error"]) for row in rows])
-    fig, ax = plt.subplots(figsize=(7.2, 4.5))
-    secondary = ax.twinx()
+    rows = read_rows(repo_root / "src/SR_analysis/results/runs/article2-steady-analysis-20260720/steady_rotation_sweep.csv")
+    magnitude = np.asarray([float(r["rear_alpha_magnitude"]) for r in rows])
+    similarity = np.asarray([float(r["legacy_dtw"]) for r in rows])
+    error = np.asarray([float(r["mean_abs_sensor_error"]) for r in rows])
+    fig, ax = plt.subplots(figsize=(7.2, 4.5)); secondary = ax.twinx()
     ax.plot(magnitude, similarity, marker="o", label="DTW similarity")
     secondary.plot(magnitude, error, marker="s", linestyle="--", label="Mean sensor error", color="tab:orange")
     ax.axvline(3.4312086805, linestyle=":", color="0.35", label="SR rear constant 3.4312")
-    ax.set_xlabel("Rear counter-rotation magnitude |α|")
-    ax.set_ylabel("Legacy DTW similarity")
-    secondary.set_ylabel("Mean absolute sensor error")
-    ax.set_ylim(0.6, 1.01)
-    ax.grid(alpha=0.25)
-    lines = ax.get_lines() + secondary.get_lines()
-    ax.legend(lines, [line.get_label() for line in lines], frameon=False, loc="center right")
-    ax.set_title("Steady-cloak constant-rotation calibration")
-    fig.tight_layout()
+    ax.set_xlabel("Rear counter-rotation magnitude |α|"); ax.set_ylabel("Legacy DTW similarity")
+    secondary.set_ylabel("Mean absolute sensor error"); ax.set_ylim(0.6, 1.01); ax.grid(alpha=0.25)
+    lines = ax.get_lines() + secondary.get_lines(); ax.legend(lines, [line.get_label() for line in lines], frameon=False, loc="center right")
+    ax.set_title("Steady-cloak constant-rotation calibration"); fig.tight_layout()
     save(fig, output, "03_steady_rotation_calibration")
 
 
-def plot_offline_residuals(root: Path, output: Path) -> None:
-    rows = read_rows(root / "article2-plotting-package-20260721/offline_predictions.csv")
-    rmse = {scene: [] for scene in TRAINING_SCENES}
+def build_offline_table(root: Path, tables: Path) -> None:
+    source = root / PACKAGE / "offline_predictions.csv"; rows = read_rows(source); output = []
     for scene in TRAINING_SCENES:
-        subset = [row for row in rows if row["scene"] == scene]
-        for action in ("front", "upper", "lower"):
-            residual = np.asarray([float(row[f"residual_alpha_{action}"]) for row in subset])
-            rmse[scene].append(float(np.sqrt(np.mean(residual**2))))
-    fig, ax = plt.subplots(figsize=(10.5, 4.5))
-    x = np.arange(len(LABELS))
-    width = 0.25
-    for index, action in enumerate(("Front", "Upper", "Lower")):
-        ax.bar(x + (index - 1) * width, [rmse[scene][index] for scene in TRAINING_SCENES], width, label=action)
-    ax.set_xticks(x, LABELS)
-    ax.set_xlabel("Training case")
-    ax.set_ylabel("Offline action RMSE in α")
-    ax.set_title("SR imitation error on PPO-visited causal states")
-    ax.grid(axis="y", alpha=0.25)
-    ax.legend(frameon=False, ncol=3)
-    fig.tight_layout()
-    save(fig, output, "04_offline_action_rmse")
+        subset = [r for r in rows if r["scene"] == scene]
+        values = {a: float(np.sqrt(np.mean([float(r[f"residual_alpha_{a}"]) ** 2 for r in subset]))) for a in ACTIONS}
+        aggregate = float(np.sqrt(np.mean([float(r[f"residual_alpha_{a}"]) ** 2 for r in subset for a in ACTIONS])))
+        output.append({"objective": subset[0]["objective"], "training_case": scene, "n_causally_aligned_states": len(subset),
+                       **{f"rmse_alpha_{a}": f"{values[a]:.8f}" for a in ACTIONS}, "aggregate_rmse_alpha": f"{aggregate:.8f}"})
+    fields = ("objective", "training_case", "n_causally_aligned_states", *(f"rmse_alpha_{a}" for a in ACTIONS), "aggregate_rmse_alpha")
+    write_rows(tables / "offline_action_rmse.csv", fields, output)
+    lines = ["# Offline SR-vs-PPO next-action RMSE", "",
+             "RMSE is evaluated on causally aligned PPO-visited states: each post-action state predicts the next PPO action. This is an offline imitation diagnostic, not closed-loop performance.", "",
+             "| Objective | Training case | States | Front RMSE | Upper RMSE | Lower RMSE | Aggregate RMSE |",
+             "|---|---|---:|---:|---:|---:|---:|"]
+    for r in output:
+        lines.append(f"| {r['objective']} | {r['training_case']} | {r['n_causally_aligned_states']} | {r['rmse_alpha_front']} | {r['rmse_alpha_upper']} | {r['rmse_alpha_lower']} | {r['aggregate_rmse_alpha']} |")
+    lines += ["", "All action errors are in dimensionless surface-speed α. Aggregate RMSE pools the three action channels and all aligned states within a case."]
+    (tables / "offline_action_rmse.md").write_text("\n".join(lines) + "\n", encoding="utf-8")
 
 
-def plot_ablation(root: Path, output: Path) -> None:
-    rows = read_rows(root / "article2-plotting-package-20260721/ablation_summary.csv")
-    karman_variants = ("k_front0", "k_rear0", "k_rear1")
-    karman_scenes = ("karman_re50", "karman_re100", "karman_re200", "karman_re400")
-    illusion_variants = ("i_front0", "i_front1", "i_rear0", "i_rear1")
-    illusion_scenes = ("illusion_0.75L", "illusion_1L", "illusion_1.5L")
-    fig, axes = plt.subplots(1, 2, figsize=(11.5, 4.5))
-    for ax, variants, scenes, title in (
-        (axes[0], karman_variants, karman_scenes, "Kármán 40-step term deletion"),
-        (axes[1], illusion_variants, illusion_scenes, "Illusion 40-step term deletion"),
-    ):
-        x = np.arange(len(scenes))
-        width = 0.8 / len(variants)
-        for index, variant in enumerate(variants):
-            values = [float(next(row["similarity"] for row in rows if row["variant"] == variant and row["scene"] == scene)) for scene in scenes]
-            ax.bar(x + (index - (len(variants) - 1) / 2) * width, values, width, label=variant)
-        ax.set_xticks(x, [scene.replace("karman_", "").replace("illusion_", "") for scene in scenes])
-        ax.set_xlabel("Case")
-        ax.set_ylabel("Legacy DTW similarity")
-        ax.set_ylim(0.6, 1.0)
-        ax.grid(axis="y", alpha=0.25)
-        ax.legend(frameon=False, fontsize=8)
-        ax.set_title(title)
-    fig.tight_layout()
-    save(fig, output, "05_term_deletion")
+def _load_json(path: Path) -> dict[str, Any]:
+    return json.loads(path.read_text(encoding="utf-8"))
 
 
-def plot_example_timeseries(root: Path, output: Path) -> None:
-    package = root / "article2-timeseries-csv-20260720/karman_re100"
-    sr = read_rows(package / "sr_wide.csv")
-    ppo = read_rows(package / "ppo_wide.csv")
-    target = read_rows(package / "target_wide.csv")
-    fig, axes = plt.subplots(2, 1, figsize=(10.5, 6.2), sharex=False)
-    axes[0].plot([float(r["t_D"]) for r in target], [float(r["sensors_center_ux"]) for r in target], label="Target reference")
-    axes[0].plot([float(r["t_D"]) for r in ppo], [float(r["sensors_center_ux"]) for r in ppo], label="PPO")
-    axes[0].plot([float(r["t_D"]) for r in sr], [float(r["sensors_center_ux"]) for r in sr], label="SR")
-    axes[0].set_ylabel("Center sensor ux / U0")
-    axes[0].set_xlabel("t U0 / D")
-    axes[0].legend(frameon=False, ncol=3)
-    axes[0].grid(alpha=0.25)
-    axes[0].set_title("Kármán Re100 observer signal")
-    for action, label in (("front", "Front"), ("upper", "Upper"), ("lower", "Lower")):
-        axes[1].plot([float(r["t_D"]) for r in sr], [float(r[f"actions_alpha_{action}"]) for r in sr], label=label)
-    axes[1].set_ylabel("SR action α")
-    axes[1].set_xlabel("t U0 / D")
-    axes[1].legend(frameon=False, ncol=3)
-    axes[1].grid(alpha=0.25)
-    fig.tight_layout()
-    save(fig, output, "06_example_timeseries_karman_re100")
+def build_term_deletion_table(repo_root: Path, root: Path, tables: Path) -> None:
+    formula_root = repo_root / "src/SR_analysis/results/runs/article-ablation-formulas-v2-20260718/formulas"
+    parent_root = repo_root / "src/SR_analysis/results/runs"
+    specs = {
+        "k_front0": ("karman", "front odd projection", "karman_front__delete_t0.json", "article-refit-karman-topology-a-20260718/formulas/joint_front.json", "article-L2-karman-20260718"),
+        "k_rear0": ("karman", "upper/lower shared-symmetry rear", "karman_rear__delete_t0.json", "article-refit-karman-topology-a-20260718/formulas/joint_rear_shared_upper.json", "article-L2-karman-20260718"),
+        "k_rear1": ("karman", "upper/lower shared-symmetry rear", "karman_rear__delete_t1.json", "article-refit-karman-topology-a-20260718/formulas/joint_rear_shared_upper.json", "article-L2-karman-20260718"),
+        "i_front0": ("illusion", "front odd projection", "illusion_front__delete_t0.json", "article-refit-illusion-topology-a-20260718/formulas/joint_front.json", "article-L2-illusionA-20260718"),
+        "i_front1": ("illusion", "front odd projection", "illusion_front__delete_t1.json", "article-refit-illusion-topology-a-20260718/formulas/joint_front.json", "article-L2-illusionA-20260718"),
+        "i_rear0": ("illusion", "upper/lower shared-symmetry rear", "illusion_rear__delete_t0.json", "article-refit-illusion-topology-a-20260718/formulas/joint_rear_shared_upper.json", "article-L2-illusionA-20260718"),
+        "i_rear1": ("illusion", "upper/lower shared-symmetry rear", "illusion_rear__delete_t1.json", "article-refit-illusion-topology-a-20260718/formulas/joint_rear_shared_upper.json", "article-L2-illusionA-20260718"),
+    }
+    ablations = read_rows(root / PACKAGE / "ablation_summary.csv"); output = []
+    for variant, (objective, head, variant_file, parent_file, parent_run) in specs.items():
+        variant_json = _load_json(formula_root / variant_file); parent_json = _load_json(parent_root / parent_file)
+        deleted = variant_json["variant_metadata"]["term_identity"]
+        variant_rows = [r for r in ablations if r["variant"] == variant]
+        similarities = [float(r["similarity"]) for r in variant_rows]
+        parent_metrics = {}
+        for path in (parent_root / parent_run / "validations").glob("*.json"):
+            record = _load_json(path); parent_metrics[record["scene"]] = float(record["metrics"]["legacy_reference_cycle_vs_last_recorded_cycle"]["similarity"])
+        for r in sorted(variant_rows, key=lambda x: TRAINING_SCENES.index(x["scene"])):
+            parent_value = parent_metrics[r["scene"]]; value = float(r["similarity"])
+            output.append({"objective": objective, "head": head, "variant": variant,
+                           "parent_formula": parent_json["deployment_expression"], "deleted_term": deleted,
+                           "retained_formula": variant_json["deployment_expression"], "case": r["scene"], "steps": r["steps"],
+                           "legacy_dtw": f"{value:.8f}", "parent_legacy_dtw_40_step": f"{parent_value:.8f}",
+                           "delta_vs_parent_40_step": f"{value-parent_value:+.8f}", "variant_mean": f"{np.mean(similarities):.8f}",
+                           "variant_min": f"{np.min(similarities):.8f}"})
+    fields = tuple(output[0]); write_rows(tables / "term_deletion.csv", fields, output)
+    lines = ["# Forty-step closed-loop term deletion", "",
+             "Every comparison uses the same 40-control-step window and the legacy DTW metric. Parent values come from the scientifically comparable parent L2 runs; deltas are deletion minus parent. Formula identities are read from the canonical parent and deletion JSON artifacts.", "",
+             "Front formulas are scalar generators deployed through the odd projection `α_F(x) = ½[f(x) − f(Gx)]`. Rear formulas define the upper action; the lower action is mapped by `α_L(x) = −α_U(Gx)`.", "",
+             "| Objective | Head | Variant | Parent formula | Deleted term | Retained formula | Case | DTW | Parent DTW | Δ | Mean | Min |",
+             "|---|---|---|---|---|---|---|---:|---:|---:|---:|---:|"]
+    for r in output:
+        fmt = lambda s: f"`{s.replace('*', '·')}`"
+        lines.append(f"| {r['objective']} | {r['head']} | {r['variant']} | {fmt(r['parent_formula'])} | {fmt(r['deleted_term'])} | {fmt(r['retained_formula'])} | {r['case']} | {r['legacy_dtw']} | {r['parent_legacy_dtw_40_step']} | {r['delta_vs_parent_40_step']} | {r['variant_mean']} | {r['variant_min']} |")
+    lines += ["", "Interpretation: Kármán deletion effects are term-dependent, with deleting the rear constant producing the largest degradation; the tested front term is weak over this short window. Illusion deletions remain stable and comparatively close to their parent, so these runs do not establish uniqueness of every term. Means and minima summarize the listed training cases only."]
+    (tables / "term_deletion.md").write_text("\n".join(lines) + "\n", encoding="utf-8")
 
 
-def build_parser() -> argparse.ArgumentParser:
-    parser = argparse.ArgumentParser(description=__doc__)
-    parser.add_argument("--repo-root", type=Path, default=Path(__file__).resolve().parents[3])
-    parser.add_argument("--output-dir", type=Path, required=True)
-    return parser
+def circular_lag(reference: np.ndarray, signal: np.ndarray) -> int:
+    """Return deterministic roll applied to signal to maximize normalized circular correlation."""
+    ref = np.asarray(reference, dtype=float) - np.mean(reference); sig = np.asarray(signal, dtype=float) - np.mean(signal)
+    if len(ref) != len(sig) or len(ref) == 0: raise ValueError("signals must have equal non-zero length")
+    scores = np.array([np.dot(ref, np.roll(sig, lag)) for lag in range(len(ref))]) / (np.linalg.norm(ref) * np.linalg.norm(sig))
+    lag = int(np.flatnonzero(np.isclose(scores, scores.max(), rtol=0, atol=1e-12))[0])
+    return lag if lag <= len(ref) // 2 else lag - len(ref)
+
+
+def plot_example_timeseries(root: Path, output: Path, package_dir: Path) -> dict[str, Any]:
+    source = root / STANDARD / "karman_re100"
+    data = {name: read_rows(source / f"{name.lower()}_wide.csv") for name in ("Target", "PPO", "SR")}
+    start, stop = 96, 146  # 50 samples, approximately three target cycles, late in the 150-row target record.
+    phase = np.linspace(0, 6 * np.pi, stop - start, endpoint=False)
+    signal_name = "sensors_center_uy"
+    reference = np.array([float(r[signal_name]) for r in data["Target"][start:stop]])
+    lags = {"Target": 0}
+    aligned: dict[str, list[dict[str, str]]] = {"Target": data["Target"][start:stop]}
+    for name in ("PPO", "SR"):
+        segment = data[name][start:stop]
+        lag = circular_lag(reference, np.array([float(r[signal_name]) for r in segment]))
+        lags[name] = lag; aligned[name] = list(np.roll(np.asarray(segment, dtype=object), lag))
+    fig, axes = plt.subplots(2, 3, figsize=(12, 7.2))
+    sensors = (("upper", "Upper sensor"), ("center", "Center sensor"), ("lower", "Lower sensor"))
+    styles = {"Target": "-", "PPO": "--", "SR": "-."}
+    for col, (sensor, title) in enumerate(sensors):
+        for name in ("Target", "PPO", "SR"):
+            rows = aligned[name]
+            axes[0, col].plot([float(r[f"sensors_{sensor}_ux"]) for r in rows], [float(r[f"sensors_{sensor}_uy"]) for r in rows],
+                              color=COLORS[name], linestyle=styles[name], linewidth=1.7, label=name)
+        axes[0, col].set_xlabel(r"$u/U_0$"); axes[0, col].set_ylabel(r"$v/U_0$"); axes[0, col].set_title(title); axes[0, col].grid(alpha=0.22)
+        axes[0, col].scatter(float(aligned["Target"][0][f"sensors_{sensor}_ux"]), float(aligned["Target"][0][f"sensors_{sensor}_uy"]), color=COLORS["Target"], marker="o", s=20, zorder=4)
+    axes[0, 0].legend(frameon=False, ncol=3, fontsize=9)
+    for col, action in enumerate(ACTIONS):
+        for name in ("PPO", "SR"):
+            axes[1, col].plot(phase / (2 * np.pi), [float(r[f"actions_alpha_{action}"]) for r in aligned[name]], color=COLORS[name], linestyle=styles[name], linewidth=1.6, label=name)
+        axes[1, col].set_xlabel("Target phase / 2π"); axes[1, col].set_ylabel("Action α"); axes[1, col].set_title(f"{action.capitalize()} action"); axes[1, col].set_xticks((0, 1, 2, 3)); axes[1, col].grid(alpha=0.22)
+    axes[1, 0].legend(frameon=False, ncol=2, fontsize=9)
+    fig.suptitle("Kármán Re100: three stable target cycles with phase-aligned controllers", y=0.995)
+    fig.tight_layout(); save(fig, output, "06_example_timeseries_karman_re100")
+    metadata = {"scene": "karman_re100", "source_package": STANDARD, "target_index_start_inclusive": start,
+                "target_index_stop_exclusive": stop, "target_t_D_start": float(data["Target"][start]["t_D"]),
+                "target_t_D_end": float(data["Target"][stop-1]["t_D"]), "samples": stop-start, "displayed_target_cycles": 3,
+                "alignment_signal": signal_name, "alignment_method": "mean-centered normalized circular cross-correlation; smallest maximizing roll selected",
+                "roll_lags_samples": lags, "control_dt_D_over_U0": 0.4, "target_actions_plotted": False}
+    (package_dir / "phase_alignment.json").write_text(json.dumps(metadata, indent=2) + "\n", encoding="utf-8")
+    return metadata
+
+
+def update_manifest(repo_root: Path, package_dir: Path) -> None:
+    old = _load_json(package_dir / "manifest.json"); summary = dict(old.get("summary", {}))
+    summary.pop("diagnostic_figures", None); summary.update({"publication_figures": 5, "publication_tables": 2, "presentation_pages": 2})
+    suffixes = {".csv", ".png", ".pdf", ".md", ".json", ".npz"}; artifacts = []
+    for path in sorted(p for p in package_dir.rglob("*") if p.is_file() and p.suffix in suffixes and p.name != "manifest.json"):
+        artifacts.append({"path": str(path.relative_to(repo_root)), "sha256": hashlib.sha256(path.read_bytes()).hexdigest()})
+    record = {"schema_version": "sr-plotting-package-v3", "source_policy": old["source_policy"], "summary": summary, "artifacts": artifacts}
+    (package_dir / "manifest.json").write_text(json.dumps(record, indent=2) + "\n", encoding="utf-8")
+
+
+def write_package_readme(package_dir: Path, metadata: Mapping[str, Any]) -> None:
+    text = f"""# SR publication plotting package
+
+This package contains five generated publication figures, two replacement tables, and two PPT-ready 16:9 summary pages. All values are derived from frozen canonical artifacts without formula refitting.
+
+- `01_training_case_performance`: panel (a) compares PPO and SR absolute closed-loop legacy DTW similarity at **400 control steps** for all seven training cases; panel (b) reports each controller's 400-minus-200 change. Here 200 and 400 are run durations in control steps, not Reynolds numbers.
+- `02_pointwise_generalization`: 200-step SR results at training conditions and sampled unseen interpolation/extrapolation conditions. Each unseen point is one realization. Categories are unconnected because the evidence does not establish a continuous parameter law or statistical robustness.
+- `03_steady_rotation_calibration`: the retained disturbance-free constant rear-rotation sweep diagnostic.
+- `06_example_timeseries_karman_re100`: late-window phase portraits and PPO/SR actions over approximately three target cycles. Target indices {metadata['target_index_start_inclusive']}–{metadata['target_index_stop_exclusive'] - 1} (`t_D={metadata['target_t_D_start']:.1f}`–`{metadata['target_t_D_end']:.1f}`) are used. PPO and SR segments at the same indices are circularly aligned to target `sensors_center_uy` by mean-centered normalized cross-correlation; deterministic roll lags are PPO={metadata['roll_lags_samples']['PPO']} and SR={metadata['roll_lags_samples']['SR']} samples. No target actions exist or are plotted.
+- `07_flow_field_comparison_karman_re100`: a deterministic stable-cycle snapshot trio. Over [96,146), downstream center-sensor `(ux,uy)` is standardized separately per trajectory; an exhaustive joint search minimizes wrapped phase errors plus a 0.001 rad/sample temporal-separation penalty with a same-direction branch check. This matches center-sensor limit-cycle phase, not six-sensor state equality or exact full-field identity.
+- `tables/offline_action_rmse.*`: offline SR-vs-PPO next-action RMSE on causally aligned PPO-visited states. This is not closed-loop performance.
+- `tables/term_deletion.*`: actual 40-step closed-loop deletion results, canonical formula identities, same-window parent comparisons, and aggregate mean/min summaries.
+
+The `presentation/` directory contains two white-background 16:9 PNG/PDF pages and a source-artifact README. PNG and PDF files share each retained figure stem. `phase_alignment.json` records the exact Figure 06 selection/alignment contract. `manifest.json` hashes every package CSV, Markdown, JSON/NPZ metadata or field artifact, and publication figure.
+"""
+    (package_dir / "README.md").write_text(text, encoding="utf-8")
 
 
 def main(argv: Sequence[str] | None = None) -> int:
-    args = build_parser().parse_args(argv)
-    repo_root = args.repo_root.resolve()
-    run_root = repo_root / "src/SR_analysis/results/runs"
-    output = args.output_dir if args.output_dir.is_absolute() else repo_root / args.output_dir
-    plt.style.use("seaborn-v0_8-whitegrid")
-    plot_performance(run_root, output)
-    plot_generalization(repo_root, output)
-    plot_steady(repo_root, output)
-    plot_offline_residuals(run_root, output)
-    plot_ablation(run_root, output)
-    plot_example_timeseries(run_root, output)
+    parser = argparse.ArgumentParser(description=__doc__); parser.add_argument("--repo-root", type=Path, default=Path(__file__).resolve().parents[3]); parser.add_argument("--output-dir", type=Path, required=True)
+    args = parser.parse_args(argv); repo_root = args.repo_root.resolve(); root = repo_root / "src/SR_analysis/results/runs"
+    output = args.output_dir if args.output_dir.is_absolute() else repo_root / args.output_dir; package_dir = root / PACKAGE
+    plt.style.use("seaborn-v0_8-whitegrid"); plt.rcParams.update({"font.size": 10, "axes.titlesize": 11, "axes.labelsize": 10, "legend.fontsize": 9, "pdf.fonttype": 42})
+    for obsolete in ("04_offline_action_rmse", "05_term_deletion"):
+        for suffix in (".png", ".pdf"): (output / f"{obsolete}{suffix}").unlink(missing_ok=True)
+    plot_performance(root, output); plot_generalization(repo_root, output); plot_steady(repo_root, output)
+    tables = package_dir / "tables"; build_offline_table(root, tables); build_term_deletion_table(repo_root, root, tables)
+    metadata = plot_example_timeseries(root, output, package_dir); write_package_readme(package_dir, metadata); update_manifest(repo_root, package_dir)
     return 0
 
 
-if __name__ == "__main__":
-    raise SystemExit(main())
+if __name__ == "__main__": raise SystemExit(main())

src/SR_analysis/tools/plot_sr_presentation.py

@@ -0,0 +1,76 @@
+#!/usr/bin/env python3
+"""Render two 16:9 presentation pages from frozen canonical SR artifacts."""
+from __future__ import annotations
+import argparse, csv, json, re
+from pathlib import Path
+from typing import Any
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+import numpy as np
+
+BLUE, ORANGE, INK, MUTED, PALE = "#1769AA", "#D65F1E", "#17202A", "#52606D", "#F3F7FA"
+
+def load_json(path: Path) -> dict[str, Any]: return json.loads(path.read_text(encoding="utf-8"))
+def rows(path: Path) -> list[dict[str,str]]:
+    with path.open(newline="", encoding="utf-8") as f: return list(csv.DictReader(f))
+def similarity(path: Path) -> float: return float(rows(path)[-1]["similarity"])
+def expression(root: Path, family: str, head: str) -> str:
+    p=root/f"article-refit-{family}-topology-a-20260718/formulas/{head}.json"
+    return load_json(p)["deployment_expression"]
+def pretty(expr: str) -> str:
+    rounded = re.sub(r"(?<![\w.])[-+]?(?:\d+\.\d*|\.\d+)(?:[eE][+-]?\d+)?", lambda match: f"{float(match.group()):.4f}", expr)
+    return rounded.replace(" * ", r"\,").replace("Cd_rear_a", r"C_{D,r}^{a}").replace("Cl_rear_s", r"C_{L,r}^{s}").replace("Cl_F", r"C_{L,F}")
+def save(fig: plt.Figure, out: Path, stem: str) -> None:
+    out.mkdir(parents=True, exist_ok=True); fig.savefig(out/f"{stem}.png", dpi=200, facecolor="white"); fig.savefig(out/f"{stem}.pdf", facecolor="white"); plt.close(fig)
+def base(title: str, subtitle: str):
+    fig=plt.figure(figsize=(16,9), facecolor="white"); fig.text(.05,.94,title,fontsize=28,weight="bold",color=INK,va="top"); fig.text(.05,.885,subtitle,fontsize=13,color=MUTED,va="top"); return fig
+
+def page_method(root: Path, out: Path) -> None:
+    fig=base("From PPO trajectories to an executable symbolic controller", "Causal alignment, symmetry, all-case refitting, then closed-loop CFD—not offline fit alone")
+    steps=[("1","PPO trajectories","states • forces • actions"),("2","Causal dataset","post-action state → next action"),("3","Topology discovery","G-symmetry-aware symbolic search"),("4","Coefficient refit","fixed topology • all family cases"),("5","CFD validation","40 / 200 / 400 steps • deletion")]
+    for i,(n,h,d) in enumerate(steps):
+        x=.05+i*.19; fig.patches.append(plt.Rectangle((x,.67),.16,.13,transform=fig.transFigure,fc=PALE,ec=BLUE,lw=1.5))
+        fig.text(x+.012,.765,n,fontsize=20,weight="bold",color=BLUE);fig.text(x+.042,.765,h,fontsize=14,weight="bold",color=INK);fig.text(x+.012,.705,d,fontsize=10.5,color=MUTED)
+        if i<4: fig.text(x+.17,.725,"→",fontsize=24,color=ORANGE,ha="center")
+    formulas={f:(pretty(expression(root,f,"joint_front")),pretty(expression(root,f,"joint_rear_shared_upper"))) for f in ("karman","illusion")}
+    for x,f,title in ((.05,"karman","Kármán shared family"),(.52,"illusion","Illusion shared numerical family")):
+        fig.patches.append(plt.Rectangle((x,.25),.43,.34,transform=fig.transFigure,fc="white",ec="#CAD5DF",lw=1.2))
+        fig.text(x+.025,.545,title,fontsize=18,weight="bold",color=BLUE if f=="karman" else ORANGE)
+        front,rear=formulas[f]
+        fig.text(x+.025,.485,r"$\alpha_F(x)=\frac{1}{2}[f(x)-f(Gx)]$",fontsize=14,color=INK)
+        fig.text(x+.025,.435,rf"$f(x)={front}$",fontsize=14,color=INK)
+        fig.text(x+.025,.375,rf"$\alpha_U(x)={rear}$",fontsize=14,color=INK)
+        fig.text(x+.025,.325,r"$\alpha_L(x)=-\alpha_U(Gx)$",fontsize=14,color=INK)
+        note=("Persistent rear counter-rotation constant dominates;\nrear lift correction is secondary." if f=="karman" else "Deployable shared family; deletion evidence does not\nestablish term uniqueness.")
+        fig.text(x+.025,.285,note,fontsize=12,color=MUTED,linespacing=1.35)
+    fig.text(.05,.12,"G mirrors upper/lower quantities and applies the signed symmetry map.  α is dimensionless cylinder surface speed.",fontsize=11.5,color=MUTED)
+    save(fig,out,"01_sr_method_and_formulas")
+
+def page_evidence(root: Path, out: Path) -> None:
+    fig=base("Quantitative evidence at a glance", "Legacy DTW similarity; deterministic frozen artifacts; control-step windows stated explicitly")
+    pkg=root/"article2-plotting-package-20260721"; standard=root/"article2-timeseries-csv-20260720"; longsr=root/"article2-long-timeseries-csv-20260720"
+    scenes={"Kármán":["karman_re50","karman_re100","karman_re200","karman_re400"],"Illusion":["illusion_0.75L","illusion_1L","illusion_1.5L"]}
+    vals={}
+    for fam,ss in scenes.items():
+        sr=np.array([similarity(longsr/s/"sr_dtw_convergence.csv") for s in ss]); ppo=np.array([similarity(pkg/"long_ppo"/s/"ppo_dtw_convergence.csv") for s in ss]); vals[fam]=(sr,ppo)
+    gen=rows(root/"article2-generalization-summary-20260720/generalization.csv")
+    deletion=rows(pkg/"tables/term_deletion.csv")
+    cards=[(.05,.61,.42,.20,"400-step training cases",[(f"Kármán SR range",f"{vals['Kármán'][0].min():.3f}–{vals['Kármán'][0].max():.3f}"),("Kármán mean SR / PPO",f"{vals['Kármán'][0].mean():.3f} / {vals['Kármán'][1].mean():.3f}"),("Illusion SR range",f"{vals['Illusion'][0].min():.3f}–{vals['Illusion'][0].max():.3f}"),("Illusion mean SR / PPO",f"{vals['Illusion'][0].mean():.3f} / {vals['Illusion'][1].mean():.3f}")]),
+    (.52,.61,.43,.20,"200-step unseen conditions",[("Kármán: 4 points",f"{min(float(r['legacy_dtw']) for r in gen if r['objective']=='karman'):.3f}–{max(float(r['legacy_dtw']) for r in gen if r['objective']=='karman'):.3f}"),("Illusion: 5 points",f"{min(float(r['legacy_dtw']) for r in gen if r['objective']=='illusion'):.3f}–{max(float(r['legacy_dtw']) for r in gen if r['objective']=='illusion'):.3f}"),("Evidence level","one realization / condition"),("Scope","sampled interpolation + extrapolation")])]
+    for x,y,w,h,title,items in cards:
+        fig.patches.append(plt.Rectangle((x,y),w,h,transform=fig.transFigure,fc=PALE,ec="#CAD5DF"));fig.text(x+.02,y+h-.045,title,fontsize=17,weight="bold",color=INK)
+        for j,(k,v) in enumerate(items): fig.text(x+.02,y+h-.09-.032*j,k,fontsize=11,color=MUTED);fig.text(x+w-.02,y+h-.09-.032*j,v,fontsize=11.5,weight="bold",color=BLUE,ha="right")
+    variants=[("K front term","k_front0"),("K rear lift","k_rear0"),("K rear constant","k_rear1"),("I worst mean deletion","i_front1")]
+    fig.text(.05,.52,"40-step term deletion: mean parent-relative change",fontsize=18,weight="bold",color=INK)
+    for i,(label,var) in enumerate(variants):
+        sub=[r for r in deletion if r["variant"]==var]; delta=np.mean([float(r["delta_vs_parent_40_step"]) for r in sub]);x=.05+i*.225
+        fig.patches.append(plt.Rectangle((x,.35),.20,.12,transform=fig.transFigure,fc="white",ec="#CAD5DF"));fig.text(x+.015,.43,label,fontsize=11.5,color=MUTED);fig.text(x+.015,.375,f"{delta:+.3f}",fontsize=22,weight="bold",color=ORANGE if delta<-.03 else BLUE)
+    fig.text(.05,.235,"Read with care",fontsize=17,weight="bold",color=INK)
+    cautions=["40-step deletion, 200-step generalization, and 400-step training summaries answer different questions.","Legacy DTW is an aligned trajectory-similarity metric—not causal proof or a physical delay estimate.","Unseen conditions are pointwise deterministic runs; no continuous parameter law or uncertainty band is claimed.","Offline action RMSE is intentionally omitted: imitation on PPO-visited states is not article performance."]
+    for i,t in enumerate(cautions): fig.text(.07,.195-.038*i,"• "+t,fontsize=11.5,color=MUTED)
+    save(fig,out,"02_sr_quantitative_evidence")
+
+def main() -> int:
+    ap=argparse.ArgumentParser();ap.add_argument("--repo-root",type=Path,default=Path(__file__).resolve().parents[3]);ap.add_argument("--output-dir",type=Path,default=None);a=ap.parse_args();root=a.repo_root.resolve()/"src/SR_analysis/results/runs";out=a.output_dir or root/"article2-plotting-package-20260721/presentation";page_method(root,out);page_evidence(root,out);return 0
+if __name__=="__main__": raise SystemExit(main())

src/SR_analysis/tools/prepare_plotting_data.py

@@ -321,8 +321,9 @@ def build_manifest(repo_root: Path, output_dir: Path, summary: Mapping[str, Any]
     artifacts = []
     for path in sorted(
         candidate
-        for pattern in ("*.csv", "*.png", "*.pdf")
+        for pattern in ("*.csv", "*.png", "*.pdf", "*.md", "*.json", "*.npz")
         for candidate in output_dir.rglob(pattern)
+        if candidate.name != "manifest.json"
     ):
         artifacts.append({"path": str(path.relative_to(repo_root)), "sha256": hash_file(path)})
     atomic_write_json(

src/steady_pinball_theory/LITERATURE.md

@@ -0,0 +1,18 @@
+# Source-to-equation ledger
+
+This ledger separates inviscid geometry, finite-Re closure, and viscous stability.
+
+- Crowdy (2006), *Analytical solutions for uniform potential flow past multiple cylinders*: exact multiply connected circular-domain uniform-flow construction. Reuse: independent reference for impermeability, circulation periods and far field; it does not impose no slip or predict separation.
+- Crowdy & Marshall (2007), *Green's functions for Laplace's equation in multiply connected domains*: Schottky–Klein Green-function machinery. Reuse: exact/reference formulation for circular multiply connected domains.
+- Kharlamov & Filip (2012), generalized method of images for several moving parallel cylinders. Reuse: iterative image cross-check and convergence logic.
+- Chan, Jameson & Smits (2011), *Vortex suppression and drag reduction in the wake of counter-rotating cylinders*: viscous doublet-like/reverse-doublet topologies and virtual-body mechanism. Reuse: mechanism observables, not a sign oracle for this code.
+- Mittal (2001), *Control of flow past bluff bodies using rotating control cylinders*: steady-wake suppression near tip-speed ratio five in a different geometry. Reuse: prior for search scale only, never a fixed optimum.
+- Watson (1996), *Slow viscous flow past two rotating cylinders*: matched/Oseen rotating-cylinder precedent. Reuse: conceptual finite-Re interface; regime and geometry differ.
+- Deng et al. (2018), fluidic-pinball bifurcations; Sierra et al. (2020), rotating-cylinder bifurcations. Reuse: require continuation/perturbation and do not assume unique symmetric steady state.
+- Marquet, Sipp & Jacquin (2008): global sensitivity of cylinder flow. Reuse: independent viscous base-flow/eigenvalue workflow after empirical stability.
+
+## Equations and claim limits
+
+The outer field is `u = U_inf e_x + sum_j grad[Q_j log|z-z_j|/(2pi)] + sum_k Gamma_k e_theta/(2pi r)` with source strengths chosen to satisfy cylinder no penetration and zero net source per body. In a strip, scalar sources use equal-sign reflections and vortices opposite-sign reflections at free-slip walls. Truncation, source radius, collocation order, precision, boundary residual, wall residual and an independent exact/image formulation must be reported.
+
+`Gamma_k` is prescribed in the inviscid problem. Mapping wall rotation to an effective circulation is a finite-Re empirical interface requiring contour plateaus and held-out outer profiles; it is not `Gamma=2pi R^2 omega` by assumption. Potential flow cannot predict viscous separation, drag, base pressure, or global stability.

src/steady_pinball_theory/README.md

@@ -0,0 +1,19 @@
+# Steady pinball theory (reset)
+
+Minimal, from-scratch implementation for the uniform-inflow/free-slip steady-cloak study.
+Nothing from the deleted implementation is promoted as evidence.
+
+## Frozen scientific scope
+
+- Canonical geometry: three equal cylinders, rear centers at `(1.3D, +/-0.75D)` from the front; `Re_D=50`; channel half-height `15D`; uniform inlet and free-slip horizontal walls.
+- Primary observer: full cross-section at `x/D=10` from the front-cylinder center.
+- Primary score: `max_y sqrt((u_ctl-u_in)^2+(v_ctl-v_in)^2)/U_inf`; passive `q_blk` is non-ranking.
+- Celeris source uses x-right/y-up lattice coordinates and `(Uw,Vw)=(-omega*ry,omega*rx)`. Both signed rear actions are retained until the paired CFD oracle is reviewed. Names such as clockwise or cloak never override numeric body IDs and actions.
+
+## Files and evidence policy
+
+The package remains below 30 files. Runtime results are written outside this tree. `contract.py` and its tests define source-level facts; `runner.py` produces paired action-reversal diagnostics; `theory.py` is the MFS outer/strip representation; `metrics.py` owns the direct-q_in objective. Literature assumptions and equation provenance are in `LITERATURE.md`.
+
+## Environments
+
+CPU theory/tests: `conda run -n pinball_math ...`. CFD: `PYTHONPATH=CelerisLab/src:src conda run -n pycuda_3_10 python -m steady_pinball_theory.runner ...`. GPU cases execute serially.

src/steady_pinball_theory/RESULTS.md

@@ -0,0 +1,46 @@
+# Reset campaign results and claim matrix
+
+All values below were generated from scratch after deletion of the previous project. Runtime arrays remain outside this source tree under `/tmp/steady-*`; these are diagnostic/production candidates, not a clean tagged archival release.
+
+## Sign oracle
+
+Celeris source defines Cartesian lattice coordinates and `(Uw,Vw)=(-omega*ry,omega*rx)`. With body order `front, rear_y_plus, rear_y_minus`, paired D20 runs at `s=5`, `tU/D=250` give:
+
+- `[0,+Omega,-Omega]`: direct `q_ctl-q_in` `E_inf_vector(x/D=10)=0.13764`; compact steady wake.
+- `[0,-Omega,+Omega]`: `E_inf_vector=0.64593`; broad deficit/wake.
+
+The accepted numerical cloak branch is therefore `[0,+Omega,-Omega]`. Its outer cardinal surfaces move upstream and its gap-facing cardinal surfaces downstream under the solver law; “outer-surface downstream boat-tail” is not an accurate description of this branch.
+
+## Rotation search
+
+At D20, `Re_D=50`, `H/D=15`, uniform inlet and free-slip horizontal walls, the accepted branch was scanned from `s=0` to `7` for `tU/D=250`. The direct profile objective falls from `1.046` at passive `s=0` to `0.0232` near `s=3.45-3.55`, then rises to `0.0604` at `s=4`, `0.1376` at `s=5`, and `0.2972` at `s=7`. The current candidate interval is `s=3.45-3.55`; resolution of the flat pointwise maximum does not justify a unique optimum more precise than this interval.
+
+## Numerical sensitivity
+
+At `s=3.55`, matched q-in references give `E_inf={0.04549,0.02317,0.02348}` for `D={10,20,30}`. D20 and D30 agree within `3.1e-4`; D10 is not converged. Lateral half-height sensitivity at D20 gives `E_inf=0.02382` for `H/D=12`, `0.02317` for `15`, and `0.02478` for `18`. Fine-grid D30 local values are `0.02353,0.02348,0.02921` at `s=3.45,3.55,3.65`.
+
+Late-window changes in the scanned steady candidates are `O(1e-5-1e-4)` in the profile diagnostic. Neighboring cold-start cases at D30 (`s=3.4,3.5,3.6,3.7`) likewise settle over `tU/D=250`. This is empirical cold-start steadiness only; no checkpoint perturbation decay or global eigenanalysis has been completed.
+
+## Analytical representation and finite-Re interface
+
+The MFS solver enforces cylinder no penetration and zero source per body for unbounded and image-strip formulations. Tests cover collocation order, source radius and image-layer sensitivity. A strip circulation scan predicts near cancellation at `|Gamma|/(U D) about 4.25` (sign depends on body/circulation convention), with boundary residual below `1e-14` in that solve. Outer-mask fits to D30 CFD at `s=3.45,3.55,3.65` yield signed `Gamma_eff/(UD)=-4.293,-4.339,-4.384`, close in magnitude to the inviscid cancellation value; residual component RMS is about `1.9-2.1% U`. This supports a circulation-dominated outer mechanism but does not close boundary layers, pressure, force or separation.
+
+## Claim matrix
+
+Supported now:
+- exact numeric action/body mapping and action-reversal ordering;
+- direct q-in profile metric at the historical `x/D=10` plane;
+- a D20/D30-converged low-error candidate interval near `s=3.45-3.55` in the tested Celeris setup;
+- numerical MFS no-penetration/image-strip representation and effective-circulation magnitude agreement.
+
+Bounded/descriptive:
+- empirical cold-start steadiness through `tU/D=250`;
+- circulation-dominated outer-flow interpretation;
+- lateral-domain sensitivity over `H/D=12-18`.
+
+Not supported:
+- global stability, unique attractor, perturbation decay, or exact optimum;
+- independent Navier-Stokes reproduction (no independent NS package is installed);
+- pressure/drag/base-bleed/separation closure;
+- JFM-ready exact Schottky-Klein derivation or clean immutable production release;
+- three-dimensional, experimental, energetic or other-Re generality.

src/steady_pinball_theory/__init__.py

@@ -0,0 +1,5 @@
+"""Minimal steady fluidic-pinball theory package."""
+from .contract import CanonicalContract, wall_velocity
+from .metrics import profile_errors
+
+__all__ = ["CanonicalContract", "wall_velocity", "profile_errors"]

src/steady_pinball_theory/contract.py

@@ -0,0 +1,70 @@
+"""Frozen geometry and sign primitives derived from Celeris source.
+
+No branch is called a cloak until the action-reversal CFD oracle is reviewed.
+Coordinates are lattice/Cartesian: x right, y up.
+"""
+from __future__ import annotations
+from dataclasses import dataclass
+import numpy as np
+
+
+def wall_velocity(omega: float, rx: float, ry: float) -> tuple[float, float]:
+    """Celeris curved-wall law: (Uw,Vw)=(-omega*ry, omega*rx)."""
+    return (-float(omega) * float(ry), float(omega) * float(rx))
+
+
+@dataclass(frozen=True)
+class CanonicalContract:
+    diameter: float = 20.0
+    u_inf: float = 0.01
+    re_d: float = 50.0
+    half_height_d: float = 15.0
+    upstream_d: float = 50.0
+    downstream_d: float = 50.0
+    rear_dx_d: float = 1.3
+    rear_dy_d: float = 0.75
+    primary_station_d: float = 10.0
+
+    @property
+    def radius(self) -> float: return self.diameter / 2.0
+    @property
+    def viscosity(self) -> float: return self.u_inf * self.diameter / self.re_d
+    @property
+    def nx(self) -> int: return int(round((self.upstream_d + self.downstream_d) * self.diameter))
+    @property
+    def ny(self) -> int: return int(round(2.0 * self.half_height_d * self.diameter))
+    @property
+    def center_y(self) -> float: return (self.ny - 1.0) / 2.0
+    @property
+    def centers(self) -> tuple[tuple[float,float], ...]:
+        x0 = self.upstream_d * self.diameter
+        return ((x0, self.center_y),
+                (x0 + self.rear_dx_d*self.diameter, self.center_y + self.rear_dy_d*self.diameter),
+                (x0 + self.rear_dx_d*self.diameter, self.center_y - self.rear_dy_d*self.diameter))
+    @property
+    def body_order(self) -> tuple[str,...]: return ("front", "rear_y_plus", "rear_y_minus")
+
+    def omega_from_s(self, s: float) -> float:
+        return float(s) * self.u_inf / self.radius
+
+    def action(self, s: float, branch: str) -> np.ndarray:
+        o = self.omega_from_s(s)
+        if branch == "plus-minus": return np.array([0.0, +o, -o])
+        if branch == "minus-plus": return np.array([0.0, -o, +o])
+        raise ValueError("branch must be 'plus-minus' or 'minus-plus'")
+
+    def oracle(self, s: float = 1.0) -> dict:
+        """Explicit outer/gap cardinal-point velocities for both branches."""
+        out = {"schema": "steady-pinball-sign-oracle/v1", "coordinates": "x-right_y-up"}
+        for branch in ("plus-minus", "minus-plus"):
+            a = self.action(s, branch)
+            # y+ rear: outer is top (ry=+R), gap-facing is bottom (ry=-R).
+            # y- rear: outer is bottom (ry=-R), gap-facing is top (ry=+R).
+            out[branch] = {
+                "action": a.tolist(),
+                "rear_y_plus_outer_top": wall_velocity(a[1], 0.0, +self.radius),
+                "rear_y_plus_gap_bottom": wall_velocity(a[1], 0.0, -self.radius),
+                "rear_y_minus_outer_bottom": wall_velocity(a[2], 0.0, -self.radius),
+                "rear_y_minus_gap_top": wall_velocity(a[2], 0.0, +self.radius),
+            }
+        return out

src/steady_pinball_theory/metrics.py

@@ -0,0 +1,21 @@
+"""Direct q_ctl-to-uniform-q_in profile objective and steady diagnostics."""
+from __future__ import annotations
+import numpy as np
+
+
+def profile_errors(ux: np.ndarray, uy: np.ndarray, *, u_inf: float) -> dict[str, float]:
+    ux=np.asarray(ux,dtype=float); uy=np.asarray(uy,dtype=float)
+    if ux.shape != uy.shape or ux.ndim != 1: raise ValueError("profiles must be aligned 1-D arrays")
+    if not np.isfinite(ux).all() or not np.isfinite(uy).all() or u_inf <= 0: raise ValueError("invalid profile")
+    du=ux/u_inf-1.0; dv=uy/u_inf
+    vec=np.hypot(du,dv)
+    return {"E_inf_x":float(np.max(np.abs(du))), "E_inf_y":float(np.max(np.abs(dv))),
+            "E_inf_vector":float(np.max(vec)), "E_L2_vector":float(np.sqrt(np.mean(vec*vec)))}
+
+
+def two_window_change(samples: np.ndarray) -> float:
+    a=np.asarray(samples,dtype=float)
+    if a.ndim < 1 or a.shape[0] < 4: return float("inf")
+    n=a.shape[0]//2
+    x=np.mean(a[:n],axis=0); y=np.mean(a[-n:],axis=0)
+    return float(np.sqrt(np.mean((y-x)**2)))

src/steady_pinball_theory/runner.py

@@ -0,0 +1,67 @@
+"""Serial Celeris diagnostic/scan runner; outputs live outside this package."""
+from __future__ import annotations
+import argparse, json, os, tempfile
+from pathlib import Path
+import numpy as np
+from .contract import CanonicalContract
+from .metrics import profile_errors
+
+
+def _cfg(c: CanonicalContract, D: float):
+    scale=D/c.diameter; nx=int(round(c.nx*scale)); ny=int(round(c.ny*scale))
+    return {"grid":{"lattice_model":"D2Q9","nx":nx,"ny":ny,"nz":1},
+      "physics":{"data_type":"FP32","viscosity":c.u_inf*D/c.re_d,"velocity":c.u_inf,"rho":1.0},
+      "method":{"collision":"MRT","streaming":"double_buffer","store_precision":"FP32","ddf_shifting":False,
+        "les":{"enabled":False,"cs":0.16,"closed_form":True},"trt":{"magic_param":0.1875},
+        "inlet":{"profile":"uniform","scheme":"regularized","trt_neq_damp":0.5,"regularized_neq_damp":0.5},
+        "outlet":{"mode":"neq_extrap","backflow_clamp":True,"blend_alpha":0.7,"srt_neq_damp":0.5},
+        "y_wall_bc":"free_slip","omega_guard":{"min":0.01,"max":1.99}},
+      "cuda":{"threads_per_block":256,"compute_capability":"auto"}}
+
+
+def run_case(out: Path, *, branch: str, s: float, D: float, t_end: float, sample_dt: float, half_height_d: float = 15.0, upstream_d: float = 50.0, downstream_d: float = 50.0):
+    from CelerisLab import Simulation
+    c=CanonicalContract(diameter=D, half_height_d=half_height_d, upstream_d=upstream_d, downstream_d=downstream_d); cfg=_cfg(c,D); centers=c.centers
+    action=np.zeros(3) if branch=="qin" else c.action(s,branch)
+    body={"objects":[]} if branch=="qin" else {"objects":[{"type":"cylinder","center":list(z),"radius":c.radius,"omega":float(w)} for z,w in zip(centers,action)]}
+    out.mkdir(parents=True,exist_ok=False)
+    (out/'contract.json').write_text(json.dumps({"branch":branch,"s":s,"D":D,"t_end":t_end,"sample_dt":sample_dt,"action":action.tolist(),"centers":centers,"oracle":c.oracle(s)},indent=2))
+    with tempfile.TemporaryDirectory(prefix='steady_cfd_') as td:
+        lp=Path(td)/'config_lbm.json'; bp=Path(td)/'config_body.json'
+        lp.write_text(json.dumps(cfg)); bp.write_text(json.dumps(body))
+        sim=Simulation(lbm_config_path=str(lp),body_config_path=str(bp)); sim.initialize()
+        steps=int(round(t_end*D/c.u_inf)); stride=max(1,int(round(sample_dt*D/c.u_inf)))
+        hist=[]
+        for stop in range(stride,steps+1,stride):
+            sim.run(min(stride,steps-(stop-stride)),zero_obs=True)
+            m=sim.get_macroscopic(); ux=np.asarray(m['ux']); uy=np.asarray(m['uy'])
+            ix=int(round(centers[0][0]+c.primary_station_d*D))
+            prof=profile_errors(ux[2:-2,ix],uy[2:-2,ix],u_inf=c.u_inf)
+            hist.append([stop,prof['E_inf_vector'],prof['E_L2_vector']])
+        m=sim.get_macroscopic(); sim.close()
+    ux=np.asarray(m['ux'],np.float32); uy=np.asarray(m['uy'],np.float32); rho=np.asarray(m['rho'],np.float32)
+    np.savez_compressed(out/'endpoint.npz',ux=ux,uy=uy,rho=rho,history=np.asarray(hist),action=action)
+    from CelerisLab.common.render import compute_vorticity,render_vorticity_field
+    render_vorticity_field(compute_vorticity(ux,uy),nx=cfg['grid']['nx'],ny=cfg['grid']['ny'],out_path=str(out/'vorticity.png'),cylinders=[] if branch=='qin' else [(z,c.radius) for z in centers])
+    result={"branch":branch,"action":action.tolist(),"steps":steps,"final_uniform_diagnostic":dict(zip(["step","E_inf_vector","E_L2_vector"],hist[-1]))}
+    (out/'result.json').write_text(json.dumps(result,indent=2)); return result
+
+
+def main(argv=None):
+    p=argparse.ArgumentParser(); p.add_argument('command',choices=['oracle','diagnostic']); p.add_argument('--out',type=Path)
+    p.add_argument('--D',type=float,default=10.0); p.add_argument('--s',type=float,default=5.0); p.add_argument('--t-end',type=float,default=40.0); p.add_argument('--sample-dt',type=float,default=2.0)
+    a=p.parse_args(argv); c=CanonicalContract(diameter=a.D)
+    if a.command=='oracle': print(json.dumps(c.oracle(a.s),indent=2)); return 0
+    if a.out is None: p.error('--out required')
+    a.out.mkdir(parents=True,exist_ok=False)
+    rows=[]
+    for b in ('qin','plus-minus','minus-plus'):
+        rows.append(run_case(a.out/b,branch=b,s=a.s,D=a.D,t_end=a.t_end,sample_dt=a.sample_dt))
+    
+    q=np.load(a.out/'qin'/'endpoint.npz'); ix=int(round(c.centers[0][0]+c.primary_station_d*a.D)); refux=q['ux'][2:-2,ix]; refuy=q['uy'][2:-2,ix]
+    for row in rows:
+        e=np.load(a.out/row['branch']/'endpoint.npz'); du=(e['ux'][2:-2,ix]-refux)/c.u_inf; dv=(e['uy'][2:-2,ix]-refuy)/c.u_inf; vec=np.hypot(du,dv)
+        row['direct_qin']={"E_inf_x":float(np.max(np.abs(du))),"E_inf_y":float(np.max(np.abs(dv))),"E_inf_vector":float(np.max(vec)),"E_L2_vector":float(np.sqrt(np.mean(vec*vec)))}
+    (a.out/'summary.json').write_text(json.dumps(rows,indent=2)); return 0
+
+if __name__=='__main__': raise SystemExit(main())

src/steady_pinball_theory/tests/test_core.py

@@ -0,0 +1,52 @@
+import numpy as np
+from steady_pinball_theory.contract import CanonicalContract,wall_velocity
+from steady_pinball_theory.metrics import profile_errors
+from steady_pinball_theory.theory import solve_mfs
+
+
+def test_canonical_geometry_and_reynolds():
+ c=CanonicalContract(); assert c.centers==((1000.0,299.5),(1026.0,314.5),(1026.0,284.5)); assert c.viscosity==0.004
+
+def test_source_level_wall_oracle():
+ c=CanonicalContract(); o=c.oracle(5.0)
+ assert o['plus-minus']['action']==[0.0,0.005,-0.005]
+ assert o['plus-minus']['rear_y_plus_outer_top'][0] < 0
+ assert o['plus-minus']['rear_y_plus_gap_bottom'][0] > 0
+ assert o['minus-plus']['rear_y_plus_outer_top'][0] > 0
+ assert wall_velocity(2,3,4)==(-8.0,6.0)
+
+def test_primary_metric_is_pointwise_vector_norm():
+ r=profile_errors(np.array([1.,.8]),np.array([0.,.1]),u_inf=1.)
+ assert np.isclose(r['E_inf_vector'],np.sqrt(.05))
+
+def test_unbounded_mfs_impermeability_and_zero_source():
+ centers=np.array([[0.,0.],[1.3,.75],[1.3,-.75]])
+ sol,diag=solve_mfs(centers,.5,[0.,2.,-2.],n_boundary=48)
+ assert diag['boundary_Linf'] < 2e-5
+ for q in sol.strengths.reshape(3,-1): assert abs(q.sum()) < 1e-10
+
+def test_strip_images_reduce_wall_normal_velocity():
+ centers=np.array([[0.,0.],[1.3,.75],[1.3,-.75]])
+ sol,diag=solve_mfs(centers,.5,[0.,2.,-2.],n_boundary=48,H=15.,image_layers=20)
+ x=np.linspace(-5,8,300); p=np.r_[np.c_[x,np.full_like(x,15.)],np.c_[x,np.full_like(x,-15.)]]
+ assert np.max(np.abs(sol.velocity(p)[:,1])) < 2e-3
+ assert diag['boundary_Linf'] < 3e-5
+
+
+def test_mfs_resolution_and_source_radius_sensitivity():
+ centers=np.array([[0.,0.],[1.3,.75],[1.3,-.75]])
+ rows=[]
+ for n in (32,48,64):
+  for ratio in (.55,.65,.75):
+   _,d=solve_mfs(centers,.5,[0.,2.,-2.],n_boundary=n,source_ratio=ratio)
+   rows.append(d['boundary_Linf'])
+ assert max(rows) < 2e-4
+
+def test_strip_image_layer_sensitivity():
+ centers=np.array([[0.,0.],[1.3,.75],[1.3,-.75]])
+ vals=[]
+ for layers in (10,20,40):
+  sol,d=solve_mfs(centers,.5,[0.,2.,-2.],n_boundary=48,H=15.,image_layers=layers)
+  vals.append(sol.velocity(np.array([[10.,0.]]))[0])
+  assert d['boundary_Linf'] < 4e-5
+ assert np.max(np.linalg.norm(np.diff(vals,axis=0),axis=1)) < 2e-4

src/steady_pinball_theory/theory.py

@@ -0,0 +1,66 @@
+"""MFS outer/strip potential-flow representation with prescribed circulation."""
+from __future__ import annotations
+from dataclasses import dataclass
+import numpy as np
+
+
+def _images(y0: float, H: float, layers: int, vortex: bool):
+    for k in range(-layers,layers+1):
+        yield y0+4*k*H, 1.0
+        yield 2*H-y0+4*k*H, -1.0 if vortex else 1.0
+
+
+def singular_velocity(points, source, *, H=None, layers=0, vortex=False):
+    p=np.asarray(points,float); x0,y0=map(float,source); vel=np.zeros_like(p)
+    imgs=[(y0,1.0)] if H is None else list(_images(y0,float(H),int(layers),vortex))
+    for yi,sgn in imgs:
+        dx=p[:,0]-x0; dy=p[:,1]-yi; r2=dx*dx+dy*dy
+        if np.any(r2 <= 1e-24): raise ValueError("evaluation at singularity")
+        if vortex: vel += sgn*np.column_stack((-dy/r2,dx/r2))/(2*np.pi)
+        else: vel += sgn*np.column_stack((dx/r2,dy/r2))/(2*np.pi)
+    return vel
+
+
+@dataclass
+class MFSSolution:
+    centers: np.ndarray; radius: float; strengths: np.ndarray; source_points: np.ndarray
+    circulations: np.ndarray; u_inf: float; H: float|None; image_layers: int
+    def velocity(self, points):
+        p=np.asarray(points,float); v=np.zeros_like(p); v[:,0]=self.u_inf
+        for q,z in zip(self.strengths,self.source_points):
+            v += q*singular_velocity(p,z,H=self.H,layers=self.image_layers)
+        for g,z in zip(self.circulations,self.centers):
+            v += g*singular_velocity(p,z,H=self.H,layers=self.image_layers,vortex=True)
+        return v
+
+
+def solve_mfs(centers, radius, circulations, *, u_inf=1.0, n_boundary=96,
+              source_ratio=0.65, H=None, image_layers=20):
+    c=np.asarray(centers,float); gam=np.asarray(circulations,float)
+    if c.shape != (len(gam),2): raise ValueError("centers/circulations mismatch")
+    th=(np.arange(n_boundary)+0.5)*2*np.pi/n_boundary
+    normals=np.column_stack((np.cos(th),np.sin(th)))
+    bpts=np.concatenate([ci+radius*normals for ci in c])
+    nrms=np.tile(normals,(len(c),1))
+    spts=np.concatenate([ci+source_ratio*radius*normals for ci in c])
+    base=np.zeros_like(bpts); base[:,0]=u_inf
+    for g,ci in zip(gam,c): base += g*singular_velocity(bpts,ci,H=H,layers=image_layers,vortex=True)
+    rhs=-np.einsum('ij,ij->i',base,nrms)
+    # Eliminate one source per body so every body's net source is exactly zero.
+    cols=[]
+    for ib in range(len(c)):
+        ref=spts[ib*n_boundary+n_boundary-1]
+        vr=singular_velocity(bpts,ref,H=H,layers=image_layers)
+        for j in range(n_boundary-1):
+            v=singular_velocity(bpts,spts[ib*n_boundary+j],H=H,layers=image_layers)-vr
+            cols.append(np.einsum('ij,ij->i',v,nrms))
+    A=np.column_stack(cols); coeff,resid,rank,sv=np.linalg.lstsq(A,rhs,rcond=1e-12)
+    strengths=np.zeros(len(spts)); k=0
+    for ib in range(len(c)):
+        q=coeff[k:k+n_boundary-1]; k+=n_boundary-1
+        strengths[ib*n_boundary:ib*n_boundary+n_boundary-1]=q
+        strengths[ib*n_boundary+n_boundary-1]=-np.sum(q)
+    sol=MFSSolution(c,float(radius),strengths,spts,gam,float(u_inf),H,int(image_layers))
+    residual=np.einsum('ij,ij->i',sol.velocity(bpts),nrms)
+    return sol,{"boundary_Linf":float(np.max(np.abs(residual))),"boundary_L2":float(np.sqrt(np.mean(residual**2))),
+                "rank":int(rank),"unknowns":int(A.shape[1]),"condition":float(sv[0]/sv[-1])}