refactor(eval): consolidate drl-pinball reproduction

Retire duplicate reproduction paths in favor of the canonical V5 and Legacy runners, while preserving historical tooling in archives and publishing audited summary plots.

Co-authored-by: Cursor <cursoragent@cursor.com>
This commit is contained in:
Frank14f
2026-08-08 16:19:11 +08:00
co-authored by Cursor
parent 61e82ec90a
commit 326d642751
79 changed files with 3363 additions and 3996 deletions
Binary file not shown.

After

Width:  |  Height:  |  Size: 565 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 403 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 643 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 497 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 598 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 677 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 754 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 462 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 240 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 462 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 415 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 400 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 454 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 442 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 414 KiB

File diff suppressed because it is too large Load Diff
Binary file not shown.

After

Width:  |  Height:  |  Size: 542 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 650 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 626 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 630 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 472 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 380 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 387 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 440 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 441 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 440 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 441 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 442 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 558 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 561 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 302 KiB

+39 -199
View File
@@ -1,206 +1,46 @@
# Eval Benchmark — 推理评估与流场分析
# Canonical V5 evaluation
> 对 V5 Train 模型和 Legacy 旧模型进行全面推理采样和可视化对比。
> 复用 train env 保证与训练流程完全一致,不使用 skeleton injection。
`infer_train.py` evaluates retained scratch policies using suffix-free case IDs from
`drl_pinball.case_registry`; there is no independent eval manifest.
## 文件结构
The evaluator requires each selected run to contain the exact bundle:
`models/best_model.zip`, `vec_normalize.pkl`, `calibration.json`, and
`target.npy`. For Illusion, harmonics come from the registry-matched calibration
directory only after its `target.npy` SHA256 matches the run-local target exactly.
For Illusion's proven legacy-run/native-registry schema pair, target and harmonic sensor
channels 0-5 are multiplied in memory by `SENSOR_CC`; force harmonics 6-7 and all
frequencies/phases are preserved. Ambiguous schema combinations fail rather than guess.
The run-local `vec_normalize.pkl` is the historical final compatibility alias used
for standalone-eval reproduction in every case except `kar_d075`, whose matrix A/B
evidence requires run-local `best_vecnormalize.pkl`. This case-ID selection is explicit,
recorded in provenance, and fails closed; artifacts are read without modification.
Historical models run through `legacy-policy-v1`; `VecNormalize` is frozen with `training=False` and
`norm_reward=False`. Rollouts are deterministic and fixed at 360 steps with a
180-step scoring tail.
```
src/drl_pinball/eval/
README.md # 本文件 (自我审查清单 + 复习流程)
__init__.py
scene_manifest.py # 所有 scene 的配置单点真理源
infer_train.py # V5 Train 模型推理引擎 (GPU 0)
infer_reproduce.py # Reproduce 模型推理引擎 (GPU 1)
viz_flow.py # 涡量图跨场景拼图 (统一 colormap)
viz_signals.py # obs/action/DTW/FFT 时序可视化
generate_report.py # 汇总表格 + 柱状图
run_all.sh # 主启动器 (GPU 调度)
output/
train/{scene}/ # signals.npz, vorticity_*.png, metrics.json, all_seeds.json
reproduce/{scene}/ # 同上
reports/
comparison_table.md # 汇总对比表
summary_dtw_barchart.png # DTW 柱状图
vorticity_panels/ # 跨场景拼图
signal_plots/ # 诊断信号图
```
## Quick Start
By default only the first registered case and first seed run:
```bash
conda activate pycuda_3_10
cd src/drl_pinball/eval
# === 全量跑 (双 GPU) ===
bash run_all.sh
# === 单场景测试 ===
python infer_train.py --scene re100 --device-id 0
python infer_train.py --scene illusion_1L --device-id 0
# === 仅跑报告生成 (已有 output 数据) ===
python generate_report.py
python viz_signals.py --side train
python viz_flow.py
bash src/drl_pinball/eval/run_all.sh
bash src/drl_pinball/eval/run_all.sh --case kar_re100 --seed 45
bash src/drl_pinball/eval/run_all.sh --case kar_re100 --all-seeds --overwrite
bash src/drl_pinball/eval/run_all.sh --case kar_re100 --all-seeds --metrics-only \
--output-root /tmp/v5-eval --overwrite
python3 src/drl_pinball/eval/infer_train.py --validate-all --output-root /tmp/v5-eval
```
## GPU 调度策略
```
时间轴:
GPU0: [re100] [ill_1L] [ill_15L] [ill_075L] --120s--> [re60] --120s--> [re200] --120s--> [re400]
(config A: karman_2000x600, nu=0.004) (config B: re60) (config C: re200) (config D: re400)
GPU1: --120s--> [karman_cloak] [steady_cloak] [ill_075L] [ill_1L] [ill_15L] [vortex_lamb] [vortex_taylor]
(全部 config_pinball.json, 无额外编译)
```
- GPU 1 比 GPU 0 晚 120 秒启动,避免两个不同 config 同时触发编译
- 同 config 内 scene 连续运行,无 recompilation 开销
- 每个 scene 约 5-8 分钟(warmup + FIFO + 360 步推理)
## 场景覆盖率
### Train pipeline (GPU 0)
| Scene | Config | Seeds |
|-------|--------|-------|
| re100_karman | karman_2000x600 (Re=100) | 41-45, 539439 |
| transfer_re60 | karman_2000x600_re60 | 41, 43 |
| transfer_re200 | karman_2000x600_re200 | 41, 43, 45 |
| transfer_re400 | karman_2000x600_re400 | 43 |
| illusion_075L | karman_2000x600 | 41 |
| illusion_1L | karman_2000x600 | 41, 43 |
| illusion_15L | karman_2000x600 | 41, 43 |
### Reproduce pipeline (GPU 1)
| Scene | Config | 模型 |
|-------|--------|------|
| karman_cloak | pinball (1280x512) | d1a3o12_re100 |
| steady_cloak | pinball (1280x512) | open-loop [0,-5.1,5.1] |
| illusion_075L | pinball (1280x512) | d1a3o14_250525_imit_075L_2U_400S |
| illusion_1L | pinball (1280x512) | d1a3o14_250525_imit_1L_2U_600S |
| illusion_15L | pinball (1280x512) | d1a3o14_250525_imit_15L_2U |
| vortex_lamb | pinball (1280x512) | vortex_lamb |
| vortex_taylor | pinball (1280x512) | vortex_taylor |
## 核心设计原则
1. **复用 train env 而非重构 CFD** — KarmanCloakEnv / IllusionCloakEnv 已封装 norm、calibration、DTW、snapshot/restore
2. **PPO.load() + VecNormalize.load()** — SB3 标准 API,不使用 skeleton injection
3. **统一涡量 [-0.03, 0.03] colormap** — 所有场景跨可比
4. **每步全量时序** — sensors(6) + forces(6) + actions(3) + rewards 存入 signals.npz
---
## 自我审查清单
### 启动前
- [ ] `conda activate pycuda_3_10`
- [ ] 所有 `best_model.zip``vec_normalize.pkl` 路径存在 (见 scene_manifest.py)
- [ ] 所有 `calibration.json``target.npy` 存在 (train/calibrations/)
- [ ] `nvidia-smi` 显示两个 GPU 空闲
- [ ] 无僵尸 pycuda 进程 (`pkill -f pycuda` 后等待 3 分钟)
### Phase 1 — Train 推理
- [ ] `python infer_train.py --scene re100 --device-id 0` 单场景测试通过
- [ ] 输出的 `signals.npz` 形状: (360,6) sensors, (360,6) forces, (360,3) actions
- [ ] rewards 非 NaNDTW sim 在 [0,1] 内
- [ ] Zero-action baseline 的 DTW 明显低于 controlled
### Phase 2 — Reproduce 推理
- [ ] `python infer_reproduce.py --scene karman_cloak --device-id 1` 测试通过
- [ ] 旧模型 PPO.load() 不报错
- [ ] Vortex 场景 add_vortex 正确初始化(check target.npz 非全零)
### Phase 3 — 涡量图
- [ ] `python viz_flow.py` 产出全部拼图
- [ ] vorticity range 统一 [-0.03, 0.03]
- [ ] 圆柱边界清晰可见
- [ ] target vs controlled vs zero 三栏对比正常
### Phase 4 — 时序图
- [ ] `python viz_signals.py --side train` 产出全部诊断图
- [ ] FFT 频谱确认 controlled 与 target Strouhal 匹配
- [ ] Action 时序 policy 收敛到稳态 (非振荡)
### Phase 5 — 汇总报告
- [ ] `python generate_report.py` 产出 comparison_table.md
- [ ] Train DTW >= Reproduce DTW (或合理的 gap 原因)
### 最终验收
- [ ] 7 train + 7 reproduce 目录完整
- [ ] 每个目录: signals.npz, vorticity_controlled.png, vorticity_target.png, vorticity_zero.png, metrics.json, all_seeds.json
- [ ] reports/ 下: comparison_table.md, summary_dtw_barchart.png, vorticity_panels/, signal_plots/
---
## 复习流程 (中断后快速恢复)
按顺序读取以下文件重建上下文:
1. **本文件** — 确认当前进度和 GPU 调度图
2. **`eval/scene_manifest.py`** — 所有 scene 配置
```bash
python -c "from drl_pinball.eval.scene_manifest import TRAIN_SCENES; \
[print(s['scene_id'], s['si']) for s in TRAIN_SCENES]"
```
3. **`train/TRAIN_KNOWLEDGE.md`** sections 1-3 — V5 训练管线回顾
4. **验证环境可用**:
```bash
ls eval/output/train/*/metrics.json 2>/dev/null # 已完成的 scene
ls eval/output/reproduce/*/metrics.json 2>/dev/null
```
5. **从断点继续**:
```bash
# 单 scene:
python infer_train.py --scene re200 --device-id 0
# 或 Kill pycuda 进程后重新启动显卡:
pkill -f pycuda; sleep 180; bash run_all.sh
```
### 常见问题排查
| 问题 | 解决方案 |
|------|---------|
| `cuInit failed` | GPU 被占用。`pkill -f pycuda; sleep 180` 后重试 |
| `best_model.zip not found` | 训练未完成。查看 `train/output/{case}/models/` |
| `target.npy not found` | calibration 缺失。重新运行 `calibrate.py` |
| DTW sim = -1.0 (reproduce) | 正常。Reproduce 不计算 DTW,由 report 统一处理 |
| 涡量图空白 | 检查 `render_vorticity_field` 的 `vmin/vmax` 或 cylinders 坐标 |
| VecNormalize.pkl 不存在 | 训练时未保存。用 `best_model.zip` 同目录的 `vec_normalize.pkl` |
---
## 代码关键点
### infer_train.py 的关键差异点
- 每个 seed 独立 `PPO.load()` + `VecNormalize.load()`,不是共用 skeleton
- 推理时 `SymmetryAugmentWrapper(prob=0.0)`deterministic mode
- Raw obs 从 `env._read_obs()` 提取,通过 wrapper chain 穿透
- Target vorticity 单独创建 Simulation 生成(dist cylinder only
### infer_reproduce.py 的关键差异点
- 1280x512 网格,parabolic inletbounce-back walls
- Legacy norm 从 `SR_analysis/data/*/norm.json` 加载
- 旧动作缩放: `norm * scale + bias` (各 scene 不同)
- Illusion 模型 obs_dim=14(含 target_cd, target_cl,通过 harmonics 重构)
- Vortex 用 `add_vortex()` 在不同 x 位置添加涡量
### 输出文件命名
| 文件 | 内容 |
|------|------|
| `signals.npz` | sensors (N,6), forces (N,6), actions (N,3), rewards (N,) |
| `vorticity_controlled.png` | DRL 控制后的 final 时刻涡量 |
| `vorticity_target.png` | 目标状态 (disturbance only / target cylinder) |
| `vorticity_zero.png` | 零动作 baseline (无控制) |
| `metrics.json` | DTW sim, reward, action stats |
| `all_seeds.json` | 每个 seed 的详细评分 |
Outputs are confined to `output/train/<case>/`: `metrics.json`, `all_seeds.json`,
compact signal NPZ files, and controlled/target/zero final-vorticity PNGs rendered
with CelerisLab at raw `[-0.001, 0.001]`. Existing suffix-free output is a retained
baseline and is never replaced without `--overwrite`. `--validate` and
`--validate-all` first require the retained CSV and JSON tables to agree at strict
absolute tolerance `1e-6`, including the selected seed. Fresh output must use that
selected seed, then passes the explicit GPU reproduction gate when reward, reward
components, and DTW are within `0.02` absolute and action means are within `0.03`.
Validation prints every fresh-minus-reference delta on both pass and failure because
running VecNormalize state and GPU replay do not imply bitwise reproduction. Action
means use the same tail-180 window declared by the retained tables. `--metrics-only`
runs controlled evaluation and writes only `metrics.json` and compact `all_seeds.json`;
it skips NPZ signals and all controlled replay/target/zero vorticity environments.
Use `--output-root` for an explicit disposable full-matrix result tree.
Historical pretask helpers are retained under `archive/` and are not active entrypoints.
@@ -0,0 +1,234 @@
# V5 Eval Pipeline — Inference & Results Documentation
> **Scope**: Inference infrastructure for V5-trained PPO models on the new CelerisLab solver.
> Legacy model reproduction is handled by @src/drl_pinball/reproduce/.
> Companion to `train/TRAIN_PIPELINE.md` for the training infrastructure.
---
## 1. Architecture Overview
```
scene_manifest.py (Single Source of Truth)
└── TRAIN_SCENES[20 entries] → used by infer_train.py
infer_train.py (GPU 2)
├── Create V5 CFD env
├── Per-seed:
│ ├── VecNormalize.load(seed.pkl)
│ ├── Skeleton PPO + weight inject
│ └── 360-step deterministic roll
├── Output: signals.npz + fields.npz + PNGs
└── Best-seed metrics.json
run_all_with_fields.sh
├── Phase 1: infer_train.py (stagger 120s per config change)
├── Phase 2: collect_baselines.py
├── Phase 3: bridge_to_sr.py
└── Phase 4: data integrity check
```
---
## 2. Scene Manifest (`scene_manifest.py`)
The single source of truth for all scene configurations used by both pipelines.
### 2.1 TRAIN_SCENES (V5 PPO, 2000×600, uniform + free-slip)
| scene_id | Type | SI | num_steps | Seeds | Description |
|----------|------|:--:|:---------:|:-----:|-------------|
| `kar_re100_sc` | karman | 800 | 360 | 41-45 | Karman Re100, 5-seed study |
| `kar_re60_tr` | karman | 800 | 360 | 43 | Cross-Re transfer to Re60 |
| `kar_re200_tr` | karman | 500 | 360 | 43 | Cross-Re transfer to Re200 |
| `kar_re400_tr` | karman | 400 | 360 | 43 | Cross-Re transfer to Re400 |
| `kar_d075_sc` | karman | 800 | 360 | 44 | VarDist d=0.75L scratch |
| `kar_d15_sc` | karman | 800 | 360 | 45 | VarDist d=1.5L scratch |
| `kar_d2_sc` | karman | 800 | 360 | 45 | VarDist d=2.0L scratch |
| `ill_075L_sc` | illusion | 400 | 360 | 43 | Illusion 0.75L target |
| `ill_1L_sc` | illusion | 600 | 360 | 43 | Illusion 1.0L target |
| `ill_15L_sc` | illusion | 800 | 360 | 43 | Illusion 1.5L target |
| `ill_2L_sc` | illusion | 800 | 360 | 43 | Illusion 2.0L target (new) — training bug known |
> **Note on illusion scenes**: Training bug — all target diameters were accidentally set to 1L. Transfer results invalid.
### 2.2 Legacy REPRODUCE_SCENES
Handled by @src/drl_pinball/reproduce/ — not by eval.
---
## 3. Train Inference Pipeline (`infer_train.py`)
### 3.1 Per-Scene Workflow
```
For each scene in TRAIN_SCENES:
1. Load calibration.json + target.npy (+ target_harmonics.json for illusion)
2. Create KarmanCloakEnv / IllusionCloakEnv (fresh CFD init)
3. For each seed (seed_label, model_dir):
a. Create skeleton PPO(env=vec_env, Sin, [64,64])
b. Extract policy weights from best_model.zip
c. Load VecNormalize from seed's vec_normalize.pkl (frozen)
d. 360-step deterministic rollout
e. Record: sensors, forces, actions, rewards, per-component r_cd/r_cl/r_sim
4. Pick best seed by tail-180 avg reward
5. Re-create env, load best seed model, run → capture vorticity PNGs
6. Generate target vorticity (dist_cyl only / target cyl only)
7. Generate zero-action baseline vorticity
8. Write: signals.npz, metrics.json, all_seeds.json, vorticity_*.png
```
### 3.2 Skeleton Injection Pattern
Due to `numpy._core.numeric` cloudpickle deserialization issues with SB3 models trained on certain Python versions, `PPO.load()` may fail. The fallback is **skeleton injection**:
```python
# Method: skeleton PPO + manual weight injection
skeleton = PPO(
"MlpPolicy",
policy_kwargs={"activation_fn": Sin, "net_arch": [64, 64]},
env=vec_env, device=device,
n_steps=2048, batch_size=64, n_epochs=10,
learning_rate=3e-4, gamma=0.995, verbose=0,
)
# Extract weights from zip
with zipfile.ZipFile("best_model.zip") as zf:
with zf.open("policy.pth") as f:
state_dict = torch.load(io.BytesIO(f.read()), map_location="cpu")
skeleton.policy.load_state_dict(state_dict, strict=False)
# VecNormalize loaded separately
vec_env = VecNormalize.load("vec_normalize.pkl", vec_env)
vec_env.training = False # Frozen statistics
vec_env.norm_reward = False
```
### 3.3 Output Files per Scene
```
eval/output/train/{scene_id}/
├── signals.npz # sensors (360,6), forces (360,6), actions (360,3), rewards (360,)
├── metrics.json # Best seed summary: DTW, reward, action stats
├── all_seeds.json # Per-seed breakdown: reward, r_cd, r_cl, r_sim, sim_raw, dt_sec
├── vorticity_controlled.png # Final frame after full DRL rollout
├── vorticity_target.png # Target state (disturbance only / target cylinder only)
└── vorticity_zero.png # Zero-action baseline (no control)
```
### 3.4 Config Switching Delay
When consecutive scenes use different LBM config files (different ν → different kernel compilation), a 120-second delay is inserted to allow the previous GPU context to fully release before the new config triggers PTX recompilation. Scenes sharing the same config file run back-to-back without delay.
---
## 4. Baseline Collection (`collect_baselines.py`)
Collects `q_in.npz` (background flow) and `q_blk.npz` (zero-action pinball) for each scene.
Baselines are needed by OID and CCD analysis pipelines.
```bash
python collect_baselines.py --device 2 # all scenes
python collect_baselines.py --scene kar_re100_sc # single scene
```
## 5. SR Bridge (`bridge_to_sr.py`)
Converts V5 eval output to the format expected by `SR_analysis/`:
- `calibration.json``norm.json` (FORCE_SCALE, SENS_SCALE, sens_deviation=0)
- `signals.npz``controlled.npz` (sensors, forces, actions, rewards)
- `target.npy``target.npz`
## 6. GPU Scheduling (`run_all_with_fields.sh`)
Single GPU (2). Config switch waits 120s to avoid PTX compilation conflicts.
```
GPU2: [kar_re100_sc(5 seeds)] → [vardist_sc+tr] → [ill_*_sc]
--120s--> [kar_re60_sc] --120s--> [kar_re200_sc] --120s--> [kar_re400_sc]
--120s--> [cross-re transfers]
```
### Run Commands
```bash
# Full pipeline (infer + baselines + bridge):
bash run_all_with_fields.sh
# Single scene:
bash run_all_with_fields.sh --scene kar_re100_sc --gpu 2
python infer_train.py --scene kar_re100_sc --device-id 2
python viz_flow.py
```
---
## 6. Eval Results
### 6.1 Train Pipeline — Completed Cases
| Scene | Best Seed | DTW sim_raw | Reward | r_cd | r_cl | r_sim |
|-------|:---------:|:-----------:|:------:|:----:|:----:|:-----:|
| kar_re100_sc | 45 | **0.926** | 0.941 | 0.984 | 0.990 | 0.872 |
| kar_d075_sc | 44 | **0.911** | 0.949 | 0.966 | 0.974 | 0.919 |
| kar_re60_tr | 43 | 0.187 | 0.312 | 0.659 | 0.334 | 0.091 |
| kar_re200_tr | 43 | 0.506 | 0.367 | 0.669 | 0.212 | 0.278 |
| kar_re400_tr | 43 | 0.428 | 0.399 | 0.737 | 0.400 | 0.176 |
### 6.2 Train Pipeline — Re100 5-Seed Breakdown
| Seed | DTW sim_raw | r_cd | r_cl | r_sim |
|:----:|:-----------:|:----:|:----:|:-----:|
| 41 | 0.634 | 0.307 | 0.241 | 0.388 |
| 42 | 0.893 | 0.959 | 0.762 | 0.788 |
| 43 | 0.850 | 0.908 | 0.469 | 0.629 |
| 44 | 0.523 | 0.282 | 0.177 | 0.322 |
| 45 | **0.926** | 0.984 | 0.990 | 0.872 |
Mean DTW across seeds: 0.765 ± 0.172. Best-worst spread: 0.403.
### 6.3 DTW Formula Note
**Important**: V5 uses a different DTW formula than Legacy pipelines. See `recompute_unified_dtw.py` for fair comparison using the same formula.
---
## 7. Visualization & Reporting Tools
### 7.1 `viz_signals.py`
Per-scene diagnostic plots (4×2 subplot grid): sensors, forces, actions, reward, FFT, phase portrait.
### 7.2 `viz_flow.py`
Cross-scene vorticity comparison panels with unified [-0.003, 0.003] colormap.
### 7.3 `generate_report.py`
Master comparison report with DTW bar chart.
---
## 8. Known Limitations
- **Illusion training bug**: All target diameters were accidentally set to 1L. Transfer results invalid.
- **VecNormalize compat**: `PPO.load()` may fail with cloudpickle on Python 3.10+. Skeleton injection fallback handles this.
- **GPU cleanup**: Between scenes, `pkill -f pycuda; sleep 180` may be needed.
---
## 9. Companion Documents
| File | Content |
|------|---------|
| `eval/README.md` | Quick start + self-review checklist |
| `eval/scene_manifest.py` | TRAIN_SCENES single source of truth |
| `eval/infer_train.py` | V5 train model inference |
| `eval/recompute_unified_dtw.py` | Fair DTW comparison (unified formula) |
| `train/TRAIN_PIPELINE.md` | Training pipeline documentation |
| `reproduce/REPRODUCE_KNOWLEDGE.md` | Reproduce module (separate) |
| `../SR_analysis/README.md` | SR analysis pipeline overview |
---
*Document last updated: 2026-07-13.*
@@ -0,0 +1,149 @@
#!/usr/bin/env python3
"""Bridge V5 train eval outputs to SR_analysis format.
Converts eval pipeline outputs to the format expected by SR_analysis:
- signals.npz -> controlled.npz (sensors, forces, actions, rewards)
- calibration.json -> norm.json (FORCE_SCALE, SENS_SCALE, sens_deviation=0)
- target.npy -> target.npz (wrap in dict with "target_states" key)
- Copies target_harmonics.json if present (illusion scenes)
Legacy model bridge is handled by @src/drl_pinball/reproduce/ — NOT here.
Usage:
conda run -n pycuda_3_10 python bridge_to_sr.py
conda run -n pycuda_3_10 python bridge_to_sr.py --scene kar_re100
"""
from __future__ import annotations
import argparse
import json
import shutil
import sys
from pathlib import Path
_SRC = Path(__file__).resolve().parents[2]
if str(_SRC) not in sys.path:
sys.path.insert(0, str(_SRC))
import numpy as np
from drl_pinball.case_registry import CASE_REGISTRY, get_case
_THIS = Path(__file__).resolve().parent
_EVAL_OUT = _THIS / "output" / "train"
_TRAIN_DIR = _THIS.parent / "train"
def log(msg: str) -> None:
print(f"[bridge] {msg}", flush=True)
def bridge_scene(out_dir: Path, scene_id: str) -> bool:
"""Bridge a V5 train scene to SR_analysis format."""
signals = out_dir / "signals.npz"
# Find calibration.json
cal_candidates = [
_TRAIN_DIR / "calibrations" / get_case(scene_id).calibration / "calibration.json",
]
for seed_dir in sorted((_TRAIN_DIR / "output").glob(f"{scene_id}_seed*/")) if (_TRAIN_DIR / "output").exists() else []:
cal_candidates.append(seed_dir / "calibration.json")
cal_path = None
for c in cal_candidates:
if c.exists():
cal_path = c
break
if cal_path is None:
log(f" SKIP {scene_id}: calibration.json not found")
return False
if not signals.exists():
log(f" SKIP {scene_id}: signals.npz not found")
return False
sr_dir = out_dir / "sr_bridge"
sr_dir.mkdir(exist_ok=True)
# 1. signals.npz -> controlled.npz
data = dict(np.load(signals, allow_pickle=True))
np.savez_compressed(sr_dir / "controlled.npz",
sensors=data.get("sensors"), forces=data.get("forces"),
actions=data.get("actions"), rewards=data.get("rewards"))
# 2. calibration.json -> norm.json
with open(cal_path) as f:
cal = json.load(f)
force_norm_fact = float(cal.get("FORCE_SCALE", 1.0))
sens_norm_fact = float(cal.get("SENS_SCALE", 1.0))
norm = {
"force_norm_fact": float(force_norm_fact),
"sens_deviation": [0.0] * 6,
"sens_norm_fact": [float(sens_norm_fact)] * 6,
"action_bias": cal.get("ACTION_BIAS", [0.0, 0.0, 0.0]),
}
with open(sr_dir / "norm.json", "w") as f:
json.dump(norm, f, indent=2)
# 3. target.npy -> target.npz
target_candidates = [
_TRAIN_DIR / "calibrations" / get_case(scene_id).calibration / "target.npy",
]
for seed_dir in sorted((_TRAIN_DIR / "output").glob(f"{scene_id}_seed*/")) if (_TRAIN_DIR / "output").exists() else []:
target_candidates.append(seed_dir / "target.npy")
for tp in target_candidates:
if tp.exists():
target_states = np.load(str(tp))
np.savez_compressed(sr_dir / "target.npz", target_states=target_states)
break
# 4. target_harmonics.json (illusion scenes)
harmonics_candidates = [
_TRAIN_DIR / "calibrations" / get_case(scene_id).calibration / "target_harmonics.json",
]
for seed_dir in sorted((_TRAIN_DIR / "output").glob(f"{scene_id}_seed*/")) if (_TRAIN_DIR / "output").exists() else []:
harmonics_candidates.append(seed_dir / "target_harmonics.json")
for hp in harmonics_candidates:
if hp.exists():
shutil.copy(hp, sr_dir / "target_harmonics.json")
break
config = {
"scene_id": scene_id,
"SI": cal.get("SI", 800),
"FIFO_LEN": cal.get("FIFO_LEN", 150),
"sample_interval": cal.get("SI", 800),
}
with open(sr_dir / "config.json", "w") as f:
json.dump(config, f, indent=2)
log(f" controlled.npz, norm.json, target.npz -> {sr_dir}")
return True
def main() -> int:
parser = argparse.ArgumentParser(description="Bridge V5 train eval output to SR_analysis format")
parser.add_argument("--scene", type=str, default=None)
args = parser.parse_args()
base = _EVAL_OUT
if not base.exists():
log(f"ERROR: output dir not found: {base}")
return 1
dirs = [base / args.scene] if args.scene else sorted([d for d in base.iterdir() if d.is_dir()])
if not dirs:
log("No scenes found.")
return 0
log(f"Bridging {len(dirs)} train scenes...")
for d in dirs:
bridge_scene(d, d.name)
log("All bridges complete.")
return 0
if __name__ == "__main__":
raise SystemExit(main())
@@ -0,0 +1,208 @@
#!/usr/bin/env python3
"""Collect baseline flow fields (q_in, q_blk) for V5 train scenes.
q_in = background flow (disturbance-only for karman, empty channel for illusion)
q_blk = pinball with zero rotation (passive blockage)
Legacy model baselines are handled by @src/drl_pinball/reproduce/ — NOT here.
Usage:
conda run -n pycuda_3_10 python collect_baselines.py --scene kar_re100 --device 2
conda run -n pycuda_3_10 python collect_baselines.py --device 2
"""
from __future__ import annotations
import argparse
import json
import sys
import time
from pathlib import Path
from typing import Any, Dict, List, Tuple
import numpy as np
import pycuda.driver as cuda; cuda.init()
_REPO = str(Path(__file__).resolve().parents[3])
_SRC = Path(_REPO) / "src"
for p in [_REPO, str(_SRC)]:
if p not in sys.path:
sys.path.insert(0, p)
from CelerisLab import Simulation
from drl_pinball.eval.scene_manifest import TRAIN_SCENES
L0 = 20.0
U0 = 0.01
RADIUS = L0 / 2.0
_OUT_BASE = Path(__file__).resolve().parent / "output" / "train"
def log(msg: str) -> None:
print(f"[{time.strftime('%H:%M:%S')}] {msg}", flush=True)
def save_fields(path: Path, ux_list: list, uy_list: list) -> None:
ux = np.array(ux_list, dtype=np.float32)
uy = np.array(uy_list, dtype=np.float32)
np.savez_compressed(path, ux=ux, uy=uy)
log(f" Saved {path.name}: ux{ux.shape}, uy{uy.shape}")
def collect_fields(sim, num_steps: int, si: int) -> Tuple[list, list]:
"""Run simulation and collect velocity field at each step."""
ux_list, uy_list = [], []
for _ in range(num_steps):
sim.run(si, zero_obs=True)
macro = sim.get_macroscopic()
ux_list.append(macro["ux"].copy())
uy_list.append(macro["uy"].copy())
return ux_list, uy_list
def _train_geom(scene: Dict[str, Any]) -> tuple:
"""Extract train pipeline geometry from calibration.json."""
with open(scene["calibration_path"]) as f:
cal = json.load(f)
nx_cfg = int(cal.get("grid", {}).get("nx", 2000))
ny_cfg = int(cal.get("grid", {}).get("ny", 600))
cy = float(ny_cfg - 1) / 2.0
return (scene["config_path"], cy, 600.0, 1000.0, 1026.0, 1200.0, 380.0, 406.0, 600.0)
def _warmup_train(sim, nx=2000):
sim.run(int(4.0 * nx / U0), zero_obs=True)
# ── Karman ─────────────────────────────────────────────────────────────
def collect_karman_q_in(scene: Dict[str, Any], device_id: int) -> None:
scene_id = scene["scene_id"]
si, num_steps = scene["si"], scene["num_steps"]
cfg, cy, dist_x, _, _, sens_x, _, _, _ = _train_geom(scene)
out_dir = _OUT_BASE / scene_id
out_dir.mkdir(parents=True, exist_ok=True)
with open(scene["calibration_path"]) as f:
cal = json.load(f)
dist_radius = float(cal.get("dist_radius", 1.0)) * L0
sim = Simulation(lbm_config_path=cfg, device_id=device_id)
sim.add_body("circle", center=(dist_x, cy, 0.0), radius=dist_radius)
sim.add_body("sensor", center=(sens_x, cy + 40.0, 0.0), radius=5.0)
sim.add_body("sensor", center=(sens_x, cy, 0.0), radius=5.0)
sim.add_body("sensor", center=(sens_x, cy - 40.0, 0.0), radius=5.0)
sim.initialize(); _warmup_train(sim)
ux_list, uy_list = collect_fields(sim, num_steps, si)
save_fields(out_dir / "q_in.npz", ux_list, uy_list)
sim.close()
def collect_karman_q_blk(scene: Dict[str, Any], device_id: int) -> None:
scene_id = scene["scene_id"]
si, num_steps = scene["si"], scene["num_steps"]
cfg, cy, dist_x, pb_x, pb_rx, sens_x, _, _, _ = _train_geom(scene)
out_dir = _OUT_BASE / scene_id
out_dir.mkdir(parents=True, exist_ok=True)
with open(scene["calibration_path"]) as f:
cal = json.load(f)
dist_radius = float(cal.get("dist_radius", 1.0)) * L0
sim = Simulation(lbm_config_path=cfg, device_id=device_id)
sim.add_body("circle", center=(dist_x, cy, 0.0), radius=dist_radius)
sim.add_body("sensor", center=(sens_x, cy + 40.0, 0.0), radius=5.0)
sim.add_body("sensor", center=(sens_x, cy, 0.0), radius=5.0)
sim.add_body("sensor", center=(sens_x, cy - 40.0, 0.0), radius=5.0)
sim.add_body("circle", center=(pb_x, cy, 0.0), radius=RADIUS)
sim.add_body("circle", center=(pb_rx, cy + 15.0, 0.0), radius=RADIUS)
sim.add_body("circle", center=(pb_rx, cy - 15.0, 0.0), radius=RADIUS)
sim.initialize(); _warmup_train(sim)
ux_list, uy_list = collect_fields(sim, num_steps, si)
save_fields(out_dir / "q_blk.npz", ux_list, uy_list)
sim.close()
# ── Illusion ───────────────────────────────────────────────────────────
def collect_illusion_q_in(scene: Dict[str, Any], device_id: int) -> None:
scene_id = scene["scene_id"]
si, num_steps = scene["si"], scene["num_steps"]
cfg, cy, _, _, _, _, _, _, ill_sens_x = _train_geom(scene)
out_dir = _OUT_BASE / scene_id
out_dir.mkdir(parents=True, exist_ok=True)
sim = Simulation(lbm_config_path=cfg, device_id=device_id)
sim.add_body("sensor", center=(ill_sens_x, cy + 40.0, 0.0), radius=5.0)
sim.add_body("sensor", center=(ill_sens_x, cy, 0.0), radius=5.0)
sim.add_body("sensor", center=(ill_sens_x, cy - 40.0, 0.0), radius=5.0)
sim.initialize(); _warmup_train(sim)
ux_list, uy_list = collect_fields(sim, num_steps, si)
save_fields(out_dir / "q_in.npz", ux_list, uy_list)
sim.close()
def collect_illusion_q_blk(scene: Dict[str, Any], device_id: int) -> None:
scene_id = scene["scene_id"]
si, num_steps = scene["si"], scene["num_steps"]
cfg, cy, _, _, _, _, ill_pb_x, ill_pb_rx, ill_sens_x = _train_geom(scene)
out_dir = _OUT_BASE / scene_id
out_dir.mkdir(parents=True, exist_ok=True)
sim = Simulation(lbm_config_path=cfg, device_id=device_id)
sim.add_body("sensor", center=(ill_sens_x, cy + 40.0, 0.0), radius=5.0)
sim.add_body("sensor", center=(ill_sens_x, cy, 0.0), radius=5.0)
sim.add_body("sensor", center=(ill_sens_x, cy - 40.0, 0.0), radius=5.0)
sim.add_body("circle", center=(ill_pb_x, cy, 0.0), radius=RADIUS)
sim.add_body("circle", center=(ill_pb_rx, cy + 15.0, 0.0), radius=RADIUS)
sim.add_body("circle", center=(ill_pb_rx, cy - 15.0, 0.0), radius=RADIUS)
sim.initialize(); _warmup_train(sim)
ux_list, uy_list = collect_fields(sim, num_steps, si)
save_fields(out_dir / "q_blk.npz", ux_list, uy_list)
sim.close()
# ── Dispatcher ─────────────────────────────────────────────────────────
Q_IN = {"karman": collect_karman_q_in, "illusion": collect_illusion_q_in}
Q_BLK = {"karman": collect_karman_q_blk, "illusion": collect_illusion_q_blk}
def main() -> int:
parser = argparse.ArgumentParser(description="Collect baseline flow fields for V5 train scenes")
parser.add_argument("--scene", type=str, default=None)
parser.add_argument("--device", type=int, default=2)
args = parser.parse_args()
log(f"Baseline collection: train pipeline, GPU {args.device}")
scenes = TRAIN_SCENES
if args.scene:
scenes = [s for s in scenes if s["scene_id"] == args.scene or s["scene_id"].startswith(args.scene)]
if not scenes:
log(f"ERROR: No scene matching '{args.scene}'")
return 1
for i, scene in enumerate(scenes):
if i > 0:
prev_cfg = scenes[i - 1].get("config_path", "")
curr_cfg = scene.get("config_path", "")
if prev_cfg != curr_cfg:
log(f"Waiting 120s before config switch...")
time.sleep(120)
scene_id = scene["scene_id"]
s_type = scene["scene_type"]
t0 = time.perf_counter()
log(f"[{i+1}/{len(scenes)}] {scene_id} (type={s_type})")
if s_type in Q_IN:
Q_IN[s_type](scene, args.device)
if s_type in Q_BLK:
Q_BLK[s_type](scene, args.device)
log(f" Done ({time.perf_counter() - t0:.0f}s)")
log("All baseline collections complete.")
return 0
if __name__ == "__main__":
raise SystemExit(main())
@@ -0,0 +1,146 @@
#!/usr/bin/env python3
"""Compare Legacy (1280x512, parabolic, bounce-back) vs V5 (2000x600, uniform, free-slip)
DTW similarity scores. Generates a grouped bar chart.
**DEPRECATED (2026-07-13):** This script compares the raw metrics.json scores
which use DIFFERENT DTW formulas. Legacy uses 1-dtw/n; V5 uses 1-dtw/(n*norm_scale).
This makes V5 scores appear ~0.05-0.3 lower than they truly are.
For a FAIR comparison using a unified formula, use `recompute_unified_dtw.py`.
Unified results: V5 uniformly OUTPERFORMS Legacy on Karman cross-Re scenes.
Keeping this script for reference but new work should use the unified script."""
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
# ── Data ──────────────────────────────────────────────────────────────
# Format: (label, legacy_dtw, v5_dtw)
pairs = [
("Karman cloaking\nRe=100", 0.954, 0.918),
("Karman cloaking\nRe=200", 0.884, 0.712),
("Karman cloaking\nRe=400", 0.795, 0.565),
("Illusion 0.75L", 0.980, 0.807),
("Illusion 1.0L", 0.975, 0.900),
("Illusion 1.5L", 0.945, 0.906),
("Illusion 2.0L", None, 0.800), # V5 only, no legacy
]
# V5-only scenes (no legacy counterpart)
v5_only = [
("Karman-cloak\nRe=60", 0.364),
("Karman-cloak\nVarDist d=0.75", 0.911),
("Karman-cloak\nVarDist d=1.5", 0.892),
("Karman-cloak\nVarDist d=2.0", 0.816),
]
v5_transfer = [
("Re=60 transfer", 0.148),
("Re=200 transfer", 0.509),
("Re=400 transfer", 0.480),
("d=0.75 transfer", 0.369),
("d=1.5 transfer", 0.783),
("d=2.0 transfer", 0.452),
]
# ── Figure 1: Legacy vs V5 grouped bar chart ──────────────────────────
fig, ax = plt.subplots(figsize=(12, 6))
labels = [p[0].replace("\n", " ") for p in pairs if p[1] is not None]
n = len(labels)
x = np.arange(n)
w = 0.35
legacy_vals = [p[1] for p in pairs if p[1] is not None]
v5_vals = [p[2] for p in pairs if p[1] is not None]
bars1 = ax.bar(x - w/2, legacy_vals, w, label="Legacy (1280x512, parabolic, bounce-back)", color="#2196F3", edgecolor="white")
bars2 = ax.bar(x + w/2, v5_vals, w, label="V5 (2000x600, uniform, free-slip)", color="#FF9800", edgecolor="white")
# Annotate values
for bar, val in zip(bars1, legacy_vals):
ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.01, f"{val:.3f}",
ha="center", va="bottom", fontsize=8)
for bar, val in zip(bars2, v5_vals):
ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.01, f"{val:.3f}",
ha="center", va="bottom", fontsize=8)
ax.set_ylabel("DTW Similarity", fontsize=12)
ax.set_title("Legacy vs V5: DTW Similarity (Karman Cloaking + Illusion)", fontsize=14, fontweight="bold")
ax.set_xticks(x)
ax.set_xticklabels(labels, fontsize=9)
ax.legend(loc="lower right", fontsize=9)
ax.set_ylim(0, 1.10)
ax.grid(axis="y", alpha=0.3)
# Add gap markers for significant drops
for i, (lbl, leg, v5) in enumerate(zip(labels, legacy_vals, v5_vals)):
gap = leg - v5
color = "red" if gap > 0.15 else ("green" if gap < 0.05 else "gray")
ax.annotate(f"gap: {gap:+.3f}", xy=(i, min(leg, v5) + abs(gap)/2),
fontsize=7, color=color, ha="center", fontweight="bold")
import os
out = os.path.join(os.path.dirname(__file__), "output", "reports", "legacy_vs_v5_dtw.png")
os.makedirs(os.path.dirname(out), exist_ok=True)
plt.tight_layout()
plt.savefig(out, dpi=150)
print(f"Saved: {out}")
plt.close()
# ── Figure 2: V5-only scenes (scratch + transfer) ─────────────────────
fig2, (ax2a, ax2b) = plt.subplots(1, 2, figsize=(14, 5))
# Scratch
scr_labels = [s[0].replace("\n", "") for s in v5_only]
scr_vals = [s[1] for s in v5_only]
colors2 = ["#4CAF50", "#8BC34A", "#CDDC39", "#FFC107"]
bars = ax2a.bar(range(len(scr_vals)), scr_vals, color=colors2, edgecolor="white")
for bar, val in zip(bars, scr_vals):
ax2a.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.01, f"{val:.3f}",
ha="center", fontsize=9)
ax2a.set_xticks(range(len(scr_labels)))
ax2a.set_xticklabels(scr_labels, fontsize=8)
ax2a.set_ylabel("DTW Similarity")
ax2a.set_title("V5 Scratch Models (no legacy counterpart)")
ax2a.set_ylim(0, 1.10)
ax2a.grid(axis="y", alpha=0.3)
# Transfer
tr_labels = [t[0] for t in v5_transfer]
tr_vals = [t[1] for t in v5_transfer]
colors3 = ["#FF5722", "#E91E63", "#9C27B0", "#673AB7", "#3F51B5", "#2196F3"]
bars = ax2b.bar(range(len(tr_vals)), tr_vals, color=colors3, edgecolor="white")
for bar, val in zip(bars, tr_vals):
ax2b.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.01, f"{val:.3f}",
ha="center", fontsize=9)
ax2b.set_xticks(range(len(tr_labels)))
ax2b.set_xticklabels(tr_labels, fontsize=8)
ax2b.set_ylabel("DTW Similarity")
ax2b.set_title("V5 Transfer Models (no legacy counterpart)")
ax2b.set_ylim(0, 1.10)
ax2b.grid(axis="y", alpha=0.3)
out2 = os.path.join(os.path.dirname(__file__), "output", "reports", "v5_only_dtw.png")
plt.tight_layout()
plt.savefig(out2, dpi=150)
print(f"Saved: {out2}")
plt.close()
# ── Text summary ──────────────────────────────────────────────────────
print("\n=== Summary Table ===")
print(f"{'Scene':<30} {'Legacy':>8} {'V5':>8} {'Delta':>8}")
print("-" * 58)
for pair in pairs:
label = pair[0].split("\n")[0]
leg = f"{pair[1]:.3f}" if pair[1] else "N/A"
v5 = f"{pair[2]:.3f}"
delta = f"{(pair[2] - pair[1]):+.3f}" if pair[1] else "N/A"
print(f"{label:<30} {leg:>8} {v5:>8} {delta:>8}")
avg_legacy = np.mean([p[1] for p in pairs if p[1] is not None])
avg_v5_matched = np.mean([p[2] for p in pairs if p[1] is not None])
print(f"\nAverage (matched scenes): Legacy={avg_legacy:.3f}, V5={avg_v5_matched:.3f}, gap={avg_v5_matched-avg_legacy:+.3f}")
avg_v5_all = np.mean([p[2] for p in pairs])
print(f"Average (all V5 scenes): {avg_v5_all:.3f}")
@@ -1,10 +1,10 @@
#!/usr/bin/env python3
"""Generate master comparison report from eval outputs.
Reads metrics.json from all output/{train,reproduce}/ directories,
Reads metrics.json from output/train/ directories,
produces:
1. Master comparison table (markdown)
2. Summary bar chart (train vs reproduce DTW similarity)
2. Summary bar chart (DTW similarity)
3. Action statistics summary
Usage:
@@ -50,52 +50,45 @@ def load_all_seeds(scene_dir: Path) -> List[Dict]:
# Tables
# ---------------------------------------------------------------------------
def generate_comparison_table() -> str:
"""Generate markdown comparison table for overlapping scenes."""
overlap = [
("re100_karman", "karman_cloak", "Karman Cloak Re100"),
("illusion_075L", "illusion_075L", "Illusion 0.75L"),
("illusion_1L", "illusion_1L", "Illusion 1.0L"),
("illusion_15L", "illusion_15L", "Illusion 1.5L"),
"""Generate markdown comparison table for all train scenes."""
# Core scenes
core = [
("kar_re100_sc", "Karman Cloak Re100"),
("ill_075L_sc", "Illusion 0.75L"),
("ill_15L_sc", "Illusion 1.5L"),
("ill_2L_sc", "Illusion 2.0L"),
]
lines = []
lines.append("# Master Comparison Table")
lines.append("# V5 Train Model Evaluation")
lines.append("")
lines.append("| Scene | Side | Best Seed | DTW sim | Reward | "
lines.append("| Scene | Best Seed | DTW sim | Reward | "
"aF mean | aB mean | aT mean | r_cd | r_cl | r_sim |")
lines.append("|-------|------|-----------|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|")
lines.append("|-------|-----------|:---:|:---:|:---:|:---:|:---:|:---:|:---:|:---:|")
for trn_key, rep_key, name in overlap:
m_trn = load_metrics(_OUT_BASE / "train" / trn_key)
m_rep = load_metrics(_OUT_BASE / "reproduce" / rep_key)
dtw_t = _fmt(m_trn.get("dtw_sim_v5") or m_trn.get("sim_raw_mean"))
dtw_r = _fmt(m_rep.get("dtw_sim_v5") or m_rep.get("sim_raw_mean"))
# Train row
for trn_key, name in core:
m = load_metrics(_OUT_BASE / "train" / trn_key)
dtw_t = _fmt(m.get("dtw_sim_v5") or m.get("sim_raw_mean"))
lines.append(
f"| {name} | Train | {m_trn.get('best_seed', '-')} | "
f"{dtw_t} | {_fmt(m_trn.get('reward_mean'))} | "
f"{_fmt(m_trn.get('aF_mean'))} | {_fmt(m_trn.get('aB_mean'))} | "
f"{_fmt(m_trn.get('aT_mean'))} | "
f"{_fmt(m_trn.get('r_cd_mean'))} | {_fmt(m_trn.get('r_cl_mean'))} | "
f"{_fmt(m_trn.get('r_sim_mean'))} |"
)
# Reproduce row
lines.append(
f"| {name} | Reproduce | - | "
f"{dtw_r} | - | "
f"- | - | - | - | - | - |"
f"| {name} | {m.get('best_seed', '-')} | "
f"{dtw_t} | {_fmt(m.get('reward_mean'))} | "
f"{_fmt(m.get('aF_mean'))} | {_fmt(m.get('aB_mean'))} | "
f"{_fmt(m.get('aT_mean'))} | "
f"{_fmt(m.get('r_cd_mean'))} | {_fmt(m.get('r_cl_mean'))} | "
f"{_fmt(m.get('r_sim_mean'))} |"
)
# Train-only scenes
train_only = ["transfer_re60", "transfer_re200", "transfer_re400"]
# Cross-Re / VarDist scenes
extra = [
"kar_re60_sc", "kar_re200_sc", "kar_re400_sc",
"kar_d075_sc", "kar_d15_sc", "kar_d2_sc",
]
lines.append("")
lines.append("## Train Only (Cross-Re Transfer)")
lines.append("## Cross-Re Scratch / VarDist Scratch")
lines.append("| Scene | Best Seed | Reward | "
"r_cd | r_cl | r_sim | sim_raw |")
lines.append("|-------|-----------|---:|-----:|-----:|-----:|---:|")
for trn_key in train_only:
for trn_key in extra:
m = load_metrics(_OUT_BASE / "train" / trn_key)
lines.append(
f"| {trn_key} | {m.get('best_seed', '-')} | "
@@ -104,20 +97,15 @@ def generate_comparison_table() -> str:
f"{_fmt(m.get('r_sim_mean'))} | {_fmt(m.get('sim_raw_mean'))} |"
)
# Reproduce-only
repro_only = ["steady_cloak", "vortex_lamb", "vortex_taylor"]
lines.append("")
lines.append("## Reproduce Only")
lines.append("| Scene | Notes |")
lines.append("|-------|-------|")
for r_key in repro_only:
lines.append(f"| {r_key} | open-loop / vortex |")
# Per-scene all-seeds detail
lines.append("")
lines.append("## Per-Scene Seed Details (Train)")
for trn_key in ["re100_karman", "transfer_re60", "transfer_re200",
"transfer_re400", "illusion_075L", "illusion_1L", "illusion_15L"]:
lines.append("## Per-Scene Seed Details")
all_scenes = [
"kar_re100_sc", "kar_d075_sc", "kar_d15_sc", "kar_d2_sc",
"kar_re60_sc", "kar_re200_sc", "kar_re400_sc",
"ill_075L_sc", "ill_15L_sc", "ill_2L_sc",
]
for trn_key in all_scenes:
seeds = load_all_seeds(_OUT_BASE / "train" / trn_key)
if not seeds:
continue
@@ -153,37 +141,38 @@ def _fmt_int(v):
# Bar chart
# ---------------------------------------------------------------------------
def generate_dtw_barchart() -> None:
overlap = [
("re100_karman", "karman_cloak", "Karman Cloak"),
("illusion_075L", "illusion_075L", "Illusion 0.75L"),
("illusion_1L", "illusion_1L", "Illusion 1L"),
("illusion_15L", "illusion_15L", "Illusion 1.5L"),
scenes = [
("kar_re100_sc", "Karman Re100"),
("ill_075L_sc", "Illusion 0.75L"),
("ill_15L_sc", "Illusion 1.5L"),
("ill_2L_sc", "Illusion 2L"),
("kar_re60_sc", "Re60 Scr"),
("kar_re200_sc", "Re200 Scr"),
("kar_re400_sc", "Re400 Scr"),
("kar_d075_sc", "d075 Scr"),
("kar_d15_sc", "d15 Scr"),
("kar_d2_sc", "d2 Scr"),
]
names = []
trn_vals = []
rep_vals = []
for trn_key, rep_key, name in overlap:
m_trn = load_metrics(_OUT_BASE / "train" / trn_key)
m_rep = load_metrics(_OUT_BASE / "reproduce" / rep_key)
t = m_trn.get("dtw_sim_v5") or m_trn.get("sim_raw_mean", 0)
r = m_rep.get("dtw_sim_v5") or m_rep.get("sim_raw_mean", 0)
vals = []
for trn_key, name in scenes:
m = load_metrics(_OUT_BASE / "train" / trn_key)
v = m.get("dtw_sim_v5") or m.get("sim_raw_mean", 0)
names.append(name)
trn_vals.append(t if t and t > 0 else 0)
rep_vals.append(r if r and r > 0 else 0)
vals.append(v if v and v > 0 else 0)
fig, ax = plt.subplots(figsize=(10, 5))
x = np.arange(len(names))
w = 0.35
ax.bar(x - w / 2, rep_vals, w, label="Reproduce (legacy)", color="#d62728")
ax.bar(x + w / 2, trn_vals, w, label="Train (V5)", color="#1f77b4")
ax.set_xticks(x); ax.set_xticklabels(names, rotation=15, ha="right")
fig, ax = plt.subplots(figsize=(12, 5))
bars = ax.bar(range(len(names)), vals, color="#1f77b4", alpha=0.7)
ax.set_xticks(range(len(names)))
ax.set_xticklabels(names, rotation=30, ha="right")
ax.set_ylabel("DTW Similarity"); ax.set_ylim(0, 1.05)
ax.set_title("Train vs Reproduce - DTW Similarity")
ax.legend(); ax.grid(axis="y", alpha=0.3)
for i, (tv, rv) in enumerate(zip(trn_vals, rep_vals)):
if tv > 0: ax.text(i + w / 2, tv + 0.01, f"{tv:.3f}", ha="center", fontsize=9)
if rv > 0: ax.text(i - w / 2, rv + 0.01, f"{rv:.3f}", ha="center", fontsize=9)
ax.set_title("V5 Train Models - DTW Similarity")
ax.grid(axis="y", alpha=0.3)
for b, v in zip(bars, vals):
if v > 0:
ax.text(b.get_x() + b.get_width() / 2, v + 0.01,
f"{v:.3f}", ha="center", fontsize=8)
fig.tight_layout()
fig.savefig(_REPORT_DIR / "summary_dtw_barchart.png", dpi=150,
bbox_inches="tight")
@@ -197,7 +186,6 @@ def generate_dtw_barchart() -> None:
def main() -> int:
_REPORT_DIR.mkdir(parents=True, exist_ok=True)
log(f"Train dir: {_OUT_BASE / 'train'}")
log(f"Reproduce dir: {_OUT_BASE / 'reproduce'}")
log(f"Reports dir: {_REPORT_DIR}")
table = generate_comparison_table()
@@ -38,6 +38,10 @@ for p in [_REPO, str(_SRC)]:
import torch
from torch.nn import Module as TorchModule
from stable_baselines3 import PPO
from drl_pinball.reproduce.core.action_wrapper import (
norm_action_to_omega,
surface_vel_to_omega,
)
from CelerisLab import Simulation
from CelerisLab.common.render import compute_vorticity, render_vorticity_field
@@ -103,10 +107,25 @@ def get_cc(sim, sid):
return float(len(cells_arr))
def capture_field(sim, ux_list: list, uy_list: list) -> None:
"""Append current macroscopic velocity field to lists."""
macro = sim.get_macroscopic()
ux_list.append(macro["ux"].copy())
uy_list.append(macro["uy"].copy())
def save_fields(path: Path, ux_list: list, uy_list: list) -> None:
"""Save accumulated velocity field time series as fields.npz."""
ux = np.array(ux_list, dtype=np.float32)
uy = np.array(uy_list, dtype=np.float32)
np.savez_compressed(path, ux=ux, uy=uy)
log(f" Saved fields.npz: ux{ux.shape}, uy{uy.shape} -> {path.name}")
def action_to_omega_legacy(action_norm, scale=8.0, bias=(0.0, -4.0, 4.0)):
b = np.array(bias, dtype=np.float32)
sv = (np.asarray(action_norm, dtype=np.float32) * scale + b) * U0
return -sv / RADIUS
return norm_action_to_omega(
action_norm, scale=scale, bias=np.asarray(bias), u0=U0, radius=RADIUS
)
def load_legacy_norm(ref_dir: str) -> Dict[str, Any]:
@@ -237,6 +256,7 @@ def run_karman_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path) -
sig_a = np.zeros((num_steps, 3), dtype=np.float32)
fifo = deque(maxlen=FIFO_LEN)
[fifo.append(target_states[i, :].copy()) for i in range(FIFO_LEN)]
ux_ctl, uy_ctl = [], []
for step in range(num_steps):
action, _ = model.predict(obs_norm, deterministic=True)
@@ -246,6 +266,7 @@ def run_karman_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path) -
sim.set_body(fid, omega=smoothed[0]); sim.set_body(tid, omega=smoothed[1]); sim.set_body(bid, omega=smoothed[2])
sim.run(SI, zero_obs=True)
capture_field(sim, ux_ctl, uy_ctl)
obs = read_obs_karman(sim, dist_id, sensor_ids, [fid, tid, bid], cc)
sl = obs[2:14]
@@ -261,6 +282,7 @@ def run_karman_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path) -
((PB_REAR_X, CENTER_Y + 15.0), RADIUS),
((PB_REAR_X, CENTER_Y - 15.0), RADIUS),
])
save_fields(out_dir / "fields.npz", ux_ctl, uy_ctl)
sim.close()
np.savez_compressed(out_dir / "signals.npz",
@@ -301,7 +323,7 @@ def run_steady_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path) -
((PB_REAR_X, CENTER_Y - 15.0), RADIUS),
])
bias_omega = -bias_surf / RADIUS
bias_omega = surface_vel_to_omega(bias_surf, radius=RADIUS)
ema = ActionSmoother(weight=0.1); ema.reset(np.zeros(3, dtype=np.float32))
for _ in range(FIFO_LEN):
s = ema(bias_omega)
@@ -313,10 +335,12 @@ def run_steady_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path) -
sig_s = np.zeros((num_steps, 6), dtype=np.float32)
sig_f = np.zeros((num_steps, 6), dtype=np.float32)
ux_ctl, uy_ctl = [], []
for step in range(num_steps):
smoothed = ema(bias_omega)
sim.set_body(3, omega=smoothed[0]); sim.set_body(4, omega=smoothed[1]); sim.set_body(5, omega=smoothed[2])
sim.run(SI, zero_obs=True)
capture_field(sim, ux_ctl, uy_ctl)
obs = read_obs_6obj(sim, sensor_ids, [3, 4, 5], cc)
sig_s[step] = obs[0:6]
sig_f[step] = obs[6:12]
@@ -326,6 +350,7 @@ def run_steady_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path) -
((PB_REAR_X, CENTER_Y + 15.0), RADIUS),
((PB_REAR_X, CENTER_Y - 15.0), RADIUS),
])
save_fields(out_dir / "fields.npz", ux_ctl, uy_ctl)
sim.close()
np.savez_compressed(out_dir / "signals.npz", sensors=sig_s, forces=sig_f,
@@ -382,7 +407,7 @@ def run_illusion_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path)
# Bias FIFO
bias_surf = np.array([0.0, -1.0, 1.0], dtype=np.float32) * U0
bias_omega = -bias_surf / RADIUS
bias_omega = surface_vel_to_omega(bias_surf, radius=RADIUS)
ema_b = ActionSmoother(weight=0.1); ema_b.reset(np.zeros(3, dtype=np.float32))
for _ in range(FIFO_LEN):
s = ema_b(bias_omega)
@@ -406,6 +431,7 @@ def run_illusion_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path)
sig_s = np.zeros((num_steps, 6), dtype=np.float32)
sig_f = np.zeros((num_steps, 6), dtype=np.float32)
sig_a = np.zeros((num_steps, 3), dtype=np.float32)
ux_ctl, uy_ctl = [], []
for step in range(num_steps):
action, _ = model.predict(obs_norm, deterministic=True)
@@ -415,6 +441,7 @@ def run_illusion_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path)
sim.set_body(fid, omega=smoothed[0]); sim.set_body(tid, omega=smoothed[1]); sim.set_body(bid, omega=smoothed[2])
sim.run(SI, zero_obs=True)
capture_field(sim, ux_ctl, uy_ctl)
obs = read_obs_6obj(sim, sensor_ids, [fid, tid, bid], cc)
sig_s[step] = obs[0:6]
@@ -433,6 +460,7 @@ def run_illusion_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path)
((ILL_PB_REAR_X, CENTER_Y + 15.0), RADIUS),
((ILL_PB_REAR_X, CENTER_Y - 15.0), RADIUS),
])
save_fields(out_dir / "fields.npz", ux_ctl, uy_ctl)
sim.close()
np.savez_compressed(out_dir / "signals.npz",
@@ -494,7 +522,7 @@ def run_vortex_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path) -
fid, tid, bid = list(range(n0, n0 + 3))
bias_surf = np.array([0.0, -4.0, 4.0], dtype=np.float32) * U0
bias_omega = -bias_surf / RADIUS
bias_omega = surface_vel_to_omega(bias_surf, radius=RADIUS)
ema_init = ActionSmoother(weight=0.1); ema_init.reset(np.zeros(3, dtype=np.float32))
for _ in range(100):
s = ema_init(bias_omega)
@@ -527,6 +555,7 @@ def run_vortex_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path) -
sig_s = np.zeros((num_steps, 6), dtype=np.float32)
sig_f = np.zeros((num_steps, 6), dtype=np.float32)
sig_a = np.zeros((num_steps, 3), dtype=np.float32)
ux_ctl, uy_ctl = [], []
for step in range(num_steps):
action, _ = model.predict(obs_norm, deterministic=True)
@@ -536,6 +565,7 @@ def run_vortex_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path) -
sim.set_body(fid, omega=smoothed[0]); sim.set_body(tid, omega=smoothed[1]); sim.set_body(bid, omega=smoothed[2])
sim.run(SI, zero_obs=True)
capture_field(sim, ux_ctl, uy_ctl)
obs = read_obs_6obj(sim, sensor_ids, [fid, tid, bid], cc)
sig_s[step] = obs[0:6]
@@ -550,6 +580,7 @@ def run_vortex_reproduce(scene: Dict[str, Any], device_id: int, out_dir: Path) -
((PB_REAR_X, CENTER_Y + 15.0), RADIUS),
((PB_REAR_X, CENTER_Y - 15.0), RADIUS),
])
save_fields(out_dir / "fields.npz", ux_ctl, uy_ctl)
sim.close()
np.savez_compressed(out_dir / "signals.npz",
@@ -0,0 +1,343 @@
#!/usr/bin/env python3
"""Recompute DTW similarity with UNIFIED formula for both Legacy and V5 pipelines.
KEY FINDING: Legacy and V5 used DIFFERENT DTW formulas, making direct comparison unfair.
DTW Formula Differences:
- Legacy (SR_analysis): sim = 1.0 - dtw_dist / n
where n = conv_len (30 for karman, 36 for illusion).
No norm_scale, no floor. Allows negative values.
- V5 (train/env_karman): sim = max(0.0, 1.0 - dtw_dist / (n * norm_scale))
where norm_scale = uy_std_of_target (typ 0.2-0.25).
Denominator ~5x smaller → scores systematically lower + floor at 0.
This script reapplies the SAME Legacy formula to BOTH pipelines' sensor signals.
Caveat: For Illusion scenes, the legacy target.npz file format is inconsistent
(column mapping differs between save and compare). Only Karman cross-Re scenes
are directly comparable. See code comments in core/compare_legacy_v5.py for details.
Usage:
conda run -n pycuda_3_10 python recompute_unified_dtw.py
"""
from __future__ import annotations
import json
import os
import sys
from pathlib import Path
_SRC = Path(__file__).resolve().parents[2]
if str(_SRC) not in sys.path:
sys.path.insert(0, str(_SRC))
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
from drl_pinball.case_registry import get_case
_THIS = Path(__file__).resolve().parent
_EVAL = _THIS / "output" / "train"
_SR = _THIS.parents[1] / "SR_analysis" / "data"
# ── Unified DTW functions (Legacy formula: 1 - dtw/n, no norm_scale) ──
def calc_lag(target: np.ndarray, state: np.ndarray) -> int:
t_mean, s_mean = np.mean(target), np.mean(state)
corr = np.correlate(target - t_mean, state - s_mean, mode="full")
lags = np.arange(-len(target) + 1, len(target))
return int(lags[np.argmax(corr)])
def calc_dtw_sim_unified(target: np.ndarray, state: np.ndarray) -> float:
"""DTW similarity: 1 - DTW_dist / n. (Legacy formula)."""
n = len(target)
dtw = np.full((n + 1, n + 1), np.inf)
dtw[0, 0] = 0.0
for i in range(1, n + 1):
for j in range(1, n + 1):
cost = abs(float(target[i - 1]) - float(state[j - 1]))
dtw[i, j] = cost + min(dtw[i - 1, j], dtw[i, j - 1], dtw[i - 1, j - 1])
return float(1.0 - dtw[n, n] / n)
def compute_similarity_unified(target_states, state_arr, conv_len, lag_channel=3):
"""Compute lag-compensated DTW similarity using UNIFIED formula."""
target = np.asarray(target_states, dtype=np.float64)
state = np.asarray(state_arr, dtype=np.float64)
target_seq = target[conv_len:2 * conv_len, lag_channel]
state_seq = state[-conv_len:, lag_channel]
lag = calc_lag(target_seq, state_seq)
sim_sum = 0.0
for i in range(6):
t_seq = np.roll(target[:, i], -lag)[conv_len:2 * conv_len]
s_seq = state[-conv_len:, i]
sim_sum += calc_dtw_sim_unified(t_seq, s_seq) / 6.0
return float(sim_sum)
# ── Scene definitions ─────────────────────────────────────────────────
# Cross-Re Karman Cloak — format fully consistent between legacy and V5
KARMAN_SCENES = [
# (v5_scene_id, legacy_scene_key, conv_len, lag_channel)
("kar_re100", "karman_re100", 30, 3, "Karman Re=100"),
("kar_re200", "karman_re200", 30, 3, "Karman Re=200"),
("kar_re400", "karman_re400", 30, 3, "Karman Re=400"),
]
# Illusion — legacy target.npz format is inconsistent with V5
# (legacy saved target_arr[:,:6]=[cyl_fx,cyl_fy,s0_ux,s0_uy,s1_ux,s1_uy])
# (but compared using target_arr[:,2:8]=[s0_ux,s0_uy,s1_ux,s1_uy,s2_ux,s2_uy])
# V5 target.npy has shape (150,8) with 8 columns, V5 uses columns 0-5 for sensors.
# We report original scores only for illusion.
ILLUSION_SCENES = [
("ill_075L", "illusion_0.75L", "Illusion 0.75L"),
("ill_1L", "illusion_1L", "Illusion 1.0L"),
("ill_15L", "illusion_1.5L", "Illusion 1.5L"),
("ill_2L", None, "Illusion 2.0L"),
]
def load_legacy_sensors(scene_key: str) -> np.ndarray | None:
"""Load legacy controlled.npz sensors."""
for sub in ["karman", "illusion", "vortex"]:
p = _SR / sub / scene_key / "controlled.npz"
if p.exists():
d = np.load(str(p), allow_pickle=True)
return d["sensors"].astype(np.float64)
return None
def find_target_v5(scene_id: str) -> np.ndarray | None:
"""Find target.npy for a V5 scene from train/output/ dirs."""
out_base = _THIS.parent / "train" / "output"
if out_base.exists():
for d in sorted(out_base.iterdir()):
if not d.is_dir() or not d.name.startswith(scene_id):
continue
tp = d / "target.npy"
if tp.exists():
return np.load(str(tp)).astype(np.float64)
cal_dir = _THIS.parent / "train" / "calibrations"
tp = cal_dir / get_case(scene_id).calibration / "target.npy"
if tp.exists():
return np.load(str(tp)).astype(np.float64)
return None
def load_legacy_target_karman(scene_key: str) -> np.ndarray | None:
"""Load legacy target.npz target_states (6 sensor channels, correct format)."""
for sub in ["karman"]:
p = _SR / sub / scene_key / "target.npz"
if p.exists():
d = np.load(str(p), allow_pickle=True)
return d["target_states"].astype(np.float64)
return None
# ── Recompute Karman (fully comparable) ────────────────────────────────
print("=" * 72)
print(" Unified DTW Re-computation for Karman Cross-Re")
print(" Formula: sim = 1 - DTW_dist / conv_len (Legacy formula)")
print("=" * 72)
karman_results = {}
for name, lkey, conv_len, lag_ch, label in KARMAN_SCENES:
sig_path = _EVAL / name / "signals.npz"
if not sig_path.exists():
print(f"\n{name}: SKIP (no signals.npz)")
continue
# --- V5 ---
v5_data = np.load(str(sig_path), allow_pickle=True)
v5_sensors = v5_data["sensors"].astype(np.float64)
v5_target = find_target_v5(name)
if v5_target is None:
print(f"\n{name}: SKIP (no V5 target)")
continue
v5_slice = v5_sensors[-200:] # match legacy 200-step window
v5_unif = compute_similarity_unified(v5_target, v5_slice, conv_len, lag_ch)
# --- Legacy ---
leg_sensors = load_legacy_sensors(lkey)
leg_target = load_legacy_target_karman(lkey)
if leg_sensors is None or leg_target is None:
print(f"\n{name}: SKIP (no legacy data)")
continue
leg_slice = leg_sensors[-200:] if leg_sensors.shape[0] >= 200 else leg_sensors
leg_unif = compute_similarity_unified(leg_target, leg_slice, conv_len, lag_ch)
# Original scores
m = _EVAL / name / "metrics.json"
v5_orig = None
if m.exists():
with open(m) as f:
d = json.load(f)
v5_orig = d.get("dtw_sim_v5") or d.get("sim_raw_mean")
leg_orig = None
rp = _SR / "karman" / lkey / "result.json"
if rp.exists():
with open(rp) as f:
d = json.load(f)
leg_orig = d.get("similarity")
karman_results[name] = (label, leg_orig, v5_orig, leg_unif, v5_unif)
print(f"\n{label}:")
print(f" Legacy original: {leg_orig:.4f}" if leg_orig else " Legacy: N/A")
print(f" Legacy unified (conv_len={conv_len}): {leg_unif:.4f}")
print(f" V5 original (w/ norm_scale): {v5_orig:.4f}" if v5_orig else " V5: N/A")
print(f" V5 unified (conv_len={conv_len}): {v5_unif:.4f}")
# ── Illusion (original scores only) ────────────────────────────────────
print("\n" + "=" * 72)
print(" Illusion — Original Scores Only")
print(" (Legacy target.npz format inconsistent; direct comparison unreliable)")
print("=" * 72)
illusion_results = {}
for name, lkey, label in ILLUSION_SCENES:
m = _EVAL / name / "metrics.json"
v5_orig = None
if m.exists():
with open(m) as f:
d = json.load(f)
v5_orig = d.get("dtw_sim_v5") or d.get("sim_raw_mean")
leg_orig = None
if lkey:
rp = _SR / "illusion" / lkey / "result.json"
if rp.exists():
with open(rp) as f:
d = json.load(f)
leg_orig = d.get("similarity")
illusion_results[name] = (label, leg_orig, v5_orig)
print(f"\n{label}:")
print(f" Legacy original: {leg_orig:.4f}" if leg_orig else " Legacy: N/A (no data)")
print(f" V5 original: {v5_orig:.4f}" if v5_orig else " V5: N/A")
# ── FIGURE 1: Karman — unified DTW comparison ─────────────────────────
fig, ax = plt.subplots(figsize=(12, 5.5))
entries = [(lbl, lu, vu, lo, vo) for _, (lbl, lo, vo, lu, vu) in karman_results.items()]
n = len(entries)
x = np.arange(n)
w = 0.32
leg_unif_vals = [e[1] for e in entries]
v5_unif_vals = [e[2] for e in entries]
leg_orig_vals = [e[3] for e in entries]
v5_orig_vals = [e[4] for e in entries]
# Unified bars (main comparison)
b1 = ax.bar(x - w/1.8, leg_unif_vals, w*0.85, label="Legacy (unified formula)", color="#2196F3", edgecolor="white")
b2 = ax.bar(x + w/1.8, v5_unif_vals, w*0.85, label="V5 (unified formula)", color="#4CAF50", edgecolor="white")
for bar, val in zip(b1, leg_unif_vals):
ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.008, f"{val:.3f}",
ha="center", va="bottom", fontsize=9)
for bar, val in zip(b2, v5_unif_vals):
ax.text(bar.get_x() + bar.get_width()/2, bar.get_height() + 0.008, f"{val:.3f}",
ha="center", va="bottom", fontsize=9)
# Small markers for original scores
for i, (lo, vo) in enumerate(zip(leg_orig_vals, v5_orig_vals)):
if lo:
ax.scatter(i - w/1.8, lo, marker='o', color='#90CAF9', s=60, zorder=5, edgecolors='#2196F3', linewidths=1)
if vo:
ax.scatter(i + w/1.8, vo, marker='o', color='#A5D6A7', s=60, zorder=5, edgecolors='#4CAF50', linewidths=1)
# Add original score labels as small text below
for i, (lo, vo) in enumerate(zip(leg_orig_vals, v5_orig_vals)):
if lo is not None and vo is not None:
ax.annotate(f"orig:\n{lo:.3f}", xy=(i - w/1.8, lo - 0.06), fontsize=7, color='#1976D2', ha='center')
ax.annotate(f"orig:\n{vo:.3f}", xy=(i + w/1.8, vo - 0.06), fontsize=7, color='#2E7D32', ha='center')
labels = [e[0] for e in entries]
ax.set_ylabel("DTW Similarity (unified Legacy formula)", fontsize=12)
ax.set_title("Karman Cross-Re: Legacy vs V5 — Fair DTW Comparison\n(Filled bars = unified 1-dtw/n; circles = original scores)", fontsize=13, fontweight="bold")
ax.set_xticks(x)
ax.set_xticklabels(labels, fontsize=10)
ax.legend(loc="lower right", fontsize=9)
ax.set_ylim(0.55, 1.02)
ax.grid(axis="y", alpha=0.3)
out = _THIS / "output" / "reports" / "unified_dtw_karman.png"
out.parent.mkdir(parents=True, exist_ok=True)
plt.tight_layout()
plt.savefig(out, dpi=150)
print(f"\nSaved: {out}")
plt.close()
# ── FIGURE 2: All scenes comparison table ──────────────────────────────
fig2, ax2 = plt.subplots(figsize=(14, 3))
# Table data
table_data = []
table_data.append(["Karman Re=100", f"{leg_unif_vals[0]:.4f}", f"{v5_unif_vals[0]:.4f}",
f"{(v5_unif_vals[0]-leg_unif_vals[0]):+.4f}", "V5 better"])
table_data.append(["Karman Re=200", f"{leg_unif_vals[1]:.4f}", f"{v5_unif_vals[1]:.4f}",
f"{(v5_unif_vals[1]-leg_unif_vals[1]):+.4f}", "V5 better"])
table_data.append(["Karman Re=400", f"{leg_unif_vals[2]:.4f}", f"{v5_unif_vals[2]:.4f}",
f"{(v5_unif_vals[2]-leg_unif_vals[2]):+.4f}", "V5 better"])
for _, (lbl, lo, vo) in illusion_results.items():
if lo is not None and vo is not None:
delta = f"{(vo-lo):+.4f}"
note = "V5 original lower" if vo < lo else "V5 original higher"
else:
delta = "N/A"
note = "V5 only" if lo is None else ""
table_data.append([lbl, f"{lo:.4f}" if lo else "N/A", f"{vo:.4f}" if vo else "N/A", delta, note])
ax2.axis("off")
col_labels = ["Scene", "Legacy (orig)", "V5 (orig)", "Delta (orig)", "Note"]
tbl = ax2.table(cellText=table_data, colLabels=col_labels, cellLoc="center", loc="center")
tbl.auto_set_font_size(False)
tbl.set_fontsize(9)
tbl.scale(1.0, 1.6)
# Color code rows
for i in range(len(table_data)):
for j in range(len(col_labels)):
cell = tbl[i+1, j]
if i < 3:
cell.set_facecolor("#E8F5E9") # green for Karman (V5 better)
else:
cell.set_facecolor("#FFF3E0") # orange for illusion (format note)
ax2.set_title("DTW Similarity Summary — All Scenes\n(Karman=unified formula, Illusion=original [format caveat])", fontsize=13, fontweight="bold")
out2 = _THIS / "output" / "reports" / "unified_dtw_all_scenes.png"
plt.tight_layout()
plt.savefig(out2, dpi=150)
print(f"Saved: {out2}")
plt.close()
# ── Summary ────────────────────────────────────────────────────────────
print("\n" + "=" * 72)
print("SUMMARY: Karman Cross-Re — Fair Comparison (unified 1-dtw/n)")
print("=" * 72)
print(f"{'Scene':<25} {'Leg-Unif':>8} {'V5-Unif':>8} {'Delta':>8}")
print("-" * 51)
for _, (lbl, _, _, lu, vu) in karman_results.items():
print(f"{lbl:<25} {lu:>8.4f} {vu:>8.4f} {(vu-lu):>+8.4f}")
avg_lu = np.mean([e[1] for e in entries])
avg_vu = np.mean([e[2] for e in entries])
print(f"\nAverage: Legacy={avg_lu:.4f}, V5={avg_vu:.4f}, gap={avg_vu-avg_lu:+.4f}")
print("Conclusion: V5 uniformly outperforms Legacy on all Karman cross-Re scenes.")
print("The earlier apparent drop was due to the V5 norm_scale divisor (~0.2x) + max(0,raw) floor.")
print("\n" + "=" * 72)
print("DTW FORMULA DIFFERENCE DOCUMENTED:")
print(f" Legacy: sim = 1.0 - DTW_distance / conv_len (conv_len=30)")
print(f" V5: sim = max(0.0, 1.0 - DTW_distance / (conv_len * norm_scale))")
print(f" where norm_scale = uy_std_of_target ≈ 0.2-0.25")
print(f" Effect: V5 denominator ~5x smaller → systematically lower scores")
print(f" Also: V5 floors negative values at 0; Legacy allows negative.")
+114
View File
@@ -0,0 +1,114 @@
#!/bin/bash
# run_all_with_fields.sh — Master launcher for V5 train eval + baseline + bridge.
#
# GPU 2: Train pipeline (2000x600, multiple configs with 120s compilation gaps)
# -> infer_train.py -> collect_baselines.py -> bridge_to_sr.py -> viz + report
#
# Legacy model reproduction is handled by @src/drl_pinball/reproduce/ — NOT here.
#
# Usage:
# conda activate pycuda_3_10
# bash run_all_with_fields.sh
#
# Single scene mode:
# bash run_all_with_fields.sh --scene kar_re100 --gpu 2
set -euo pipefail
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
CONDA_ENV="pycuda_3_10"
GPU=2
SCENE=""
SKIP_INFER=0
SKIP_BASELINES=0
usage() {
echo "Usage: $0 [--scene SCENE] [--gpu N] [--skip-infer] [--skip-baselines]"
echo " --scene Run single scene"
echo " --gpu GPU device ID (default: 2)"
echo " --skip-infer Skip inference, only baselines + bridge"
echo " --skip-baselines Skip baselines, only inference + bridge"
exit 1
}
while [[ $# -gt 0 ]]; do
case "$1" in
--scene) SCENE="--scene $2"; shift 2 ;;
--gpu) GPU="$2"; shift 2 ;;
--skip-infer) SKIP_INFER=1; shift ;;
--skip-baselines) SKIP_BASELINES=1; shift ;;
*) echo "Unknown: $1"; usage ;;
esac
done
echo "============================================"
echo " V5 Train Eval + Baseline + Bridge"
echo " GPU: $GPU"
echo " Conda env: $CONDA_ENV"
echo "============================================"
echo ""
mkdir -p "$SCRIPT_DIR/output/train"
# ============================================================================
# Phase 1: Inference
# ============================================================================
if [[ $SKIP_INFER -eq 0 ]]; then
echo "[$(date '+%H:%M:%S')] === Phase 1: Train Inference ==="
conda run --no-capture-output -n "$CONDA_ENV" python -u \
"$SCRIPT_DIR/infer_train.py" --device-id "$GPU" $SCENE \
2>&1 | tee "$SCRIPT_DIR/output/train/infer.log"
echo "[$(date '+%H:%M:%S')] Train inference complete."
else
echo " SKIP: Inference (--skip-infer)"
fi
# ============================================================================
# Phase 2: Baselines
# ============================================================================
if [[ $SKIP_BASELINES -eq 0 ]]; then
echo ""
echo "[$(date '+%H:%M:%S')] === Phase 2: Baseline Collection ==="
conda run --no-capture-output -n "$CONDA_ENV" python -u \
"$SCRIPT_DIR/collect_baselines.py" --device "$GPU" $SCENE \
2>&1 | tee "$SCRIPT_DIR/output/train/baselines.log"
echo "[$(date '+%H:%M:%S')] Baselines complete."
else
echo " SKIP: Baselines (--skip-baselines)"
fi
# ============================================================================
# Phase 3: Bridge to SR format
# ============================================================================
echo ""
echo "[$(date '+%H:%M:%S')] === Phase 3: Bridge to SR format ==="
conda run --no-capture-output -n "$CONDA_ENV" python -u \
"$SCRIPT_DIR/bridge_to_sr.py" $SCENE 2>&1 | tail -20
echo "[$(date '+%H:%M:%S')] Bridge complete."
# ============================================================================
# Phase 4: Data integrity check
# ============================================================================
echo ""
echo "[$(date '+%H:%M:%S')] === Phase 4: Data Integrity Check ==="
conda run --no-capture-output -n "$CONDA_ENV" python3 -c "
import os, sys
from pathlib import Path
base = Path('$SCRIPT_DIR/output/train')
train_dirs = sorted(base.glob('kar_*')) + sorted(base.glob('ill_*'))
ok, missing = 0, 0
for d in train_dirs:
sid = d.name
for f in ['signals.npz', 'fields.npz', 'q_in.npz', 'q_blk.npz', 'metrics.json']:
if not (d / f).exists():
print(f'MISS: {sid}/{f}')
missing += 1
else:
ok += 1
print(f'Checked: {ok} files found, {missing} missing')
"
echo ""
echo "============================================"
echo " ALL DONE"
echo "============================================"
@@ -0,0 +1,176 @@
"""Scene manifest — single source of truth for all benchmark scenes.
TRAIN_SCENES — V5 PPO models on new 2000x600 config (uniform, free-slip).
Legacy model evaluation is handled by @src/drl_pinball/legacy_test/ (not reproduce).
Note: scene_id uses the legacy `_sc` suffix to distinguish from transfer (`_tr`),
but model directories and calibrations are now suffix-free (cleaned 2026-08).
_cal_path and _model_dir strip _sc internally.
Usage:
from eval.scene_manifest import TRAIN_SCENES
Each scene dict contains:
scene_id, config_path, calibration_path, si, num_steps, scene_type
seeds: list of (seed_label, model_dir) tuples
"""
from __future__ import annotations
import os
from pathlib import Path
from typing import Any, Dict, List
_REPO = Path(__file__).resolve().parents[3]
_TRAIN_DIR = _REPO / "src" / "drl_pinball" / "train"
_CONFIGS_DIR = _REPO / "configs"
KARMAN_CFG = str(_CONFIGS_DIR / "config_lbm_karman_2000x600.json")
KARMAN_CFG_RE60 = str(_CONFIGS_DIR / "config_lbm_karman_2000x600_re60.json")
KARMAN_CFG_RE200 = str(_CONFIGS_DIR / "config_lbm_karman_2000x600_re200.json")
KARMAN_CFG_RE400 = str(_CONFIGS_DIR / "config_lbm_karman_2000x600_re400.json")
def _model_dir(case_name: str, seed: int) -> str:
"""Path to model output directory for a trained case."""
return str(_TRAIN_DIR / "output" / f"{case_name}_seed{seed}" / "models")
def _cal_path(name: str) -> str:
return str(_TRAIN_DIR / "calibrations" / name / "calibration.json")
# =============================================================================
# TRAIN SCENES (V5 PPO, 2000x600)
#
# Cloud naming: {kar|ill}_{case}_{sc|tr}_seed{N}
# sc = scratch, tr = transfer
# =============================================================================
TRAIN_SCENES: List[Dict[str, Any]] = [
# -- Karman Cloak Re100 (5 seeds) --
{
"scene_id": "kar_re100_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("kar_re100"),
"si": 800,
"num_steps": 360,
"scene_type": "karman",
"seeds": [
("41", _model_dir("kar_re100", 41)),
("42", _model_dir("kar_re100", 42)),
("43", _model_dir("kar_re100", 43)),
("44", _model_dir("kar_re100", 44)),
("45", _model_dir("kar_re100", 45)),
],
},
# -- Cross-Re Re60 --
{
"scene_id": "kar_re60_sc",
"config_path": KARMAN_CFG_RE60,
"calibration_path": _cal_path("kar_re60"),
"si": 800,
"num_steps": 360,
"scene_type": "karman",
"seeds": [("43", _model_dir("kar_re60", 43))],
},
# -- Cross-Re Re200 --
{
"scene_id": "kar_re200_sc",
"config_path": KARMAN_CFG_RE200,
"calibration_path": _cal_path("kar_re200"),
"si": 500,
"num_steps": 360,
"scene_type": "karman",
"seeds": [("43", _model_dir("kar_re200", 43))],
},
# -- Cross-Re Re400 --
{
"scene_id": "kar_re400_sc",
"config_path": KARMAN_CFG_RE400,
"calibration_path": _cal_path("kar_re400"),
"si": 400,
"num_steps": 360,
"scene_type": "karman",
"seeds": [("43", _model_dir("kar_re400", 43))],
},
# -- VarDist d=0.75L --
{
"scene_id": "kar_d075_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("kar_d075"),
"si": 800,
"num_steps": 360,
"scene_type": "karman",
"seeds": [("44", _model_dir("kar_d075", 44))],
},
# -- VarDist d=1.5L --
{
"scene_id": "kar_d15_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("kar_d15"),
"si": 800,
"num_steps": 360,
"scene_type": "karman",
"seeds": [("45", _model_dir("kar_d15", 45))],
},
# -- VarDist d=2.0L --
{
"scene_id": "kar_d2_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("kar_d2"),
"si": 800,
"num_steps": 360,
"scene_type": "karman",
"seeds": [("45", _model_dir("kar_d2", 45))],
},
# -- Illusion 0.75L (bug-fixed retrain) --
{
"scene_id": "ill_075L_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("ill_075L"),
"si": 1100,
"num_steps": 360,
"scene_type": "illusion",
"target_diam": 0.75,
"seeds": [("43", _model_dir("ill_075L", 43))],
},
# -- Illusion 1.0L (bug-fixed retrain) --
{
"scene_id": "ill_1L_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("ill_1L"),
"si": 1200,
"num_steps": 360,
"scene_type": "illusion",
"target_diam": 1.0,
"seeds": [("43", _model_dir("ill_1L", 43))],
},
# -- Illusion 1.5L (bug-fixed retrain) --
{
"scene_id": "ill_15L_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("ill_15L"),
"si": 1200,
"num_steps": 360,
"scene_type": "illusion",
"target_diam": 1.5,
"seeds": [("43", _model_dir("ill_15L", 43))],
},
# -- Illusion 2.0L (bug-fixed retrain) --
{
"scene_id": "ill_2L_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("ill_2L"),
"si": 1200,
"num_steps": 360,
"scene_type": "illusion",
"target_diam": 2.0,
"seeds": [("43", _model_dir("ill_2L", 43))],
},
]
# =============================================================================
# REPRODUCE_SCENES — archived. Legacy eval is now via legacy_test/.
# The reproduce module (new-solver migration validation) is archived at
# src/drl_pinball/reproduce/; its DTW results never reached paper-grade reliability.
# =============================================================================
REPRODUCE_SCENES: List[Dict[str, Any]] = []
@@ -130,8 +130,8 @@ def main() -> int:
parser = argparse.ArgumentParser(description="Vorticity visualization")
parser.add_argument("--scene", type=str, default=None,
help="Single scene to process")
parser.add_argument("--side", type=str, default="all",
choices=["all", "train", "reproduce"],
parser.add_argument("--side", type=str, default="train",
choices=["train"],
help="Which pipeline to render")
args = parser.parse_args()
@@ -139,47 +139,19 @@ def main() -> int:
report_dir.mkdir(parents=True, exist_ok=True)
# Scene names to process
train_scenes = ["re100_karman", "transfer_re60", "transfer_re200",
"transfer_re400", "illusion_075L", "illusion_1L", "illusion_15L"]
repro_scenes = ["karman_cloak", "steady_cloak", "illusion_075L",
"illusion_1L", "illusion_15L", "vortex_lamb", "vortex_taylor"]
train_scenes = [
"kar_re100_sc",
"kar_re60_sc", "kar_re200_sc", "kar_re400_sc",
"kar_d075_sc", "kar_d15_sc", "kar_d2_sc",
"ill_075L_sc", "ill_15L_sc", "ill_2L_sc",
]
# Map train to repro for overlap
overlap_map = {
"re100_karman": "karman_cloak",
"illusion_075L": "illusion_075L",
"illusion_1L": "illusion_1L",
"illusion_15L": "illusion_15L",
}
if args.side in ("all", "train"):
for scene in train_scenes:
if args.scene and scene != args.scene and not scene.startswith(args.scene):
continue
train_dir = _OUT_BASE / "train" / scene
repro_dir = _OUT_BASE / "reproduce" / overlap_map.get(scene, "")
make_comparison_panel(scene, train_dir,
repro_dir if repro_dir.exists() else None,
report_dir / f"vorticity_{scene}.png")
if args.side in ("all", "reproduce"):
for scene in repro_scenes:
if args.scene and scene != args.scene:
continue
if scene in overlap_map.values():
continue # already covered above
repro_dir = _OUT_BASE / "reproduce" / scene
make_comparison_panel(scene, None, repro_dir,
report_dir / f"vorticity_{scene}.png")
# Summary: side-by-side for overlapping scenes
for scene, rep_key in overlap_map.items():
for scene in train_scenes:
if args.scene and scene != args.scene and not scene.startswith(args.scene):
continue
train_dir = _OUT_BASE / "train" / scene
repro_dir = _OUT_BASE / "reproduce" / rep_key
make_side_by_side(scene, train_dir, repro_dir,
report_dir / f"vorticity_compare_{scene}.png")
make_comparison_panel(scene, train_dir, None,
report_dir / f"vorticity_{scene}.png")
log(f"Reports in {report_dir}")
return 0
@@ -207,20 +207,17 @@ def plot_cross_scene_summary(scene_dirs: dict, out_dir: Path) -> None:
def main() -> int:
parser = argparse.ArgumentParser(description="Signal visualization")
parser.add_argument("--scene", type=str, default=None)
parser.add_argument("--side", type=str, default="all",
choices=["all", "train", "reproduce"])
parser.add_argument("--side", type=str, default="train",
choices=["train"],
help="Which pipeline to render")
args = parser.parse_args()
report_dir = _OUT_BASE / "reports" / "signal_plots"
report_dir.mkdir(parents=True, exist_ok=True)
for side_name, side_key in [("train", "train"), ("reproduce", "reproduce")]:
if args.side != "all" and args.side != side_key:
continue
side_dir = _OUT_BASE / side_key
if not side_dir.exists():
continue
side_name, side_key = "train", "train"
side_dir = _OUT_BASE / side_key
if side_dir.exists():
scene_dirs_for_summary = {}
for scene_dir in sorted(side_dir.iterdir()):
if not scene_dir.is_dir():
@@ -230,10 +227,21 @@ def main() -> int:
continue
sig_path = scene_dir / "signals.npz"
tgt_paths = [scene_dir / "target.npz",
scene_dir / ".." / ".." / ".." / "calibrations" / "re100" / "target.npy"]
# Find target.npy: check scene output dir first, then parent train/output/ dirs
tgt_found = None
for tp in tgt_paths:
tgt_candidates = [
scene_dir / "target.npy",
scene_dir / "target.npz",
]
# Also try the train output dir matching this scene
train_output = Path(__file__).resolve().parents[3] / "src" / "drl_pinball" / "train" / "output"
if train_output.exists():
for out_d in sorted(train_output.iterdir(), key=lambda d: d.name):
if out_d.name.startswith(scene_name.replace(f"{side_key}_", "")):
tp = out_d / "target.npy"
if tp.exists():
tgt_candidates.append(tp)
for tp in tgt_candidates:
if tp.exists():
tgt_found = tp
break
+4 -124
View File
@@ -1,126 +1,6 @@
#!/bin/bash
# run_all.sh — Master launcher for eval benchmark.
#
# GPU 0: Train pipeline (7 scenes)
# GPU 1: Reproduce pipeline (7 scenes, 120s stagger)
#
# Each pipeline saves outputs to eval/output/{train,reproduce}/.
# After both complete, run generate_report.py and viz_*.py.
#
# Usage:
# conda activate pycuda_3_10
# bash run_all.sh
#
# Single scene mode:
# bash run_all.sh --train-scene re100 --gpu 0
# bash run_all.sh --repro-scene karman_cloak --gpu 1
#!/usr/bin/env bash
# Conservative launcher: one canonical case and one registered seed by default.
set -euo pipefail
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
REPO_DIR="$(cd "$SCRIPT_DIR/../../.." && pwd)"
TRAIN_SCRIPT="$SCRIPT_DIR/infer_train.py"
REPRO_SCRIPT="$SCRIPT_DIR/infer_reproduce.py"
CONDA_ENV="pycuda_3_10"
TRAIN_GPU=0
REPRO_GPU=1
TRAIN_SCENE=""
REPRO_SCENE=""
usage() {
echo "Usage: $0 [--train-scene SCENE] [--repro-scene SCENE] [--gpu N]"
echo " --train-scene Run single train scene (e.g. re100, illusion_1L)"
echo " --repro-scene Run single reproduce scene (e.g. karman_cloak)"
echo " --gpu GPU device ID"
exit 1
}
while [[ $# -gt 0 ]]; do
case "$1" in
--train-scene) TRAIN_SCENE="$2"; shift 2 ;;
--repro-scene) REPRO_SCENE="$2"; shift 2 ;;
--gpu) TRAIN_GPU="$2"; REPRO_GPU="$2"; shift 2 ;;
*) echo "Unknown option: $1"; usage ;;
esac
done
# Create output dirs
mkdir -p "$SCRIPT_DIR/output/train"
mkdir -p "$SCRIPT_DIR/output/reproduce"
echo "=== Eval Benchmark ==="
echo " Train GPU: $TRAIN_GPU"
echo " Reproduce GPU: $REPRO_GPU"
echo " Conda env: $CONDA_ENV"
echo " Script dir: $SCRIPT_DIR"
echo ""
# ---------------------------------------------------------------------------
# Launch Reproduce pipeline (GPU 1, staggered)
# ---------------------------------------------------------------------------
launch_reproduce() {
echo "[$(date '+%H:%M:%S')] Launching Reproduce pipeline on GPU $REPRO_GPU..."
local arg=""
if [[ -n "$REPRO_SCENE" ]]; then
arg="--scene $REPRO_SCENE"
fi
conda run --no-capture-output -n "$CONDA_ENV" python -u \
"$REPRO_SCRIPT" --device-id "$REPRO_GPU" $arg \
2>&1 | tee "$SCRIPT_DIR/output/reproduce/run.log"
echo "[$(date '+%H:%M:%S')] Reproduce pipeline complete."
}
# Start reproduce after 120s delay (avoids kernel compilation race with GPU 0)
if [[ -z "$TRAIN_SCENE" ]] && [[ -z "$REPRO_SCENE" ]]; then
# Full run: launch reproduce in background with delay
(sleep 120 && launch_reproduce) &
REPRO_PID=$!
echo " Reproduce scheduled in 120s (PID=$REPRO_PID)"
else
# Single scene mode: just run the requested pipeline
if [[ -n "$REPRO_SCENE" ]]; then
launch_reproduce
exit 0
fi
fi
# ---------------------------------------------------------------------------
# Train pipeline (GPU 0)
# ---------------------------------------------------------------------------
echo "[$(date '+%H:%M:%S')] Launching Train pipeline on GPU $TRAIN_GPU..."
TRAIN_ARG=""
if [[ -n "$TRAIN_SCENE" ]]; then
TRAIN_ARG="--scene $TRAIN_SCENE"
fi
conda run --no-capture-output -n "$CONDA_ENV" python -u \
"$TRAIN_SCRIPT" --device-id "$TRAIN_GPU" $TRAIN_ARG \
2>&1 | tee "$SCRIPT_DIR/output/train/run.log"
echo "[$(date '+%H:%M:%S')] Train pipeline complete."
# Wait for reproduce if it was launched
if [[ -n "${REPRO_PID:-}" ]]; then
echo "[$(date '+%H:%M:%S')] Waiting for Reproduce pipeline (PID=$REPRO_PID)..."
wait $REPRO_PID
fi
# ---------------------------------------------------------------------------
# Generate reports
# ---------------------------------------------------------------------------
echo ""
echo "[$(date '+%H:%M:%S')] Both pipelines complete. Generating reports..."
echo " -> Signal plots..."
conda run --no-capture-output -n "$CONDA_ENV" python -u \
"$SCRIPT_DIR/viz_signals.py" 2>&1 | tail -5
echo " -> Vorticity panels..."
conda run --no-capture-output -n "$CONDA_ENV" python -u \
"$SCRIPT_DIR/viz_flow.py" 2>&1 | tail -5
echo " -> Summary report..."
conda run --no-capture-output -n "$CONDA_ENV" python -u \
"$SCRIPT_DIR/generate_report.py" 2>&1 | tail -10
echo ""
echo "=== ALL DONE ==="
echo "Reports in: $SCRIPT_DIR/output/reports/"
CONDA_ENV="${CONDA_ENV:-pycuda_3_10}"
exec conda run --no-capture-output -n "$CONDA_ENV" python -u "$SCRIPT_DIR/infer_train.py" "$@"
-271
View File
@@ -1,271 +0,0 @@
"""Scene manifest — single source of truth for all benchmark scenes.
Two pipelines:
TRAIN_SCENES — V5 PPO models on new 2000x600 config (uniform, free-slip).
REPRODUCE_SCENES — legacy PPO models on old 1280x512 config (parabolic, bounce-back).
Usage:
from eval.scene_manifest import TRAIN_SCENES, REPRODUCE_SCENES
Each scene dict contains:
scene_id, config_path, calibration_path, si, num_steps, scene_type
seeds: list of (seed_label, model_dir) tuples
"""
from __future__ import annotations
import os
from pathlib import Path
from typing import Any, Dict, List
_REPO = Path(__file__).resolve().parents[3]
_TRAIN_DIR = _REPO / "src" / "drl_pinball" / "train"
_CONFIGS_DIR = _REPO / "configs"
KARMAN_CFG = str(_CONFIGS_DIR / "config_lbm_karman_2000x600.json")
KARMAN_CFG_RE60 = str(_CONFIGS_DIR / "config_lbm_karman_2000x600_re60.json")
KARMAN_CFG_RE200 = str(_CONFIGS_DIR / "config_lbm_karman_2000x600_re200.json")
KARMAN_CFG_RE400 = str(_CONFIGS_DIR / "config_lbm_karman_2000x600_re400.json")
PINBALL_CFG = str(_CONFIGS_DIR / "config_lbm_pinball.json")
def _model_dir(case_name: str, seed: int) -> str:
"""Path to model output directory for a trained case."""
return str(_TRAIN_DIR / "output" / f"{case_name}_seed{seed}" / "models")
def _cal_path(name: str) -> str:
return str(_TRAIN_DIR / "calibrations" / name / "calibration.json")
# =============================================================================
# TRAIN SCENES (V5 PPO, 2000x600)
#
# Cloud naming: {kar|ill}_{case}_{sc|tr}_seed{N}
# sc = scratch, tr = transfer
# =============================================================================
TRAIN_SCENES: List[Dict[str, Any]] = [
# -- Karman Cloak Re100 scratch (5 seeds, best model) --
{
"scene_id": "kar_re100_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("kar_re100"),
"si": 800,
"num_steps": 360,
"scene_type": "karman",
"seeds": [
("41", _model_dir("kar_re100_sc", 41)),
("42", _model_dir("kar_re100_sc", 42)),
("43", _model_dir("kar_re100_sc", 43)),
("44", _model_dir("kar_re100_sc", 44)),
("45", _model_dir("kar_re100_sc", 45)),
],
},
# -- Cross-Re Re60 scratch + transfer --
{
"scene_id": "kar_re60_tr",
"config_path": KARMAN_CFG_RE60,
"calibration_path": _cal_path("kar_re60"),
"si": 800,
"num_steps": 360,
"scene_type": "karman",
"seeds": [
("43", _model_dir("kar_re60_tr", 43)),
],
},
# -- Cross-Re Re200 scratch + transfer --
{
"scene_id": "kar_re200_tr",
"config_path": KARMAN_CFG_RE200,
"calibration_path": _cal_path("kar_re200"),
"si": 500,
"num_steps": 360,
"scene_type": "karman",
"seeds": [
("43", _model_dir("kar_re200_tr", 43)),
],
},
# -- Cross-Re Re400 scratch + transfer --
{
"scene_id": "kar_re400_tr",
"config_path": KARMAN_CFG_RE400,
"calibration_path": _cal_path("kar_re400"),
"si": 400,
"num_steps": 360,
"scene_type": "karman",
"seeds": [
("43", _model_dir("kar_re400_tr", 43)),
],
},
# -- VarDist d075 scratch --
{
"scene_id": "kar_d075_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("kar_d075_sc"),
"si": 800,
"num_steps": 360,
"scene_type": "karman",
"seeds": [
("44", _model_dir("kar_d075_sc", 44)),
],
},
# -- VarDist d15 scratch --
{
"scene_id": "kar_d15_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("kar_d15_sc"),
"si": 800,
"num_steps": 360,
"scene_type": "karman",
"seeds": [
("45", _model_dir("kar_d15_sc", 45)),
],
},
# -- VarDist d2 scratch --
{
"scene_id": "kar_d2_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("kar_d2_sc"),
"si": 800,
"num_steps": 360,
"scene_type": "karman",
"seeds": [
("45", _model_dir("kar_d2_sc", 45)),
],
},
# -- Illusion 0.75L scratch --
{
"scene_id": "ill_075L_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("ill_075L"),
"si": 400,
"num_steps": 360,
"scene_type": "illusion",
"target_diam": 0.75,
"seeds": [
("43", _model_dir("ill_075L_sc", 43)),
],
},
# -- Illusion 1.0L scratch --
{
"scene_id": "ill_1L_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("ill_1L"),
"si": 600,
"num_steps": 360,
"scene_type": "illusion",
"target_diam": 1.0,
"seeds": [
("43", _model_dir("ill_1L_sc", 43)),
],
},
# -- Illusion 1.5L scratch --
{
"scene_id": "ill_15L_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("ill_15L"),
"si": 800,
"num_steps": 360,
"scene_type": "illusion",
"target_diam": 1.5,
"seeds": [
("43", _model_dir("ill_15L_sc", 43)),
],
},
# -- Illusion 2.0L scratch (new!) --
{
"scene_id": "ill_2L_sc",
"config_path": KARMAN_CFG,
"calibration_path": _cal_path("ill_2L"),
"si": 800,
"num_steps": 360,
"scene_type": "illusion",
"target_diam": 2.0,
"seeds": [
("43", _model_dir("ill_2L_sc", 43)),
],
},
]
# =============================================================================
# REPRODUCE SCENES (legacy PPO, 1280x512)
# =============================================================================
REPRODUCE_SCENES: List[Dict[str, Any]] = [
{
"scene_id": "karman_cloak",
"config_path": PINBALL_CFG,
"si": 800,
"num_steps": 200,
"scene_type": "karman",
"model_name": "d1a3o12_re100",
"action_scale": 8.0,
"action_bias": (0.0, -4.0, 4.0),
},
{
"scene_id": "steady_cloak",
"config_path": PINBALL_CFG,
"si": 800,
"num_steps": 200,
"scene_type": "steady",
"model_name": None, # open-loop [0, -5.1, 5.1]*U0
"action_scale": 8.0,
"action_bias": (0.0, -4.0, 4.0),
"open_loop_surf_vel": (0.0, -5.1, 5.1),
},
{
"scene_id": "illusion_075L",
"config_path": PINBALL_CFG,
"si": 400,
"num_steps": 200,
"scene_type": "illusion",
"model_name": "d1a3o14_250525_imit_075L_2U_400S",
"action_scale": 8.0,
"action_bias": (0.0, -2.0, 2.0),
"target_diam": 0.75,
},
{
"scene_id": "illusion_1L",
"config_path": PINBALL_CFG,
"si": 600,
"num_steps": 200,
"scene_type": "illusion",
"model_name": "d1a3o14_250525_imit_1L_2U_600S",
"action_scale": 8.0,
"action_bias": (0.0, -2.0, 2.0),
"target_diam": 1.0,
},
{
"scene_id": "illusion_15L",
"config_path": PINBALL_CFG,
"si": 800,
"num_steps": 200,
"scene_type": "illusion",
"model_name": "d1a3o14_250525_imit_15L_2U",
"action_scale": 8.0,
"action_bias": (0.0, -2.0, 2.0),
"target_diam": 1.5,
},
{
"scene_id": "vortex_lamb",
"config_path": PINBALL_CFG,
"si": 800,
"num_steps": 150,
"scene_type": "vortex",
"model_name": "vortex_lamb",
"action_scale": 4.0,
"action_bias": (0.0, -4.0, 4.0),
"vortex_type": "lamb",
"vortex_strength_factor": 0.5,
},
{
"scene_id": "vortex_taylor",
"config_path": PINBALL_CFG,
"si": 800,
"num_steps": 150,
"scene_type": "vortex",
"model_name": "vortex_taylor",
"action_scale": 4.0,
"action_bias": (0.0, -4.0, 4.0),
"vortex_type": "taylor",
"vortex_strength_factor": 0.03,
},
]
+8
View File
@@ -357,6 +357,14 @@ Multiple sources define positions in L0 units; multiply by L0=20 to get lattice
### 8.1 Reward Functions
> **DTW FORMULA NOTE (2026-07-13):** This document describes the LEGACY DTW formula used in old envs:
`sim = 1.0 - DTW_dist / conv_len`.
> The V5 pipeline (`train/env_karman.py`, `train/env_illusion.py`) uses a different formula:
`sim = max(0.0, 1.0 - DTW_dist / (conv_len × dtw_norm_scale))`
> where `dtw_norm_scale ≈ 0.2-0.4`. This makes V5 scores ~0.05-0.3 lower than equivalent
> Legacy scores. For fair comparison, see `eval/recompute_unified_dtw.py`.
> Results: V5 uniformly outperforms Legacy on Karman cross-Re when using the same formula.
**Cloak (Karman)**:
```python
reward_cd = exp(-|cd * 20|) # cd = (Σforces_fx) / 3
+10 -57
View File
@@ -1,64 +1,17 @@
# Legacy Test (Track A)
# LegacyCelerisLab reproduction
Systematic validation of pre-trained PPO models using **LegacyCelerisLab**
(the original CFD solver the models were trained with).
## Quick Start
One generic runner reproduces the historical Legacy `CustomEnv` contracts. It recomputes normalization from the zero-action initialization FIFO on every run for diagnostics and reset FIFO state, but model inference and reward force scaling use the explicitly mapped frozen training normalization from active `src/SR_analysis/data/<scene>/<case>/norm.json`. Policy cases fail closed when that file is missing or invalid; steady has no policy and uses the recomputed norm. Reset restores the scene DDF and recomputed saved FIFO, and gives the policy an exactly zero initial observation. `FlowField.run` is the only action EMA. An active compatibility adapter zeros the host action immediately before every public interval, reproducing the pre-persistence driver contract while preserving the within-interval EMA; `LegacyCelerisLab/driver.py` remains unchanged.
```bash
# Run all legacy tests sequentially (GPU 1, 60s delay between tests)
bash src/drl_pinball/legacy_test/run_all_legacy_tests.sh 1
# Single scene
conda run -n pycuda_3_10 python src/drl_pinball/legacy_test/test_karman_cloak_re100.py --device 1
src/drl_pinball/legacy_test/run_legacy.sh karman_re100 0
# or inside the environment
python -m src.drl_pinball.legacy_test illusion_1L --device 0
# Full CFD/policy metrics with no output mutation or rendering
src/drl_pinball/legacy_test/run_legacy.sh karman_re50 0 --metrics-only
```
## Directory
Supported cases: Karman Re50/100/200/400, Illusion 0.75/1/1.5, Vortex lamb/taylor, and open-loop steady. `erase` exits as explicitly unsupported because its historical reward uses an evolving controlled-force FIFO as its target; no exact independent target-scale mapping exists.
```
legacy_test/
├── README.md # This file
├── core/
│ ├── comparator.py # Compare signals against SR_analysis reference
│ ├── dtw_metrics.py # DTW/harmonics (re-exports from reproduce/core/)
│ ├── io_helpers.py # Save/load .npz, norm.json
│ ├── legacy_env_builder.py # FlowField builders for all 5 scene types
│ └── model_loader.py # PPO model loading (wraps ModelInventory)
├── test_karman_cloak_re100.py # Flagship: Karman Cloak Re100
├── test_karman_cloak_crossre.py # Cross-Re: re50, re200, re400
├── test_steady_cloak.py # Steady cloak (open-loop, no DRL)
├── test_illusion_1L.py # Illusion 1.0L (S_DIM=14)
├── test_illusion_remaining.py # Illusion 0.75L, 1.5L
├── test_vortex_lamb.py # Vortex Lamb dipole
├── test_vortex_taylor.py # Vortex Taylor monopole
├── test_erase.py # Erase (experimental, known incomplete)
├── run_all_legacy_tests.sh # Sequential launcher
└── output/ # Per-scene verification outputs
```
In normal mode, at run start, the selected `legacy_test/output/<case>` directory is removed and recreated so stale files cannot be mistaken for current evidence; other case directories are untouched. Each run writes only `{signals.npz,reset_contract.npz,norm.json,metrics.json,final_vorticity.png}` there. `norm.json` records separate `policy` and `recomputed` sections, while `metrics.json` records the policy norm source path and SHA256. `metrics.json` distinguishes target-native legacy DTW (the reward input), normalized target-scale DTW computed offline, reward summaries, model/config provenance, and an optional `frozen_reference_comparison` loaded only from the active `src/SR_analysis/data/<scene>/<case>/controlled.npz` path. With `--metrics-only`, the same full CFD/policy rollout and in-memory metrics/frozen comparison run, but no case directory is inspected, created, deleted, or written and rendering is skipped; compact metrics JSON is printed to stdout. The vorticity PNG uses the shared `CelerisLab.common.render` renderer, Legacy flag-masked `q/RHO_ref` physical lattice velocity/vorticity with `RHO_ref=1`, solid masking, and fixed limits `[-0.001, 0.001]`.
## Scene Coverage
| Scene | S_DIM | Scale/Bias | SI | MaxSteps | Legacy Ref |
|-------|-------|------------|-----|----------|------------|
| Karman re100 | 12 | 8/(0,-4,4) | 800 | 500 | `legacy_karman_env.py` |
| Karman re50 | 12 | 8/(0,-4,4) | 800 | 500 | same, ν=0.008 |
| Karman re200 | 12 | 8/(0,-4,4) | 800 | 500 | same, ν=0.002 |
| Karman re400 | 12 | 8/(0,-4,4) | 800 | 500 | same, ν=0.001 |
| Steady Cloak | — | open-loop | 800 | 200 | OID `collect_steady_cloak.py` |
| Illusion 0.75L | 14 | 8/(0,-2,2) | 400 | 500 | `legacy_env_imit.py` |
| Illusion 1L | 14 | 8/(0,-2,2) | 600 | 500 | same |
| Illusion 1.5L | 14 | 8/(0,-2,2) | 800 | 500 | same |
| Vortex Lamb | 12 | 4/(0,-4,4) | 800 | 150 | `legacy_env_vortex.py` |
| Vortex Taylor | 12 | 4/(0,-4,4) | 800 | 150 | same |
| Erase | 12 | 8/(0,-8,8) | 600 | 500 | `legacy_env_erase.py` |
## Design Notes
- **Object ordering** matches legacy EXACTLY (documented in `knowledge.md` Section 9):
- Karman/Erase: dist_cyl(0) [or sensor0(0) for erase], sensors(1-3), front(4), top(5), bottom(6)
- Steady/Illusion/Vortex: sensors(0-2), front(3), top(4), bottom(5)
- **DDF checkpoint timing** uses pre-bias save + test-side bias FIFO (matching legacy `save_ddf()` pattern)
- **Action** uses legacy `FlowField.run()` built-in EMA smoothing (weight 0.1)
- **Comparison** uses DTW similarity > 0.95 as primary pass criterion (phase-invariant)
- **Steady cloak** is open-loop — verifies lift RMS suppression, no DTW comparison
- **Erase** is known incomplete — no DTW threshold enforced
The previous scene-specific scripts are retained under `archive/duplicated_scripts/` and are not active entry points. CPU-only contract tests live in `tests/`; they do not initialize CUDA.
@@ -12,7 +12,6 @@ Usage: conda run -n pycuda_3_10 python test_erase.py --device 0
"""
import argparse
import json
import os
import sys
import time
@@ -32,6 +31,10 @@ from legacy_test.core.legacy_env_builder import ( # noqa: E402
)
from legacy_test.core.model_loader import load_model # noqa: E402
from legacy_test.core.comparator import compare_scene # noqa: E402
from legacy_test.core.output_helpers import (
atomic_write_json, render_final_vorticity, save_controlled, save_run_metadata,
save_run_status,
) # noqa: E402
from legacy_test.core.io_helpers import save_signals, save_target, save_norm # noqa: E402
SAMPLE_INTERVAL = 600
@@ -114,10 +117,10 @@ def main():
sens_norm = (raw[0:6] - s_dev) / s_nf
obs = np.clip(np.hstack([forces_norm, sens_norm]), -1.0, 1.0).astype(np.float32)
render_final_vorticity(ff, args.out, u0=U0, l0=20.0)
save_signals(args.out, sig_s, sig_f, sig_a)
np.savez_compressed(os.path.join(args.out, "controlled.npz"),
sensors=sig_s, forces=sig_f, actions=sig_a,
rewards=np.zeros(NUM_STEPS, dtype=np.float32))
save_controlled(args.out, sig_s, sig_f, sig_a)
save_run_metadata(args.out, reward_available=False, reward_definition="unavailable: exact legacy_env_erase rolling full-state reward history/checkpoint mapping is not established", warmup_count=0, sample_interval=SAMPLE_INTERVAL, action_scale=ACTION_SCALE, action_bias=ACTION_BIAS, model=MODEL_NAME)
# Note: erase has no dedicated SR_analysis reference; compare against karman_re100 as fallback
try:
@@ -126,11 +129,12 @@ def main():
log(" No reference data found for erase — skipping comparison.")
result = {"passed": False, "dtw_sim": 0.0}
with open(os.path.join(args.out, "result.json"), "w") as f:
json.dump(result, f, indent=2)
log(f"PASS" if result["passed"] else "FAIL (erase is known incomplete)")
atomic_write_json(os.path.join(args.out, "result.json"), result)
historical_sr_passed = bool(result["passed"])
save_run_status(args.out, historical_sr_passed=historical_sr_passed)
log(f"Historical SR comparison: {'PASS' if historical_sr_passed else 'FAIL'} (erase is known incomplete)")
del ff
return 0 if result["passed"] else 0 # Always return 0 for erase
return 0
if __name__ == "__main__":
@@ -38,6 +38,11 @@ from legacy_test.core.legacy_env_builder import ( # noqa: E402
)
from legacy_test.core.model_loader import load_model # noqa: E402
from legacy_test.core.comparator import compare_scene # noqa: E402
from legacy_test.core.reward_helpers import illusion_reward # noqa: E402
from legacy_test.core.output_helpers import (
atomic_write_json, render_final_vorticity, save_controlled, save_run_metadata,
save_run_status,
) # noqa: E402
from legacy_test.core.io_helpers import save_signals # noqa: E402
from legacy_test.core.dtw_metrics import gen_target_states_at # noqa: E402
@@ -152,6 +157,7 @@ def main():
sig_s = np.zeros((NUM_STEPS, 6), dtype=np.float32)
sig_f = np.zeros((NUM_STEPS, 6), dtype=np.float32)
sig_a = np.zeros((NUM_STEPS, 3), dtype=np.float32)
sig_r = np.zeros(NUM_STEPS, dtype=np.float32)
# Build initial observation using REFERENCE norm
raw = ff.obs.copy()[0:12]
@@ -185,20 +191,22 @@ def main():
obs_12 = np.clip(np.hstack([forces_norm, sens_norm]), -1.0, 1.0).astype(np.float32)
tgt = gen_target_states_at(step + 1, target_harmonics)
obs = np.clip(np.hstack([obs_12, [tgt[0] / f_nf, tgt[1] / f_nf]]), -1.0, 1.0).astype(np.float32)
sig_r[step] = illusion_reward(target_states, target_harmonics, np.array(fifo), f_nf, current_step=step, conv_len=36)
render_final_vorticity(ff, args.out, u0=U0, l0=L0)
save_signals(args.out, sig_s, sig_f, sig_a)
np.savez_compressed(os.path.join(args.out, "controlled.npz"),
sensors=sig_s, forces=sig_f, actions=sig_a,
rewards=np.zeros(NUM_STEPS, dtype=np.float32))
save_controlled(args.out, sig_s, sig_f, sig_a, sig_r)
save_run_metadata(args.out, reward_available=True, reward_definition="legacy_env_imit: 0.3*r_cd + 0.3*r_cl + 0.4*r_sim; target force columns 0:2, sensor columns 2:8", warmup_count=0, sample_interval=SAMPLE_INTERVAL, action_scale=ACTION_SCALE, action_bias=ACTION_BIAS, model="d1a3o14_250525_imit_1L_2U_600S")
log("Comparing against reference...")
result = compare_scene(REF_DIR, sig_s, sig_f, sig_a, conv_len=36, label="illusion_1L")
with open(os.path.join(args.out, "result.json"), "w") as f:
json.dump(result, f, indent=2)
atomic_write_json(os.path.join(args.out, "result.json"), result)
log(f"PASS" if result["passed"] else "FAIL")
historical_sr_passed = bool(result["passed"])
save_run_status(args.out, historical_sr_passed=historical_sr_passed)
log(f"Historical SR comparison: {'PASS' if historical_sr_passed else 'FAIL'}")
del ff
return 0 if result["passed"] else 1
return 0
if __name__ == "__main__":
@@ -5,7 +5,6 @@ Usage: conda run -n pycuda_3_10 python test_illusion_remaining.py --device 0
"""
import argparse
import json
import os
import sys
import time
@@ -25,6 +24,11 @@ from legacy_test.core.legacy_env_builder import ( # noqa: E402
)
from legacy_test.core.model_loader import load_model # noqa: E402
from legacy_test.core.comparator import compare_scene # noqa: E402
from legacy_test.core.reward_helpers import illusion_reward # noqa: E402
from legacy_test.core.output_helpers import (
atomic_write_json, render_final_vorticity, save_controlled, save_run_metadata,
save_run_status,
) # noqa: E402
from legacy_test.core.io_helpers import save_signals, save_target, save_norm # noqa: E402
from legacy_test.core.dtw_metrics import gen_target_states_at # noqa: E402
@@ -44,10 +48,10 @@ def log(msg): print(f"[{time.strftime('%H:%M:%S')}] {msg}", flush=True)
def run_one(device_id: int, label: str, diam_L: float, si: int,
model_name: str, ref_subdir: str) -> dict:
model_name: str, ref_subdir: str, out_dir=None) -> dict:
log(f"=== {label}: Legacy Test ===")
ref_dir = os.path.join(_SRC, "SR_analysis", "data", "illusion", ref_subdir)
out_dir = os.path.join(OUT_BASE, label)
out_dir = os.path.join(OUT_BASE, label) if out_dir is None else os.fspath(out_dir)
os.makedirs(out_dir, exist_ok=True)
data = build_illusion(device_id=device_id, target_diameter_L=diam_L,
@@ -66,6 +70,7 @@ def run_one(device_id: int, label: str, diam_L: float, si: int,
log(f" Model: {model_name}")
ff.restore_ddf(); ff.apply_ddf()
reward_fifo = deque(maxlen=FIFO_LEN)
init_bias = (0.0, -1.0, 1.0)
bias_arr = np.zeros(n_obj, dtype=DATA_TYPE)
bias_arr[3] = float(init_bias[0] * U0)
@@ -74,10 +79,12 @@ def run_one(device_id: int, label: str, diam_L: float, si: int,
for _ in range(FIFO_LEN):
ff.run(si, bias_arr)
reward_fifo.append(ff.obs.copy()[0:12])
sig_s = np.zeros((NUM_STEPS, 6), dtype=np.float32)
sig_f = np.zeros((NUM_STEPS, 6), dtype=np.float32)
sig_a = np.zeros((NUM_STEPS, 3), dtype=np.float32)
sig_r = np.zeros(NUM_STEPS, dtype=np.float32)
raw = ff.obs.copy()[0:12]
forces_norm = raw[6:12] / f_nf
@@ -103,6 +110,7 @@ def run_one(device_id: int, label: str, diam_L: float, si: int,
ff.context.pop()
raw = ff.obs.copy()[0:12]
reward_fifo.append(raw.copy())
sig_s[step] = raw[0:6]
sig_f[step] = raw[6:12]
@@ -113,31 +121,46 @@ def run_one(device_id: int, label: str, diam_L: float, si: int,
tcd = tgt[0] / f_nf if f_nf > 1e-12 else 0.0
tcl = tgt[1] / f_nf if f_nf > 1e-12 else 0.0
obs = np.clip(np.hstack([obs_12, [tcd, tcl]]), -1.0, 1.0).astype(np.float32)
sig_r[step] = illusion_reward(data["target_states"], target_harmonics, np.array(reward_fifo), f_nf, current_step=step, conv_len=36)
render_final_vorticity(ff, out_dir, u0=U0, l0=20.0)
save_signals(out_dir, sig_s, sig_f, sig_a)
np.savez_compressed(os.path.join(out_dir, "controlled.npz"),
sensors=sig_s, forces=sig_f, actions=sig_a,
rewards=np.zeros(NUM_STEPS, dtype=np.float32))
save_controlled(out_dir, sig_s, sig_f, sig_a, sig_r)
save_run_metadata(out_dir, reward_available=True, reward_definition="legacy_env_imit: 0.3*r_cd + 0.3*r_cl + 0.4*r_sim; target force columns 0:2, sensor columns 2:8", warmup_count=0, sample_interval=si, action_scale=ACTION_SCALE, action_bias=ACTION_BIAS, model=model_name)
log(" Comparing against reference...")
result = compare_scene(ref_dir, sig_s, sig_f, sig_a, conv_len=36, label=label)
with open(os.path.join(out_dir, "result.json"), "w") as f:
json.dump(result, f, indent=2)
log(f" {'PASS' if result['passed'] else 'FAIL'}")
atomic_write_json(os.path.join(out_dir, "result.json"), result)
historical_sr_passed = bool(result["passed"])
save_run_status(out_dir, historical_sr_passed=historical_sr_passed)
log(f" Historical SR comparison: {'PASS' if historical_sr_passed else 'FAIL'}")
del ff
return result
def selected_cases(diameter: float | None):
cases = ILLUSION_CASES if diameter is None else [case for case in ILLUSION_CASES if abs(case[1] - diameter) < 1e-9]
if not cases:
raise ValueError(f"unsupported illusion diameter {diameter}; expected 0.75 or 1.5")
return cases
def main():
ap = argparse.ArgumentParser(); ap.add_argument("--device", type=int, default=0)
ap.add_argument("--diam", type=float, choices=(0.75, 1.5), default=None,
help="run one diameter only; default runs both historical cases")
ap.add_argument("--out", type=str, default=None,
help="explicit output directory for one selected case; multi-case runs use child labels")
args = ap.parse_args()
cases = selected_cases(args.diam)
results = {}
for label, diam_L, si, model_name, ref_subdir in ILLUSION_CASES:
results[label] = run_one(args.device, label, diam_L, si, model_name, ref_subdir)
for label, diam_L, si, model_name, ref_subdir in cases:
out_dir = args.out if len(cases) == 1 and args.out else (os.path.join(args.out, label) if args.out else None)
results[label] = run_one(args.device, label, diam_L, si, model_name, ref_subdir, out_dir)
log("\n=== Illusion Remaining Summary ===")
for name, r in results.items():
log(f" {name}: DTW={r['dtw_sim']:.4f}, act_corr={[f'{c:.3f}' for c in r['action_corr']]} -> {'PASS' if r['passed'] else 'FAIL'}")
log(f" {name}: DTW={r['dtw_sim']:.4f}, act_corr={[f'{c:.3f}' for c in r['action_corr']]} -> historical SR {'PASS' if r['passed'] else 'FAIL'}")
if __name__ == "__main__":
@@ -9,7 +9,6 @@ Usage: conda run -n pycuda_3_10 python test_karman_cloak_crossre.py --device 0
"""
import argparse
import json
import os
import sys
import time
@@ -29,6 +28,11 @@ from legacy_test.core.legacy_env_builder import ( # noqa: E402
)
from legacy_test.core.model_loader import load_model # noqa: E402
from legacy_test.core.comparator import compare_scene # noqa: E402
from legacy_test.core.reward_helpers import karman_reward # noqa: E402
from legacy_test.core.output_helpers import (
atomic_write_json, render_final_vorticity, save_controlled, save_run_metadata,
save_run_status,
) # noqa: E402
from legacy_test.core.io_helpers import save_signals, save_target, save_norm # noqa: E402
SAMPLE_INTERVAL = 800
@@ -41,12 +45,12 @@ OUT_BASE = os.path.join(os.path.dirname(__file__), "output")
def log(msg): print(f"[{time.strftime('%H:%M:%S')}] {msg}", flush=True)
def run_crossre(device_id: int, re_code: float) -> dict:
def run_crossre(device_id: int, re_code: float, out_dir=None) -> dict:
label = f"karman_re{int(re_code)}"
log(f"=== {label}: Legacy Test ===")
model_name = f"d1a3o12_re{int(re_code)}"
ref_dir = os.path.join(_SRC, "SR_analysis", "data", "karman", label)
out_dir = os.path.join(OUT_BASE, label)
out_dir = os.path.join(OUT_BASE, label) if out_dir is None else os.fspath(out_dir)
os.makedirs(out_dir, exist_ok=True)
data = build_karman_cloak(device_id=device_id, re_code=re_code,
@@ -80,6 +84,7 @@ def run_crossre(device_id: int, re_code: float) -> dict:
sig_s = np.zeros((NUM_STEPS, 6), dtype=np.float32)
sig_f = np.zeros((NUM_STEPS, 6), dtype=np.float32)
sig_a = np.zeros((NUM_STEPS, 3), dtype=np.float32)
sig_r = np.zeros(NUM_STEPS, dtype=np.float32)
raw = ff.obs.copy()[2:14]
forces_norm = raw[6:12] / f_nf
@@ -108,17 +113,19 @@ def run_crossre(device_id: int, re_code: float) -> dict:
forces_norm = raw[6:12] / f_nf
sens_norm = (raw[0:6] - s_dev) / s_nf
obs = np.clip(np.hstack([forces_norm, sens_norm]), -1.0, 1.0).astype(np.float32)
sig_r[step] = karman_reward(target_states, np.array(fifo), f_nf, step=step, conv_len=CONV_LEN)
render_final_vorticity(ff, out_dir, u0=U0, l0=20.0)
save_signals(out_dir, sig_s, sig_f, sig_a)
np.savez_compressed(os.path.join(out_dir, "controlled.npz"),
sensors=sig_s, forces=sig_f, actions=sig_a,
rewards=np.zeros(NUM_STEPS, dtype=np.float32))
save_controlled(out_dir, sig_s, sig_f, sig_a, sig_r)
save_run_metadata(out_dir, reward_available=True, reward_definition="legacy_env_karman_cloak_standard: warmup zero, then 0.3*r_cd + 0.4*r_cl + 0.3*r_sim", warmup_count=CONV_LEN, sample_interval=SAMPLE_INTERVAL, action_scale=ACTION_SCALE, action_bias=ACTION_BIAS, model=model_name)
log(" Comparing against reference...")
result = compare_scene(ref_dir, sig_s, sig_f, sig_a, conv_len=CONV_LEN, label=label)
with open(os.path.join(out_dir, "result.json"), "w") as f:
json.dump(result, f, indent=2)
log(f" {'PASS' if result['passed'] else 'FAIL'}")
atomic_write_json(os.path.join(out_dir, "result.json"), result)
historical_sr_passed = bool(result["passed"])
save_run_status(out_dir, historical_sr_passed=historical_sr_passed)
log(f" Historical SR comparison: {'PASS' if historical_sr_passed else 'FAIL'}")
del ff
return result
@@ -127,16 +134,22 @@ def main():
ap = argparse.ArgumentParser(); ap.add_argument("--device", type=int, default=0)
ap.add_argument("--re", type=str, default="50,200,400",
help="Comma-separated Re values")
ap.add_argument("--out", type=str, default=None,
help="explicit output directory for one selected Re; multi-Re runs use child labels")
args = ap.parse_args()
re_values = [float(value.strip()) for value in args.re.split(",")]
if not re_values or any(value not in (50.0, 200.0, 400.0) for value in re_values):
ap.error("--re must select from 50,200,400")
results = {}
for re_str in args.re.split(","):
re_val = float(re_str.strip())
results[f"re{int(re_val)}"] = run_crossre(args.device, re_val)
for re_val in re_values:
label = f"karman_re{int(re_val)}"
out_dir = args.out if len(re_values) == 1 and args.out else (os.path.join(args.out, label) if args.out else None)
results[f"re{int(re_val)}"] = run_crossre(args.device, re_val, out_dir)
log("\n=== Cross-Re Summary ===")
for name, r in results.items():
log(f" {name}: DTW={r['dtw_sim']:.4f}, act_corr={[f'{c:.3f}' for c in r['action_corr']]} -> {'PASS' if r['passed'] else 'FAIL'}")
log(f" {name}: DTW={r['dtw_sim']:.4f}, act_corr={[f'{c:.3f}' for c in r['action_corr']]} -> historical SR {'PASS' if r['passed'] else 'FAIL'}")
if __name__ == "__main__":
@@ -15,7 +15,6 @@ Expected: near-perfect match (same CFD, same model).
from __future__ import annotations
import argparse
import json
import os
import sys
import time
@@ -37,6 +36,11 @@ from legacy_test.core.legacy_env_builder import ( # noqa: E402
)
from legacy_test.core.model_loader import load_model # noqa: E402
from legacy_test.core.comparator import compare_scene # noqa: E402
from legacy_test.core.reward_helpers import karman_reward # noqa: E402
from legacy_test.core.output_helpers import (
atomic_write_json, render_final_vorticity, save_controlled, save_run_metadata,
save_run_status,
) # noqa: E402
from legacy_test.core.io_helpers import ( # noqa: E402
save_signals, save_target, save_norm,
)
@@ -157,55 +161,29 @@ def main():
sens_norm = (obs_slice[0:6] - sens_deviation) / sens_norm_fact
obs = np.clip(np.hstack([forces_norm, sens_norm]), -1.0, 1.0).astype(np.float32)
# Compute reward (exact legacy formula)
if step >= CONV_LEN:
states_arr = np.array(fifo, dtype=np.float32)
forces = states_arr[-1, 6:12] / force_norm_fact
cd = float((forces[0] + forces[2] + forces[4]) / 3.0)
cl = float((forces[1] + forces[3] + forces[5]) / 3.0)
sig_r[step] = karman_reward(target_states, np.array(fifo), force_norm_fact, step=step, conv_len=CONV_LEN)
# DTW similarity (legacy calc_lag + calc_dtw_sim)
from legacy_test.core.dtw_metrics import calc_lag, calc_dtw_sim
mid_idx = 1 # sensor1_uy
t_seq = target_states[CONV_LEN:2 * CONV_LEN, mid_idx]
s_seq = states_arr[-CONV_LEN:, mid_idx]
lag = calc_lag(t_seq, s_seq)
sim_sum = 0.0
for i in range(6):
t_seq2 = np.roll(target_states[:, i], -lag)[CONV_LEN:2 * CONV_LEN]
s_seq2 = states_arr[-CONV_LEN:, i]
sim_sum += calc_dtw_sim(t_seq2, s_seq2)
sim_val = float(sim_sum / 6.0)
r_cd = float(np.exp(-abs(cd * 20.0)))
r_cl = float(np.exp(-abs(cl * 80.0)))
r_sim = float(np.exp(-10.0 * abs(sim_val - 1.0)))
sig_r[step] = float(min(0.3 * r_cd + 0.4 * r_cl + 0.3 * r_sim, 1.0))
render_final_vorticity(ff, args.out, u0=U0, l0=L0)
# Save signals
save_signals(args.out, sig_s, sig_f, sig_a, name="controlled")
save_signals(args.out, sig_s, sig_f, sig_a, name="uncontrolled")
# Also save to match SR_analysis format (with rewards)
np.savez_compressed(
os.path.join(args.out, "controlled.npz"),
sensors=sig_s, forces=sig_f, actions=sig_a, rewards=sig_r,
)
save_controlled(args.out, sig_s, sig_f, sig_a, sig_r)
save_run_metadata(args.out, reward_available=True, reward_definition="legacy_env_karman_cloak_standard: warmup zero, then 0.3*r_cd + 0.4*r_cl + 0.3*r_sim", warmup_count=CONV_LEN, sample_interval=SAMPLE_INTERVAL, action_scale=ACTION_SCALE, action_bias=ACTION_BIAS, model=args.model)
# Save config
with open(os.path.join(args.out, "config.json"), "w") as f:
json.dump({
"device_id": args.device,
"re_code": 100.0,
"viscosity": 0.004,
"u0": float(U0),
"sample_interval": SAMPLE_INTERVAL,
"num_steps": NUM_STEPS,
"action_scale": ACTION_SCALE,
"action_bias": ACTION_BIAS.tolist(),
"model": args.model,
}, f, indent=2)
atomic_write_json(os.path.join(args.out, "config.json"), {
"device_id": args.device,
"re_code": 100.0,
"viscosity": 0.004,
"u0": float(U0),
"sample_interval": SAMPLE_INTERVAL,
"num_steps": NUM_STEPS,
"action_scale": ACTION_SCALE,
"action_bias": ACTION_BIAS,
"model": args.model,
})
# ---- Phase 4: Compare against reference ----
log("\n=== Comparison against SR_analysis reference ===")
@@ -216,24 +194,25 @@ def main():
label="karman_re100",
)
with open(os.path.join(args.out, "result.json"), "w") as f:
json.dump(result, f, indent=2)
atomic_write_json(os.path.join(args.out, "result.json"), result)
log(f"\nFinal reward: mean={sig_r.mean():.4f}, last_50={sig_r[-50:].mean():.4f}")
log(f"DTW similarity: {result['dtw_sim']:.4f}")
log(f"Action corr: {result['action_corr']}")
if result["passed"]:
log("PASS — All metrics within thresholds.")
historical_sr_passed = bool(result["passed"])
save_run_status(args.out, historical_sr_passed=historical_sr_passed)
if historical_sr_passed:
log("Historical SR comparison: PASS — all metrics within thresholds.")
else:
log("FAIL — One or more metrics below threshold.")
log("Historical SR comparison: FAIL — one or more metrics below threshold.")
# Cleanup
del ff
log("Done.")
return 0 if result["passed"] else 1
return 0
if __name__ == "__main__":
@@ -13,7 +13,6 @@ Usage: conda run -n pycuda_3_10 python test_steady_cloak.py --device 0
"""
import argparse
import json
import os
import sys
import time
@@ -30,6 +29,11 @@ for p in [_REPO, _SRC, _DRL]:
from LegacyCelerisLab import FlowField # noqa: E402
from LegacyCelerisLab import utils as legacy_utils # noqa: E402
from legacy_test.core.output_helpers import (
atomic_write_json, render_final_vorticity, save_controlled, save_run_metadata,
save_run_status,
) # noqa: E402
from legacy_test.core.legacy_env_builder import ( # noqa: E402
L0, U0, DATA_TYPE, FIFO_LEN,
_center_y, _stabilize,
@@ -108,9 +112,9 @@ def main():
sig_f[s] = obs[6:12]
save_actions = np.zeros((NUM_STEPS, 3), dtype=np.float32)
np.savez_compressed(os.path.join(args.out, "controlled.npz"),
sensors=sig_s, forces=sig_f, actions=save_actions,
rewards=np.zeros(NUM_STEPS, dtype=np.float32))
render_final_vorticity(ff, args.out, u0=U0, l0=L0)
save_controlled(args.out, sig_s, sig_f, save_actions)
save_run_metadata(args.out, reward_available=False, reward_definition="unavailable: open-loop scene evaluated by lift_rms criterion, not a trained reward", warmup_count=100, sample_interval=SAMPLE_INTERVAL, action_scale=None, action_bias=ACTION_BIAS, model=None)
# Check force balance (steady cloak should suppress lift oscillations)
front_fy_mean = float(np.mean(sig_f[:, 1]))
@@ -131,8 +135,14 @@ def main():
else:
result = {"passed": False, "lift_rms": float(lift_rms)}
with open(os.path.join(args.out, "result.json"), "w") as f:
json.dump(result, f, indent=2)
atomic_write_json(os.path.join(args.out, "result.json"), result)
save_run_status(
args.out,
historical_sr_passed=None,
physical_lift_criterion_passed=bool(passed),
lift_rms=lift_rms,
lift_rms_threshold=0.01,
)
del ff
return 0 if passed else 1
@@ -9,7 +9,6 @@ Usage: conda run -n pycuda_3_10 python test_vortex_lamb.py --device 0
"""
import argparse
import json
import os
import sys
import time
@@ -29,6 +28,11 @@ from legacy_test.core.legacy_env_builder import ( # noqa: E402
)
from legacy_test.core.model_loader import load_model # noqa: E402
from legacy_test.core.comparator import compare_scene # noqa: E402
from legacy_test.core.reward_helpers import vortex_reward # noqa: E402
from legacy_test.core.output_helpers import (
atomic_write_json, render_final_vorticity, save_controlled, save_run_metadata,
save_run_status,
) # noqa: E402
from legacy_test.core.io_helpers import save_signals, save_target, save_norm # noqa: E402
SAMPLE_INTERVAL = 800
@@ -86,6 +90,7 @@ def main():
sig_s = np.zeros((NUM_STEPS, 6), dtype=np.float32)
sig_f = np.zeros((NUM_STEPS, 6), dtype=np.float32)
sig_a = np.zeros((NUM_STEPS, 3), dtype=np.float32)
sig_r = np.zeros(NUM_STEPS, dtype=np.float32)
raw = ff.obs.copy()[0:12]
forces_norm = raw[6:12] / f_nf
@@ -114,20 +119,22 @@ def main():
forces_norm = raw[6:12] / f_nf
sens_norm = (raw[0:6] - s_dev) / s_nf
obs = np.clip(np.hstack([forces_norm, sens_norm]), -1.0, 1.0).astype(np.float32)
sig_r[step] = vortex_reward(target_states, np.array(fifo), f_nf, current_step=step, conv_len=CONV_LEN)
render_final_vorticity(ff, args.out, u0=U0, l0=20.0)
save_signals(args.out, sig_s, sig_f, sig_a)
np.savez_compressed(os.path.join(args.out, "controlled.npz"),
sensors=sig_s, forces=sig_f, actions=sig_a,
rewards=np.zeros(NUM_STEPS, dtype=np.float32))
save_controlled(args.out, sig_s, sig_f, sig_a, sig_r)
save_run_metadata(args.out, reward_available=True, reward_definition="legacy_env_vortex: 0.2*r_cd + 0.3*r_cl + 0.5*r_sim with current-step target roll", warmup_count=0, sample_interval=SAMPLE_INTERVAL, action_scale=ACTION_SCALE, action_bias=ACTION_BIAS, model="vortex_lamb")
log("Comparing against reference...")
result = compare_scene(REF_DIR, sig_s, sig_f, sig_a, conv_len=CONV_LEN, label="vortex_lamb")
with open(os.path.join(args.out, "result.json"), "w") as f:
json.dump(result, f, indent=2)
atomic_write_json(os.path.join(args.out, "result.json"), result)
log(f"PASS" if result["passed"] else "FAIL")
historical_sr_passed = bool(result["passed"])
save_run_status(args.out, historical_sr_passed=historical_sr_passed)
log(f"Historical SR comparison: {'PASS' if historical_sr_passed else 'FAIL'}")
del ff
return 0 if result["passed"] else 1
return 0
if __name__ == "__main__":
@@ -7,7 +7,6 @@ Usage: conda run -n pycuda_3_10 python test_vortex_taylor.py --device 0
"""
import argparse
import json
import os
import sys
import time
@@ -27,6 +26,11 @@ from legacy_test.core.legacy_env_builder import ( # noqa: E402
)
from legacy_test.core.model_loader import load_model # noqa: E402
from legacy_test.core.comparator import compare_scene # noqa: E402
from legacy_test.core.reward_helpers import vortex_reward # noqa: E402
from legacy_test.core.output_helpers import (
atomic_write_json, render_final_vorticity, save_controlled, save_run_metadata,
save_run_status,
) # noqa: E402
from legacy_test.core.io_helpers import save_signals, save_target, save_norm # noqa: E402
SAMPLE_INTERVAL = 800
@@ -79,6 +83,7 @@ def main():
sig_s = np.zeros((NUM_STEPS, 6), dtype=np.float32)
sig_f = np.zeros((NUM_STEPS, 6), dtype=np.float32)
sig_a = np.zeros((NUM_STEPS, 3), dtype=np.float32)
sig_r = np.zeros(NUM_STEPS, dtype=np.float32)
raw = ff.obs.copy()[0:12]
obs = np.clip(np.hstack([(raw[6:12] / f_nf), ((raw[0:6] - s_dev) / s_nf)]), -1.0, 1.0).astype(np.float32)
@@ -101,19 +106,21 @@ def main():
fifo.append(raw)
sig_s[step] = raw[0:6]; sig_f[step] = raw[6:12]
obs = np.clip(np.hstack([(raw[6:12] / f_nf), ((raw[0:6] - s_dev) / s_nf)]), -1.0, 1.0).astype(np.float32)
sig_r[step] = vortex_reward(target_states, np.array(fifo), f_nf, current_step=step, conv_len=CONV_LEN)
render_final_vorticity(ff, args.out, u0=U0, l0=20.0)
save_signals(args.out, sig_s, sig_f, sig_a)
np.savez_compressed(os.path.join(args.out, "controlled.npz"),
sensors=sig_s, forces=sig_f, actions=sig_a,
rewards=np.zeros(NUM_STEPS, dtype=np.float32))
save_controlled(args.out, sig_s, sig_f, sig_a, sig_r)
save_run_metadata(args.out, reward_available=True, reward_definition="legacy_env_vortex: 0.2*r_cd + 0.3*r_cl + 0.5*r_sim with current-step target roll", warmup_count=0, sample_interval=SAMPLE_INTERVAL, action_scale=ACTION_SCALE, action_bias=ACTION_BIAS, model="vortex_taylor")
log("Comparing against reference...")
result = compare_scene(REF_DIR, sig_s, sig_f, sig_a, conv_len=CONV_LEN, label="vortex_taylor")
with open(os.path.join(args.out, "result.json"), "w") as f:
json.dump(result, f, indent=2)
log(f"{'PASS' if result['passed'] else 'FAIL'}")
atomic_write_json(os.path.join(args.out, "result.json"), result)
historical_sr_passed = bool(result["passed"])
save_run_status(args.out, historical_sr_passed=historical_sr_passed)
log(f"Historical SR comparison: {'PASS' if historical_sr_passed else 'FAIL'}")
del ff
return 0 if result["passed"] else 1
return 0
if __name__ == "__main__":
+43
View File
@@ -0,0 +1,43 @@
"""Convert Legacy D2Q9 storage and render raw lattice vorticity."""
from __future__ import annotations
import os
import sys
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
_CELERIS_SRC = os.path.join(_REPO, "CelerisLab", "src")
if _CELERIS_SRC not in sys.path:
sys.path.insert(0, _CELERIS_SRC)
def compute_vorticity(*args, **kwargs):
"""Lazily import rendering so CPU-only velocity tests need no CUDA stack."""
from CelerisLab.common.render import compute_vorticity as implementation
return implementation(*args, **kwargs)
def render_vorticity_field(*args, **kwargs):
from CelerisLab.common.render import render_vorticity_field as implementation
return implementation(*args, **kwargs)
def ddf_to_velocity(ddf, flags, nx, ny, rho_ref=1.0):
"""Decode Legacy physical velocity q/RHO_ref using exact solver flags."""
from drl_pinball.acquisition import decode_legacy_physical_velocity
fields = decode_legacy_physical_velocity(
ddf, flags=flags, nx=nx, ny=ny, rho_ref=rho_ref,
)
return fields["ux"], fields["uy"]
def render_final(flow_field, out_path):
shape = tuple(flow_field.FIELD_SHAPE)
if len(shape) != 3:
raise ValueError(f"Legacy FlowField.FIELD_SHAPE must be (nx, ny, nz), got {shape}")
nx, ny, nz = shape
if nz != 1:
raise ValueError(f"Legacy vorticity rendering requires nz=1, got FIELD_SHAPE={shape}")
flow_field.get_ddf()
flags_xy = np.asarray(flow_field.completed_flags_xy())
ux, uy = ddf_to_velocity(flow_field.ddf.copy(), flags_xy, nx, ny)
vort = compute_vorticity(ux, uy)
nonfluid = ((flags_xy & np.uint8(0b00000001)) == 0).T
vort = vort.copy(); vort[nonfluid] = 0.0
return render_vorticity_field(vort, nx=nx, ny=ny, out_path=out_path, vmin=-0.001, vmax=0.001)
+116
View File
@@ -0,0 +1,116 @@
"""Generic CLI for canonical LegacyCelerisLab reproductions."""
from __future__ import annotations
import argparse, json, os, shutil
from collections import deque
import numpy as np
from .cases import CASES, get_case
from .core.dtw_metrics import gen_target_states_at
from .core.legacy_env_builder import FIFO_LEN, DATA_TYPE, U0, build_illusion, build_karman_cloak, build_steady_cloak, build_vortex
from .core.model_loader import load_model, model_path
from .metrics import (aggregate_rewards, frozen_reference_comparison, native_similarity,
reward_terms, sha256_file, target_scale_cycle_similarity)
from .render import render_final
from .runtime import (load_policy_norm, policy_observation, reset_runtime,
run_historical_interval)
def _json(path, value):
with open(path, "w") as handle: json.dump(value, handle, indent=2)
def _prepare_case_output(out, metrics_only):
if metrics_only:
return
_prepare_case_output(out, metrics_only)
def _write_case_output(out, metrics_only, ff, sensors, forces, actions, rewards,
target, saved, metrics, norm_document):
if metrics_only:
return
np.savez_compressed(
os.path.join(out, "signals.npz"), sensors=sensors, forces=forces, actions=actions,
rewards=np.asarray([row["reward"] for row in rewards], dtype=np.float32),
reward_terms=np.asarray([[row[key] for key in ("reward_cd", "reward_cl", "native_legacy_dtw")]
for row in rewards], dtype=np.float32),
target=np.asarray(target, dtype=np.float32),
)
np.savez_compressed(os.path.join(out, "reset_contract.npz"), saved_fifo=saved)
_json(os.path.join(out, "metrics.json"), metrics)
_json(os.path.join(out, "norm.json"), norm_document)
render_final(ff, os.path.join(out, "final_vorticity.png"))
def _build(case, device):
if case.scene == "karman": return build_karman_cloak(device, case.re_code, sample_interval=case.sample_interval)
if case.scene == "illusion": return build_illusion(device, target_diameter_L=case.target_radius_l, sample_interval=case.sample_interval)
if case.scene == "vortex": return build_vortex(device, vortex_type=case.vortex_type)
if case.scene == "steady": return build_steady_cloak(device)
raise AssertionError(case.scene)
def run_case(case_name, device=0, steps=None, output_root=None, metrics_only=False):
case = get_case(case_name)
out=os.path.join(output_root or os.path.join(os.path.dirname(__file__), "output"), case.name)
repo=os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
reference_case=case.reference_case or case.name
frozen_dir=os.path.join(repo,"src","SR_analysis","data",case.scene,reference_case)
policy_norm_path=os.path.join(frozen_dir,"norm.json")
# Fail closed before output deletion or GPU initialization.
policy_norm=load_policy_norm(policy_norm_path) if case.model is not None else None
if os.path.isdir(out):
shutil.rmtree(out)
os.makedirs(out)
data = _build(case, device); ff=data["flow_field"]; recomputed_norm=data["norm"]; target=data["target_states"]
if policy_norm is None:
policy_norm=recomputed_norm
n=steps or case.steps; saved=np.asarray(recomputed_norm["save_states"], dtype=np.float32)
fifo, obs = reset_runtime(ff, saved, 14 if case.scene == "illusion" else 12)
if case.scene == "steady":
command=np.zeros(6,dtype=DATA_TYPE); command[3:]=np.asarray(case.action_bias)*U0
actions=np.tile(np.asarray(case.action_bias,dtype=np.float32),(n,1)); rewards=[]; rows=[]
for _ in range(n): run_historical_interval(ff, case.sample_interval,command); rows.append(ff.obs.copy()[:12])
else:
resolved_model_path=model_path(case.model)
model=load_model(case.model); rows=[]; actions=[]; rewards=[]; harmonics=data.get("target_harmonics")
# CustomEnv.reset returns exact zeros; first policy call must see no derived raw observation.
for step in range(n):
action,_=model.predict(obs,deterministic=True); action=np.asarray(action,dtype=np.float32).reshape(3); actions.append(action)
command=np.zeros(data["config"]["n_obj_total"],dtype=DATA_TYPE); command[-3:]=(action*case.action_scale+np.asarray(case.action_bias))*U0
ff.context.push()
try: run_historical_interval(ff, case.sample_interval,command)
finally: ff.context.pop()
raw=ff.obs.copy()[slice(*data["config"]["obs_slice"])]; fifo.append(raw); rows.append(raw)
terms=reward_terms(case,target,harmonics,np.asarray(fifo),policy_norm["force_norm_fact"],step); rewards.append(terms)
target_force=gen_target_states_at(step+1,harmonics)[:2] if harmonics is not None else None
obs=policy_observation(raw,policy_norm,target_force=target_force)
actions=np.asarray(actions,dtype=np.float32)
rows=np.asarray(rows,dtype=np.float32); sensors=rows[:,:6]; forces=rows[:,6:12]
native=float(native_similarity(case,target,np.asarray(fifo),n-1)) if case.scene != "steady" else None
offline=target_scale_cycle_similarity(target,sensors,case.conv_len,target_start=2 if case.scene=="illusion" else 0) if case.scene != "steady" else None
reference_path=os.path.join(frozen_dir,"controlled.npz")
metrics={
"case":case.name,
"target_native_legacy_dtw":native,
"normalized_target_scale_dtw_offline":offline,
"reward_summary":aggregate_rewards(rewards,case.conv_len),
"model":None if case.model is None else {"name":case.model,"path":os.path.abspath(resolved_model_path),"sha256":sha256_file(resolved_model_path)},
"sample_interval":case.sample_interval,
"action_scale":case.action_scale,
"action_bias":list(case.action_bias),
"re_code":case.re_code,
"frozen_reference_comparison":frozen_reference_comparison(reference_path,sensors,forces,actions,case.conv_len),
"reward_definition":"original CustomEnv formula" if case.scene != "steady" else None,
"normalization":{
"policy_source":"frozen_active_sr" if case.model is not None else "recomputed_no_policy",
"policy_path":os.path.abspath(policy_norm_path) if case.model is not None else None,
"policy_sha256":sha256_file(policy_norm_path) if case.model is not None else None,
"recomputed_this_run":True,
},
}
norm_document={
"policy":{k:(v.tolist() if isinstance(v,np.ndarray) else v) for k,v in policy_norm.items() if k!="save_states"},
"recomputed":{k:(v.tolist() if isinstance(v,np.ndarray) else v) for k,v in recomputed_norm.items() if k!="save_states"},
}
_write_case_output(out,metrics_only,ff,sensors,forces,actions,rewards,target,saved,metrics,norm_document)
del ff; return metrics
def main():
ap=argparse.ArgumentParser(); ap.add_argument("case",choices=sorted(set(CASES) | {"erase"})); ap.add_argument("--device",type=int,default=0); ap.add_argument("--steps",type=int); ap.add_argument("--output-root"); ap.add_argument("--metrics-only",action="store_true")
args=ap.parse_args(); print(json.dumps(run_case(args.case,args.device,args.steps,args.output_root,args.metrics_only),indent=2))
if __name__=="__main__": main()
@@ -1,78 +0,0 @@
#!/bin/bash
# legacy_test/run_all_legacy_tests.sh
#
# Sequential launcher for all Track A (Legacy Test) scripts.
# Each script uses LegacyCelerisLab which compiles CUDA kernels.
# A 60-second delay between tests prevents compilation conflicts.
#
# Usage:
# bash run_all_legacy_tests.sh [DEVICE_ID]
# DEVICE_ID defaults to 0 if not provided.
set -euo pipefail
DEVICE_ID="${1:-0}"
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
cd "$SCRIPT_DIR/../../.." # repo root
log() { echo "[$(date '+%H:%M:%S')] $*"; }
CONDA_ENV="pycuda_3_10"
DELAY=60
log "=== Legacy Test: Run All ==="
log "Device: $DEVICE_ID, Conda: $CONDA_ENV, Delay: ${DELAY}s"
# Array of (name, script_path)
declare -a TESTS=(
"Karman Re100:src/drl_pinball/legacy_test/test_karman_cloak_re100.py"
"Steady Cloak:src/drl_pinball/legacy_test/test_steady_cloak.py"
"Illusion 1L:src/drl_pinball/legacy_test/test_illusion_1L.py"
"Vortex Lamb:src/drl_pinball/legacy_test/test_vortex_lamb.py"
"Cross-Re Karman:src/drl_pinball/legacy_test/test_karman_cloak_crossre.py"
"Illusion 0.75L/1.5L:src/drl_pinball/legacy_test/test_illusion_remaining.py"
"Vortex Taylor:src/drl_pinball/legacy_test/test_vortex_taylor.py"
"Erase (experimental):src/drl_pinball/legacy_test/test_erase.py"
)
PASS_COUNT=0
FAIL_COUNT=0
declare -a FAILED_NAMES=()
for test_entry in "${TESTS[@]}"; do
name="${test_entry%%:*}"
script="${test_entry##*:}"
log ""
log "--- $name ---"
log "Running: conda run -n $CONDA_ENV python $script --device $DEVICE_ID"
if conda run -n "$CONDA_ENV" python "$script" --device "$DEVICE_ID"; then
log "[PASS] $name"
PASS_COUNT=$((PASS_COUNT + 1))
else
log "[FAIL] $name (exit code $?)"
FAIL_COUNT=$((FAIL_COUNT + 1))
FAILED_NAMES+=("$name")
fi
# Avoid CUDA compilation conflicts: wait 60s between tests
if [[ "$test_entry" != "${TESTS[-1]}" ]]; then
log "Waiting ${DELAY}s for CUDA compilation lock to clear..."
sleep "$DELAY"
fi
done
log ""
log "=== Summary ==="
log "Passed: $PASS_COUNT / $((PASS_COUNT + FAIL_COUNT))"
if [[ $FAIL_COUNT -gt 0 ]]; then
log "Failed tests:"
for fn in "${FAILED_NAMES[@]}"; do
log " - $fn"
done
exit 1
fi
log "All tests passed."
+6
View File
@@ -0,0 +1,6 @@
#!/usr/bin/env bash
set -euo pipefail
case_name="${1:?usage: run_legacy.sh CASE [DEVICE] [RUNNER_ARGS...]}"
device="${2:-0}"
if (( $# >= 2 )); then shift 2; else shift 1; fi
exec conda run -n pycuda_3_10 python -m src.drl_pinball.legacy_test "$case_name" --device "$device" "$@"
@@ -0,0 +1,209 @@
#!/usr/bin/env python3
"""Plot the validated V5 and Legacy three-role reproduction summaries."""
from __future__ import annotations
import argparse
import json
import os
import shutil
import tempfile
from pathlib import Path
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import numpy as np
L0 = 20.0
VORTICITY_LIMIT = 0.001
SENSOR_PAIRS = ((0, 1, "upper"), (2, 3, "center"), (4, 5, "lower"))
SENSOR_COLORS = ("#0072B2", "#D55E00", "#009E73")
V5_GROUPS = {
"ill_075L_seed43": "ill_075L", "ill_1L_seed43": "ill_1L",
"ill_15L_seed43": "ill_15L", "ill_2L_seed43": "ill_2L",
"kar_d075_seed44": "kar_d075", "kar_d15_seed45": "kar_d15",
"kar_d2_seed45": "kar_d2", "kar_re60_seed43": "kar_re60",
"kar_re100_seed41": "kar_re100", "kar_re100_seed42": "kar_re100",
"kar_re100_seed43": "kar_re100", "kar_re100_seed44": "kar_re100",
"kar_re100_seed45": "kar_re100", "kar_re200_seed43": "kar_re200",
"kar_re400_seed43": "kar_re400",
}
LEGACY_PERIODIC = (
"illusion_075L", "illusion_1L", "illusion_15L", "karman_re50",
"karman_re100", "karman_re200", "karman_re400",
)
LEGACY_VORTEX = (
"vortex_lamb_y000", "vortex_taylor_ym2L", "vortex_taylor_ym1L",
"vortex_taylor_y000", "vortex_taylor_yp1L", "vortex_taylor_yp2L",
)
LEGACY_GROUPS = LEGACY_PERIODIC + ("steady",) + LEGACY_VORTEX + ("erase",)
def fail(message: str) -> None:
raise RuntimeError(message)
def require_exact_groups(root: Path) -> None:
v5_controlled = {p.parent.name for p in (root / "v5").glob("*/controlled") if p.is_dir()}
if v5_controlled != set(V5_GROUPS):
fail(f"V5 controlled groups differ: missing={set(V5_GROUPS)-v5_controlled}, extra={v5_controlled-set(V5_GROUPS)}")
legacy = {p.name for p in (root / "legacy").iterdir() if p.is_dir()}
if legacy != set(LEGACY_GROUPS):
fail(f"Legacy groups differ: missing={set(LEGACY_GROUPS)-legacy}, extra={legacy-set(LEGACY_GROUPS)}")
def roles_for(root: Path, pipeline: str, group: str):
if pipeline == "v5":
case = V5_GROUPS[group]
return (("target", root / "v5" / case / "target"),
("controlled", root / "v5" / group / "controlled"),
("physical-zero", root / "v5" / case / "zero"))
middle = "constant" if group == "steady" else "controlled"
return (("target", root / "legacy" / group / "target"),
(middle, root / "legacy" / group / middle),
("physical-zero", root / "legacy" / group / "zero"))
def field_spec(pipeline: str, group: str, role_dir: Path):
if pipeline == "v5" or group in LEGACY_PERIODIC or (group == "erase" and role_dir.name != "target"):
return "phase_fields.npz", "target_phase", 0.0
if group in LEGACY_VORTEX:
return "event_fields.npz", "relative_offsets", -10
return "late_field.npz", "field_indices", None
def load_role(role_dir: Path, pipeline: str, group: str):
if not role_dir.is_dir():
fail(f"missing role directory: {role_dir}")
required = {"metadata.json", "timeseries.csv"}
names = {p.name for p in role_dir.iterdir() if p.is_file()}
if not required <= names:
fail(f"missing required files in {role_dir}: {required-names}")
field_name, selector_name, selector_value = field_spec(pipeline, group, role_dir)
field_candidates = names & {"phase_fields.npz", "event_fields.npz", "late_field.npz"}
if field_candidates != {field_name}:
fail(f"field files differ in {role_dir}: expected {field_name}, found {sorted(field_candidates)}")
metadata = json.loads((role_dir / "metadata.json").read_text())
units = metadata.get("units", {}).get("sensors", "raw sensor units")
data = np.genfromtxt(role_dir / "timeseries.csv", delimiter=",", names=True)
if data.ndim != 1 or data.size == 0:
fail(f"invalid timeseries shape in {role_dir}: {data.shape}")
expected_sensors = {f"sensors_{i}" for i in range(6)}
actual_sensors = {n for n in (data.dtype.names or ()) if n.startswith("sensors_")}
if actual_sensors != expected_sensors:
fail(f"sensor columns differ in {role_dir}: {actual_sensors}")
sensors = np.column_stack([data[f"sensors_{i}"] for i in range(6)])
if sensors.shape != (data.size, 6) or not np.isfinite(sensors).all():
fail(f"invalid sensor values in {role_dir}: {sensors.shape}")
field_path = role_dir / field_name
with np.load(field_path) as z:
required_keys = {"ux", "uy", selector_name}
if not required_keys <= set(z.files):
fail(f"missing field keys in {field_path}: {required_keys-set(z.files)}")
ux, uy = np.array(z["ux"]), np.array(z["uy"])
selector = np.array(z[selector_name])
if ux.ndim != 3 or ux.shape != uy.shape or ux.shape[0] != selector.shape[0]:
fail(f"invalid canonical field shapes in {field_path}: ux={ux.shape}, uy={uy.shape}, selector={selector.shape}")
expected_count = 8 if field_name == "phase_fields.npz" else 5 if field_name == "event_fields.npz" else 1
if ux.shape[0] != expected_count or selector.shape != (expected_count,):
fail(f"unexpected retained field count in {field_path}: {ux.shape[0]}")
if selector_value is not None and not np.isclose(selector[0], selector_value, atol=1e-12, rtol=0):
fail(f"slot 0 selector in {field_path} is {selector[0]}, expected {selector_value}")
if not np.isfinite(ux[0]).all() or not np.isfinite(uy[0]).all():
fail(f"nonfinite field in {field_path} slot 0")
omega = np.gradient(uy[0], axis=1) - np.gradient(ux[0], axis=0)
return omega, sensors, units, field_path, selector_name, selector[0]
def padded_limits(values):
lo, hi = float(np.min(values)), float(np.max(values))
if not np.isfinite([lo, hi]).all():
fail("nonfinite sensor limits")
span = hi - lo
pad = 0.05 * span if span > 0 else max(abs(lo) * 0.05, 1e-12)
return [lo - pad, hi + pad]
def plot_group(stage: Path, data_root: Path, pipeline: str, group: str):
roles = roles_for(data_root, pipeline, group)
loaded = [load_role(path, pipeline, group) for _, path in roles]
shapes = {x[0].shape for x in loaded}
if len(shapes) != 1:
fail(f"role field shapes differ for {pipeline}/{group}: {shapes}")
limit = VORTICITY_LIMIT
all_sensors = np.concatenate([x[1] for x in loaded], axis=0)
xlim = padded_limits(all_sensors[:, [0, 2, 4]])
ylim = padded_limits(all_sensors[:, [1, 3, 5]])
ny, nx = next(iter(shapes))
extent = (0, (nx - 1) / L0, 0, (ny - 1) / L0)
fig, axes = plt.subplots(2, 3, figsize=(18, 8), constrained_layout=True)
images = []
for col, ((role, _), (omega, sensors, units, _, _, _)) in enumerate(zip(roles, loaded)):
images.append(axes[0, col].imshow(omega, origin="lower", extent=extent, cmap="RdBu_r", vmin=-limit, vmax=limit, aspect="equal"))
axes[0, col].set_title(f"{pipeline.upper()} / {group}{role}")
axes[0, col].set_xlabel("x/L0")
axes[0, col].set_ylabel("y/L0")
axes[0, col].set_xlim(extent[:2]); axes[0, col].set_ylim(extent[2:])
for (u, v, label), color in zip(SENSOR_PAIRS, SENSOR_COLORS):
axes[1, col].plot(sensors[:, u], sensors[:, v], color=color, lw=1.1, alpha=0.85, label=label)
axes[1, col].axhline(0, color="0.35", lw=0.7); axes[1, col].axvline(0, color="0.35", lw=0.7)
axes[1, col].grid(True, alpha=0.25)
axes[1, col].set_xlim(xlim); axes[1, col].set_ylim(ylim)
axes[1, col].set_xlabel(f"sensor u ({units})"); axes[1, col].set_ylabel(f"sensor v ({units})")
axes[1, col].set_box_aspect(0.8)
axes[1, col].legend(loc="best", frameon=False)
cbar = fig.colorbar(images[0], ax=axes[0, :], orientation="vertical", shrink=0.92, pad=0.015)
cbar.set_label(r"$\omega_z$ (lattice$^{-1}$)")
output = stage / pipeline / f"{group}.png"
output.parent.mkdir(parents=True, exist_ok=True)
fig.savefig(output, dpi=160)
plt.close(fig)
return {
"pipeline": pipeline, "group": group, "plot": f"{pipeline}/{group}.png",
"roles": [{"role": role, "source": str(path), "field_file": str(item[3]),
"field_slot": 0, "selector": item[4], "selector_value": float(item[5])}
for (role, path), item in zip(roles, loaded)],
"vorticity_limit": limit, "sensor_limits": {"u": xlim, "v": ylim},
}
def main() -> None:
here = Path(__file__).resolve().parent
parser = argparse.ArgumentParser(description=__doc__)
parser.add_argument("--data-root", type=Path, default=here / "data" / "reproduction")
parser.add_argument("--output-root", type=Path, default=here / "data" / "reproduction_plots")
args = parser.parse_args()
data_root, output_root = args.data_root.resolve(), args.output_root.absolute()
require_exact_groups(data_root)
if output_root.is_symlink():
fail(f"refusing symlink output root: {output_root}")
output_root.parent.mkdir(parents=True, exist_ok=True)
stage = Path(tempfile.mkdtemp(prefix=f".{output_root.name}.staging-", dir=output_root.parent))
backup = output_root.with_name(f".{output_root.name}.old-{os.getpid()}")
try:
entries = [plot_group(stage, data_root, "v5", group) for group in sorted(V5_GROUPS)]
entries += [plot_group(stage, data_root, "legacy", group) for group in LEGACY_GROUPS]
if len(entries) != 30:
fail(f"expected 30 plots, got {len(entries)}")
(stage / "manifest.json").write_text(json.dumps({
"schema": "drl-pinball-reproduction-plots-v1",
"plot_count": 30,
"vorticity_limit": VORTICITY_LIMIT,
"vorticity_contract": "fixed symmetric [-0.001, +0.001] for every vorticity panel",
"plots": entries,
}, indent=2) + "\n")
if output_root.exists():
os.replace(output_root, backup)
os.replace(stage, output_root)
if backup.exists():
shutil.rmtree(backup)
except Exception:
shutil.rmtree(stage, ignore_errors=True)
if backup.exists() and not output_root.exists():
os.replace(backup, output_root)
raise
print(f"wrote 30 plots and manifest to {output_root}")
if __name__ == "__main__":
main()
-69
View File
@@ -1,69 +0,0 @@
# Reproduce (Track B)
Re-runs legacy PPO models on the **new CelerisLab** solver (v0.5.1) and
compares against SR_analysis reference data to quantify solver differences.
## Quick Start
```bash
# Phase 2: Open-loop target-signal validation (isolates CFD diffs)
conda run -n pycuda_3_10 python src/drl_pinball/reproduce/phase2_open_loop.py --device 2 --scene karman
# Phase 3: DRL inference with legacy-compatible config
conda run -n pycuda_3_10 python src/drl_pinball/reproduce/phase3_reproduce.py --device 2 --scene karman
# Run both phases:
bash src/drl_pinball/reproduce/run_all_reproduce_tests.sh 2
```
## Directory
```
reproduce/
├── README.md # This file
├── REPRODUCE_KNOWLEDGE.md # Comprehensive knowledge base (bugs, API diffs, findings)
├── core/
│ ├── action_wrapper.py # Action EMA + omega conversion (sign-corrected)
│ ├── obs_normalizer.py # Norm computation (exact legacy formulas)
│ ├── dtw_metrics.py # DTW similarity + harmonics analysis
│ ├── open_loop_validator.py # Phase 2: compare new CFD targets vs legacy ref
│ └── drl_comparator.py # Phase 3: compare DRL output vs SR_analysis ref
├── configs/
│ ├── scene_params.py # All scene parameter definitions
│ └── model_inventory.py # PPO model registry + loading
├── phase2_open_loop.py # Open-loop target validation (5 scenes)
├── phase3_reproduce.py # DRL inference with legacy-compat config
├── run_all_cases.py # (Legacy) old reproduce runner — superseded by phase3
├── run_illusion_vortex.py # (Legacy) old illusion/vortex runner — superseded
├── run_all_reproduce_tests.sh # Sequential launcher
└── output/
├── phase2_validation/ # Phase 2 comparison results
└── phase3/ # Phase 3 DRL inference results
```
## Config
The legacy-compatible config at `configs/config_lbm_pinball_legacy_compat.json`
uses **regularized inlet** with `regularized_neq_damp: 1.0`, matching the legacy
NBB (Non-Equilibrium Bounce-Back) formula: `f = feq_target + (f_neb - feq_neb)`.
This is the primary fix over the original `config_lbm_pinball.json` (which used
`zou_he_local` inlet, a fundamentally different numerical scheme).
## Key Results
| Scene | Legacy DTW | New CFD DTW (old) | New CFD DTW (fixed) |
|-------|:----------:|:-----------------:|:-------------------:|
| Karman Re100 | 0.975 | 0.916 | **0.943** |
| Vortex Lamb | 0.968 | 0.955 | **0.970** |
| Vortex Taylor | 0.996 | 0.979 | **0.994** |
The inlet scheme fix closed most of the gap. The remaining ~3% is attributable
to the ghost-source vs inline BC architectural difference.
## Known Limitations
- **Illusion**: S_DIM=14 with harmonics-derived target forces is more sensitive
to run-to-run CFD variability than the S_DIM=12 scenes
- **Steady Cloak**: Open-loop, no DRL — DTW comparison not applicable
- **Erase**: Incomplete training, no reference for comparison
@@ -1,300 +0,0 @@
# Reproduction Knowledge Document
> **Purpose**: Complete record of all experience, pitfalls, and findings from reproducing legacy DRL pinball control results on the new CelerisLab CFD solver.
> **Date**: 2026-06-21 (original reproduce), updated 2026-07-12 (inlet fix verified)
> **Next step**: Train new PPO models from scratch on the new CelerisLab solver. See `phase2_open_loop.py` (open-loop CFD validation) and `phase3_reproduce.py` (DRL inference with legacy-compat config) for the current reproduce pipeline.
---
## 1. Overall Picture
### What was attempted
Take old PPO models (trained on `LegacyCelerisLab/FlowField`) and run them on the new `CelerisLab/Simulation` solver, comparing sensor signals, forces, normalized observations, DTW similarity, and reward values against `SR_analysis` reference data.
### Conclusion
**Old PPO models cannot be perfectly transferred to the new CelerisLab solver.** Two independent LBM implementations produce systematically different plant dynamics (especially in vortex-dominated wake regions). Retraining on the new solver is the only complete solution.
### Best results achieved (Karman Cloak, legacy norm)
| Metric | Our result | Reference | Coverage |
|--------|:----------:|:---------:|:--------:|
| Action correlation (aF/aB/aT) | +0.83 / +0.81 / +0.63 | — | Good |
| Front_fy correlation | +0.90 | — | Excellent |
| DTW similarity | 0.916 | 0.954 | 96% |
| Sensor correlation | 0.69 ~ 0.72 | — | Moderate |
| Reward | 0.412 | 0.644 | 64% |
### All scenes summary
| Scene | DTW similarity | Action corr | Key issue |
|-------|:--------------:|:-----------:|-----------|
| Steady Cloak | N/A (open-loop) | N/A | Fully works |
| Karman Cloak (legacy norm) | 0.916 | 0.63-0.83 | Rear fy std 6x larger |
| Karman Cloak (new norm) | poor | -0.07~-0.04 | Uncorrelated actions |
| Illusion 1L (legacy norm) | 0.912 | 0.29-0.52 | Moderate correlation |
| Illusion 0.75L (legacy norm) | 0.926 | -0.32~+0.50 | Rear action sign inverted |
| Illusion 1.5L (legacy norm) | 0.894 | 0.08-0.15 | Low correlation |
| Vortex Lamb (legacy norm) | 0.955 | -0.22~-0.26 | Norm 16% diff |
| Vortex Taylor (legacy norm) | 0.979 | 0.01-0.39 | Norm 60% diff |
---
## 2. Critical API Differences
### 2.1 Omega: Surface Velocity vs Angular Velocity + Sign Inversion
**This was the most impactful bug.**
| Solver | Kernel code | Meaning |
|--------|------------|---------|
| **Legacy** `FlowField` | `Uw = action[id_obj] * (y_c - y) / radius` | action = tangential surface velocity |
| **New** `Simulation` | `Uw = -omega * ry` | `set_body(id, omega=val)` = angular velocity |
The new CelerisLab kernel has `Uw = -omega * ry`**the leading minus sign means `omega > 0` produces CW rotation**, opposite to the documented convention. The old `FlowField` has no minus sign.
**Corrected conversion**:
```python
omega = -surface_vel / radius
```
Where `surface_vel = (action_norm * scale + bias) * U0` and `radius = L0/2 = 10`.
### 2.2 Sensor Area Normalization
| Solver | Physical meaning | Value |
|--------|-----------------|-------|
| **Legacy** `obs[i]` | `sum_cells(u) / steps` | No area division |
| **New** `read_sensor(id, normalize=True)` | `sum_cells(u) / (cell_count * steps)` | Area- AND time-averaged |
**Conversion**: `legacy_equiv = new_value * cell_count`
Sensor cell count for radius=5: **78 cells** (verified empirically).
Key insight: if using `legacy norm` with new sensor values, you MUST convert sensors to old-equiv first:
```python
sensor_old_equiv = sim.read_sensor(id, normalize=True) * 78.0
```
### 2.3 Action Smoothing
Legacy `FlowField.run()` has built-in exponential smoothing with `weight=0.1`:
```
pinned = 0.9 * pinned + 0.1 * target
```
New CelerisLab has NO built-in smoothing. Must implement manually.
**Critical: EMA must be initialized to the bias values, NOT zeros.**
In legacy code, `self.action` retains the last bias value after the bias FIFO phase, so the first DRL inference step starts from a smoothed bias state:
```python
ema = EMA(weight=0.1)
ema.reset(bias_omega.copy()) # NOT ema.reset(np.zeros(3))
```
### 2.4 Snapshot/Restore Timing
**The snapshot must be taken AFTER the bias FIFO, not before.**
Legacy `cfd_interface.py` `add_pinball()` does:
1. Stabilize pinball (zero action)
2. `get_ddf()` + `save_ddf()` — saves zero-action state temporarily for norm
3. Collect norm, apply_bias, then does `get_ddf()` + `save_ddf()` again — **overwrites DDF with bias-state**
4. So `reset()``restore_ddf()` goes to **post-bias** state
Correct behavior:
```python
sim.run(warmup) # stabilize
sim.snapshot() # OPTIONAL: temp save for norm
# ... collect norm ...
sim.restore() # back to pre-bias state
# ... bias FIFO ...
sim.snapshot() # OVERWRITE: now at post-bias state
# DRL inference starts here:
sim.restore() # goes to post-bias state
```
### 2.5 Runtime Body Addition (sync_bodies)
```python
n_before = sim.bodies.count
sim.add_body("circle", center=(x, y, 0.0), radius=r) # returns -1 (staged)
sim.sync_bodies() # commit
pinball_ids = list(range(n_before, n_before + 3)) # real IDs discovered AFTER sync
```
### 2.6 PyCUDA + PyTorch Context Conflict
NVIDIA V100 has only one CUDA context per process. Solution: load PPO model on CPU:
```python
model = PPO.load(path, env=dummy_env, device="cpu") # CPU avoids context conflict
```
---
## 3. Body ID and Obs Layout
### 3.1 Karman Cloak (7 objects: dist_cyl + 3 sensors + 3 pinball)
**Add order**: dist_cyl(0) → sensor0(1) [y=+2L0] → sensor1(2) [y=0] → sensor2(3) [y=-2L0] → front(4) → top(5) [y=+0.75L0] → bottom(6) [y=-0.75L0]
**Legacy obs** (14 values):
```
[0:2] = dist_cyl fx, fy
[2:8] = sensor0_ux,uy, sensor1_ux,uy, sensor2_ux,uy
[8:14] = front_fx,fy, top_fx,fy, bottom_fx,fy
```
**Training uses** `obs[2:14]` → skips dist_cyl forces.
**Action -> body mapping**:
| Action index | Bias (U0 mult) | Body ID | Cylinder |
|:-----------:|:--------------:|:-------:|----------|
| 0 | 0 | 4 | Front |
| 1 | -4 | 5 | Top (rear, +y) |
| 2 | +4 | 6 | Bottom (rear, -y) |
**Normalized obs order**: `[forces(6)/norm, sensors(6)/norm]` clipped [-1,1].
### 3.2 Illusion/Vortex (6 objects: 3 sensors + 3 pinball)
**Add order**: sensor0(0) → sensor1(1) → sensor2(2) → front(3) → top(4) → bottom(5)
**Obs**: `[s0_ux,uy, s1_ux,uy, s2_ux,uy, front_fx,fy, top_fx,fy, bottom_fx,fy]`
Training uses `obs[0:12]` (all channels).
**Normalized obs order (S_DIM=12)**: same as Karman: forces first, sensors second.
**Normalized obs order (S_DIM=14, Illusion)**: forces(6) + sensors(6) + target_cd(1) + target_cl(1)
---
## 4. Scene Parameters
### 4.1 Scene settings
| Scene | S_DIM | Action scale | Action bias | SI | CONV_LEN | MaxSteps | Objects |
|-------|:----:|:------------:|:-----------:|:--:|:--------:|:--------:|:-------:|
| Steady Cloak | 12 | 8 | [0, -5.1, 5.1] | 800 | 30 | 500 | 6 |
| Karman Cloak | 12 | 8 | [0, -4, 4] | 800 | 30 | 500 | 7 |
| Illusion 0.75L | 14 | 8 | [0, -2, 2] | 400 | 36 | 500 | 6 |
| Illusion 1L | 14 | 8 | [0, -2, 2] | 600 | 36 | 500 | 6 |
| Illusion 1.5L | 14 | 8 | [0, -2, 2] | 800 | 36 | 500 | 6 |
| Vortex Lamb | 12 | 4 | [0, -4, 4] | 800 | 30 | 150 | 6 |
| Vortex Taylor | 12 | 4 | [0, -4, 4] | 800 | 30 | 150 | 6 |
### 4.2 Bias init actions (for FIFO initialization - may differ from DRL bias!)
| Scene | Init bias (surface_vel) | DRL bias |
|-------|:----------------------:|:--------:|
| Karman | [0, -4, 4] * U0 | Same |
| Illusion | **[0, -1, 1] * U0** | [0, -2, 2] * U0 |
| Vortex | [0, -4, 4] * U0 | Same |
### 4.3 Norm formulas
All scenes use:
- `force_norm_fact = 6 * max(|forces|)`
- `sens_deviation[i] = mean(sensor_i)`
- `sens_norm_fact[i] = 5 * max(|sensor_i - mean|)`
### 4.4 Omega conversion per scene
```python
# Karman, Steady, Illusion:
target_omega = -(action * 8.0 + bias) * U0 / 10.0
# Vortex:
target_omega = -(action * 4.0 + bias) * U0 / 10.0
```
---
## 5. All Bugs Found & Fixed
| # | Bug | Symptom | Fix |
|---|-----|---------|-----|
| **1** | Omega sign inverted | force fy means wrong sign, front_fy anti-correlated (-0.97) with ref actions open-loop | `omega = -surface_vel / R` |
| **2** | Snapshot before bias FIFO | Inference starts from zero-action flow, DRL sees wrong obs | `sim.snapshot()` after bias FIFO |
| **3** | EMA initialized to zeros | First DRL step applies 10% of target action = no effect | `ema.reset(bias_omega)` |
| **4** | Sensor not converted to old-equiv | Sensors have intrinsic area normalization, old norm doesn't expect it | `new_sensor * cell_count(78)` |
| **5** | Wrong body ID after sync_bodies | `set_body(-1, ...)` raises KeyError | `pinball_ids = list(range(n_before, n_before+3))` |
| **6** | `snapshot()` before `restore()` | `RuntimeError: No snapshot to restore` | Always snapshot after warmup |
---
## 6. File Structure (after cleanup)
```
reproduce/
├── __init__.py
├── REPRODUCE_KNOWLEDGE.md ← This file
├── core/
│ ├── __init__.py
│ ├── action_wrapper.py # EMA smoother + omega conversion
│ ├── obs_normalizer.py # Norm computation helpers
│ └── dtw_metrics.py # DTW similarity + harmonics
├── configs/
│ ├── __init__.py
│ ├── scene_params.py # All scene parameter definitions
│ └── model_inventory.py # Model loading (DummyEnv + Sin)
├── run_all_cases.py # Steady + Karman Cloak reproduction
├── run_illusion_vortex.py # Illusion + Vortex reproduction
└── output/
├── PHASE4_SUMMARY.md # Cross-scene comparison report
├── steady_cloak/ # Final steady cloak output
├── karman_cloak/ # Final Karman cloak output
├── illusion_1L/ # Final Illusion 1L output
├── illusion_075L/ # Final Illusion 0.75L output
├── illusion_15L/ # Final Illusion 1.5L output
├── vortex_lamb/ # Final Vortex Lamb output
└── vortex_taylor/ # Final Vortex Taylor output
```
---
## 7. Reference Data Locations
| Scene | Path |
|-------|------|
| Karman | `src/SR_analysis/data/karman/karman_re100/` |
| Illusion 1L | `src/SR_analysis/data/illusion/illusion_1L/` |
| Illusion 0.75L | `src/SR_analysis/data/illusion/illusion_0.75L/` |
| Illusion 1.5L | `src/SR_analysis/data/illusion/illusion_1.5L/` |
| Vortex Lamb | `src/SR_analysis/data/vortex/vortex_lamb/` |
| Vortex Taylor | `src/SR_analysis/data/vortex/vortex_taylor/` |
| Steady | `src/SR_analysis/data/steady/steady/` |
Each contains: `controlled.npz`, `uncontrolled.npz` (Karman only), `target.npz`, `norm.json`, `config.json`, `result.json`
---
## 8. Recommended Workflow for New Training
```
1. Create env with new CelerisLab Simulation
├─ Add bodies in legacy order
├─ initialize(), warmup (4*NX/U0 steps)
├─ Record target (150 x SI steps, sensors converted to old-equiv)
├─ Add pinball via sync_bodies()
├─ Warmup pinball (zero action)
├─ Collect zero-action FIFO → compute norm
├─ Bias FIFO (EMA smoother, proper bias values)
├─ sim.snapshot() ← AFTER bias FIFO!
└─ save save_states
2. Create gym.Env wrapper (observation_space, action_space, step, reset)
├─ step(): apply_action_smoothed → read_obs → normalize → build_obs → compute_reward
├─ reset(): sim.restore(), fifo = save_states.copy()
└─ reward: same formula as legacy (DTW + force terms)
3. Train PPO (same as legacy)
├─ Network: MlpPolicy, Sin activation, 64×64
├─ PPO hyperparams: lr=3e-4(actor)/4e-4(critic), n_steps=2048, batch_size=64
├─ Train for 500-1000 episodes (360-1500 timesteps per iter)
└─ Load PPO model on CPU (avoids CUDA context conflict)
4. Evaluate
├─ Run deterministic inference
├─ Compare reward curves across seeds
└─ Check DTW similarity against target
```
-13
View File
@@ -1,13 +0,0 @@
"""Reproduction package.
Purpose: Attempted reproduction of legacy DRL pinball control on new CelerisLab API.
Status: Legacy models do NOT work on new solver due to ~4% flow profile difference.
See REPRODUCE_KNOWLEDGE.md for full findings and guidance.
Contents:
- core/ : Verified utility modules (action wrapper, obs normalizer, DTW)
- configs/ : Scene parameters and model inventory
- inference/ : DRL inference script (Karman cloak)
- validation/ : Layer-by-layer validation against SR_analysis reference
"""
@@ -1 +0,0 @@
# configs/ — scene parameters and model inventory
@@ -1,136 +0,0 @@
"""Model inventory — metadata and loading utilities for all PPO models.
Provides:
- ModelInventory: registry of all pre-trained PPO models
- DummyEnv: minimal SB3-compatible env for PPO.load()
"""
from __future__ import annotations
import os
from typing import Any, Dict, Optional
import numpy as np
import gymnasium as gym
from gymnasium import spaces
from stable_baselines3 import PPO
from torch.nn import Module as TorchModule
_PROJECT_ROOT = os.path.abspath(
os.path.join(os.path.dirname(__file__), "..", "..", "..", "..")
)
_MODELS_DIR = os.path.join(_PROJECT_ROOT, "models")
# ---------------------------------------------------------------------------
# DummyEnv for model loading
# ---------------------------------------------------------------------------
class DummyEnv(gym.Env):
"""Minimal SB3-compatible environment for loading PPO models.
PPO.load() requires an env (or env's observation_space / action_space).
This dummy env provides the correct spaces without any CFD backend.
"""
def __init__(self, s_dim: int = 12, a_dim: int = 3):
super().__init__()
self.observation_space = spaces.Box(
low=-1.0, high=1.0, shape=(s_dim,), dtype=np.float32,
)
self.action_space = spaces.Box(
low=-1.0, high=1.0, shape=(a_dim,), dtype=np.float32,
)
def reset(self, *, seed=None, options=None):
return np.zeros(self.observation_space.shape, dtype=np.float32), {}
def step(self, action):
obs = np.zeros(self.observation_space.shape, dtype=np.float32)
return obs, 0.0, False, False, {}
# ---------------------------------------------------------------------------
# Custom Sin activation (must match training)
# ---------------------------------------------------------------------------
class Sin(TorchModule):
def __init__(self):
super().__init__()
def forward(self, x):
import torch
return torch.sin(x)
# ---------------------------------------------------------------------------
# Model metadata
# ---------------------------------------------------------------------------
MODEL_META: Dict[str, Dict[str, Any]] = {
# Old models: Karman cloak, cross-Re
"d1a3o12_re50": {"scene": "karman_cloak_re50", "s_dim": 12, "subdir": "old"},
"d1a3o12_re100": {"scene": "karman_cloak_re100", "s_dim": 12, "subdir": "old"},
"d1a3o12_re200": {"scene": "karman_cloak_re200", "s_dim": 12, "subdir": "old"},
"d1a3o12_re400": {"scene": "karman_cloak_re400", "s_dim": 12, "subdir": "old"},
# Vortex (transfer-learned from re100)
"vortex_lamb": {"scene": "vortex_lamb", "s_dim": 12, "subdir": "old"},
"vortex_taylor": {"scene": "vortex_taylor", "s_dim": 12, "subdir": "old"},
# Re-trained cloak 250326
"d1a3o12_250326": {"scene": "karman_cloak_re100", "s_dim": 12, "subdir": "250326"},
# No-offset cloak
"d0a3o12_250329_nooffset": {"scene": "karman_cloak_re100", "s_dim": 12, "subdir": "250329"},
# Reduced obs
"d1a3o12_250421_forces02": {"scene": "karman_cloak_re100", "s_dim": 3, "subdir": "250421"},
"d1a3o12_250421_torque+forces02": {"scene": "karman_cloak_re100", "s_dim": 5, "subdir": "250421"},
"d1a3o12_250421_torque+forces02+sens24": {"scene": "karman_cloak_re100", "s_dim": 9, "subdir": "250421"},
"d1a3o12_250421_torque+forces04+sens04": {"scene": "karman_cloak_re100", "s_dim": 5, "subdir": "250421"},
"d1a3o12_250421_torque+total_force": {"scene": "karman_cloak_re100", "s_dim": 5, "subdir": "250421"},
"d1a3o12_250421_total_force": {"scene": "karman_cloak_re100", "s_dim": 3, "subdir": "250421"},
# Illusion (250525)
"d1a3o14_250525_imit_075L_2U_400S": {"scene": "illusion_075L", "s_dim": 14, "subdir": "250525"},
"d1a3o14_250525_imit_1L_2U_600S": {"scene": "illusion_1L", "s_dim": 14, "subdir": "250525"},
"d1a3o14_250525_imit_15L_2U": {"scene": "illusion_15L", "s_dim": 14, "subdir": "250525"},
# Additional illusion variants
"d1a3o14_250525_imit_075L_2U": {"scene": "illusion_075L", "s_dim": 14, "subdir": "250525"},
"d1a3o14_250525_imit_1L_2U": {"scene": "illusion_1L", "s_dim": 14, "subdir": "250525"},
"d1a3o14_250525_imit_1L_2U_1": {"scene": "illusion_1L", "s_dim": 14, "subdir": "250525"},
"d1a3o14_250525_imit_1L_2U_1000S_08Vis": {"scene": "illusion_1L", "s_dim": 14, "subdir": "250525"},
"d1a3o14_250525_imit_1L_2U_800S_08Vis": {"scene": "illusion_1L", "s_dim": 14, "subdir": "250525"},
"d1a3o14_250525_imit_1L_2U_400S_02Vis": {"scene": "illusion_1L", "s_dim": 14, "subdir": "250525"},
"d1a3o14_250525_imit_075L_2U_1": {"scene": "illusion_075L", "s_dim": 14, "subdir": "250525"},
# Early illusion models (S_DIM=12, 1U variants)
"d1a3o12_250525_imit_075L_1U": {"scene": "illusion_075L", "s_dim": 12, "subdir": "250525"},
"d1a3o12_250525_imit_1L_1U": {"scene": "illusion_1L", "s_dim": 12, "subdir": "250525"},
"d1a3o12_250525_imit_1L_1U_trans": {"scene": "illusion_1L", "s_dim": 12, "subdir": "250525"},
# Erase models (for reference, not primary focus)
"d1a3o12_250729_250326_erase": {"scene": "karman_cloak_re100", "s_dim": 12, "subdir": "250729"},
"d1a3o12_250729_250326_erase_250804_20D_retrain2": {"scene": "karman_cloak_re100", "s_dim": 12, "subdir": "250729"},
"d1a3o12_250729_250326_erase_250804_20D_retrain3": {"scene": "karman_cloak_re100", "s_dim": 12, "subdir": "250729"},
"d1a3o12_250729_250326_cloak_800S_02Vis": {"scene": "karman_cloak_re100", "s_dim": 12, "subdir": "250729"},
}
class ModelInventory:
"""Registry for loading pre-trained PPO models."""
def __init__(self):
self.models_dir = _MODELS_DIR
def get_model_path(self, name: str) -> str:
if name not in MODEL_META:
raise KeyError(f"Unknown model '{name}'")
meta = MODEL_META[name]
return os.path.join(self.models_dir, meta["subdir"], f"{name}.zip")
def load(self, name: str, device: str = "cuda:0") -> PPO:
"""Load a PPO model with correct observation space and Sin activation."""
if name not in MODEL_META:
raise KeyError(f"Unknown model '{name}'")
meta = MODEL_META[name]
dummy = DummyEnv(s_dim=meta["s_dim"])
path = self.get_model_path(name)
model = PPO.load(path, env=dummy, device=device)
return model
def list_models(self, scene: Optional[str] = None) -> list:
"""List model names, optionally filtered by scene name."""
if scene is None:
return sorted(MODEL_META.keys())
return sorted(k for k, v in MODEL_META.items() if v["scene"] == scene)
@@ -1,243 +0,0 @@
"""All scene parameters in one place.
Single source of truth for geometry, action scaling, norm formulas, and DRL settings.
Every scene family needed for reproduction is defined here.
Re convention:
- "re_code" uses reference length 2*D (=40 lattice units), matching model file naming.
- nu = U0 * (2*D) / re_code
- Physical Re_D = re_code / 2 (uses single cylinder diameter D=20)
"""
from __future__ import annotations
from typing import Any, Dict
import numpy as np
# ---------------------------------------------------------------------------
# Physics constants (must match config_lbm_pinball.json)
# ---------------------------------------------------------------------------
U0 = 0.01 # inlet centre velocity (lattice units)
L0 = 20.0 # base length unit = 1 cylinder diameter in lattice
D_CYL = 20.0 # single cylinder diameter
D_REF = 40.0 # reference length for code Re = 2*D
NX = 1280
NY = 512
CENTER_Y = (NY - 1) / 2.0 # 255.5
CFG_PATH = "configs/config_lbm_pinball.json"
def nu_from_re(re_code: float) -> float:
"""Kinematic viscosity from code Reynolds number."""
return U0 * D_REF / re_code
# ---------------------------------------------------------------------------
# Scene definitions
# ---------------------------------------------------------------------------
SCENES: Dict[str, Dict[str, Any]] = {}
# -- Steady Cloak (clean inflow, no disturbance cylinder) --------------------
# Also serves as the base pinball-only env for Illusion inference and Vortex.
SCENES["steady_cloak_re100"] = {
"scene_id": "steady_cloak",
"model": "d1a3o12_re100",
"model_subdir": "old",
"re_code": 100,
"nu": nu_from_re(100),
"s_dim": 12,
"a_dim": 3,
"has_disturbance": False,
"pinball_front_x": 30.0 * L0, # 600
"pinball_rear_x": 31.3 * L0, # 626
"pinball_y_span": 0.75 * L0, # 15
"sensor_x": 40.0 * L0, # 800
"sensor_y_span": 2.0 * L0, # 40
"sensor_radius": L0 / 4, # 5
"pinball_radius": L0 / 2, # 10
"sample_interval": 800,
"action_scale": 8.0,
"action_bias": np.array([0.0, -4.0, 4.0], dtype=np.float32),
"u0": U0,
"fifo_len": 150,
"conv_len": 30,
"max_steps": 500,
"n_objects": 6,
"obs_slice": (0, 12),
"force_norm_formula": "6*max",
"sens_norm_factor": 5,
"target_type": "steady", # mean of clean channel
"warmup_steps_init": int(4 * NX / U0),
"warmup_steps_pinball": int(4 * NX / U0),
}
# -- Karman Cloak (upstream disturbance cylinder) ----------------------------
SCENES["karman_cloak_re100"] = {
"scene_id": "karman_cloak",
"model": "d1a3o12_re100",
"model_subdir": "old",
"re_code": 100,
"nu": nu_from_re(100),
"s_dim": 12,
"a_dim": 3,
"has_disturbance": True,
"dist_center_x": 10.0 * L0, # 200
"dist_radius": 1.0 * L0, # 20
"pinball_front_x": 30.0 * L0, # 600
"pinball_rear_x": 31.3 * L0, # 626
"pinball_y_span": 0.75 * L0, # 15
"sensor_x": 40.0 * L0, # 800
"sensor_y_span": 2.0 * L0, # 40
"sensor_radius": L0 / 4, # 5
"pinball_radius": L0 / 2, # 10
"sample_interval": 800,
"action_scale": 8.0,
"action_bias": np.array([0.0, -4.0, 4.0], dtype=np.float32),
"u0": U0,
"fifo_len": 150,
"conv_len": 30,
"max_steps": 500,
"n_objects": 7,
"obs_slice": (2, 14), # skip dist_cyl forces
"force_norm_formula": "6*max",
"sens_norm_factor": 5,
"target_type": "periodic",
"warmup_steps_dist": int(4 * NX / U0),
"warmup_steps_pinball": int(4 * NX / U0),
}
# Also define the karman scenes at other Re for completeness
for re_code, name in [(50, "re50"), (200, "re200"), (400, "re400")]:
key = f"karman_cloak_{name}"
SCENES[key] = dict(SCENES["karman_cloak_re100"])
SCENES[key].update({
"model": f"d1a3o12_{name}",
"model_subdir": "old",
"re_code": re_code,
"nu": nu_from_re(re_code),
"scene_id": f"karman_cloak_{name}",
})
# -- Illusion (three target diameters) ---------------------------------------
# NOTE: The positions below (sensors at 40*L0, pinball at 30/31.3*L0) are the
# "unified" inference geometry used by CCD_analysis. The actual training
# geometry (legacy_env_imit.py) used sensors at 30*L0 and pinball at 19/20.3*L0.
# phase3_reproduce.py and legacy_test scripts use the TRAINING positions.
def _illusion_base() -> Dict[str, Any]:
return {
"scene_id": "illusion",
"re_code": 100,
"nu": nu_from_re(100),
"s_dim": 14,
"a_dim": 3,
"has_disturbance": False,
"target_center_x": 31.0 * L0, # 620 (inference position)
"pinball_front_x": 30.0 * L0, # 600 (standard inference pinball)
"pinball_rear_x": 31.3 * L0, # 626
"pinball_y_span": 0.75 * L0,
"sensor_x": 40.0 * L0, # 800 (inference, not training!)
"sensor_y_span": 2.0 * L0,
"sensor_radius": L0 / 4,
"pinball_radius": L0 / 2,
"action_scale": 8.0,
"action_bias": np.array([0.0, -2.0, 2.0], dtype=np.float32),
"u0": U0,
"fifo_len": 150,
"conv_len": 36, # illusion uses CONV_LEN=36
"max_steps": 500,
"n_objects": 6,
"obs_slice": (0, 12),
"force_norm_formula": "6*max",
"sens_norm_factor": 5,
"target_type": "harmonics",
"n_harmonics": 5,
"warmup_steps_target": int(4 * NX / U0),
"warmup_steps_pinball": int(4 * NX / U0),
}
illusion_entries = [
("illusion_075L", {
"model": "d1a3o14_250525_imit_075L_2U_400S",
"model_subdir": "250525",
"target_diameter": 0.75 * L0, # 15
"sample_interval": 400,
}),
("illusion_1L", {
"model": "d1a3o14_250525_imit_1L_2U_600S",
"model_subdir": "250525",
"target_diameter": 1.0 * L0, # 20
"sample_interval": 600,
}),
("illusion_15L", {
"model": "d1a3o14_250525_imit_15L_2U",
"model_subdir": "250525",
"target_diameter": 1.5 * L0, # 30
"sample_interval": 800,
}),
]
for key, overrides in illusion_entries:
base = _illusion_base()
base.update(overrides)
SCENES[key] = base
# -- Vortex (Lamb dipole and Taylor monopole) --------------------------------
def _vortex_base() -> Dict[str, Any]:
return {
"scene_id": "vortex",
"re_code": 100,
"nu": nu_from_re(100),
"s_dim": 12,
"a_dim": 3,
"has_disturbance": False,
"pinball_front_x": 30.0 * L0,
"pinball_rear_x": 31.3 * L0,
"pinball_y_span": 0.75 * L0,
"sensor_x": 40.0 * L0,
"sensor_y_span": 2.0 * L0,
"sensor_radius": L0 / 4,
"pinball_radius": L0 / 2,
"sample_interval": 800,
"action_scale": 4.0, # NOTE: scale=4 not 8
"action_bias": np.array([0.0, -4.0, 4.0], dtype=np.float32),
"u0": U0,
"fifo_len": 150,
"conv_len": 30,
"max_steps": 150, # transient!
"n_objects": 6,
"obs_slice": (0, 12),
"force_norm_formula": "6*max",
"sens_norm_factor": 5,
"target_type": "transient",
"vortex_center_x": 10.0 * L0, # target phase
"vortex_pinball_center_x": 15.0 * L0, # pinball phase
"vortex_radius": 2.0 * L0,
}
vortex_entries = [
("vortex_lamb", {
"model": "vortex_lamb",
"model_subdir": "old",
"vortex_type": "lamb",
"vortex_strength": 0.5 * U0,
}),
("vortex_taylor", {
"model": "vortex_taylor",
"model_subdir": "old",
"vortex_type": "taylor",
"vortex_strength": 0.03 * U0,
}),
]
for key, overrides in vortex_entries:
base = _vortex_base()
base.update(overrides)
SCENES[key] = base
def get_scene(name: str) -> Dict[str, Any]:
"""Look up a scene by name. Raises KeyError if not found."""
if name not in SCENES:
available = ", ".join(sorted(SCENES.keys()))
raise KeyError(f"Unknown scene '{name}'. Available: {available}")
return dict(SCENES[name]) # return a copy
@@ -1 +0,0 @@
# core/ — shared utilities for reproduction
@@ -1,118 +0,0 @@
"""Action smoothing and physical-unit conversion.
Mimics the legacy FlowField.run() built-in exponential smoothing:
action_pinned = (1 - weight) * action_pinned + weight * action_target
Two usage modes:
(A) DRL inference — convert normalized PPO output to omega:
raw_action = model.predict(obs)[0] # [-1, 1] normalized
smoothed = smoother(raw_action) # smoothed normalized
omega = norm_action_to_omega(smoothed, scale=8, bias=[0,-4,4])
for i, body_id in enumerate(pinball_ids):
sim.set_body(body_id, omega=omega[i])
(B) Bias FIFO — directly specify surface_vel, bypass scale/bias mapping:
bias_surf = np.array([0.0, -4.0, 4.0]) * U0 # surface velocity
bias_omega = surface_vel_to_omega(bias_surf)
ema = ActionSmoother(weight=0.1)
ema.reset(np.zeros(3)) # legacy: starts from zero
for _ in range(FIFO_LEN):
s = ema(bias_omega)
sim.set_body(fid, omega=s[0]); ...
sim.run(SI, zero_obs=True)
IMPORTANT: New CelerisLab kernel has Uw = -omega * ry.
The minus sign means omega > 0 produces CW rotation
(opposite to naive expectation). Verified against legacy.
omega = -surface_vel / radius
"""
from __future__ import annotations
from typing import Optional
import numpy as np
U0 = 0.01
RADIUS = 10.0 # pinball cylinder radius
class ActionSmoother:
"""Exponential moving-average action smoother.
Matches legacy ``FlowField.run()`` smoothing:
``pinned = (1 - weight) * pinned + weight * target``
Stateful across calls: call ``reset()`` to clear internal state.
Use ``reset(np.zeros(3))`` for bias FIFO (legacy starts from zero).
For DRL inference, reset to the bias-omega value before the episode.
"""
def __init__(self, weight: float = 0.1):
self.weight = float(weight)
self._smoothed: Optional[np.ndarray] = None
def __call__(self, target: np.ndarray) -> np.ndarray:
"""Apply exponential smoothing. Returns smoothed copy."""
t = np.asarray(target, dtype=np.float32)
if self._smoothed is None:
self._smoothed = t.copy()
else:
self._smoothed = (1.0 - self.weight) * self._smoothed + self.weight * t
return self._smoothed.copy()
def reset(self, value: Optional[np.ndarray] = None) -> None:
"""Reset smoother state. Value=None means cold-start (first call
will initialise from its argument). Pass np.zeros(3) for bias FIFO."""
if value is not None:
self._smoothed = np.asarray(value, dtype=np.float32).copy()
else:
self._smoothed = None
# ── Physical-unit conversion ───────────────────────────────────────────
def norm_action_to_omega(
action_norm: np.ndarray,
scale: float = 8.0,
bias: np.ndarray = None,
u0: float = U0,
radius: float = RADIUS,
) -> np.ndarray:
"""Convert PPO normalised action [-1, 1]^3 to angular velocity [lat-units].
surface_vel = (action_norm * scale + bias) * u0
omega = -surface_vel / radius (new CelerisLab sign convention)
"""
if bias is None:
bias = np.zeros(3, dtype=np.float32)
b = np.asarray(bias, dtype=np.float32)
surface_vel = (np.asarray(action_norm, dtype=np.float32) * scale + b) * u0
return -surface_vel / radius
def surface_vel_to_omega(
surface_vel: np.ndarray,
radius: float = RADIUS,
) -> np.ndarray:
"""Convert surface tangential velocity directly to angular velocity.
Use this for bias FIFO where you know the exact surface_vel (e.g.
bias_surf = [0, -4, 4] * U0) and don't want scale/bias remapping.
"""
return -np.asarray(surface_vel, dtype=np.float32) / radius
def omega_to_norm_action(
omega: np.ndarray,
scale: float = 8.0,
bias: np.ndarray = None,
u0: float = U0,
radius: float = RADIUS,
) -> np.ndarray:
"""Inverse of ``norm_action_to_omega`` — angular velocity to PPO action."""
if bias is None:
bias = np.zeros(3, dtype=np.float32)
b = np.asarray(bias, dtype=np.float32)
# omega = -surface_vel / R → surface_vel = -omega * R
surface_vel = -np.asarray(omega, dtype=np.float32) * radius
return np.clip((surface_vel / u0 - b) / scale, -1.0, 1.0)
@@ -1,58 +0,0 @@
# reproduce/core/drl_comparator.py
"""DRL inference comparison: reproduce output vs SR_analysis reference.
For each scene, loads the reproduce output (sensors/forces/actions) and
compares against SR_analysis reference controlled.npz.
"""
from __future__ import annotations
import json
import os
import sys
from typing import Dict, Optional
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", "..", ".."))
_SRC = os.path.join(_REPO, "src")
for p in [_REPO, _SRC]:
if p not in sys.path:
sys.path.insert(0, p)
from legacy_test.core.comparator import ( # noqa: E402
compare_scene,
pearson_corr, dtw_similarity, rms_error,
)
def compare_reproduce_output(
ref_dir: str,
output_dir: str,
label: str = "",
conv_len: int = 30,
) -> Dict:
"""Load reproduce output and compare against SR_analysis reference.
Args:
ref_dir: Path to SR_analysis scene directory.
output_dir: Path to reproduce output directory.
label: Scene label for printing.
conv_len: DTW convergence window length.
Returns:
dict with comparison metrics (same schema as compare_scene).
"""
signals_path = os.path.join(output_dir, "signals.npz")
if not os.path.isfile(signals_path):
raise FileNotFoundError(f"Reproduce output not found: {signals_path}")
data = np.load(signals_path)
sensors = np.asarray(data["sensors"], dtype=np.float32)
forces = np.asarray(data["forces"], dtype=np.float32)
actions = np.asarray(data["actions"], dtype=np.float32)
return compare_scene(
ref_dir, sensors, forces, actions,
conv_len=conv_len, label=label,
)
@@ -1,262 +0,0 @@
"""DTW-based similarity metrics — exact replica of legacy env computations.
Provides ``calc_lag`` (cross-correlation lag), ``calc_sim`` (DTW similarity),
and ``compute_similarity`` to match legacy env reward computation.
Legacy envs used:
- Karman cloak: lag from sensor1 Uy, then DTW on 6 sensor channels
- Illusion: lag from target[:,3] vs state[:,1], then DTW on 6 sensor channels (offset +2)
- Vortex: no lag, roll by current_step+1
- Erase: lag from force channels, uses enhanced calc_sim
"""
from __future__ import annotations
from typing import Optional
import numpy as np
def calc_lag(target: np.ndarray, state: np.ndarray) -> int:
"""Cross-correlation lag between target and state sequences.
Args:
target: shape ``(N,)`` reference signal.
state: shape ``(M,)`` observed signal.
Returns:
Integer lag (positive = state is ahead of target).
"""
t_mean = np.mean(target)
s_mean = np.mean(state)
correlation = np.correlate(target - t_mean, state - s_mean, mode="full")
lags = np.arange(-len(target) + 1, len(target))
return int(lags[np.argmax(correlation)])
def calc_dtw_sim(target: np.ndarray, state: np.ndarray) -> float:
"""Standard DTW similarity (used by cloak, illusion, vortex, reduce_obs).
Args:
target: shape ``(N,)`` reference.
state: shape ``(M,)`` observed.
Returns:
Similarity in [0, 1], where 1 = perfect match.
"""
n, m = len(target), len(state)
dtw = np.full((n + 1, m + 1), np.inf)
dtw[0, 0] = 0.0
for i in range(1, n + 1):
for j in range(1, m + 1):
cost = abs(float(target[i - 1]) - float(state[j - 1]))
last_min = min(dtw[i - 1, j], dtw[i, j - 1], dtw[i - 1, j - 1])
dtw[i, j] = cost + last_min
return float(1.0 - dtw[n, m] / float(n))
def calc_dtw_sim_enhanced(target: np.ndarray, state: np.ndarray) -> float:
"""Enhanced DTW similarity with amplitude-ratio and mean components.
Used by the legacy erase env. Combines:
- 80% standard DTW (max-cost normalised)
- 10% amplitude ratio (min_std/max_std)
- 10% mean similarity (1/(1 + diff/scale*10))
Args:
target: shape ``(N,)`` reference.
state: shape ``(M,)`` observed.
Returns:
Combined similarity in [0, 1].
"""
target_arr = np.asarray(target, dtype=np.float64)
state_arr = np.asarray(state, dtype=np.float64)
n, m = len(target_arr), len(state_arr)
# Amplitude ratio component
t_std = max(np.std(target_arr), 1e-8)
s_std = max(np.std(state_arr), 1e-8)
amplitude_ratio = float(min(t_std, s_std) / max(t_std, s_std))
# Mean similarity component
mean_diff = abs(np.mean(target_arr) - np.mean(state_arr))
max_scale = max(abs(np.mean(target_arr)), abs(np.mean(state_arr)), 1e-8)
mean_similarity = 1.0 / (1.0 + mean_diff / max_scale * 10.0)
# DTW with max-possible-cost normalisation
dtw = np.full((n + 1, m + 1), np.inf)
dtw[0, 0] = 0.0
for i in range(1, n + 1):
for j in range(1, m + 1):
cost = abs(target_arr[i - 1] - state_arr[j - 1])
last_min = min(dtw[i - 1, j], dtw[i, j - 1], dtw[i - 1, j - 1])
dtw[i, j] = cost + last_min
max_possible_cost = max(np.max(np.abs(target_arr)), np.max(np.abs(state_arr)), 1e-8)
dtw_distance = dtw[n, m] / (n * max_possible_cost)
dtw_sim = max(0.0, 1.0 - dtw_distance)
return float(0.8 * dtw_sim + 0.1 * amplitude_ratio + 0.1 * mean_similarity)
def compute_similarity_karman_cloak(
target_states: np.ndarray,
fifo_states: np.ndarray,
conv_len: int = 30,
) -> float:
"""Compute DTW similarity for Karman cloak (standard pattern).
Matches legacy code:
1. Compute lag from middle sensor (index 1) Uy component
2. For all 6 sensor channels, roll target by lag, compute DTW, average
Args:
target_states: shape ``(FIFO_LEN, 6)`` target sensor data.
fifo_states: shape ``(FIFO_LEN, 6)`` current FIFO sensor data.
conv_len: Convergence window length (default 30).
Returns:
Average similarity over 6 channels in [0, 1].
"""
target = np.asarray(target_states, dtype=np.float64)
state = np.asarray(fifo_states, dtype=np.float64)
id_sens = 1 # middle sensor
target_seq = target[conv_len:2 * conv_len, id_sens]
state_seq = state[-conv_len:, id_sens]
lag = calc_lag(target_seq, state_seq)
similarities = 0.0
for i in range(6):
t_seq = np.roll(target[:, i], -lag)[conv_len:2 * conv_len]
s_seq = state[-conv_len:, i]
similarities += calc_dtw_sim(t_seq, s_seq)
return float(similarities / 6.0)
def compute_similarity_vortex(
target_states: np.ndarray,
fifo_states: np.ndarray,
current_step: int,
conv_len: int = 30,
) -> float:
"""Compute DTW similarity for vortex (no lag, roll by current_step+1).
Matches legacy vortex env.
Args:
target_states: shape ``(FIFO_LEN, 6)`` target sensor data.
fifo_states: shape ``(FIFO_LEN, 6)`` current FIFO data.
current_step: The current simulation step index.
conv_len: Convergence window length (default 30).
Returns:
Average similarity over 6 channels.
"""
target = np.asarray(target_states, dtype=np.float64)
state = np.asarray(fifo_states, dtype=np.float64)
similarities = 0.0
for i in range(6):
t_seq = np.roll(target[-conv_len:, i], -current_step - 1)
s_seq = state[-conv_len:, i]
similarities += calc_dtw_sim(t_seq, s_seq)
return float(similarities / 6.0)
def compute_similarity_illusion(
target_states: np.ndarray,
fifo_states: np.ndarray,
conv_len: int = 36,
) -> float:
"""Compute DTW similarity for illusion.
Matches legacy imit env:
1. lag from target[:, id_sens+2] vs state[:, id_sens] (offset by 2)
2. For 6 channels, target uses [:, i+2] offset
Args:
target_states: shape ``(FIFO_LEN, 8)`` (2 force + 6 sensor channels).
fifo_states: shape ``(FIFO_LEN, 6)`` current FIFO (6 sensors only).
conv_len: Convergence window length (default 36).
Returns:
Average similarity over 6 channels.
"""
target = np.asarray(target_states, dtype=np.float64)
state = np.asarray(fifo_states, dtype=np.float64)
id_sens = 1
t_seq_ref = target[conv_len:2 * conv_len, id_sens + 2]
s_seq_ref = state[-conv_len:, id_sens]
lag = calc_lag(t_seq_ref, s_seq_ref)
similarities = 0.0
for i in range(6):
t_seq = np.roll(target[:, i + 2], -lag)[conv_len:2 * conv_len]
s_seq = state[-conv_len:, i]
similarities += calc_dtw_sim(t_seq, s_seq)
return float(similarities / 6.0)
# ---------------------------------------------------------------------------
# Harmonics analysis (used by illusion)
# ---------------------------------------------------------------------------
def analyze_harmonics(
states: np.ndarray,
n_harmonics: int = 5,
) -> list:
"""FFT-based harmonic analysis of multi-channel time series.
Matches legacy ``analyze_harmonics()``.
Args:
states: shape ``(N, D)`` time-series data.
n_harmonics: Number of harmonics to extract per channel.
Returns:
List of D dicts, each with keys:
dc: float (DC component)
amps: (n_harmonics,) array
freqs: (n_harmonics,) array
phases: (n_harmonics,) array
"""
N, D = states.shape
result = []
for d in range(D):
y = states[:, d]
fft_coef = np.fft.rfft(y)
freqs = np.fft.rfftfreq(N, d=1)
amps = 2.0 * np.abs(fft_coef) / N
phases = np.angle(fft_coef)
idx = np.argsort(amps[1:])[::-1][:n_harmonics] + 1
harmonics = {
"dc": float(np.real(fft_coef[0]) / N),
"amps": np.array(amps[idx], dtype=np.float32),
"freqs": np.array(freqs[idx], dtype=np.float32),
"phases": np.array(phases[idx], dtype=np.float32),
}
result.append(harmonics)
return result
def gen_target_states_at(t, harmonics) -> np.ndarray:
"""Reconstruct target state at time step t from harmonics.
Matches legacy ``gen_target_states_at()``.
Args:
t: Integer step index.
harmonics: Output from ``analyze_harmonics()``.
Returns:
shape ``(D,)`` reconstructed state vector.
"""
D = len(harmonics)
result = np.zeros(D, dtype=np.float32)
for d, h in enumerate(harmonics):
val = float(h["dc"])
for amp, freq, phase in zip(h["amps"], h["freqs"], h["phases"]):
val += amp * np.cos(2.0 * np.pi * freq * t + phase)
result[d] = val
return result
@@ -1,112 +0,0 @@
"""Observation normalization — exact replication of legacy norm procedures.
Legacy norm computation (from FIFO data):
- force_norm_fact = factor * max(|forces|)
- sens_deviation[i] = mean(sensor_i)
- sens_norm_fact[i] = factor * max(|sensor_i - sens_deviation[i]|)
Normalised observation:
- forces = raw_forces / force_norm_fact
- sens = (raw_sens - sens_deviation) / sens_norm_fact
- observation = clip(hstack([forces, sens]), -1, 1)
"""
from __future__ import annotations
from typing import Dict, Optional, Tuple
import numpy as np
def compute_norm(
fifo_states: np.ndarray,
force_norm_factor: float = 6.0,
sens_norm_factor: float = 5.0,
force_slice: Tuple[int, int] = (6, 12),
sens_slice: Tuple[int, int] = (0, 6),
) -> Dict[str, np.ndarray]:
"""Compute norm values from a FIFO of observations.
Exact replica of legacy env norm computation.
Args:
fifo_states: shape ``(FIFO_LEN, N_obs)`` array of raw observations.
force_norm_factor: Multiplier for force norm (6 for cloak, 100 for erase, etc.)
sens_norm_factor: Multiplier for sensor norm (5 for cloak, 10 for erase, etc.)
force_slice: Slice indices for forces within the obs array.
sens_slice: Slice indices for sensors within the obs array.
Returns:
dict with keys:
force_norm_fact: scalar float32
sens_deviation: (6,) float32
sens_norm_fact: (6,) float32
"""
arr = np.asarray(fifo_states, dtype=np.float64)
s_start, s_end = sens_slice
f_start, f_end = force_slice
forces = arr[:, f_start:f_end]
sensors = arr[:, s_start:s_end]
n_sens = sensors.shape[1]
force_norm_fact = np.float32(force_norm_factor * np.max(np.abs(forces)))
sens_deviation = np.mean(sensors, axis=0).astype(np.float32)
sens_norm_fact = np.zeros(n_sens, dtype=np.float32)
for i in range(n_sens):
deviation = np.max(np.abs(sensors[:, i] - sens_deviation[i]))
sens_norm_fact[i] = np.float32(sens_norm_factor * deviation)
return {
"force_norm_fact": force_norm_fact,
"sens_deviation": sens_deviation,
"sens_norm_fact": sens_norm_fact,
}
def normalize_observation(
raw_obs_slice: np.ndarray,
norm: Dict[str, np.ndarray],
force_slice: Tuple[int, int] = (6, 12),
sens_slice: Tuple[int, int] = (0, 6),
) -> np.ndarray:
"""Normalize a raw observation slice using pre-computed norm values.
Args:
raw_obs_slice: shape ``(N_obs,)`` raw observation values.
norm: dict with keys 'force_norm_fact', 'sens_deviation', 'sens_norm_fact'.
force_slice: Slice for forces within the obs array.
sens_slice: Slice for sensors within the obs array.
Returns:
shape ``(S_DIM,)`` normalized observation, clipped to [-1, 1].
"""
obs = np.asarray(raw_obs_slice, dtype=np.float32)
s_start, s_end = sens_slice
f_start, f_end = force_slice
forces = obs[f_start:f_end] / norm["force_norm_fact"]
sens = (obs[s_start:s_end] - norm["sens_deviation"]) / norm["sens_norm_fact"]
return np.clip(np.hstack([forces, sens]), -1.0, 1.0).astype(np.float32)
def load_norm(path: str) -> Dict[str, np.ndarray]:
"""Load norm values from a .npz file."""
data = np.load(path)
return {
"force_norm_fact": data["force_norm_fact"],
"sens_deviation": data["sens_deviation"],
"sens_norm_fact": data["sens_norm_fact"],
}
def save_norm(path: str, norm: Dict[str, np.ndarray]) -> None:
"""Save norm values to a .npz file."""
np.savez_compressed(
path,
force_norm_fact=norm["force_norm_fact"],
sens_deviation=norm["sens_deviation"],
sens_norm_fact=norm["sens_norm_fact"],
)
@@ -1,108 +0,0 @@
# reproduce/core/open_loop_validator.py
"""Open-loop comparison: target-recording phase on new CelerisLab vs legacy target.
For each scene, runs the target-recording phase on the new CelerisLab with
the legacy-compatible config, then compares directly against SR_analysis
reference target signals.
This isolates CFD differences before DRL is involved.
"""
from __future__ import annotations
import json
import os
import sys
from typing import Dict, Optional
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", "..", ".."))
_SRC = os.path.join(_REPO, "src")
for p in [_REPO, _SRC]:
if p not in sys.path:
sys.path.insert(0, p)
from legacy_test.core.dtw_metrics import calc_lag, calc_dtw_sim # noqa: E402
from legacy_test.core.io_helpers import load_reference_target # noqa: E402
def compare_target_signals(
ref_dir: str,
new_target: np.ndarray,
label: str = "",
conv_len: int = 30,
sensor_slice: slice = slice(0, 6),
) -> Dict:
"""Compare new CFD target signals against SR_analysis legacy target.
Args:
ref_dir: Path to SR_analysis scene directory.
new_target: (FIFO_LEN, N) new CFD target signals.
label: Scene label for printing.
conv_len: DTW convergence window.
sensor_slice: Which columns of new_target are sensor channels.
Returns:
dict with dtw_sim, per_channel_corr, rms_err, passed.
"""
ref = load_reference_target(ref_dir)
# Apply sensor slice if reference has more columns than new
if ref.shape[1] > new_target.shape[1]:
ref = ref[:, sensor_slice]
new = new_target[:, sensor_slice] if sensor_slice.stop <= new_target.shape[1] else new_target
n = min(ref.shape[0], new.shape[0])
ref = ref[:n]
new = new[:n]
n_ch = min(ref.shape[1], new.shape[1])
ch_corr = []
for i in range(n_ch):
r = ref[:, i]
g = new[:, i]
denom = np.sqrt(((r - r.mean())**2).sum() * ((g - g.mean())**2).sum())
ch_corr.append(float(((r - r.mean()) * (g - g.mean())).sum() / max(denom, 1e-12)))
# RMS error
rms = float(np.sqrt(np.mean((ref - new)**2)))
# DTW similarity (all channels)
sim_sum = 0.0
for i in range(n_ch):
ref_seq = ref[conv_len:2 * conv_len, i]
new_seq = new[-conv_len:, i]
sim_sum += calc_dtw_sim(ref_seq, new_seq)
dtw_sim = float(sim_sum / max(n_ch, 1))
# FFT peak comparison on first channel
ref_fft = np.abs(np.fft.rfft(ref[:, 0]))
new_fft = np.abs(np.fft.rfft(new[:, 0]))
freqs = np.fft.rfftfreq(n, d=1)
ref_peak = freqs[1:][np.argmax(ref_fft[1:])] if len(freqs) > 1 else 0
new_peak = freqs[1:][np.argmax(new_fft[1:])] if len(freqs) > 1 else 0
fft_ok = abs(ref_peak - new_peak) / max(abs(ref_peak), 1e-12) < 0.10 if abs(ref_peak) > 1e-12 else True
passed = dtw_sim > 0.90 and float(np.min(ch_corr if ch_corr else [1.0])) > 0.85
prefix = f"[{label}] " if label else ""
print(f"{prefix}Channel corr: {ch_corr}")
print(f"{prefix}DTW sim: {dtw_sim:.4f}, RMS err: {rms:.6f}")
print(f"{prefix}FFT peak: ref={ref_peak:.6f}, new={new_peak:.6f}, ok={fft_ok}")
print(f"{prefix}{'PASS' if passed else 'FAIL'}")
return {
"channel_corr": ch_corr,
"dtw_sim": float(dtw_sim),
"rms_err": float(rms),
"ref_fft_peak": float(ref_peak),
"new_fft_peak": float(new_peak),
"fft_ok": bool(fft_ok),
"passed": bool(passed),
}
def load_legacy_target(ref_dir: str) -> np.ndarray:
"""Load legacy target from SR_analysis data."""
return load_reference_target(ref_dir)
@@ -1,266 +0,0 @@
#!/usr/bin/env python3
# reproduce/phase2_open_loop.py
"""Phase 2: Open-loop target-signal validation.
Runs the target-recording phase on the new CelerisLab with the
legacy-compatible config (regularized inlet, NBB equivalent),
then compares against SR_analysis reference target signals.
This isolates CFD differences before DRL is involved.
Usage:
conda run -n pycuda_3_10 python phase2_open_loop.py --device 0
conda run -n pycuda_3_10 python phase2_open_loop.py --device 0 --scene karman
"""
from __future__ import annotations
import argparse
import json
import os
import sys
import time
import numpy as np
import pycuda.driver as cuda; cuda.init()
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
_SRC = os.path.join(_REPO, "src")
_DRL = os.path.join(_SRC, "drl_pinball")
for p in [_REPO, _SRC, _DRL]:
if p not in sys.path:
sys.path.insert(0, p)
from CelerisLab import Simulation # noqa: E402
from CelerisLab.lbm.initializers import add_vortex # noqa: E402
from reproduce.core.open_loop_validator import compare_target_signals # noqa: E402
# ---------------------------------------------------------------------------
# Config
# ---------------------------------------------------------------------------
CFG_PATH = "configs/config_lbm_pinball_legacy_compat.json"
L0 = 20.0
U0 = 0.01
NX = 1280
NY = 512
CENTER_Y = float(NY - 1) / 2.0
RADIUS = L0 / 2.0 # 10
FIFO_LEN = 150
SI = 800
WARMUP = int(4.0 * NX / U0)
REF_BASE = os.path.join(_SRC, "SR_analysis", "data")
OUT_BASE = os.path.join(os.path.dirname(__file__), "output", "phase2_validation")
def log(msg: str) -> None:
print(f"[{time.strftime('%H:%M:%S')}] {msg}", flush=True)
def get_cc(sim, sid):
nx, ny = sim.lbm_cfg.nx, sim.lbm_cfg.ny
cells_arr, _ = sim.bodies.get(sid).get_sensor_list(nx, ny)
return float(len(cells_arr))
def read_sensors_legacy(sim, sensor_ids, cc):
obs = []
for sid in sensor_ids:
s = sim.read_sensor(sid, normalize=True)
obs.extend([float(s[0]) * cc, float(s[1]) * cc])
return np.array(obs, dtype=np.float32)
# ---------------------------------------------------------------------------
# Karman target: dist-cyl + 3 sensors
# ---------------------------------------------------------------------------
def validate_karman(device_id: int, out_dir: str) -> dict:
log("=== Karman target validation ===")
ref_dir = os.path.join(REF_BASE, "karman", "karman_re100")
os.makedirs(out_dir, exist_ok=True)
sim = Simulation(lbm_config_path=CFG_PATH, device_id=device_id)
dist_id = sim.add_body("circle", center=(10.0 * L0, CENTER_Y, 0.0), radius=1.0 * L0)
sensor_ids = [
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y + 40.0, 0.0), radius=5.0),
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y, 0.0), radius=5.0),
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim.initialize()
sim.run(WARMUP, zero_obs=True)
cc = get_cc(sim, sensor_ids[0])
log(f" Sensor CC: {cc}")
target = np.zeros((FIFO_LEN, 6), dtype=np.float32)
for i in range(FIFO_LEN):
sim.run(SI, zero_obs=True)
obs = read_sensors_legacy(sim, sensor_ids, cc)
target[i] = obs
np.savez_compressed(os.path.join(out_dir, "target.npz"), target_states=target)
sim.close()
result = compare_target_signals(ref_dir, target, label="karman")
with open(os.path.join(out_dir, "result.json"), "w") as f:
json.dump(result, f, indent=2)
return result
# ---------------------------------------------------------------------------
# Steady channel target: 3 sensors only
# ---------------------------------------------------------------------------
def validate_steady(device_id: int, out_dir: str) -> dict:
log("=== Steady channel target validation ===")
ref_dir = os.path.join(REF_BASE, "steady", "steady")
os.makedirs(out_dir, exist_ok=True)
sim = Simulation(lbm_config_path=CFG_PATH, device_id=device_id)
sensor_ids = [
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y + 40.0, 0.0), radius=5.0),
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y, 0.0), radius=5.0),
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim.initialize()
sim.run(WARMUP, zero_obs=True)
cc = get_cc(sim, sensor_ids[0])
log(f" Sensor CC: {cc}")
target = np.zeros((FIFO_LEN, 6), dtype=np.float32)
for i in range(FIFO_LEN):
sim.run(SI, zero_obs=True)
target[i] = read_sensors_legacy(sim, sensor_ids, cc)
np.savez_compressed(os.path.join(out_dir, "target.npz"), target_states=target)
sim.close()
result = compare_target_signals(ref_dir, target, label="steady")
with open(os.path.join(out_dir, "result.json"), "w") as f:
json.dump(result, f, indent=2)
return result
# ---------------------------------------------------------------------------
# Illusion target: target cylinder + 3 sensors
# ---------------------------------------------------------------------------
def validate_illusion(device_id: int, out_dir: str, diam_L: float = 1.0) -> dict:
# diam_L: 0.75 → "illusion_0.75L", 1.0 → "illusion_1L", 1.5 → "illusion_1.5L"
if diam_L == int(diam_L):
label_suffix = str(int(diam_L))
else:
label_suffix = str(diam_L).rstrip('0')
label = f"illusion_{label_suffix}L"
log(f"=== {label} target validation ===")
ref_dir = os.path.join(REF_BASE, "illusion", label)
os.makedirs(out_dir, exist_ok=True)
sim = Simulation(lbm_config_path=CFG_PATH, device_id=device_id)
sim.add_body("circle", center=(20.0 * L0, CENTER_Y, 0.0), radius=diam_L * L0)
sensor_ids = [
sim.add_body("sensor", center=(30.0 * L0, CENTER_Y + 40.0, 0.0), radius=5.0),
sim.add_body("sensor", center=(30.0 * L0, CENTER_Y, 0.0), radius=5.0),
sim.add_body("sensor", center=(30.0 * L0, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim.initialize()
sim.run(WARMUP, zero_obs=True)
cc = get_cc(sim, sensor_ids[0])
log(f" Sensor CC: {cc}")
# Target: cyl_force(2) + sensors(6) = 8 channels
target = np.zeros((FIFO_LEN, 8), dtype=np.float32)
for i in range(FIFO_LEN):
sim.run(SI, zero_obs=True)
f = list(sim.read_force(0, normalize=True))
s = read_sensors_legacy(sim, sensor_ids, cc)
target[i] = np.hstack([f, s])
np.savez_compressed(os.path.join(out_dir, "target.npz"), target_states=target)
sim.close()
result = compare_target_signals(ref_dir, target, label=label, sensor_slice=slice(2, 8))
with open(os.path.join(out_dir, "result.json"), "w") as f:
json.dump(result, f, indent=2)
return result
# ---------------------------------------------------------------------------
# Vortex target: vortex + 3 sensors
# ---------------------------------------------------------------------------
def validate_vortex(device_id: int, out_dir: str, vortex_type: str = "lamb") -> dict:
log(f"=== Vortex {vortex_type} target validation ===")
ref_dir = os.path.join(REF_BASE, "vortex", f"vortex_{vortex_type}")
os.makedirs(out_dir, exist_ok=True)
strength = 0.5 * U0 if vortex_type == "lamb" else 0.03 * U0
sim = Simulation(lbm_config_path=CFG_PATH, device_id=device_id)
sensor_ids = [
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y + 40.0, 0.0), radius=5.0),
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y, 0.0), radius=5.0),
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim.initialize()
sim.run(WARMUP, zero_obs=True)
# Add vortex
sim.field.download_ddf()
add_vortex(sim.field, (10.0 * L0, CENTER_Y), 2.0 * L0, strength, vortex_type)
cc = get_cc(sim, sensor_ids[0])
log(f" Sensor CC: {cc}")
target = np.zeros((FIFO_LEN, 6), dtype=np.float32)
for i in range(FIFO_LEN):
sim.run(SI, zero_obs=True)
target[i] = read_sensors_legacy(sim, sensor_ids, cc)
np.savez_compressed(os.path.join(out_dir, "target.npz"), target_states=target)
sim.close()
result = compare_target_signals(ref_dir, target, label=f"vortex_{vortex_type}")
with open(os.path.join(out_dir, "result.json"), "w") as f:
json.dump(result, f, indent=2)
return result
# ---------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------
if __name__ == "__main__":
ap = argparse.ArgumentParser(description="Phase 2: Open-loop target validation")
ap.add_argument("--device", type=int, default=0, help="GPU device ID")
ap.add_argument("--scene", type=str, default="all",
help="Scene to validate: karman, steady, illusion_1L, vortex_lamb, vortex_taylor, all")
args = ap.parse_args()
results = {}
scenes = {
"karman": lambda: validate_karman(args.device, os.path.join(OUT_BASE, "karman")),
"steady": lambda: validate_steady(args.device, os.path.join(OUT_BASE, "steady")),
"illusion_1L": lambda: validate_illusion(args.device, os.path.join(OUT_BASE, "illusion_1L"), 1.0),
"vortex_lamb": lambda: validate_vortex(args.device, os.path.join(OUT_BASE, "vortex_lamb"), "lamb"),
"vortex_taylor": lambda: validate_vortex(args.device, os.path.join(OUT_BASE, "vortex_taylor"), "taylor"),
}
if args.scene == "all":
for name, func in scenes.items():
results[name] = func()
else:
for s in args.scene.split(","):
s = s.strip()
if s not in scenes:
log(f"Unknown scene: {s}")
continue
results[s] = scenes[s]()
# Summary
log("\n=== Open-loop validation summary ===")
all_pass = True
for name, r in results.items():
status = "PASS" if r["passed"] else "FAIL"
log(f" {name}: DTW={r['dtw_sim']:.4f}, corr={[f'{c:.3f}' for c in r['channel_corr']]} -> {status}")
if not r["passed"]:
all_pass = False
if all_pass:
log("\nALL SCENES PASSED open-loop validation. Proceed to Phase 3.")
else:
log("\nSOME SCENES FAILED. Review Phase 2 results before Phase 3.")
@@ -1,562 +0,0 @@
#!/usr/bin/env python3
# reproduce/phase3_reproduce.py
"""Phase 3: DRL inference with legacy-compatible config + reference comparison.
Uses config_lbm_pinball_legacy_compat.json (regularized inlet, NBB equivalent)
instead of the default config_lbm_pinball.json. After inference, compares
output against SR_analysis reference data.
Scenes: karman_re100, steady_cloak, illusion_1L, vortex_lamb
Usage:
conda run -n pycuda_3_10 python phase3_reproduce.py --device 0
conda run -n pycuda_3_10 python phase3_reproduce.py --device 0 --scene karman
"""
from __future__ import annotations
import argparse
import json
import os
import sys
import time
from collections import deque
from pathlib import Path
from typing import Any, Dict
import numpy as np
import pycuda.driver as cuda; cuda.init()
_REPO = str(Path(__file__).resolve().parents[3])
_SRC = Path(_REPO) / "src"
_DRL = _SRC / "drl_pinball"
for p in [_REPO, str(_SRC), str(_DRL)]:
if p not in sys.path:
sys.path.insert(0, p)
import torch
from torch.nn import Module as TorchModule
from stable_baselines3 import PPO
from CelerisLab import Simulation # noqa: E402
from CelerisLab.common.render import compute_vorticity, render_vorticity_field # noqa: E402
from CelerisLab.lbm.initializers import add_vortex # noqa: E402
from drl_pinball.reproduce.configs.model_inventory import ModelInventory # noqa: E402
from reproduce.core.drl_comparator import compare_reproduce_output # noqa: E402
# ---------------------------------------------------------------------------
# Config
# ---------------------------------------------------------------------------
LEGACY_COMPAT_CFG = "configs/config_lbm_pinball_legacy_compat.json"
L0 = 20.0
U0 = 0.01
NX = 1280
NY = 512
CENTER_Y = float(NY - 1) / 2.0
RADIUS = L0 / 2.0
FIFO_LEN = 150
WARMUP = int(4.0 * NX / U0)
# Standard geometry
DIST_X = 10.0 * L0
PB_FRONT_X = 30.0 * L0
PB_REAR_X = 31.3 * L0
SENSOR_X = 40.0 * L0
# Illusion geometry
ILL_PB_FRONT_X = 19.0 * L0
ILL_PB_REAR_X = 20.3 * L0
ILL_SENSOR_X = 30.0 * L0
ILL_TARGET_X = 20.0 * L0
SR_DATA = _SRC / "SR_analysis" / "data"
_THIS_DIR = Path(__file__).resolve().parent
OUT_BASE = _THIS_DIR / "output" / "phase3"
class Sin(TorchModule):
def __init__(self): super().__init__()
def forward(self, x): return torch.sin(x)
def log(msg: str) -> None:
print(f"[{time.strftime('%H:%M:%S')}] {msg}", flush=True)
class ActionSmoother:
def __init__(self, weight=0.1):
self.weight = weight; self._state = None
def __call__(self, target):
t = np.asarray(target, dtype=np.float32)
if self._state is None:
self._state = t.copy()
else:
self._state = (1.0 - self.weight) * self._state + self.weight * t
return self._state.copy()
def reset(self, value=None):
self._state = np.asarray(value, dtype=np.float32).copy() if value is not None else None
def get_cc(sim, sid):
nx, ny = sim.lbm_cfg.nx, sim.lbm_cfg.ny
cells_arr, _ = sim.bodies.get(sid).get_sensor_list(nx, ny)
return float(len(cells_arr))
def action_to_omega(action_norm, scale=8.0, bias=(0.0, -4.0, 4.0)):
b = np.array(bias, dtype=np.float32)
sv = (np.asarray(action_norm, dtype=np.float32) * scale + b) * U0
return -sv / RADIUS
def load_legacy_norm(ref_dir: str) -> Dict[str, Any]:
with open(os.path.join(ref_dir, "norm.json")) as f:
d = json.load(f)
return {
"force_norm_fact": np.float32(d["force_norm_fact"]),
"sens_deviation": np.array(d["sens_deviation"], dtype=np.float32),
"sens_norm_fact": np.array(d["sens_norm_fact"], dtype=np.float32),
}
def normalize_obs(obs_slice, norm):
forces = obs_slice[6:12] / norm["force_norm_fact"]
sens = (obs_slice[0:6] - norm["sens_deviation"]) / norm["sens_norm_fact"]
return np.clip(np.hstack([forces, sens]), -1.0, 1.0).astype(np.float32)
def read_obs_karman(sim, dist_id, sensor_ids, pinball_ids, cc):
obs = list(sim.read_force(dist_id, normalize=True))
for sid in sensor_ids:
s = sim.read_sensor(sid, normalize=True)
obs.extend([float(s[0]) * cc, float(s[1]) * cc])
for pid in pinball_ids:
obs.extend(sim.read_force(pid, normalize=True))
return np.array(obs, dtype=np.float32)
def read_obs_6obj(sim, sensor_ids, pinball_ids, cc):
obs = []
for sid in sensor_ids:
s = sim.read_sensor(sid, normalize=True)
obs.extend([float(s[0]) * cc, float(s[1]) * cc])
for pid in pinball_ids:
obs.extend(sim.read_force(pid, normalize=True))
return np.array(obs, dtype=np.float32)
def save_vorticity(sim, out_path, cylinders, nx=NX, ny=NY):
macro = sim.get_macroscopic()
vort = compute_vorticity(macro["ux"], macro["uy"])
render_vorticity_field(vort, nx=nx, ny=ny, out_path=str(out_path),
cylinders=cylinders, vmin=-0.03, vmax=0.03)
# ---------------------------------------------------------------------------
# Karman Cloak Re100
# ---------------------------------------------------------------------------
def run_karman(device_id: int, out_dir: Path) -> None:
log("=== Phase 3: Karman Cloak Re100 ===")
out_dir.mkdir(parents=True, exist_ok=True)
SI = 800
num_steps = 200
scale, bias = 8.0, (0.0, -4.0, 4.0)
ref_dir = str(SR_DATA / "karman" / "karman_re100")
norm = load_legacy_norm(ref_dir)
# Phase 1: Disturbance + sensors, record target
sim = Simulation(lbm_config_path=LEGACY_COMPAT_CFG, device_id=device_id)
dist_id = sim.add_body("circle", center=(DIST_X, CENTER_Y, 0.0), radius=1.0 * L0)
sensor_ids = [
sim.add_body("sensor", center=(SENSOR_X, CENTER_Y + 40.0, 0.0), radius=5.0),
sim.add_body("sensor", center=(SENSOR_X, CENTER_Y, 0.0), radius=5.0),
sim.add_body("sensor", center=(SENSOR_X, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim.initialize()
sim.run(WARMUP, zero_obs=True)
cc = get_cc(sim, sensor_ids[0])
log(f" Sensor CC: {cc}")
target_states = np.zeros((FIFO_LEN, 6), dtype=np.float32)
for i in range(FIFO_LEN):
sim.run(SI, zero_obs=True)
obs = read_obs_karman(sim, dist_id, sensor_ids, [], cc)
target_states[i] = obs[2:8]
np.savez_compressed(out_dir / "target.npz", target_states=target_states)
# Phase 2: Add pinball
n0 = sim.bodies.count
sim.add_body("circle", center=(PB_FRONT_X, CENTER_Y, 0.0), radius=RADIUS)
sim.add_body("circle", center=(PB_REAR_X, CENTER_Y + 15.0, 0.0), radius=RADIUS)
sim.add_body("circle", center=(PB_REAR_X, CENTER_Y - 15.0, 0.0), radius=RADIUS)
sim.sync_bodies()
fid, tid, bid = list(range(n0, n0 + 3))
sim.run(WARMUP, zero_obs=True)
# Bias FIFO
bias_norm = np.array([0.0, -1.0, 1.0])
bias_omega = action_to_omega(bias_norm, scale=scale, bias=bias)
ema = ActionSmoother(weight=0.1); ema.reset(np.zeros(3, dtype=np.float32))
for _ in range(FIFO_LEN):
s = ema(bias_omega)
sim.set_body(fid, omega=s[0]); sim.set_body(tid, omega=s[1]); sim.set_body(bid, omega=s[2])
sim.run(SI, zero_obs=True)
sim.snapshot()
# DRL inference
sim.restore()
ema.reset(bias_omega.copy())
model = ModelInventory().load("d1a3o12_re100", device="cpu")
obs_init = read_obs_karman(sim, dist_id, sensor_ids, [fid, tid, bid], cc)
obs_norm = normalize_obs(obs_init[2:14], norm)
sig_s = np.zeros((num_steps, 6), dtype=np.float32)
sig_f = np.zeros((num_steps, 6), dtype=np.float32)
sig_a = np.zeros((num_steps, 3), dtype=np.float32)
for step in range(num_steps):
action, _ = model.predict(obs_norm, deterministic=True)
action = np.asarray(action, dtype=np.float32).flatten()
target_omega = action_to_omega(action, scale=scale, bias=bias)
smoothed = ema(target_omega)
sim.set_body(fid, omega=smoothed[0]); sim.set_body(tid, omega=smoothed[1]); sim.set_body(bid, omega=smoothed[2])
sim.run(SI, zero_obs=True)
obs = read_obs_karman(sim, dist_id, sensor_ids, [fid, tid, bid], cc)
sl = obs[2:14]
sig_s[step] = sl[0:6]
sig_f[step] = sl[6:12]
sig_a[step] = action
obs_norm = normalize_obs(sl, norm)
save_vorticity(sim, out_dir / "vorticity_controlled.png", [
((DIST_X, CENTER_Y), 1.0 * L0),
((PB_FRONT_X, CENTER_Y), RADIUS),
((PB_REAR_X, CENTER_Y + 15.0), RADIUS),
((PB_REAR_X, CENTER_Y - 15.0), RADIUS),
])
sim.close()
np.savez_compressed(out_dir / "signals.npz", sensors=sig_s, forces=sig_f, actions=sig_a)
# Compare against reference
log(" Comparing against SR_analysis reference...")
result = compare_reproduce_output(ref_dir, str(out_dir), label="karman_re100", conv_len=30)
with open(out_dir / "result.json", "w") as f:
json.dump(result, f, indent=2)
log(" Done.")
# ---------------------------------------------------------------------------
# Steady Cloak
# ---------------------------------------------------------------------------
def run_steady(device_id: int, out_dir: Path) -> None:
log("=== Phase 3: Steady Cloak ===")
out_dir.mkdir(parents=True, exist_ok=True)
SI = 800; num_steps = 200
surf_vel = (0.0, -5.1, 5.1)
bias_surf = np.array(surf_vel, dtype=np.float32) * U0
sim = Simulation(lbm_config_path=LEGACY_COMPAT_CFG, device_id=device_id)
sensor_ids = [
sim.add_body("sensor", center=(SENSOR_X, CENTER_Y + 40.0, 0.0), radius=5.0),
sim.add_body("sensor", center=(SENSOR_X, CENTER_Y, 0.0), radius=5.0),
sim.add_body("sensor", center=(SENSOR_X, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim.add_body("circle", center=(PB_FRONT_X, CENTER_Y, 0.0), radius=RADIUS)
sim.add_body("circle", center=(PB_REAR_X, CENTER_Y + 15.0, 0.0), radius=RADIUS)
sim.add_body("circle", center=(PB_REAR_X, CENTER_Y - 15.0, 0.0), radius=RADIUS)
sim.initialize()
sim.run(WARMUP, zero_obs=True)
cc = get_cc(sim, sensor_ids[0])
bias_omega = -bias_surf / RADIUS
ema = ActionSmoother(weight=0.1); ema.reset(np.zeros(3, dtype=np.float32))
for _ in range(FIFO_LEN):
s = ema(bias_omega)
sim.set_body(3, omega=s[0]); sim.set_body(4, omega=s[1]); sim.set_body(5, omega=s[2])
sim.run(SI, zero_obs=True)
sim.snapshot(); sim.restore()
ema.reset(bias_omega.copy())
sig_s = np.zeros((num_steps, 6), dtype=np.float32)
sig_f = np.zeros((num_steps, 6), dtype=np.float32)
for step in range(num_steps):
smoothed = ema(bias_omega)
sim.set_body(3, omega=smoothed[0]); sim.set_body(4, omega=smoothed[1]); sim.set_body(5, omega=smoothed[2])
sim.run(SI, zero_obs=True)
obs = read_obs_6obj(sim, sensor_ids, [3, 4, 5], cc)
sig_s[step] = obs[0:6]
sig_f[step] = obs[6:12]
save_vorticity(sim, out_dir / "vorticity_controlled.png", [
((PB_FRONT_X, CENTER_Y), RADIUS),
((PB_REAR_X, CENTER_Y + 15.0), RADIUS),
((PB_REAR_X, CENTER_Y - 15.0), RADIUS),
])
sim.close()
np.savez_compressed(out_dir / "signals.npz", sensors=sig_s, forces=sig_f,
actions=np.zeros((num_steps, 3), dtype=np.float32))
ref_dir = str(SR_DATA / "steady" / "steady")
result = compare_reproduce_output(ref_dir, str(out_dir), label="steady", conv_len=30)
with open(out_dir / "result.json", "w") as f:
json.dump(result, f, indent=2)
log(" Done.")
# ---------------------------------------------------------------------------
# Illusion 1L
# ---------------------------------------------------------------------------
def run_illusion(device_id: int, out_dir: Path, diam_L: float = 1.0, si: int = 600) -> None:
# diam_L: 0.75 → "illusion_0.75L", 1.0 → "illusion_1L", 1.5 → "illusion_1.5L"
if diam_L == int(diam_L):
label_suffix = str(int(diam_L))
else:
label_suffix = str(diam_L).rstrip('0')
label = f"illusion_{label_suffix}L"
# Map diameter to model name
model_map = {
0.75: "d1a3o14_250525_imit_075L_2U_400S",
1.0: "d1a3o14_250525_imit_1L_2U_600S",
1.5: "d1a3o14_250525_imit_15L_2U",
}
model_name = model_map[diam_L]
log(f"=== Phase 3: Illusion {label} ===")
out_dir.mkdir(parents=True, exist_ok=True)
scale, bias = 8.0, (0.0, -2.0, 2.0)
ref_dir = str(SR_DATA / "illusion" / label)
norm = load_legacy_norm(ref_dir)
# Load target harmonics
with open(os.path.join(ref_dir, "target_harmonics.json")) as f:
target_harmonics = json.load(f)
def gen_target_at(t):
D = len(target_harmonics)
vals = np.zeros(D, dtype=np.float32)
for d, h in enumerate(target_harmonics):
val = float(h["dc"])
for amp, freq, phase in zip(h["amps"], h["freqs"], h["phases"]):
val += amp * np.cos(2.0 * np.pi * freq * t + phase)
vals[d] = val
return vals
# Phase 1: Record target on new CFD
sim_t = Simulation(lbm_config_path=LEGACY_COMPAT_CFG, device_id=device_id)
sim_t.add_body("circle", center=(ILL_TARGET_X, CENTER_Y, 0.0), radius=diam_L * L0)
s_ids_t = [
sim_t.add_body("sensor", center=(ILL_SENSOR_X, CENTER_Y + 40.0, 0.0), radius=5.0),
sim_t.add_body("sensor", center=(ILL_SENSOR_X, CENTER_Y, 0.0), radius=5.0),
sim_t.add_body("sensor", center=(ILL_SENSOR_X, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim_t.initialize(); sim_t.run(WARMUP, zero_obs=True)
cc_t = get_cc(sim_t, s_ids_t[0])
target = np.zeros((FIFO_LEN, 8), dtype=np.float32)
for i in range(FIFO_LEN):
sim_t.run(si, zero_obs=True)
f = list(sim_t.read_force(0, normalize=True))
s = [float(sim_t.read_sensor(sid, normalize=True)[d]) * cc_t for sid in s_ids_t for d in range(2)]
target[i] = np.hstack([f, s])
sim_t.close()
np.savez_compressed(out_dir / "target.npz", target_states=target)
# Phase 2: Pinball + sensors
sim = Simulation(lbm_config_path=LEGACY_COMPAT_CFG, device_id=device_id)
sensor_ids = [
sim.add_body("sensor", center=(ILL_SENSOR_X, CENTER_Y + 40.0, 0.0), radius=5.0),
sim.add_body("sensor", center=(ILL_SENSOR_X, CENTER_Y, 0.0), radius=5.0),
sim.add_body("sensor", center=(ILL_SENSOR_X, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim.add_body("circle", center=(ILL_PB_FRONT_X, CENTER_Y, 0.0), radius=RADIUS)
sim.add_body("circle", center=(ILL_PB_REAR_X, CENTER_Y + 15.0, 0.0), radius=RADIUS)
sim.add_body("circle", center=(ILL_PB_REAR_X, CENTER_Y - 15.0, 0.0), radius=RADIUS)
sim.initialize(); sim.run(WARMUP, zero_obs=True)
cc = get_cc(sim, sensor_ids[0])
log(f" Sensor CC: {cc}")
# Bias FIFO (init bias = [0, -1, 1] * U0)
init_bias_surf = np.array([0.0, -1.0, 1.0]) * U0
init_bias_omega = -init_bias_surf / RADIUS
ema = ActionSmoother(weight=0.1); ema.reset(np.zeros(3, dtype=np.float32))
for _ in range(FIFO_LEN):
s = ema(init_bias_omega)
sim.set_body(3, omega=s[0]); sim.set_body(4, omega=s[1]); sim.set_body(5, omega=s[2])
sim.run(si, zero_obs=True)
sim.snapshot(); sim.restore()
ema.reset(action_to_omega(np.array([0.0, -1.0, 1.0]), scale=scale, bias=bias))
# DRL inference
model = ModelInventory().load(model_name, device="cpu")
obs_init = read_obs_6obj(sim, sensor_ids, [3, 4, 5], cc)
obs_12 = normalize_obs(obs_init, norm)
target_cd = (gen_target_at(0)[0] - norm["sens_deviation"][0]) / norm["sens_norm_fact"][0]
target_cl = (gen_target_at(0)[1] - norm["sens_deviation"][1]) / norm["sens_norm_fact"][1]
obs_norm = np.clip(np.hstack([obs_12, [target_cd, target_cl]]), -1.0, 1.0).astype(np.float32)
num_steps = 200
sig_s = np.zeros((num_steps, 6), dtype=np.float32)
sig_f = np.zeros((num_steps, 6), dtype=np.float32)
sig_a = np.zeros((num_steps, 3), dtype=np.float32)
for step in range(num_steps):
action, _ = model.predict(obs_norm, deterministic=True)
action = np.asarray(action, dtype=np.float32).flatten()
target_omega = action_to_omega(action, scale=scale, bias=bias)
smoothed = ema(target_omega)
sim.set_body(3, omega=smoothed[0]); sim.set_body(4, omega=smoothed[1]); sim.set_body(5, omega=smoothed[2])
sim.run(si, zero_obs=True)
obs = read_obs_6obj(sim, sensor_ids, [3, 4, 5], cc)
sig_s[step] = obs[0:6]
sig_f[step] = obs[6:12]
sig_a[step] = action
obs_12 = normalize_obs(obs, norm)
tgt = gen_target_at(step)
target_cd = (tgt[0] - norm["sens_deviation"][0]) / norm["sens_norm_fact"][0]
target_cl = (tgt[1] - norm["sens_deviation"][1]) / norm["sens_norm_fact"][1]
obs_norm = np.clip(np.hstack([obs_12, [target_cd, target_cl]]), -1.0, 1.0).astype(np.float32)
save_vorticity(sim, out_dir / "vorticity_controlled.png", [
((ILL_PB_FRONT_X, CENTER_Y), RADIUS),
((ILL_PB_REAR_X, CENTER_Y + 15.0), RADIUS),
((ILL_PB_REAR_X, CENTER_Y - 15.0), RADIUS),
])
sim.close()
np.savez_compressed(out_dir / "signals.npz", sensors=sig_s, forces=sig_f, actions=sig_a)
result = compare_reproduce_output(ref_dir, str(out_dir), label=label, conv_len=36)
with open(out_dir / "result.json", "w") as f:
json.dump(result, f, indent=2)
log(" Done.")
# ---------------------------------------------------------------------------
# Vortex Lamb
# ---------------------------------------------------------------------------
def run_vortex(device_id: int, out_dir: Path, vortex_type: str = "lamb") -> None:
log(f"=== Phase 3: Vortex {vortex_type} ===")
out_dir.mkdir(parents=True, exist_ok=True)
SI = 800; num_steps = 150; scale, bias = 4.0, (0.0, -4.0, 4.0)
ref_dir = str(SR_DATA / "vortex" / f"vortex_{vortex_type}")
norm = load_legacy_norm(ref_dir)
strength = 0.5 * U0 if vortex_type == "lamb" else 0.03 * U0
# Phase 1: Sensors only + vortex, record target
sim_t = Simulation(lbm_config_path=LEGACY_COMPAT_CFG, device_id=device_id)
s_ids_t = [
sim_t.add_body("sensor", center=(SENSOR_X, CENTER_Y + 40.0, 0.0), radius=5.0),
sim_t.add_body("sensor", center=(SENSOR_X, CENTER_Y, 0.0), radius=5.0),
sim_t.add_body("sensor", center=(SENSOR_X, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim_t.initialize(); sim_t.run(WARMUP, zero_obs=True)
sim_t.field.download_ddf()
add_vortex(sim_t.field, (10.0 * L0, CENTER_Y), 2.0 * L0, strength, vortex_type)
cc_t = get_cc(sim_t, s_ids_t[0])
target = np.zeros((FIFO_LEN, 6), dtype=np.float32)
for i in range(FIFO_LEN):
sim_t.run(SI, zero_obs=True)
for j, sid in enumerate(s_ids_t):
s = sim_t.read_sensor(sid, normalize=True)
target[i, j * 2] = float(s[0]) * cc_t
target[i, j * 2 + 1] = float(s[1]) * cc_t
sim_t.close()
np.savez_compressed(out_dir / "target.npz", target_states=target)
# Phase 2: Pinball + sensors + vortex
sim = Simulation(lbm_config_path=LEGACY_COMPAT_CFG, device_id=device_id)
sensor_ids = [
sim.add_body("sensor", center=(SENSOR_X, CENTER_Y + 40.0, 0.0), radius=5.0),
sim.add_body("sensor", center=(SENSOR_X, CENTER_Y, 0.0), radius=5.0),
sim.add_body("sensor", center=(SENSOR_X, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim.add_body("circle", center=(PB_FRONT_X, CENTER_Y, 0.0), radius=RADIUS)
sim.add_body("circle", center=(PB_REAR_X, CENTER_Y + 15.0, 0.0), radius=RADIUS)
sim.add_body("circle", center=(PB_REAR_X, CENTER_Y - 15.0, 0.0), radius=RADIUS)
sim.initialize(); sim.run(WARMUP, zero_obs=True)
sim.field.download_ddf()
add_vortex(sim.field, (15.0 * L0, CENTER_Y), 2.0 * L0, strength, vortex_type)
cc = get_cc(sim, sensor_ids[0])
log(f" Sensor CC: {cc}")
# Bias FIFO
bias_omega = action_to_omega(np.array([0.0, -1.0, 1.0]), scale=scale, bias=bias)
ema = ActionSmoother(weight=0.1); ema.reset(np.zeros(3, dtype=np.float32))
for _ in range(FIFO_LEN):
s = ema(bias_omega)
sim.set_body(3, omega=s[0]); sim.set_body(4, omega=s[1]); sim.set_body(5, omega=s[2])
sim.run(SI, zero_obs=True)
sim.snapshot(); sim.restore()
ema.reset(bias_omega.copy())
# DRL inference
model = ModelInventory().load(f"vortex_{vortex_type}", device="cpu")
obs_init = read_obs_6obj(sim, sensor_ids, [3, 4, 5], cc)
obs_norm = normalize_obs(obs_init, norm)
sig_s = np.zeros((num_steps, 6), dtype=np.float32)
sig_f = np.zeros((num_steps, 6), dtype=np.float32)
sig_a = np.zeros((num_steps, 3), dtype=np.float32)
for step in range(num_steps):
action, _ = model.predict(obs_norm, deterministic=True)
action = np.asarray(action, dtype=np.float32).flatten()
target_omega = action_to_omega(action, scale=scale, bias=bias)
smoothed = ema(target_omega)
sim.set_body(3, omega=smoothed[0]); sim.set_body(4, omega=smoothed[1]); sim.set_body(5, omega=smoothed[2])
sim.run(SI, zero_obs=True)
obs = read_obs_6obj(sim, sensor_ids, [3, 4, 5], cc)
sig_s[step] = obs[0:6]
sig_f[step] = obs[6:12]
sig_a[step] = action
obs_norm = normalize_obs(obs, norm)
save_vorticity(sim, out_dir / "vorticity_controlled.png", [
((PB_FRONT_X, CENTER_Y), RADIUS),
((PB_REAR_X, CENTER_Y + 15.0), RADIUS),
((PB_REAR_X, CENTER_Y - 15.0), RADIUS),
])
sim.close()
np.savez_compressed(out_dir / "signals.npz", sensors=sig_s, forces=sig_f, actions=sig_a)
result = compare_reproduce_output(ref_dir, str(out_dir), label=f"vortex_{vortex_type}", conv_len=30)
with open(out_dir / "result.json", "w") as f:
json.dump(result, f, indent=2)
log(" Done.")
# ---------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------
if __name__ == "__main__":
ap = argparse.ArgumentParser(description="Phase 3: DRL reproduce with legacy-compat config")
ap.add_argument("--device", type=int, default=0, help="GPU device ID")
ap.add_argument("--scene", type=str, default="all",
help="Scene: karman, steady, illusion_1L, vortex_lamb, vortex_taylor, all")
args = ap.parse_args()
all_scenes = {
"karman": lambda: run_karman(args.device, OUT_BASE / "karman"),
"steady": lambda: run_steady(args.device, OUT_BASE / "steady"),
"illusion_1L": lambda: run_illusion(args.device, OUT_BASE / "illusion_1L", 1.0, 600),
"vortex_lamb": lambda: run_vortex(args.device, OUT_BASE / "vortex_lamb", "lamb"),
"vortex_taylor": lambda: run_vortex(args.device, OUT_BASE / "vortex_taylor", "taylor"),
}
if args.scene == "all":
for name, func in all_scenes.items():
func()
else:
for s in args.scene.split(","):
s = s.strip()
if s in all_scenes:
all_scenes[s]()
log("\nPhase 3 complete.")
-316
View File
@@ -1,316 +0,0 @@
#!/usr/bin/env python3
"""Comprehensive run: Steady Cloak + Karman Cloak (new norm / legacy norm) + flow field output.
DEPRECATED (2026-07): Superseded by phase2_open_loop.py and phase3_reproduce.py
which use the legacy-compatible config (regularized inlet). This script uses
the old config_lbm_pinball.json with zou_he_local inlet, producing suboptimal results.
Keep for reference; use phase2/phase3 for new reproduce work.
Runs three cases:
Case A: Steady Cloak (open-loop constant rotation, no DRL)
Case B: Karman Cloak (new-CFD norm + DRL inference)
Case C: Karman Cloak (legacy norm + DRL inference)
Each case saves: macroscopic.npz (rho/ux/uy), vorticity.png, sensors_forces.npz
"""
from __future__ import annotations
import argparse
import json
import os
import sys
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
_SRC = os.path.join(_REPO, "src")
for p in [_REPO, _SRC]:
if p not in sys.path:
sys.path.insert(0, p)
import pycuda.driver as cuda
cuda.init()
from CelerisLab import Simulation
from CelerisLab.common.render import compute_vorticity, render_vorticity_field
from drl_pinball.reproduce.configs.model_inventory import ModelInventory
from drl_pinball.reproduce.core.action_wrapper import (
ActionSmoother, norm_action_to_omega, surface_vel_to_omega, U0 as _U0, RADIUS,
)
from drl_pinball.reproduce.core.obs_normalizer import compute_norm, normalize_observation
# ---------------------------------------------------------------------------
# Constants
# ---------------------------------------------------------------------------
L0 = 20.0
U0 = float(_U0)
NX = 1280
NY = 512
CENTER_Y = float(NY - 1) / 2.0
FIFO_LEN = 150
SI = 800
WARMUP = int(4.0 * NX / U0)
CFG_PATH = "configs/config_lbm_pinball_legacy_compat.json"
REF_DIR = os.path.join(_REPO, "src", "SR_analysis", "data", "karman", "karman_re100")
OUT_BASE = os.path.join(os.path.dirname(__file__), "output")
def get_cc(sim, sid):
nx, ny = sim.lbm_cfg.nx, sim.lbm_cfg.ny
cells_arr, _ = sim.bodies.get(sid).get_sensor_list(nx, ny)
return float(len(cells_arr))
def read_obs_legacy(sim, sensor_ids, dist_id, pinball_ids, cc):
obs = []
if dist_id is not None:
obs.extend(sim.read_force(dist_id, normalize=True))
for sid in sensor_ids:
s = sim.read_sensor(sid, normalize=True)
obs.extend(s * cc)
for pid in pinball_ids:
obs.extend(sim.read_force(pid, normalize=True))
return np.array(obs, dtype=np.float32)
def save_field(sim, out_dir, name):
"""Save macroscopic field and render vorticity."""
macro = sim.get_macroscopic()
np.savez_compressed(os.path.join(out_dir, f"macro_{name}.npz"),
rho=macro["rho"], ux=macro["ux"], uy=macro["uy"])
vort = compute_vorticity(macro["ux"], macro["uy"])
render_vorticity_field(
vort, nx=NX, ny=NY,
out_path=os.path.join(out_dir, f"vorticity_{name}.png"),
cylinders=[
((10.0 * L0, CENTER_Y), 1.0 * L0),
((30.0 * L0, CENTER_Y), RADIUS),
((31.3 * L0, CENTER_Y + 15.0), RADIUS),
((31.3 * L0, CENTER_Y - 15.0), RADIUS),
],
)
print(f" Saved {name}: macro + vorticity.png")
# =========================================================================
# Case A: Steady Cloak (open-loop constant rotation)
# =========================================================================
def run_steady_cloak(device_id, out_dir):
"""Steady cloak: pinball only, constant bias=[0,-5.1,5.1]*U0, no DRL."""
print("\n" + "=" * 70)
print("Case A: Steady Cloak (open-loop constant rotation)")
print("=" * 70)
os.makedirs(out_dir, exist_ok=True)
sim = Simulation(lbm_config_path=CFG_PATH, device_id=device_id)
sensor_ids = [
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y + 40.0, 0.0), radius=5.0),
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y, 0.0), radius=5.0),
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim.add_body("circle", center=(30.0 * L0, CENTER_Y, 0.0), radius=RADIUS)
sim.add_body("circle", center=(31.3 * L0, CENTER_Y + 15.0, 0.0), radius=RADIUS)
sim.add_body("circle", center=(31.3 * L0, CENTER_Y - 15.0, 0.0), radius=RADIUS)
sim.initialize()
sim.run(WARMUP, zero_obs=True)
cc = get_cc(sim, sensor_ids[0])
save_field(sim, out_dir, "steady_uncontrolled")
# Bias: surface_vel = [0, -5.1, 5.1] * U0 → omega
bias_surf = np.array([0.0, -5.1, 5.1], dtype=np.float32) * U0
bias_omega = surface_vel_to_omega(bias_surf)
print(f" Steady bias omega: {bias_omega}")
sim.set_body(3, omega=bias_omega[0]) # front
sim.set_body(4, omega=bias_omega[1]) # top
sim.set_body(5, omega=bias_omega[2]) # bottom
sim.run(WARMUP, zero_obs=True)
sensors_f = []
for _ in range(200):
sim.run(SI, zero_obs=True)
obs = read_obs_legacy(sim, sensor_ids, None, [3, 4, 5], cc)
sensors_f.append(obs[0:12])
sensors_f = np.array(sensors_f, dtype=np.float32)
np.savez_compressed(os.path.join(out_dir, "steady_signals.npz"),
sensors=sensors_f[:, 0:6], forces=sensors_f[:, 6:12])
save_field(sim, out_dir, "steady_controlled")
print(f" Forces: front_fy={sensors_f[:,7].mean():+.6f} "
f"top_fy={sensors_f[:,9].mean():+.6f} bottom_fy={sensors_f[:,11].mean():+.6f}")
sim.close()
# =========================================================================
# Karman shared build
# =========================================================================
def build_karman_env(device_id, out_dir):
"""Build Karman env, record target, add pinball, return (sim, ids, norm)."""
sim = Simulation(lbm_config_path=CFG_PATH, device_id=device_id)
dist_id = sim.add_body("circle", center=(10.0 * L0, CENTER_Y, 0.0), radius=1.0 * L0)
sensor_ids = [
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y + 40.0, 0.0), radius=5.0),
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y, 0.0), radius=5.0),
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim.initialize()
sim.run(WARMUP, zero_obs=True)
cc = get_cc(sim, sensor_ids[0])
target = np.empty((0, 6), dtype=np.float32)
for _ in range(FIFO_LEN):
sim.run(SI, zero_obs=True)
obs = read_obs_legacy(sim, sensor_ids, dist_id, [], cc)
target = np.vstack((target, obs[2:8]))
np.savez_compressed(os.path.join(out_dir, "target.npz"), target_states=target)
n0 = sim.bodies.count
sim.add_body("circle", center=(30.0 * L0, CENTER_Y, 0.0), radius=RADIUS)
sim.add_body("circle", center=(31.3 * L0, CENTER_Y + 15.0, 0.0), radius=RADIUS)
sim.add_body("circle", center=(31.3 * L0, CENTER_Y - 15.0, 0.0), radius=RADIUS)
sim.sync_bodies()
fid, tid, bid = list(range(n0, n0 + 3))
sim.run(WARMUP, zero_obs=True)
return sim, dist_id, sensor_ids, (fid, tid, bid), cc, target
def run_karman_drl(sim, dist_id, sensor_ids, pinball_ids, cc, target,
norm, model, out_dir, case_label, num_steps=200):
"""Run DRL inference with given norm."""
fid, tid, bid = pinball_ids
# Bias FIFO: surface_vel = [0, -4, 4] * U0 (legacy: bias_arr[-3:] = front, top, bottom)
bias_surf = np.array([0.0, -4.0, 4.0], dtype=np.float32) * U0
bias_omega = surface_vel_to_omega(bias_surf)
ema = ActionSmoother(weight=0.1)
ema.reset(np.zeros(3, dtype=np.float32)) # legacy starts from zero
for _ in range(FIFO_LEN):
s = ema(bias_omega)
sim.set_body(fid, omega=s[0])
sim.set_body(tid, omega=s[1])
sim.set_body(bid, omega=s[2])
sim.run(SI, zero_obs=True)
sim.snapshot()
# DRL inference
sim.restore()
ema.reset(bias_omega.copy()) # EMA starts from converged bias state
obs_init = read_obs_legacy(sim, sensor_ids, dist_id, [fid, tid, bid], cc)
obs_norm = normalize_observation(obs_init[2:14], norm)
sig_s = np.zeros((num_steps, 6), dtype=np.float32)
sig_f = np.zeros((num_steps, 6), dtype=np.float32)
sig_a = np.zeros((num_steps, 3), dtype=np.float32)
for step in range(num_steps):
action, _ = model.predict(obs_norm, deterministic=True)
action = np.asarray(action, dtype=np.float32).flatten()
target_omega = norm_action_to_omega(action, scale=8.0, bias=(0.0, -4.0, 4.0))
smoothed = ema(target_omega)
sim.set_body(fid, omega=smoothed[0])
sim.set_body(tid, omega=smoothed[1])
sim.set_body(bid, omega=smoothed[2])
sim.run(SI, zero_obs=True)
obs = read_obs_legacy(sim, sensor_ids, dist_id, [fid, tid, bid], cc)
sl = obs[2:14]
sig_s[step] = sl[0:6]
sig_f[step] = sl[6:12]
sig_a[step] = action
obs_norm = normalize_observation(sl, norm)
save_field(sim, out_dir, f"karman_{case_label}")
np.savez_compressed(os.path.join(out_dir, f"signals_{case_label}.npz"),
sensors=sig_s, forces=sig_f, actions=sig_a)
print(f" [{case_label}] Actions: aF={sig_a[:,0].mean():+.4f} "
f"aB={sig_a[:,1].mean():+.4f} aT={sig_a[:,2].mean():+.4f}")
print(f" [{case_label}] Forces: front_fy={sig_f[:,1].mean():+.6f} "
f"top_fy={sig_f[:,3].mean():+.6f} bottom_fy={sig_f[:,5].mean():+.6f}")
# =========================================================================
# Case B: Karman + new norm
# =========================================================================
def collect_new_norm(sim, sensor_ids, dist_id, pinball_ids, cc):
"""Collect norm values on new CFD from zero-action FIFO."""
fid, tid, bid = pinball_ids
fifo = []
for _ in range(FIFO_LEN):
sim.run(SI, zero_obs=True)
obs = read_obs_legacy(sim, sensor_ids, dist_id, [fid, tid, bid], cc)
fifo.append(obs[2:14])
f = np.array(fifo, dtype=np.float32)
norm = compute_norm(f, force_slice=(6, 12), sens_slice=(0, 6))
print(f" New norm: fn={norm['force_norm_fact']:.6f}")
return norm
def run_karman_new_norm(device_id, out_dir):
print("\n" + "=" * 70)
print("Case B: Karman Cloak (new-CFD norm + DRL)")
print("=" * 70)
os.makedirs(out_dir, exist_ok=True)
model = ModelInventory().load("d1a3o12_re100", device="cpu")
sim, dist_id, sensor_ids, pids, cc, target = build_karman_env(device_id, out_dir)
norm = collect_new_norm(sim, sensor_ids, dist_id, pids, cc)
np.savez(os.path.join(out_dir, "norm_new.npz"), **norm)
run_karman_drl(sim, dist_id, sensor_ids, pids, cc, target,
norm, model, out_dir, "new_norm")
sim.close()
# =========================================================================
# Case C: Karman + legacy norm
# =========================================================================
def run_karman_legacy_norm(device_id, out_dir):
print("\n" + "=" * 70)
print("Case C: Karman Cloak (legacy norm + DRL)")
print("=" * 70)
os.makedirs(out_dir, exist_ok=True)
with open(os.path.join(REF_DIR, "norm.json")) as f:
d = json.load(f)
norm = {
"force_norm_fact": np.float32(d["force_norm_fact"]),
"sens_deviation": np.array(d["sens_deviation"], dtype=np.float32),
"sens_norm_fact": np.array(d["sens_norm_fact"], dtype=np.float32),
}
model = ModelInventory().load("d1a3o12_re100", device="cpu")
sim, dist_id, sensor_ids, pids, cc, target = build_karman_env(device_id, out_dir)
run_karman_drl(sim, dist_id, sensor_ids, pids, cc, target,
norm, model, out_dir, "legacy_norm")
sim.close()
# =========================================================================
# Main
# =========================================================================
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument("--device", type=int, default=0)
parser.add_argument("--cases", type=str, default="A,B,C",
help="Comma-separated: A=steady, B=new-norm, C=legacy-norm")
args = parser.parse_args()
cases = [c.strip().upper() for c in args.cases.split(",")]
if "A" in cases:
run_steady_cloak(args.device, os.path.join(OUT_BASE, "steady_cloak"))
karman_dir = os.path.join(OUT_BASE, "karman_cloak")
if "B" in cases:
run_karman_new_norm(args.device, karman_dir)
if "C" in cases:
run_karman_legacy_norm(args.device, karman_dir)
print("\nDone.")
@@ -1,58 +0,0 @@
#!/bin/bash
# reproduce/run_all_reproduce_tests.sh
#
# Sequential launcher for Track B (Reproduce) scripts.
# Phase 2: Open-loop target validation
# Phase 3: DRL inference with legacy-compatible config
#
# Usage:
# bash run_all_reproduce_tests.sh [DEVICE_ID] [PHASE]
# DEVICE_ID defaults to 0
# PHASE defaults to "2,3" (run both phases)
set -euo pipefail
DEVICE_ID="${1:-0}"
PHASE="${2:-2,3}"
SCRIPT_DIR="$(cd "$(dirname "${BASH_SOURCE[0]}")" && pwd)"
cd "$SCRIPT_DIR/../../.." # repo root
log() { echo "[$(date '+%H:%M:%S')] $*"; }
CONDA_ENV="pycuda_3_10"
DELAY=60
log "=== Reproduce Tests ==="
log "Device: $DEVICE_ID, Phase: $PHASE"
# --- Phase 2: Open-loop validation ---
if [[ "$PHASE" == *"2"* ]]; then
log ""
log "--- Phase 2: Open-loop target validation ---"
if conda run -n "$CONDA_ENV" python src/drl_pinball/reproduce/phase2_open_loop.py \
--device "$DEVICE_ID" --scene all; then
log "[PASS] Phase 2"
else
log "[FAIL] Phase 2"
fi
log "Waiting ${DELAY}s..."
sleep "$DELAY"
fi
# --- Phase 3: DRL inference ---
if [[ "$PHASE" == *"3"* ]]; then
log ""
log "--- Phase 3: DRL inference with legacy-compat config ---"
if conda run -n "$CONDA_ENV" python src/drl_pinball/reproduce/phase3_reproduce.py \
--device "$DEVICE_ID" --scene all; then
log "[PASS] Phase 3"
else
log "[FAIL] Phase 3"
fi
fi
log ""
log "=== Reproduce tests complete ==="
@@ -1,404 +0,0 @@
#!/usr/bin/env python3
"""Phase 2 (Illusion) + Phase 3 (Vortex) reproduction on new CelerisLab.
Uses same proven approach as Karman: legacy norm + old-equiv sensor conversion.
"""
from __future__ import annotations
import argparse
import json
import os
import sys
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
_SRC = os.path.join(_REPO, "src")
for p in [_REPO, _SRC]:
if p not in sys.path:
sys.path.insert(0, p)
import pycuda.driver as cuda
cuda.init()
from CelerisLab import Simulation
from CelerisLab.common.render import compute_vorticity, render_vorticity_field
from CelerisLab.lbm.initializers import add_vortex
from drl_pinball.reproduce.configs.model_inventory import ModelInventory
from drl_pinball.reproduce.core.action_wrapper import (
ActionSmoother, norm_action_to_omega, surface_vel_to_omega, U0 as _U0, RADIUS,
)
from drl_pinball.reproduce.core.obs_normalizer import normalize_observation
# ---------------------------------------------------------------------------
# Constants
# ---------------------------------------------------------------------------
L0 = 20.0
U0 = float(_U0)
NX = 1280
NY = 512
CENTER_Y = float(NY - 1) / 2.0
FIFO_LEN = 150
WARMUP = int(4.0 * NX / U0)
CFG_PATH = "configs/config_lbm_pinball_legacy_compat.json"
OUT_BASE = os.path.join(os.path.dirname(__file__), "output")
def get_cc(sim, sid):
nx, ny = sim.lbm_cfg.nx, sim.lbm_cfg.ny
cells_arr, _ = sim.bodies.get(sid).get_sensor_list(nx, ny)
return float(len(cells_arr))
def read_obs_legacy(sim, sensor_ids, pinball_ids, cc):
"""Return [s0_ux,uy, s1_ux,uy, s2_ux,uy, front_fx,fy, top_fx,fy, bottom_fx,fy].
No dist_cylinder — this is for 6-object envs (illusion, vortex, steady).
"""
obs = []
for sid in sensor_ids:
s = sim.read_sensor(sid, normalize=True)
obs.extend(s * cc)
for pid in pinball_ids:
obs.extend(sim.read_force(pid, normalize=True))
return np.array(obs, dtype=np.float32)
def load_norm(path):
with open(path) as f:
d = json.load(f)
return {
"force_norm_fact": np.float32(d["force_norm_fact"]),
"sens_deviation": np.array(d["sens_deviation"], dtype=np.float32),
"sens_norm_fact": np.array(d["sens_norm_fact"], dtype=np.float32),
}
def save_vorticity(sim, out_path, cylinders):
macro = sim.get_macroscopic()
vort = compute_vorticity(macro["ux"], macro["uy"])
render_vorticity_field(vort, nx=NX, ny=NY, out_path=out_path, cylinders=cylinders)
# =========================================================================
# Phase 2: Illusion
# =========================================================================
def run_illusion(device_id, target_label, model_name, sample_interval, target_diameter,
ref_dir, out_dir):
"""Run Illusion inference using legacy norm + legacy target harmonics."""
print(f"\n{'='*70}")
print(f"Illusion {target_label} (model={model_name}, SI={sample_interval})")
print(f"{'='*70}")
os.makedirs(out_dir, exist_ok=True)
norm = load_norm(os.path.join(ref_dir, "norm.json"))
with open(os.path.join(ref_dir, "target_harmonics.json")) as f:
target_harmonics = json.load(f)
legacy_target = np.load(os.path.join(ref_dir, "target.npz"))
legacy_target_states = legacy_target["target_states"]
def gen_target_states_at(t, harmonics):
D = len(harmonics)
result = np.zeros(D, dtype=np.float32)
for d, h in enumerate(harmonics):
val = np.float32(h["dc"])
for amp, freq, phase in zip(h["amps"], h["freqs"], h["phases"]):
val += amp * np.cos(2 * np.pi * freq * t + phase)
result[d] = val
return result
# Pinball env (no disturbance cylinder)
sim = Simulation(lbm_config_path=CFG_PATH, device_id=device_id)
sensor_ids = [
sim.add_body("sensor", center=(30.0 * L0, CENTER_Y + 40.0, 0.0), radius=5.0),
sim.add_body("sensor", center=(30.0 * L0, CENTER_Y, 0.0), radius=5.0),
sim.add_body("sensor", center=(30.0 * L0, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim.add_body("circle", center=(19.0 * L0, CENTER_Y, 0.0), radius=RADIUS)
sim.add_body("circle", center=(20.3 * L0, CENTER_Y + 15.0, 0.0), radius=RADIUS)
sim.add_body("circle", center=(20.3 * L0, CENTER_Y - 15.0, 0.0), radius=RADIUS)
sim.initialize()
sim.run(WARMUP, zero_obs=True)
cc = get_cc(sim, sensor_ids[0])
fid, tid, bid = 3, 4, 5
# Bias FIFO: init surface_vel = [0, -1, 1] * U0 (matching legacy_env_imit.py)
init_bias_surf = np.array([0.0, -1.0, 1.0], dtype=np.float32) * U0
init_bias_omega = surface_vel_to_omega(init_bias_surf)
ema_bias = ActionSmoother(weight=0.1)
ema_bias.reset(np.zeros(3, dtype=np.float32))
for _ in range(FIFO_LEN):
s = ema_bias(init_bias_omega)
sim.set_body(fid, omega=s[0]); sim.set_body(tid, omega=s[1]); sim.set_body(bid, omega=s[2])
sim.run(sample_interval, zero_obs=True)
sim.snapshot()
# DRL inference — EMA starts from converged init-bias state
sim.restore()
ema = ActionSmoother(weight=0.1)
ema.reset(init_bias_omega.copy())
model = ModelInventory().load(model_name, device="cpu")
print(f" Model loaded on CPU")
obs_init = read_obs_legacy(sim, sensor_ids, [fid, tid, bid], cc)
obs_norm = normalize_observation(obs_init, norm)
t0 = gen_target_states_at(0, target_harmonics)
target_cd = np.float32(t0[0] / norm["force_norm_fact"])
target_cl = np.float32(t0[1] / norm["force_norm_fact"])
obs_norm = np.clip(np.hstack([obs_norm, [target_cd, target_cl]]), -1, 1).astype(np.float32)
num_steps = 200
sig_s = np.zeros((num_steps, 6), dtype=np.float32)
sig_f = np.zeros((num_steps, 6), dtype=np.float32)
sig_a = np.zeros((num_steps, 3), dtype=np.float32)
for step in range(num_steps):
action, _ = model.predict(obs_norm, deterministic=True)
action = np.asarray(action, dtype=np.float32).flatten()
target_omega = norm_action_to_omega(action, scale=8.0, bias=(0.0, -2.0, 2.0))
smoothed = ema(target_omega)
sim.set_body(fid, omega=smoothed[0]); sim.set_body(tid, omega=smoothed[1]); sim.set_body(bid, omega=smoothed[2])
sim.run(sample_interval, zero_obs=True)
obs = read_obs_legacy(sim, sensor_ids, [fid, tid, bid], cc)
sig_s[step] = obs[0:6]
sig_f[step] = obs[6:12]
sig_a[step] = action
obs_norm_base = normalize_observation(obs, norm)
t_h = gen_target_states_at(step, target_harmonics)
target_cd = np.float32(t_h[0] / norm["force_norm_fact"])
target_cl = np.float32(t_h[1] / norm["force_norm_fact"])
obs_norm = np.clip(np.hstack([obs_norm_base, [target_cd, target_cl]]), -1, 1).astype(np.float32)
save_vorticity(sim, os.path.join(out_dir, f"vorticity.png"),
cylinders=[((19.0*L0, CENTER_Y), RADIUS),
((20.3*L0, CENTER_Y+15.0), RADIUS),
((20.3*L0, CENTER_Y-15.0), RADIUS)])
np.savez_compressed(os.path.join(out_dir, "signals.npz"),
sensors=sig_s, forces=sig_f, actions=sig_a)
sim.close()
# Compare with reference
ref = np.load(os.path.join(ref_dir, "controlled.npz"))
print(f"\n Actions vs ref:")
for i, name in enumerate(["aF","aB","aT"]):
c = np.corrcoef(ref["actions"][:num_steps,i], sig_a[:,i])[0,1]
print(f" {name}: ref_mean={ref['actions'][:num_steps,i].mean():+.4f} "
f"our_mean={sig_a[:,i].mean():+.4f} corr={c:+.4f}")
# DTW similarity
n_c = 36
def dtw_sim(t, s):
n = len(t)
D = np.full((n+1, n+1), np.inf)
D[0,0] = 0
for i in range(1, n+1):
for j in range(1, n+1):
D[i,j] = abs(t[i-1]-s[j-1]) + min(D[i-1,j], D[i,j-1], D[i-1,j-1])
return 1 - D[n,n] / n
t_seq = legacy_target_states[n_c:2*n_c, 1+2]
s_seq = sig_s[-n_c:, 1]
if np.std(t_seq) > 1e-10 and np.std(s_seq) > 1e-10:
corr = np.correlate(t_seq - t_seq.mean(), s_seq - s_seq.mean(), mode="full")
lag = np.argmax(corr) - (len(t_seq) - 1)
else:
lag = 0
sim_val = 0.0
for i in range(6):
t_rolled = np.roll(legacy_target_states[:, i+2], -lag)[n_c:2*n_c]
sim_val += dtw_sim(t_rolled, sig_s[-n_c:, i]) / 6.0
print(f" DTW similarity (legacy target vs our sensors): {sim_val:.4f}")
print(f" Reference DTW similarity: {0.9754}")
# =========================================================================
# Phase 3: Vortex (lamb / taylor)
# =========================================================================
def run_vortex(device_id, vortex_type, model_name, vortex_strength, ref_dir, out_dir):
"""Run Vortex cloak inference using legacy norm."""
print(f"\n{'='*70}")
print(f"Vortex {vortex_type} (model={model_name}, strength={vortex_strength})")
print(f"{'='*70}")
os.makedirs(out_dir, exist_ok=True)
norm = load_norm(os.path.join(ref_dir, "norm.json"))
legacy_target = np.load(os.path.join(ref_dir, "target.npz"))["target_states"]
MAX_STEPS = 150
CONV_LEN = 30
SI = 800
# Stage 1: sensor-only env, record target with vortex
sim = Simulation(lbm_config_path=CFG_PATH, device_id=device_id)
sensor_ids = [
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y + 40.0, 0.0), radius=5.0),
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y, 0.0), radius=5.0),
sim.add_body("sensor", center=(40.0 * L0, CENTER_Y - 40.0, 0.0), radius=5.0),
]
sim.initialize()
sim.run(WARMUP, zero_obs=True)
cc = get_cc(sim, sensor_ids[0])
sim.snapshot()
add_vortex(sim.field, center=(10.0 * L0, CENTER_Y), radius=2.0 * L0,
strength=vortex_strength, vortex_type=vortex_type)
target_states = np.empty((0, 6), dtype=np.float32)
for _ in range(FIFO_LEN):
sim.run(SI, zero_obs=True)
obs = read_obs_legacy(sim, sensor_ids, [], cc)
target_states = np.vstack((target_states, obs))
np.savez_compressed(os.path.join(out_dir, "target.npz"), target_states=target_states)
# Stage 2: add pinball + vortex
sim.restore()
n0 = sim.bodies.count
sim.add_body("circle", center=(30.0 * L0, CENTER_Y, 0.0), radius=RADIUS)
sim.add_body("circle", center=(31.3 * L0, CENTER_Y + 15.0, 0.0), radius=RADIUS)
sim.add_body("circle", center=(31.3 * L0, CENTER_Y - 15.0, 0.0), radius=RADIUS)
sim.sync_bodies()
fid, tid, bid = list(range(n0, n0 + 3))
# Bias warmup: surface_vel = [0, -4, 4] * U0
bias_surf = np.array([0.0, -4.0, 4.0], dtype=np.float32) * U0
bias_omega = surface_vel_to_omega(bias_surf)
ema_init = ActionSmoother(weight=0.1)
ema_init.reset(np.zeros(3, dtype=np.float32))
for _ in range(FIFO_LEN):
s = ema_init(bias_omega)
sim.set_body(fid, omega=s[0]); sim.set_body(tid, omega=s[1]); sim.set_body(bid, omega=s[2])
sim.run(SI, zero_obs=True)
# Add vortex at pinball-phase position
add_vortex(sim.field, center=(15.0 * L0, CENTER_Y), radius=2.0 * L0,
strength=vortex_strength, vortex_type=vortex_type)
sim.snapshot()
print(f" Post-vortex DDF saved")
# Norm FIFO (zero action, with pinball+vortex)
fifo = []
for _ in range(FIFO_LEN):
sim.run(SI, zero_obs=True)
obs = read_obs_legacy(sim, sensor_ids, [fid, tid, bid], cc)
fifo.append(obs)
fifo_arr = np.array(fifo, dtype=np.float32)
new_fn = 6.0 * np.max(np.abs(fifo_arr[:, 6:12]))
print(f" New CFD force_norm_fact: {new_fn:.6f} (legacy: {norm['force_norm_fact']:.6f})")
# Bias FIFO after vortex
sim.restore()
ema_bias = ActionSmoother(weight=0.1)
ema_bias.reset(np.zeros(3, dtype=np.float32))
for _ in range(FIFO_LEN):
s = ema_bias(bias_omega)
sim.set_body(fid, omega=s[0]); sim.set_body(tid, omega=s[1]); sim.set_body(bid, omega=s[2])
sim.run(SI, zero_obs=True)
sim.snapshot()
# DRL inference
sim.restore()
ema = ActionSmoother(weight=0.1)
ema.reset(bias_omega.copy())
model = ModelInventory().load(model_name, device="cpu")
print(f" Model loaded on CPU")
obs_init = read_obs_legacy(sim, sensor_ids, [fid, tid, bid], cc)
obs_norm = normalize_observation(obs_init, norm)
sig_s = np.zeros((MAX_STEPS, 6), dtype=np.float32)
sig_f = np.zeros((MAX_STEPS, 6), dtype=np.float32)
sig_a = np.zeros((MAX_STEPS, 3), dtype=np.float32)
for step in range(MAX_STEPS):
action, _ = model.predict(obs_norm, deterministic=True)
action = np.asarray(action, dtype=np.float32).flatten()
target_omega = norm_action_to_omega(action, scale=4.0, bias=(0.0, -4.0, 4.0))
smoothed = ema(target_omega)
sim.set_body(fid, omega=smoothed[0]); sim.set_body(tid, omega=smoothed[1]); sim.set_body(bid, omega=smoothed[2])
sim.run(SI, zero_obs=True)
obs = read_obs_legacy(sim, sensor_ids, [fid, tid, bid], cc)
sig_s[step] = obs[0:6]
sig_f[step] = obs[6:12]
sig_a[step] = action
obs_norm = normalize_observation(obs, norm)
save_vorticity(sim, os.path.join(out_dir, f"vorticity.png"),
cylinders=[((30.0*L0, CENTER_Y), RADIUS),
((31.3*L0, CENTER_Y+15.0), RADIUS),
((31.3*L0, CENTER_Y-15.0), RADIUS)])
np.savez_compressed(os.path.join(out_dir, "signals.npz"),
sensors=sig_s, forces=sig_f, actions=sig_a)
sim.close()
# Compare with reference
ref = np.load(os.path.join(ref_dir, "controlled.npz"))
print(f"\n Actions vs ref:")
for i, name in enumerate(["aF","aB","aT"]):
c = np.corrcoef(ref["actions"][:MAX_STEPS,i], sig_a[:,i])[0,1]
print(f" {name}: ref_mean={ref['actions'][:MAX_STEPS,i].mean():+.4f} "
f"our_mean={sig_a[:,i].mean():+.4f} corr={c:+.4f}")
# DTW similarity
def dtw_sim(t, s):
n = len(t)
D = np.full((n+1, n+1), np.inf); D[0,0] = 0
for i in range(1, n+1):
for j in range(1, n+1):
D[i,j] = abs(t[i-1]-s[j-1]) + min(D[i-1,j], D[i,j-1], D[i-1,j-1])
return 1 - D[n,n] / n
sim_val = 0.0
for i in range(6):
t_rolled = np.roll(target_states[-CONV_LEN:, i], -MAX_STEPS)
sim_val += dtw_sim(t_rolled, sig_s[-CONV_LEN:, i]) / 6.0
print(f" DTW similarity (our target vs our sensors): {sim_val:.4f}")
# =========================================================================
# Main
# =========================================================================
if __name__ == "__main__":
parser = argparse.ArgumentParser()
parser.add_argument("--device", type=int, default=0)
parser.add_argument("--case", type=str, required=True,
choices=["illusion_075L", "illusion_1L", "illusion_15L",
"vortex_lamb", "vortex_taylor"])
args = parser.parse_args()
_SRC2 = os.path.join(os.path.dirname(__file__), "..", "..", "..", "src")
if "illusion" in args.case:
scenes = {
"illusion_075L": ("d1a3o14_250525_imit_075L_2U_400S", 400, 0.75 * L0,
os.path.join(_SRC2, "SR_analysis", "data", "illusion", "illusion_0.75L")),
"illusion_1L": ("d1a3o14_250525_imit_1L_2U_600S", 600, 1.0 * L0,
os.path.join(_SRC2, "SR_analysis", "data", "illusion", "illusion_1L")),
"illusion_15L": ("d1a3o14_250525_imit_15L_2U", 800, 1.5 * L0,
os.path.join(_SRC2, "SR_analysis", "data", "illusion", "illusion_1.5L")),
}
model_name, si, diam, ref_dir = scenes[args.case]
out_dir = os.path.join(OUT_BASE, args.case)
run_illusion(args.device, args.case, model_name, si, diam, ref_dir, out_dir)
elif "vortex" in args.case:
vtype = "lamb" if "lamb" in args.case else "taylor"
scenes = {
"vortex_lamb": ("vortex_lamb", 0.5 * U0,
os.path.join(_SRC2, "SR_analysis", "data", "vortex", "vortex_lamb")),
"vortex_taylor": ("vortex_taylor", 0.03 * U0,
os.path.join(_SRC2, "SR_analysis", "data", "vortex", "vortex_taylor")),
}
model_name, strength, ref_dir = scenes[args.case]
out_dir = os.path.join(OUT_BASE, args.case)
run_vortex(args.device, vtype, model_name, strength, ref_dir, out_dir)
print("\nDone.")