SR Analysis: Phase-state SINDy + ablation study + documentation

Core changes:
- New phase-state features (PHASE_STATE_KEYS, ILLUSION_PHASE_KEYS) with obs dynamics
- Derivative and absolute output modes (output_mode="deriv"|"absolute")
- predict_v23_deriv() with integration support for closed-loop
- Offline multi-step rollout evaluator (eval_rollout.py)

Key results:
- Illusion 0.75L/1L: phase-state+error-state+abs achieves 0.974/0.958 closed-loop
  with zero action history features — proving the new route works
- Karman re100: phase-state+abs reaches 0.699 (vs 0.901 with action history)
- 1.5L confirmed as bang-bang regime (R2=0.12 for linear SINDy)
- Feature ablation: 6-dim phase-state outperforms 16-dim full-lag in closed-loop

Documentation:
- docs/SR_analysis_results.md: comprehensive analysis report
- docs/HANDOVER_SR_ANALYSIS.md: handover notes for next coder
- 6 figures in docs/figures/SR_analysis/
- Updated README.md, sindy_sr_notes.md, sindy_sr_knowledge.md
- Updated configs.py with generalization scenes

Co-authored-by: Cursor <cursoragent@cursor.com>
This commit is contained in:
Frank14f 2026-06-22 16:55:03 +08:00
parent dd71af901d
commit 8e62716ce4
26 changed files with 3591 additions and 630 deletions

View File

@ -0,0 +1,81 @@
# SR Analysis Pipeline — Handover Notes (2026-06-15)
## 交接人 → 接手人
### 当前管线状态
| 模块 | 状态 | 说明 |
|------|:----:|------|
| Illusion 0.75L + 1L phase-state + abs | **已验证** | 闭环 0.974 / 0.958,无动作历史,可进 PySR |
| Illusion 1.5L | **边界 case** | bang-bang 机制,线性 SINDy 不适用 |
| Karman phase-state + abs | **0.699** | 优于 deriv 模式,但低于旧 v23 的 0.901 |
| Karman old v23 (a_lag) | **0.901** | 保留作基线对照 |
| Karman 泛化测试 | **已完成** | Re70/150/300/25 约 0.54-0.60 |
| PySR 符号回归 | **有 shell** | 需要修复 run_pysr.py 后重新跑 |
| Vortex 偏移扩展 | **未做** | 低优先级 |
### 核心文件改动2026-06-14~15
| 文件 | 改动类型 |
|------|----------|
| `utils/feature_builder.py` | **新增** PHASE_STATE_KEYS, ILLUSION_PHASE_KEYS, KARMAN_EXPANDED_KEYS, obs dynamics, error-state, mu_Cl_tot |
| `utils/sindy_fitter.py` | **新增** compute_action_deriv, get_feature_matrix_deriv(output_mode) |
| `sindy/run_all_v2.py` | **新增** --deriv, --phase, --karman-expand, --karman-mu, --output-mode, --augment-level CLI |
| `validate/run_closed_loop.py` | **新增** predict_v23_deriv, mode="abs", load_sindy_coefs返回mode |
| `validate/run_closed_loop_illusion.py` | **修改** 支持 predict_v23_deriv, 自动检测模式 |
| `validate/eval_rollout.py` | **新文件** 离线多步 rollout 评估 |
| `scripts/plot_sr_results.py` | **新文件** 结果可视化图表 |
| `docs/SR_analysis_results.md` | **新文件** 完整分析报告 |
| `docs/figures/SR_analysis/fig*.png` | **新文件** 6 张图表 |
### 关键设计决策(接手前必读)
1. **phase-state 特征** = `u_a, du_a/dt, Cl_tot, dCl_tot/dt, Cd_tot, Cd_rear` (6维)
2. **ILLUSION_PHASE_KEYS** = phase-state + `Cd_err, Cl_err, dCd_err/dt, dCl_err/dt` (10维)
3. **绝对动作输出** `output_mode="absolute"` 优于导数 `"deriv"`,无积分累积
4. **v23 结构**始终默认front no-bias, rear shared-head
5. **时间一阶导**统一除以 `dt_c = SAMPLE_INTERVAL/2000`,跨场景可比较
6. `controlled.npz` 中新加了 `target_forces` 字段illusion 场景必须有
7. **FIFO bias ≠ DRL action bias**1U vs 2U 不要混淆
### 常用命令速查
```bash
# 拟合 + 验证 Illusion phase-state + absolute (完整流程)
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes illusion_0.75L,illusion_1L --deriv --phase --output-mode absolute
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop_illusion.py \
--scene illusion_0.75L --device 0 --steps 320 \
--sindy-results src/SR_analysis/sindy/illusion/sindy_results_deriv.json
# 拟合 + 验证 Karman phase-state + absolute
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re100 --deriv --phase --output-mode absolute
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \
--scene karman_re100 --device 0 --steps 200 --mode abs \
--sindy-results src/SR_analysis/sindy/karman/sindy_results_deriv.json
# 离线 rollout 评估
python3 src/SR_analysis/validate/eval_rollout.py \
--sindy-results src/SR_analysis/sindy/karman/sindy_results_deriv.json \
--scene karman_re100
# PySR (需要先修复滞后的 bug)
conda run -n sr_env python src/SR_analysis/sindy/run_pysr.py --scene illusion_1L
# 画图
python3 scripts/plot_sr_results.py
```
### 目前最适合推进的方向
1. **Illusion PySR 符号回归**0.75L + 1L separate → 公式比较)
2. **Karman 状态补强**(配合 CCD/OID 分析找出缺失的状态量,再回 SR
3. **Karman 跨 Re 联合**(在 phase-state + mu 基础上做跨 Re 联合拟合 + 泛化)
### 注意重新运行 run_pysr.py
当前 `run_pysr.py` 有路径/导入问题,接手后需先确认:
- `env_sr``sr_env` 环境内的 PySR 可用性
- whitelist 特征与 `controlled.npz` 中的字段匹配(特别是 illusion 的 `target_forces`
- 滞后构造正确(`a_prev[1:] = actions_phys[:-1]`

225
docs/SR_analysis_results.md Normal file
View File

@ -0,0 +1,225 @@
# SR Analysis: Phase-State SINDy Results Report
> Date: 2026-06-15
> Project: DynamisLab — Active hydrodynamic cloaking and illusion using DRL on a fluidic pinball.
> Analysis pipeline: SINDy (STLSQ) + feature engineering for interpretable control law extraction.
---
## Executive Summary
Three key findings from the SR analysis pipeline:
1. **Illusion 0.75L + 1L: New route validated.** Phase-state features + error-state + absolute action output achieves closed-loop similarity of **0.96+** (97% of PPO) with **zero action history features**. This proves that physically meaningful control laws can be extracted without relying on action memory.
2. **Illusion 1.5L: Regime shift identified.** The 1.5L target exhibits bang-bang/saturated control (alpha range [-8, 8], autocorrelation r=0.07). Linear SINDy is fundamentally inadequate (R2=0.12). This is a regime boundary, not a modelling failure.
3. **Karman: State representation still incomplete.** The same phase-state approach reaches 0.699 (vs 0.901 baseline with action history). The problem is not the output form (derivative vs absolute) but insufficient input state information for the Karman scene.
---
## 1. Methodology
### 1.1 Feature Architecture (Final)
Three feature levels were tested for fitting `alpha = f(state)` or `d(alpha)/dt = g(state)`:
| Level | Features | Dim | Description |
|-------|----------|:---:|-------------|
| Static | u_m, u_a, u_c, v_a, Cd_tot, Cd_rear, Cl_tot, Cl_diff | 8 | Current-step physics only, no memory |
| Phase-state | u_a, du_a/dt, Cl_tot, dCl_tot/dt, Cd_tot, Cd_rear (+error terms) | 6+4 | Oscillation phase + rate + drag feedback |
| Full-lag | Static + lag-1 of all 8 | 16 | Brute-force temporal context (baseline) |
### 1.2 Output Modes
Two output targets were compared:
- **Derivative mode**: predict d(alpha)/dt, then integrate `alpha(t) = alpha(t-1) + dt_c * dalpha/dt`
- **Absolute mode**: predict alpha directly, no integration needed
### 1.3 Evaluation Metrics
Models are evaluated on:
1. **One-step R2**: fit quality on training data
2. **Offline multi-step rollout**: 1/5/20/50 step recursive prediction on held-out PPO data
3. **CFD closed-loop**: full CFD environment with SINDy control law, DTW similarity vs target
---
## 2. Illusion Results
### 2.1 Three-Scene Summary
| Scene | S | Old v23 (a_lag) | New phase+error+abs | PPO baseline | % of PPO | Features | Action history? |
|------|:--:|:----------------:|:-------------------:|:------------:|:--------:|----------|:---------------:|
| 0.75L | 400 | 0.908 | **0.974** | 0.972 | 100.2% | 10-dim (ILLUSION_PHASE) | **No** |
| 1L | 600 | 0.962 | **0.958** | 0.973 | 98.5% | 10-dim (ILLUSION_PHASE) | **No** |
| 1.5L | 800 | 0.926 | **N/A** | 0.945 | — | — | Bang-bang regime |
**Key insight**: The 0.75L scene's new route **outperforms** the old v23 (0.974 vs 0.908), while 1L matches it within 1.5%. This definitively proves that "action history is not necessary."
### 2.2 Illusion 1L Front Coefficients (Phase-State, Absolute Output)
```
alpha_F = f(u_a, du_a/dt, Cl_tot, dCl_tot/dt, Cd_tot, Cd_rear, Cd_err, Cl_err, dCd_err/dt, dCl_err/dt)
R2 = 0.987
```
Key contributors (sorted by |coef|): Cd_rear > Cd_tot > Cd_err > Cl_tot > Cl_err > dCl_err/dt > dCd_err/dt
The formula is dominated by **drag-based feedback** (Cd_tot, Cd_rear, Cd_err), with the oscillation phase (u_a, Cl_tot) providing secondary modulation.
### 2.3 Illusion 1.5L — Regime Shift Evidence
| Metric | 1L | 1.5L | Ratio |
|--------|:--:|:----:|:-----:|
| Alpha range (front) | [-0.36, 1.11] | **[-8.0, 8.0]** | 10x |
| 1-step autocorrelation | 0.957 | **0.065** | — |
| d(alpha)/dt std | 0.12 | **6.07** | 50x |
| Linear SINDy R2 | 0.987 | **0.124** | — |
The 1.5L controller operates at saturation limits (maximum rotation speed), flipping between extremes. This is consistent with the frequency-doubling control strategy reported in the confirmation report.
---
## 3. Karman Results
### 3.1 Ablation Study (re100)
| Configuration | Feat | Output | R2 | Closed-loop | % of PPO |
|--------------|:----:|:-----:|:--:|:----------:|:--------:|
| old v23 (a_lag1 dominant) | 14+3 | alpha | 0.996 | **0.901** | 94.4% |
| **Phase-state -> abs (best new)** | **6** | **alpha** | **0.965** | **0.699** | 73.3% |
| Phase-state -> deriv | 6 | dalpha/dt | 0.837 | 0.656 | 68.8% |
| Phase-state + mu -> abs | 9 | alpha | 0.979 | 0.700 | 73.4% |
| Expanded 10-dim -> abs | 10 | alpha | 0.980 | **0.580** | 60.8% |
| Full-lag -> deriv | 16 | dalpha/dt | 0.939 | 0.619 | 64.9% |
| Static -> deriv | 8 | dalpha/dt | 0.321 | 0.745 | 78.1% |
**Key insight from ablation**:
1. **Phase-state + absolute output is the best new route** (0.699), but still well below the action-history baseline.
2. Adding extra static features (expanded 10-dim) **hurts** closed-loop despite higher R2 — classic covariate shift.
3. The static->deriv paradox: low R2 (0.321) but good closed-loop (0.745), because no-memory models are naturally robust to rollout divergence.
4. Mu modulation doesn't help at single-Re; its value will appear in cross-Re fitting.
### 3.2 Cross-Re Generalization (old v23 model)
| Re | Type | Closed-loop | Note |
|:--:|:----:|:----------:|------|
| 50 | Training | 0.582 | Low-frequency shedding |
| 100 | Training | **0.901** | Default / best |
| 200 | Training | 0.793 | Moderate degradation |
| 400 | Training | 0.664 | High-Re challenge |
| 25 | Unseen (subcritical) | 0.567 | Below Hopf bifurcation |
| 70 | Unseen | 0.577 | Between Re50 and Re100 |
| 150 | Unseen | 0.595 | Between Re100 and Re200 |
| 300 | Unseen | 0.541 | Outer extrapolation |
### 3.3 One-Step R2 vs Closed-Loop Paradox
A key phenomenon discovered during this work: **high one-step R2 does not predict good closed-loop performance** when temporal features (lags, derivatives) are present in the input. This is because:
- **Training**: features use ground-truth PPO observations
- **Deployment**: features use SINDY-controlled observations (distribution shift)
- Temporal features amplify this shift recursively
---
## 4. Roadmap to Next Steps
### 4.1 Illusion: Proceed to PySR
The 0.75L and 1L scenes are ready for symbolic regression. Recommended input:
```python
ILLUSION_PHASE_KEYS = [
"u_a", "du_a_dt", # oscillation phase
"Cl_tot", "dCl_tot_dt", # lift dynamics
"Cd_tot", "Cd_rear", # drag feedback
"Cd_err", "Cl_err", # force error
"dCd_err_dt", "dCl_err_dt", # error dynamics
]
```
Output: absolute alpha (non-dimensional action), no integration.
First do separate PySR for each scene, then compare formula structure for shared backbone.
### 4.2 Karman: CCD/OID for state completion
The 0.699 ceiling suggests missing state variables. Candidate directions:
- Recirculation zone length / reattachment point
- Wake centerline deflection
- POD mode coefficients for phase completion
- CCD modes correlated with control action
### 4.3 SR / CCD / OID Integration Framework
```
┌─────────────────────────────────────┐
│ Control Objective │
│ (stealth / illusion / erase) │
└──────────┬──────────────────────────┘
┌──────────▼──────────┐
│ DRL Policy │
│ (PPO + Sin act) │
└──────────┬──────────┘
┌────────────────┼────────────────┐
│ │ │
┌────────▼──────┐ ┌─────▼──────┐ ┌──────▼─────────┐
│ SINDy / SR │ │ CCD / OID │ │ Validation │
│ obs -> act │ │ obs -> z │ │ CFD closed- │
│ white-box │ │ structure │ │ loop + DTW │
│ control law │ │ analysis │ │ │
└────────┬──────┘ └─────┬──────┘ └──────┬─────────┘
│ │ │
└────────────────┼─────────────────┘
┌──────────▼──────────┐
│ Interpretable │
│ Control Mechanics │
│ obs -> z -> act │
│ -> structure -> │
│ -> signature │
└─────────────────────┘
```
---
## 5. Code Changes Summary
### Files Modified (8 core files)
| File | Changes |
|------|---------|
| `src/SR_analysis/utils/feature_builder.py` | PHYSICS_FEAT_KEYS, ILLUSION_ERR_KEYS, PHASE_STATE_KEYS, ILLUSION_PHASE_KEYS, KARMAN_EXPANDED_KEYS; obs dynamics (lag1 + derivative); error-state computation; 1D target_forces support |
| `src/SR_analysis/utils/sindy_fitter.py` | compute_action_deriv(); get_feature_matrix_deriv() with output_mode="deriv"|"absolute"; include_mu support |
| `src/SR_analysis/sindy/run_all_v2.py` | run_single_scene_deriv(); run_joint_karman_deriv(); --deriv, --phase, --karman-expand, --karman-mu, --output-mode, --augment-level CLI |
| `src/SR_analysis/validate/run_closed_loop.py` | predict_v23_deriv() with output_mode; mode="abs" branch; load_sindy_coefs returns "mode" |
| `src/SR_analysis/validate/run_closed_loop_illusion.py` | Support predict_v23_deriv; auto-detect SINDy mode from coefs |
| `src/SR_analysis/utils/__init__.py` | Export new constants and functions |
| `src/SR_analysis/sindy/wrap_joint.py` | Parameterized for karman/illusion |
| `src/SR_analysis/validate/eval_rollout.py` | **New**: offline multi-step rollout evaluation |
### Key Design Principles
1. **No action history in features** for all new routes (PHASE_STATE_KEYS, ILLUSION_PHASE_KEYS)
2. **Time-normalized derivatives** for cross-scene compatibility: `dx/dt = (x(t) - x(t-1)) / dt_c`
3. **Error-state** for Illusion: Cd_err, Cl_err encode current-to-target deviation
4. **Output flexibility**: both derivative and absolute modes supported
---
## 6. Figures
All figures in `docs/figures/SR_analysis/`:
- **fig1_illusion_comparison.png**: Illusion 3-scene comparison bar chart
- **fig2_karman_ablation.png**: Karman re100 ablation across 7 configurations
- **fig3_karman_generalization.png**: Karman cross-Re generalization (training vs unseen)
- **fig4_r2_vs_closedloop.png**: One-step R2 vs closed-loop paradox scatter
- **fig5_illusion_coefficients.png**: Illusion 1L front feature coefficients
- **fig6_roadmap.png**: Research progress roadmap

Binary file not shown.

After

Width:  |  Height:  |  Size: 35 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 63 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 55 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 48 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 31 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 34 KiB

254
scripts/plot_sr_results.py Normal file
View File

@ -0,0 +1,254 @@
#!/usr/bin/env python3
"""
Generate SR analysis result charts for reporting.
Output: docs/figures/SR_analysis/ (PNG files)
"""
import json, os, sys, glob
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
import matplotlib.ticker as ticker
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
_OUT = os.path.join(_REPO, "docs", "figures", "SR_analysis")
os.makedirs(_OUT, exist_ok=True)
# ---------------------------------------------------------------------------
# Data
# ---------------------------------------------------------------------------
# Illusion three-scenario comparison (new SR route: phase-state + error-state + absolute)
illusion_results = {
"0.75L": {"old_v23": 0.908, "new_phase": 0.974, "ppo": 0.972, "model": "d1a3o14_250525_imit_075L_2U_400S", "S": 400},
"1L": {"old_v23": 0.962, "new_phase": 0.958, "ppo": 0.973, "model": "d1a3o14_250525_imit_1L_2U_600S", "S": 600},
"1.5L": {"old_v23": 0.926, "new_phase": "N/A", "ppo": 0.945, "model": "d1a3o14_250525_imit_15L_2U", "S": 800},
}
# Karman ablation
karman_ablation = {
"old v23\n(a_lag+da)": 0.901,
"static->deriv\n(8dim)": 0.745,
"full-lag->deriv\n(16dim)": 0.619,
"phase->deriv\n(6dim)": 0.656,
"phase->abs\n(6dim)": 0.699,
"phase+mu->abs\n(9dim)": 0.700,
"expanded->abs\n(10dim)": 0.580,
}
# Karman generalization
karman_gen = {
"Re25": 0.567, "Re50": 0.582, "Re70": 0.577,
"Re100": 0.901, "Re150": 0.595, "Re200": 0.793,
"Re300": 0.541, "Re400": 0.664
}
# Ablation one-step R2
ablation_r2 = {
"static\n0": 0.321,
"phase\n6": 0.837,
"phase+abs\n6": 0.965,
"full-lag\n16": 0.939,
"expanded\n10": 0.980,
"phase+mu\n9": 0.979,
}
# ---------------------------------------------------------------------------
# Color scheme
# ---------------------------------------------------------------------------
C_OLD = "#d62728" # red
C_NEW_PHASE = "#2ca02c" # green
C_PPO = "#1f77b4" # blue
C_DERIV = "#ff7f0e" # orange
C_ABS = "#2ca02c" # green
C_GEN = "#9467bd" # purple
C_BG = "#f0f0f0"
# ---------------------------------------------------------------------------
# Fig 1: Illusion comparison bar chart
# ---------------------------------------------------------------------------
fig, ax = plt.subplots(figsize=(10, 5))
labels = list(illusion_results.keys())
x = np.arange(len(labels))
w = 0.25
old_vals = [illusion_results[k]["old_v23"] for k in labels]
new_vals = [illusion_results[k]["new_phase"] for k in labels]
ppo_vals = [illusion_results[k]["ppo"] for k in labels]
# Convert N/A to NaN for plotting
new_vals_plot = [v if isinstance(v, (int, float)) else 0 for v in new_vals]
ax.bar(x - w, old_vals, w, label="Old v23 (动作历史)", color=C_OLD, alpha=0.8)
ax.bar(x, new_vals_plot, w, label="New phase-state (无动作历史)", color=C_NEW_PHASE, alpha=0.8)
ax.bar(x + w, ppo_vals, w, label="PPO 基线", color=C_PPO, alpha=0.6)
# 1.5L label
ax.text(x[2], 0.05, "N/A\n(bang-bang)", ha="center", va="bottom", fontsize=9, color="gray")
ax.set_xticks(x)
ax.set_xticklabels(labels)
ax.set_ylabel("闭环相似度 (DTW)")
ax.set_title("Figure 1: Illusion 三场景 — 新路线 vs 旧版 vs PPO 基线", fontsize=12)
ax.legend(fontsize=9)
ax.set_ylim(0, 1.05)
ax.grid(axis="y", alpha=0.3)
fig.tight_layout()
fig.savefig(os.path.join(_OUT, "fig1_illusion_comparison.png"), dpi=150)
print("Saved fig1_illusion_comparison.png")
plt.close(fig)
# ---------------------------------------------------------------------------
# Fig 2: Karman ablation bar chart
# ---------------------------------------------------------------------------
fig, ax = plt.subplots(figsize=(12, 4.5))
keys = list(karman_ablation.keys())
vals = list(karman_ablation.values())
colors = []
for k in keys:
if "old" in k.lower(): colors.append(C_OLD)
elif "abs" in k: colors.append(C_ABS)
elif "deriv" in k: colors.append(C_DERIV)
else: colors.append("#7f7f7f")
bars = ax.bar(range(len(keys)), vals, color=colors, alpha=0.85)
ax.axhline(y=0.901, color=C_OLD, linestyle="--", alpha=0.5, label="Old v23 baseline (0.901)")
ax.set_xticks(range(len(keys)))
ax.set_xticklabels(keys, fontsize=8, rotation=20, ha="right")
ax.set_ylabel("闭环相似度")
ax.set_title("Figure 2: Karman re100 消融实验 — 输入/输出形式对比", fontsize=12)
ax.legend(fontsize=9)
ax.set_ylim(0, 1.0)
ax.grid(axis="y", alpha=0.3)
# Add red dashed line at 0.699 highlighting best new route
ax.axhline(y=0.699, color=C_ABS, linestyle=":", alpha=0.5)
ax.text(5.5, 0.705, "phase+abs 最佳: 0.699", fontsize=8, color=C_ABS)
fig.tight_layout()
fig.savefig(os.path.join(_OUT, "fig2_karman_ablation.png"), dpi=150)
print("Saved fig2_karman_ablation.png")
plt.close(fig)
# ---------------------------------------------------------------------------
# Fig 3: Karman generalization across Re
# ---------------------------------------------------------------------------
fig, ax = plt.subplots(figsize=(10, 4.5))
re_vals = sorted(karman_gen.keys(), key=lambda s: int(s.replace("Re","")))
sims = [karman_gen[k] for k in re_vals]
ax.plot(range(len(re_vals)), sims, "o-", color=C_GEN, linewidth=2, markersize=8)
# Mark training Re
train_re = [0, 3, 5, 7] # indices of Re50/100/200/400
for i in train_re:
ax.plot(i, sims[i], "o", color=C_OLD, markersize=12, markeredgecolor="black", markeredgewidth=1.5)
ax.set_xticks(range(len(re_vals)))
ax.set_xticklabels(re_vals)
ax.set_ylabel("闭环相似度")
ax.set_title("Figure 3: Karman 跨 Re 泛化 (旧 v23 模型)", fontsize=12)
ax.set_xlabel("Reynolds Number (2D reference)")
ax.axhline(y=0.5, color="gray", linestyle=":", alpha=0.5)
ax.grid(axis="y", alpha=0.3)
ax.set_ylim(0, 1.0)
# Annotations
ax.annotate("训练 Re", xy=(1.8, 0.92), fontsize=9, color=C_OLD)
ax.annotate("泛化 Re\n(未见过的)", xy=(4.5, 0.55), fontsize=9, color=C_GEN)
fig.tight_layout()
fig.savefig(os.path.join(_OUT, "fig3_karman_generalization.png"), dpi=150)
print("Saved fig3_karman_generalization.png")
plt.close(fig)
# ---------------------------------------------------------------------------
# Fig 4: One-step R2 vs closed-loop scatter (diagnostic)
# ---------------------------------------------------------------------------
fig, ax = plt.subplots(figsize=(7, 5.5))
points = [
("old v23", 0.996, 0.901, C_OLD),
("static→deriv", 0.321, 0.745, C_DERIV),
("full-lag→deriv", 0.939, 0.619, "#7f7f7f"),
("phase→deriv", 0.837, 0.656, C_DERIV),
("phase→abs", 0.965, 0.699, C_ABS),
("expanded→abs", 0.980, 0.580, "#7f7f7f"),
("phase+mu→abs", 0.979, 0.700, C_ABS),
]
for name, r2, sim, color in points:
ax.scatter(r2, sim, s=100, color=color, zorder=5)
ax.annotate(name, (r2, sim), textcoords="offset points", xytext=(5, 5), fontsize=8)
ax.set_xlabel("One-step R²")
ax.set_ylabel("CFD 闭环相似度")
ax.set_title("Figure 4: Karman re100 — One-step R² 与闭环不一致性", fontsize=12)
ax.grid(alpha=0.3)
# Upper-left region = good closed-loop, bad one-step (static-deriv)
# Upper-right region = good both (old v23)
# Lower-right region = good one-step, bad closed-loop (most new methods)
ax.annotate("稳健欠拟合", xy=(0.15, 0.85), fontsize=9, color="gray", fontstyle="italic")
ax.annotate("分布偏移\n(训练好, 闭环差)", xy=(0.75, 0.45), fontsize=9, color="gray", fontstyle="italic")
fig.tight_layout()
fig.savefig(os.path.join(_OUT, "fig4_r2_vs_closedloop.png"), dpi=150)
print("Saved fig4_r2_vs_closedloop.png")
plt.close(fig)
# ---------------------------------------------------------------------------
# Fig 5: Phase-state feature coefficients (Illusion 1L)
# ---------------------------------------------------------------------------
# Load the SINDy results for illusion 1L phase-state
try:
with open(os.path.join(_REPO, "src/SR_analysis/sindy/illusion/sindy_results_deriv.json")) as f:
sr = json.load(f)
per = sr["per_scene"]["illusion_1L"]
fn_f = per["feature_names_front"]
coef_f = per["front"]["best_coef"][:len(fn_f)]
fig, ax = plt.subplots(figsize=(8, 4.5))
# Sort by |coef|
pairs = sorted(zip(fn_f, coef_f), key=lambda p: -abs(p[1]))
names = [p[0] for p in pairs]
vals = [p[1] for p in pairs]
colors_bar = [C_NEW_PHASE if v > 0 else C_OLD for v in vals]
ax.barh(range(len(names)), vals, color=colors_bar, alpha=0.8)
ax.set_yticks(range(len(names)))
ax.set_yticklabels(names)
ax.axvline(x=0, color="black", linewidth=0.5)
ax.set_xlabel("系数值")
ax.set_title("Figure 5: Illusion 1L Front — Phase-state 特征系数", fontsize=12)
ax.grid(axis="x", alpha=0.3)
fig.tight_layout()
fig.savefig(os.path.join(_OUT, "fig5_illusion_coefficients.png"), dpi=150)
print("Saved fig5_illusion_coefficients.png")
plt.close(fig)
except Exception as e:
print(f"fig5 skipped: {e}")
# ---------------------------------------------------------------------------
# Fig 6: Summary timeline / roadmap
# ---------------------------------------------------------------------------
fig, ax = plt.subplots(figsize=(10, 3.5))
phases = [
("Phase 0\nBug Audit", "2026-06-12\n12 bugs\nfixed", 0.8),
("Phase 1\nTarget fix", "2026-06-13\nIllusion target\ninfo added", 0.85),
("Phase 2\nAblation", "2026-06-14\nPhase-state\nvalidated", 0.90),
("Phase 2b\nOutput mode", "2026-06-15\nPhase+abs\nIllusion 0.97", 0.95),
("Phase 3\nIllusion SR", "Next\nPySR on\n0.75L/1L", 0.7),
]
y_pos = 1
for i, (label, desc, conf) in enumerate(phases):
color = plt.cm.RdYlGn(conf)
ax.barh(y_pos, 1, left=i, height=0.5, color=color, alpha=0.8)
ax.text(i + 0.5, y_pos, label, ha="center", va="center", fontsize=8, fontweight="bold")
ax.text(i + 0.5, y_pos - 0.3, desc, ha="center", va="top", fontsize=6, color="gray")
ax.set_xlim(0, len(phases))
ax.set_ylim(0, 2)
ax.axis("off")
ax.set_title("Figure 6: 研究进展路线图", fontsize=12)
fig.tight_layout()
fig.savefig(os.path.join(_OUT, "fig6_roadmap.png"), dpi=150)
print("Saved fig6_roadmap.png")
plt.close(fig)
print(f"\nAll figures saved to {_OUT}/")

View File

@ -12,278 +12,174 @@ cloak and illusion scenes, using dimensionless physical features,
G-equivariant structural constraints, and STLSQ threshold grids. G-equivariant structural constraints, and STLSQ threshold grids.
For background, see: For background, see:
- `src/sindy_sr_notes.md` -- execution plan - `sindy_sr_notes.md` -- execution plan and task tracking
- `src/sindy_sr_knoeledge.md` -- confirmed facts and known pitfalls - `sindy_sr_knowledge.md` -- confirmed facts and known pitfalls
- `../../docs/SR_analysis_results.md` -- comprehensive results report
---
## Directory Structure ## Directory Structure
``` ```
SR_analysis/ SR_analysis/
configs.py # Unified scene metadata (all 10 scenes) configs.py # Unified scene metadata (all 10+ scenes)
configs/ configs/
legacy/ # Legacy CFD configs (config_cuda.json, config_flowfield.json) legacy/ # Legacy CFD configs
utils/ utils/
__init__.py # Selective exports (no pycuda dependency) __init__.py # Selective exports (no pycuda dependency)
feature_builder.py # Dimensionless features + G-operator (from analysis_cloak) feature_builder.py # Dimensionless features + G-operator + phase-state features
sindy_fitter.py # STLSQ threshold grid, feature matrix builder sindy_fitter.py # STLSQ + feature matrices + derivative/absolute modes
cfd_interface.py # LegacyCelerisLab wrapper (requires pycuda_3_10) cfd_interface.py # LegacyCelerisLab wrapper (requires pycuda_3_10)
g_operator.py # Equivariance diagnostics g_operator.py # Equivariance diagnostics
data/ data/
karman/ # Karman cloak: karman_re50, re100, re200, re400 karman/ # Karman cloak: karman_re50/100/200/400
steady/ # Steady cloak: steady_data.npz steady/ # Steady cloak
illusion/ # Illusion: illusion_0.75L, illusion_1L, illusion_1.5L illusion/ # Illusion: illusion_0.75L/1L/1.5L
vortex/ # Vortex cloak: vortex_lamb, vortex_taylor vortex/ # Vortex cloak
scripts/ scripts/
infer_karman.py # Inference: LegacyCFD + PPO -> controlled.npz infer_karman.py # Inference: LegacyCFD + PPO -> controlled.npz
infer_illusion.py # Inference: for 0.75L, 1L, 1.5L diameters infer_illusion.py # Inference for illusion scenes
infer_vortex.py # Inference: for Lamb dipole + Taylor monopole infer_vortex.py # Inference for vortex scenes
sindy/ sindy/
run_karman.py # SINDy fitting for Karman scenes run_all_v2.py # Unified SINDy fitting (supports --deriv, --phase, --output-mode etc.)
run_illusion.py # SINDy fitting for Illusion scenes run_pysr.py # Restricted PySR symbolic regression
run_vortex.py # SINDy fitting for Vortex scenes wrap_joint.py # Joint model -> wrapped format for validator
run_pareto.py # Pareto-front analysis from SINDy results compare_v2.py # Cross-scene comparison report
karman/ # Output: sindy_results.json, pareto_*.json karman/illusion/vortex/ # SINDy output JSONs
illusion/ # Output: sindy_results.json, pareto_*.json
vortex/ # Output: sindy_results.json, pareto_*.json
validate/ validate/
run_closed_loop.py # Unified closed-loop validator (v23 + unstructured modes) run_closed_loop.py # Karman closed-loop validator (v23/deriv/abs modes)
run_closed_loop_illusion.py # Illusion closed-loop validator
eval_rollout.py # Offline multi-step rollout evaluation
results/ # Validation result JSONs
compare/ compare/
support_overlap.py # Pairwise support set comparison support_overlap.py # Support set comparison
shared_core.py # Multi-scene shared-core detection shared_core.py # Shared core detection
``` ```
---
## Key Design Decisions ## Key Design Decisions
### 1. Scene Metadata Driven ### 1. Scene Metadata Driven
All scene parameters (Re, action scaling, geometry, model paths) are defined All scene parameters defined once in `configs.py`.
once in `configs.py`, not hard-coded in scripts. Adding a new scene means
adding one dict to `configs.py`.
### 2. Data / Features / Models Separation ### 2. Feature Levels
- `data/` -- raw sensor/force/action arrays (.npz), one-time generation | Level | Features | Dim | Description |
- `sindy/` -- SINDy fitting results (JSON), reusable for comparison |-------|----------|:---:|-------------|
- `scripts/` -- inference pipelines that produce `data/` | Static | u_m, u_a, u_c, v_a, Cd_tot, Cd_rear, Cl_tot, Cl_diff | 8 | Current-step only |
| **Phase-state** | u_a, du_a/dt, Cl_tot, dCl_tot/dt, Cd_tot, Cd_rear | **6** | Oscillation phase + rate |
| Illusion-phase | Phase-state + Cd_err, Cl_err, dCd_err/dt, dCl_err/dt | **10** | Phase + error-state |
| Karman-expanded | Phase-state + u_m, u_c, v_a, Cl_diff | **10** | Phase + supplementary |
| Full-lag | Static + lag-1 | 16 | Full temporal context |
### 3. Unified Feature Builder ### 3. Output Modes
Every scene uses the same `utils/feature_builder.py`, which produces - **deriv**: predict `d(alpha)/dt`, then `alpha(t) = alpha(t-1) + dt_c * dalpha/dt`
21 dimensionless features from raw lattice-unit sensor/force data: - **absolute**: predict `alpha(t)` directly (no integration drift)
**Sensor features (nondim):**
- `u_m`, `u_a`, `u_c` -- streamwise: mean, antisymmetric, centre
- `v_a` -- antisymmetric cross-stream
- `sin_ua`, `cos_ua` -- phase encoding via u_a
**Force features (Cd/Cl):**
- `Cd_tot`, `Cd_rear` -- total and rear-cylinder drag
- `Cl_tot`, `Cl_diff` -- total and differential lift
**Memory features (nondim alpha):**
- `aF_lag1`, `aB_lag1`, `aT_lag1` -- lagged actions (t-1)
- `daF`, `daB`, `daT` -- action increments (t-1)-(t-2)
**Reynolds modulation:**
- `mu` (= 1/Re_D), `mu_u_a`, `mu_v_a`, `mu_Cd_tot`, `mu_Cl_diff`
### 4. G-Equivariant Structure (v23) ### 4. G-Equivariant Structure (v23)
Default control law structure (confirmed as the best v23 model):
``` ```
Front(t) = f_front(x(t)) # no bias, odd under G Front(t) = f_front(x(t)) # no bias, odd under G
Top(t) = f_rear(x(t)) # with bias Top(t) = f_rear(x(t)) # with bias
Bottom(t) = -f_rear(G[x(t)]) # shared-head: bottom = -top(Gx) Bottom(t) = -f_rear(G[x(t)]) # shared-head
``` ```
Where G is the mirror operator (y -> -y) with corrected sign rules: ---
- `[aF, aT, aB] -> [-aF, -aB, -aT]`
- Sensor swap: top <-> bottom, negate v
- Force swap: front unchanged, bottom <-> top, negate Cl
### 5. STLSQ Threshold Grid ## Current Best Results (2026-06-15)
Default thresholds: `[0, 0.001, 0.002, 0.005, 0.01, 0.015, 0.02, 0.03, 0.05, 0.1]` ### Illusion — New Route: Phase-state + Error-state + Absolute Action
Per-channel: front (no bias), top (shared-head), bottom (independent, for comparison)
## Scene Inventory | Scene | Closed-loop | % of PPO | Action history? | Features |
|-------|:----------:|:--------:|:---------------:|----------|
| 0.75L | **0.974** | 100.2% | **No** | ILLUSION_PHASE (10dim) |
| 1L | **0.958** | 98.5% | **No** | ILLUSION_PHASE (10dim) |
| 1.5L | N/A | — | **No** | Bang-bang regime |
| Scene Name | Description | Re_code | Sample Interval | Action | U0 | ### Karman re100 — Ablation
|---|---|---|---|---|---|
| karman_re50 | Karman cloak at low Re | 50 | 800 | 8x + [0,-4,4] | 0.01 |
| karman_re100 | Karman cloak (default) | 100 | 800 | 8x + [0,-4,4] | 0.01 |
| karman_re200 | Karman cloak at high Re | 200 | 800 | 8x + [0,-4,4] | 0.01 |
| karman_re400 | Karman cloak at highest Re | 400 | 800 | 8x + [0,-4,4] | 0.01 |
| steady | Open-loop constant rotation | 100 | 800 | 8x + [0,-5.1,5.1] | 0.01 |
| illusion_0.75L | Imitate 0.75D cylinder | 100 | 600 | 8x + [0,-2,2] | 0.01 |
| illusion_1L | Imitate 1.0D cylinder | 100 | 600 | 8x + [0,-2,2] | 0.01 |
| illusion_1.5L | Imitate 1.5D cylinder | 100 | 600 | 8x + [0,-2,2] | 0.02 |
| vortex_lamb | Cloak Lamb dipole | 100 | 800 | 4x + [0,-4,4] | 0.01 |
| vortex_taylor | Cloak Taylor monopole | 100 | 800 | 4x + [0,-4,4] | 0.01 |
Note: "Re_code" uses reference length 2*D (code convention). | Config | Feat | Output | R2 | Closed-loop | Note |
Physical Re_D = Re_code / 2. E.g. Re_code=100 -> Re_D=50. |--------|:----:|:-----:|:--:|:----------:|------|
| old v23 (a_lag) | 14+3 | alpha | 0.996 | **0.901** | Baseline |
| **Phase->abs** | **6** | **alpha** | **0.965** | **0.699** | Best new route |
| Phase->deriv | 6 | dalpha/dt | 0.837 | 0.656 | |
| Phase+mu->abs | 9 | alpha | 0.979 | 0.700 | mu helps cross-Re |
| Expanded->abs | 10 | alpha | 0.980 | 0.580 | Overfitting |
## Re-generation Commands ---
All commands run from repo root (`/home/frank14f/DynamisLab`). ## Commands
### Data Generation (requires GPU, pycuda_3_10 env) All from repo root (`/home/frank14f/DynamisLab`).
### SINDy Fitting
```bash ```bash
# Karman cloak -- all 4 training Re # Illusion phase-state + absolute (recommended for 0.75L/1L)
conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_karman.py --re all --device 0 conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes illusion_0.75L,illusion_1L --deriv --phase --output-mode absolute
# Karman cloak -- single Re # Karman phase-state + absolute
conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_karman.py --re 100 --device 0 --steps 200 conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re100 --deriv --phase --output-mode absolute
# Illusion -- all 3 diameters # Karman expanded (10 dim)
conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_illusion.py --diameter all --device 0 conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re100 --deriv --karman-expand --output-mode absolute
# Vortex -- both types # Karman with mu modulation
conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_vortex.py --type all --device 0 conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re100 --deriv --karman-mu --output-mode absolute
# Old-style (v2, with action history)
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re50,karman_re100 --joint
``` ```
### SINDy Fitting (no GPU needed, pycuda_3_10 env for pysindy) ### Closed-loop Validation
```bash
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_karman.py
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_illusion.py
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_vortex.py
```
### Pareto Analysis (no GPU, no conda needed)
```bash
python3 src/SR_analysis/sindy/run_pareto.py --scene karman_re100
python3 src/SR_analysis/sindy/run_pareto.py --scene illusion_1L
```
### Closed-loop Validation (requires GPU)
```bash ```bash
# Karman with absolute action
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \ conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \
--scene karman_re70 --device 2 \ --scene karman_re100 --device 0 --steps 200 --mode abs \
--sindy-results src/SR_analysis/sindy/karman/sindy_results.json --sindy-results src/SR_analysis/sindy/karman/sindy_results_deriv.json
# With custom mode # Karman old v23
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \ conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \
--scene karman_re70 --device 2 --mode unstructured --scene karman_re100 --device 0 --steps 200 --mode v23 \
--sindy-results src/SR_analysis/sindy/karman/sindy_joint_wrapped.json
# Illusion with absolute action
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop_illusion.py \
--scene illusion_1L --device 0 --steps 320 \
--sindy-results src/SR_analysis/sindy/illusion/sindy_results_deriv.json
``` ```
### Cross-scene Comparison (no GPU) ### PySR Symbolic Regression
```bash ```bash
# Pairwise support overlap conda run -n sr_env python src/SR_analysis/sindy/run_pysr.py --scene illusion_1L
python3 src/SR_analysis/compare/support_overlap.py \
--sindy-results src/SR_analysis/sindy/karman/sindy_results.json \
--scenes karman_re100 illusion_1L
# Multi-scene shared core
python3 src/SR_analysis/compare/shared_core.py \
--sindy-results src/SR_analysis/sindy/karman/sindy_results.json \
--scenes karman_re50 karman_re100 karman_re200 karman_re400
``` ```
## Key Results Summary ### Offline Rollout Evaluation
### Data Quality (similarity scores) ```bash
python3 src/SR_analysis/validate/eval_rollout.py \
--sindy-results src/SR_analysis/sindy/karman/sindy_results_deriv.json \
--scene karman_re100
```
| Scene | PPO Similarity | ---
|---|---|
| karman_re50 | 0.962 |
| karman_re100 | 0.954 |
| karman_re200 | 0.884 |
| karman_re400 | 0.795 (inferred, not verified) |
| vortex_lamb | 0.942 |
| vortex_taylor | 0.916 |
| illusion_1L | ~0.55 (metric not directly comparable) |
### SINDy Fit Quality (R2 scores for one-step prediction) ## Important Reminders
| Scene | Front | Top (shared) | Bottom | - `controlled.npz` actions are **normalized [-1,1]** — must convert via `(norm * scale + bias) * u0`
|---|---|---|---| - **FIFO bias ≠ DRL action bias** for Illusion: FIFO=[0, -0.01, 0.01], decode=[0, -0.02, 0.02]
| karman_re50 | 0.998 | 0.989 | 0.996 | - "2U" in model name = S_DIM=14 (not 2x velocity), u0 always 0.01
| karman_re100 | 0.995 | 0.993 | 0.997 | - SAMPLE_INTERVAL: 0.75L=400, 1L=600, 1.5L/Karman=800
| karman_re200 | 0.957 | 0.914 | 0.918 | - Closed-loop steps auto-set: S=400→320, S=600→214, S=800→160
| karman_re400 | 0.991 | 0.979 | 0.969 | - One-step R² high ≠ closed-loop good — always validate
| illusion_0.75L | 0.991 | 0.989 | 0.990 | - For phase-state features, always pass `sensors_raw`/`forces_raw` to enable derivative computation
| illusion_1L | 0.979 | 0.984 | 0.984 |
| illusion_1.5L | 0.959 | 0.928 | 0.932 |
| vortex_lamb | 0.904 | 0.980 | 0.933 |
| vortex_taylor | 0.960 | 0.810 | 0.643 |
### Shared Core Features
**Karman cross-Re (active in all re50/100/200):**
- Front core: `mu`, `mu_Cd_tot`, `mu_Cl_diff`, `mu_v_a` (mu-modulated terms dominate)
- Top core: `Cl_tot`, `bias`, `mu_Cd_tot`, `mu_Cl_diff`, `mu_u_a`, `mu_v_a`
- Scene-specific: lower-Re scenes have additional `Cd_tot`, `Cl_diff`, `aT_lag1` etc.
**Illusion cross-diameter (active in all 0.75L/1L/1.5L):**
- Front core: `mu`, `mu_Cd_tot`, `mu_Cl_diff` (same structure as Karman front!)
- Top core: `Cd_rear`, `Cl_tot`, `bias`, `mu_Cd_tot`, `mu_Cl_diff`
- This suggests a **shared mu-modulated feedback structure** exists across both scenes
## Known Issues and Caveats
1. **Vortex Taylor rear channels** have low R2 (0.64-0.81). The weak monopole
produces near-zero rear action, making SINDy fitting noisy. Use Lamb as the
primary vortex reference.
2. **Closed-loop validator** (`validate/run_closed_loop.py`) has been ported but
NOT yet tested end-to-end. The original `validate_v23.py` verified Karman
but the new unified version has not been run.
3. **Illusion similarity scores** use the Karman CONV_LEN=30 metric, giving
lower raw numbers. The controlled.npz data itself is valid for SINDy.
4. **Steady cloak** is open-loop constant rotation, not PPO-derived. It serves
as a physical consistency check, not a primary comparison scene.
5. **SINDy one-step R2 is not sufficient** -- a high R2 does not guarantee good
closed-loop performance. Always validate via `validate/run_closed_loop.py`.
6. **Scene key naming**: keys like `illusion_1L`, `illusion_1.5L` use the short
float format from Python (1.0 -> "1L", 1.5 -> "1.5L", 0.75 -> "0.75L").
## Next Steps (Future Work)
1. **PySR symbolic regression** -- Run PySR on the SINDy-identified active
features (in `sr_env` conda env) to find closed-form formulas. Essential
reading: `src/pysr.md`.
2. **Closed-loop validation of all new scenes** -- Run
`validate/run_closed_loop.py` for illusion and vortex scenes using their
SINDy coefficients.
3. **Cross-scene shared backbone test** -- Fit a single SINDy model on merged
Karman + Illusion data, test if it performs on both.
4. **Time-scale explicit formulation** -- Make the sample interval an explicit
feature to compare control laws across different frequencies.
5. **Steady as consistency check** -- Validate that Karman-derived control laws
can reproduce the steady cloak result as a sanity check.
## File Reference
| File | Lines | Purpose |
|---|---|---|
| configs.py | ~205 | Unified scene metadata |
| utils/feature_builder.py | ~212 | Dimensionless features + G-op |
| utils/sindy_fitter.py | ~175 | STLSQ fitting, feature matrix builder |
| utils/cfd_interface.py | ~370 | LegacyCelerisLab wrapper |
| utils/g_operator.py | ~170 | Equivariance diagnostics |
| utils/__init__.py | ~10 | Selective exports |
| scripts/infer_karman.py | ~250 | Karman inference pipeline |
| scripts/infer_illusion.py | ~270 | Illusion inference pipeline |
| scripts/infer_vortex.py | ~280 | Vortex inference pipeline |
| sindy/run_karman.py | ~160 | Karman SINDy fitting |
| sindy/run_illusion.py | ~110 | Illusion SINDy fitting |
| sindy/run_vortex.py | ~110 | Vortex SINDy fitting |
| sindy/run_pareto.py | ~140 | Pareto analysis |
| validate/run_closed_loop.py | ~270 | Closed-loop validator |
| compare/support_overlap.py | ~150 | Pairwise support comparison |
| compare/shared_core.py | ~140 | Multi-scene shared core detection |

View File

@ -27,7 +27,7 @@ NX = 1280
NY = 512 NY = 512
CENTER_Y = (NY - 1) / 2.0 CENTER_Y = (NY - 1) / 2.0
FIFO_LEN = 150 FIFO_LEN = 150
CONV_LEN = 30 CONV_LEN = 30 # default; per-scene conv_len overrides this (Illusion=36)
def nu_from_re(re_code: float, u0: float = U0) -> float: def nu_from_re(re_code: float, u0: float = U0) -> float:
@ -72,7 +72,30 @@ for rc, mn in [(50, "d1a3o12_re50"), (100, "d1a3o12_re100"),
"u0": U0, "u0": U0,
} }
# -- Steady Cloak (open-loop constant rotation) ----------------------------- # -- Karman Cloak generalization (unseen Re, no PPO model) --------------------
for rc in [25, 70, 150, 300]:
key = f"karman_re{rc}"
SCENES[key] = {
"scene_id": "karman",
"re_code": rc,
"mu": 2.0 / rc,
"nu": nu_from_re(rc),
"has_disturbance": True,
"sample_interval": 800,
"action_scale": 8.0,
"action_bias": (0.0, -4.0, 4.0),
"source": "generalization",
"model_name": None,
"n_objects_env": 7,
"obs_slice": (2, 14),
"sensor_x": 40.0,
"pinball_front_x": 30.0,
"pinball_rear_x": 31.3,
"target_type": "periodic",
"s_dim": 12,
"u0": U0,
}
SCENES["steady"] = { SCENES["steady"] = {
"scene_id": "steady", "scene_id": "steady",
"re_code": 100, "re_code": 100,
@ -94,26 +117,29 @@ SCENES["steady"] = {
"u0": U0, "u0": U0,
} }
# -- Illusion (cylinder imitation, 3 diameters, 1U=0.01) -------------------- # -- Illusion (cylinder imitation, 3 diameters) ----------------------------
# "2U" in model name means S_DIM=14 (2 extra target force channels), NOT 2x velocity.
# u0 is always 0.01, nu is always 0.004 (default) unless Vis suffix present.
def _illusion_key(diam: float) -> str: def _illusion_key(diam: float) -> str:
"""Generate clean illusion scene key.""" """Generate clean illusion scene key."""
s = f"{diam:.3f}".rstrip("0").rstrip(".") s = f"{diam:.3f}".rstrip("0").rstrip(".")
return f"illusion_{s}L" return f"illusion_{s}L"
_ILLUSION_1U = [ _ILLUSION_BEST = [
(0.75, "d1a3o12_250525_imit_075L_1U"), (0.75, "d1a3o14_250525_imit_075L_2U_400S", 400, 14),
(1.0, "d1a3o12_250525_imit_1L_1U"), (1.0, "d1a3o14_250525_imit_1L_2U_600S", 600, 14),
] ]
for diam, mn in _ILLUSION_1U: for diam, mn, si, sd in _ILLUSION_BEST:
key = _illusion_key(diam) key = _illusion_key(diam)
SCENES[key] = { SCENES[key] = {
"scene_id": "illusion", "scene_id": "illusion",
"target_diameter": diam, "target_diameter": diam,
"re_code": 100, "re_code": 100,
"mu": 2.0 / 100, "mu": 2.0 / 100,
"nu": nu_from_re(100), "nu": 0.004, # default, no Vis suffix in model name
"has_disturbance": False, "has_disturbance": False,
"sample_interval": 600, "sample_interval": si,
"conv_len": 36, # Illusion uses 36, not default 30
"action_scale": 8.0, "action_scale": 8.0,
"action_bias": (0.0, -2.0, 2.0), "action_bias": (0.0, -2.0, 2.0),
"source": "PPO_inference", "source": "PPO_inference",
@ -124,19 +150,20 @@ for diam, mn in _ILLUSION_1U:
"pinball_front_x": 19.0, "pinball_front_x": 19.0,
"pinball_rear_x": 20.3, "pinball_rear_x": 20.3,
"target_type": "periodic", "target_type": "periodic",
"s_dim": 12, "s_dim": sd,
"u0": U0, "u0": U0, # always 0.01
} }
# 1.5L Illusion (2U=0.02 model) # 1.5L Illusion (SAMPLE_INTERVAL=800, default, no S suffix in name)
SCENES[_illusion_key(1.5)] = { SCENES[_illusion_key(1.5)] = {
"scene_id": "illusion", "scene_id": "illusion",
"target_diameter": 1.5, "target_diameter": 1.5,
"re_code": 100, "re_code": 100,
"mu": 2.0 / 100, "mu": 2.0 / 100,
"nu": nu_from_re(100, u0=0.02), "nu": 0.004, # default, no Vis suffix
"has_disturbance": False, "has_disturbance": False,
"sample_interval": 600, "sample_interval": 800, # no S suffix in name = default 800
"conv_len": 36, # Illusion uses 36
"action_scale": 8.0, "action_scale": 8.0,
"action_bias": (0.0, -2.0, 2.0), "action_bias": (0.0, -2.0, 2.0),
"source": "PPO_inference", "source": "PPO_inference",
@ -148,7 +175,7 @@ SCENES[_illusion_key(1.5)] = {
"pinball_rear_x": 20.3, "pinball_rear_x": 20.3,
"target_type": "periodic", "target_type": "periodic",
"s_dim": 14, "s_dim": 14,
"u0": 0.02, "u0": U0, # always 0.01
} }
# -- Vortex Cloak (Lamb dipole + Taylor monopole) -------------------------- # -- Vortex Cloak (Lamb dipole + Taylor monopole) --------------------------
@ -183,6 +210,37 @@ for vtype, mn, strength in _SCENES_VORTEX:
} }
# -- Illusion generalization (unseen diameters, no PPO model) ----------------
for diam, mn, si in [
(0.5, None, 400),
(1.2, None, 600),
(2.0, None, 800),
]:
key = _illusion_key(diam)
SCENES[key] = {
"scene_id": "illusion",
"target_diameter": diam,
"re_code": 100,
"mu": 2.0 / 100,
"nu": 0.004,
"has_disturbance": False,
"sample_interval": si,
"conv_len": 36,
"action_scale": 8.0,
"action_bias": (0.0, -2.0, 2.0),
"source": "generalization",
"model_name": mn,
"n_objects_env": 6,
"obs_slice": (0, 12),
"sensor_x": 30.0,
"pinball_front_x": 19.0,
"pinball_rear_x": 20.3,
"target_type": "periodic",
"s_dim": 14,
"u0": U0,
}
# -- Utility helpers --------------------------------------------------------- # -- Utility helpers ---------------------------------------------------------
def get_scene(name: str) -> dict: def get_scene(name: str) -> dict:

View File

@ -4,9 +4,9 @@ Generates controlled data for a given target cylinder diameter using
LegacyCelerisLab + trained PPO model. LegacyCelerisLab + trained PPO model.
Usage: Usage:
conda run -n pycuda_3_10 python scripts/infer_illusion.py \\ conda run -n pycuda_3_10 python scripts/infer_illusion.py \
--diameter 1.0 --device 0 --diameter 1.0 --device 0
conda run -n pycuda_3_10 python scripts/infer_illusion.py \\ conda run -n pycuda_3_10 python scripts/infer_illusion.py \
--diameter all --device 2 --diameter all --device 2
""" """
from __future__ import annotations from __future__ import annotations
@ -37,12 +37,61 @@ from SR_analysis.utils.cfd_interface import (
) )
from SR_analysis.configs import ( from SR_analysis.configs import (
get_scene, get_scene_list, model_path_for_scene, get_scene, get_scene_list, model_path_for_scene,
LEGACY_CFG_DIR, FIFO_LEN, CONV_LEN, LEGACY_CFG_DIR, FIFO_LEN,
) )
DATA_TYPE = np.float32 DATA_TYPE = np.float32
def analyze_harmonics(states: np.ndarray, n_harmonics: int = 5) -> list:
"""FFT-based harmonics analysis matching legacy_env_imit.py exactly.
Args:
states: (T, D) time series
n_harmonics: number of harmonics (excluding DC) to keep
Returns:
list of dict per channel: {dc, amps, freqs, phases}
"""
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 * 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': amps[idx].tolist(),
'freqs': freqs[idx].tolist(),
'phases': phases[idx].tolist(),
}
result.append(harmonics)
return result
def gen_target_states_at(t: int, harmonics: list) -> np.ndarray:
"""Reconstruct target forces at step t from harmonics.
Args:
t: current step index (0-based)
harmonics: from analyze_harmonics()
Returns:
ndarray of shape (D,) with reconstructed values
"""
D = len(harmonics)
result = np.zeros(D, dtype=np.float32)
for d, h in enumerate(harmonics):
val = 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] = float(val)
return result
def run_single_illusion( def run_single_illusion(
scene_name: str, scene_name: str,
device_id: int, device_id: int,
@ -55,13 +104,15 @@ def run_single_illusion(
u0 = cfg["u0"] u0 = cfg["u0"]
l0 = 20.0 l0 = 20.0
sample_interval = cfg["sample_interval"] sample_interval = cfg["sample_interval"]
conv_len = cfg.get("conv_len", 30) # Illusion uses 36
action_scale = cfg["action_scale"] action_scale = cfg["action_scale"]
action_bias = cfg["action_bias"] action_bias = cfg["action_bias"] # (front, bottom, top) for DRL decoding
n_obj_total = cfg["n_objects_env"] n_obj_total = cfg["n_objects_env"]
sensor_x = cfg["sensor_x"] # 30.0 for illusion sensor_x = cfg["sensor_x"] # 30.0 for illusion
front_x = cfg["pinball_front_x"] # 19.0 front_x = cfg["pinball_front_x"] # 19.0
rear_x = cfg["pinball_rear_x"] # 20.3 rear_x = cfg["pinball_rear_x"] # 20.3
target_diam = cfg["target_diameter"] target_diam = cfg["target_diameter"]
s_dim = cfg["s_dim"] # 14 for all best illusion models
os.makedirs(output_root, exist_ok=True) os.makedirs(output_root, exist_ok=True)
@ -74,17 +125,15 @@ def run_single_illusion(
json.dump({k: str(v) if not isinstance(v, (int, float, list, bool)) json.dump({k: str(v) if not isinstance(v, (int, float, list, bool))
else v for k, v in cfg.items()}, f, indent=2) else v for k, v in cfg.items()}, f, indent=2)
# Load legacy CFD configs with overridden viscosity and velocity # Load legacy CFD configs with overridden viscosity
cuda_cfg, field_cfg = load_legacy_configs(LEGACY_CFG_DIR) cuda_cfg, field_cfg = load_legacy_configs(LEGACY_CFG_DIR)
field_cfg = field_cfg._replace(viscosity=float(nu)) field_cfg = field_cfg._replace(viscosity=float(nu))
if u0 != 0.01:
field_cfg = field_cfg._replace(velocity=float(u0))
# -- Phase 1: Target recording (target cylinder + 3 sensors) ------------ # -- Phase 1: Target recording (target cylinder + 3 sensors) ------------
ff = FlowField(field_cfg, cuda_cfg, device_id=device_id) ff = FlowField(field_cfg, cuda_cfg, device_id=device_id)
ny = ff.FIELD_SHAPE[1] ny = ff.FIELD_SHAPE[1]
# Add target cylinder at x=20*L0, with radius = target_diam * L0 # Add target cylinder at x=20*L0
print(f" Adding target cylinder: diam={target_diam}L, pos=({20*l0:.0f}, {ny/2:.0f})") print(f" Adding target cylinder: diam={target_diam}L, pos=({20*l0:.0f}, {ny/2:.0f})")
ff.add_cylinder((20.0 * l0, (ny - 1) / 2, 0.0), target_diam * l0) ff.add_cylinder((20.0 * l0, (ny - 1) / 2, 0.0), target_diam * l0)
@ -101,16 +150,30 @@ def run_single_illusion(
print(f" Stabilising ({stabilize_steps} steps)...") print(f" Stabilising ({stabilize_steps} steps)...")
ff.run(stabilize_steps, np.zeros(n_obj_phase1, dtype=DATA_TYPE)) ff.run(stabilize_steps, np.zeros(n_obj_phase1, dtype=DATA_TYPE))
# Record target # Record target: 8 channels = [cylinder_fx, cylinder_fy, sensor0_ux,uy, sensor1_ux,uy, sensor2_ux,uy]
fifo_len = FIFO_LEN # 150
target_states = np.empty((0, 8), dtype=DATA_TYPE) target_states = np.empty((0, 8), dtype=DATA_TYPE)
for _ in range(FIFO_LEN): for _ in range(fifo_len):
ff.run(sample_interval, np.zeros(n_obj_phase1, dtype=DATA_TYPE)) ff.run(sample_interval, np.zeros(n_obj_phase1, dtype=DATA_TYPE))
new_state = ff.obs.copy()[0:8] # sensor[6] + cylinder force[2] new_state = ff.obs.copy()[0:8] # cylinder[2] + sensor[6] = 8 channels
target_states = np.vstack((target_states, new_state)) target_states = np.vstack((target_states, new_state))
print(f" Target recorded: {target_states.shape}") print(f" Target recorded: {target_states.shape}")
# Save target # Analyze harmonics of force channels (channel 0,1 = force; 2..7 = sensors)
np.savez(os.path.join(output_root, "target.npz"), target_states=target_states) # Legacy env uses target_states[:, 2:8] (sensors only) for DTW comparison,
# and target_states[:, 0:2] (cylinder force) for target_cd/target_cl via harmonics.
target_sensors = target_states[:, 2:8].copy() # (150, 6) for similarity
target_forces = target_states[:, 0:2].copy() # (150, 2) for harmonics
target_harmonics = analyze_harmonics(target_forces, n_harmonics=5)
print(f" Target harmonics computed: {len(target_harmonics)} channels")
print(f" DC: {target_harmonics[0]['dc']:.6f}, {target_harmonics[1]['dc']:.6f}")
# Save target data
np.savez(os.path.join(output_root, "target.npz"),
target_states=target_states,
target_sensors=target_sensors)
with open(os.path.join(output_root, "target_harmonics.json"), "w") as f:
json.dump(target_harmonics, f, indent=2)
# Clean up and create pinball env # Clean up and create pinball env
del ff del ff
@ -123,7 +186,6 @@ def run_single_illusion(
ff.add_sensor(sc, l0 / 4.0) ff.add_sensor(sc, l0 / 4.0)
# Add 3 pinball cylinders (illusion positions) # Add 3 pinball cylinders (illusion positions)
# Front at x=front_x*L0, rear at x=rear_x*L0
ff.add_cylinder((front_x * l0, (ny - 1) / 2, 0.0), l0 / 2.0) ff.add_cylinder((front_x * l0, (ny - 1) / 2, 0.0), l0 / 2.0)
ff.add_cylinder((rear_x * l0, (ny - 1) / 2 + 0.75 * l0, 0.0), l0 / 2.0) ff.add_cylinder((rear_x * l0, (ny - 1) / 2 + 0.75 * l0, 0.0), l0 / 2.0)
ff.add_cylinder((rear_x * l0, (ny - 1) / 2 - 0.75 * l0, 0.0), l0 / 2.0) ff.add_cylinder((rear_x * l0, (ny - 1) / 2 - 0.75 * l0, 0.0), l0 / 2.0)
@ -140,9 +202,9 @@ def run_single_illusion(
ff.get_ddf() ff.get_ddf()
ff.save_ddf() ff.save_ddf()
# Norm collection (zero action) # Norm collection (zero action) — matches legacy_env_imit.py
fifo = deque(maxlen=FIFO_LEN) fifo = deque(maxlen=fifo_len)
for _ in range(FIFO_LEN): for _ in range(fifo_len):
ff.run(sample_interval, np.zeros(n_obj, dtype=DATA_TYPE)) ff.run(sample_interval, np.zeros(n_obj, dtype=DATA_TYPE))
fifo.append(ff.obs.copy()[0:12]) fifo.append(ff.obs.copy()[0:12])
@ -161,21 +223,28 @@ def run_single_illusion(
} }
print(f" norm: force_norm_fact={force_norm_fact:.6f}") print(f" norm: force_norm_fact={force_norm_fact:.6f}")
# Bias-action rollout # Bias-action rollout — match legacy_env_imit.py EXACTLY:
# legacy line 142-143: np.array([0.0, 0.0, 0.0, 0.0, -1*U0, 1*U0])
# This is NOT action_bias (which is [0, -2, 2]*U0).
# It's a specific weaker bias: [0, -U0, U0] = [0, -0.01, 0.01]
ff.apply_ddf() ff.apply_ddf()
bias_arr = np.zeros(n_obj, dtype=DATA_TYPE) bias_arr = np.zeros(n_obj, dtype=DATA_TYPE)
bias_arr[3] = float(action_bias[0] * u0) bias_arr[3] = 0.0 # front = 0
bias_arr[4] = float(action_bias[1] * u0) bias_arr[4] = -1.0 * u0 # bottom = -0.01 (cf. action_bias bottom = -2*u0 = -0.02)
bias_arr[5] = float(action_bias[2] * u0) bias_arr[5] = 1.0 * u0 # top = 0.01 (cf. action_bias top = 2*u0 = 0.02)
print(f" bias action: {bias_arr}") print(f" bias action: {bias_arr}")
fifo.clear() fifo.clear()
for _ in range(FIFO_LEN): for _ in range(fifo_len):
ff.run(sample_interval, bias_arr) ff.run(sample_interval, bias_arr)
fifo.append(ff.obs.copy()[0:12]) fifo.append(ff.obs.copy()[0:12])
save_states = np.array(list(fifo), dtype=DATA_TYPE) save_states = np.array(list(fifo), dtype=DATA_TYPE)
norm["save_states"] = save_states norm["save_states"] = save_states
ff.apply_ddf()
# Save DDF AFTER bias rollout (so reset restores post-bias state)
ff.get_ddf()
ff.save_ddf()
ff.apply_ddf() # restore to checkpoint for inference start
# Save norm # Save norm
norm_json = {k: v for k, v in norm.items() if not isinstance(v, np.ndarray)} norm_json = {k: v for k, v in norm.items() if not isinstance(v, np.ndarray)}
@ -183,11 +252,10 @@ def run_single_illusion(
json.dump(norm_json, f, indent=2) json.dump(norm_json, f, indent=2)
# -- Phase 3: Controlled inference --------------------------------------- # -- Phase 3: Controlled inference ---------------------------------------
result = {"scene": scene_name, "controlled": False} result = {"scene": scene_name, "controlled": False, "similarity": 0.0}
model_path = model_path_for_scene(scene_name) model_path = model_path_for_scene(scene_name)
if model_path is not None: if model_path is not None:
s_dim = cfg.get("s_dim", 12)
print(f" loading model: {model_path} (s_dim={s_dim})") print(f" loading model: {model_path} (s_dim={s_dim})")
model = load_ppo_model(model_path, device=f"cuda:{device_id}", s_dim=s_dim) model = load_ppo_model(model_path, device=f"cuda:{device_id}", s_dim=s_dim)
model.set_random_seed(0) model.set_random_seed(0)
@ -196,15 +264,16 @@ def run_single_illusion(
ff.restore_ddf() ff.restore_ddf()
ff.apply_ddf() ff.apply_ddf()
# Re-bias FIFO # Re-bias FIFO (using same bias action as Phase 2)
fifo = deque(maxlen=FIFO_LEN) fifo = deque(maxlen=fifo_len)
for _ in range(FIFO_LEN): for _ in range(fifo_len):
ff.context.push() ff.context.push()
ff.run(sample_interval, bias_arr) ff.run(sample_interval, bias_arr)
ff.context.pop() ff.context.pop()
fifo.append(ff.obs.copy()[0:12]) fifo.append(ff.obs.copy()[0:12])
sens_list, forc_list, action_list = [], [], [] sens_list, forc_list, action_list = [], [], []
target_forces_list = []
obs = np.zeros(s_dim, dtype=np.float32) obs = np.zeros(s_dim, dtype=np.float32)
for step in range(n_infer_steps): for step in range(n_infer_steps):
@ -227,28 +296,34 @@ def run_single_illusion(
sens_list.append(obs_slice[0:6]) sens_list.append(obs_slice[0:6])
forc_list.append(obs_slice[6:12]) forc_list.append(obs_slice[6:12])
# Build normalized obs (just forces + sens for S_DIM=12) # Build normalized obs
forces_norm = obs_slice[6:12] / force_norm_fact forces_norm = obs_slice[6:12] / force_norm_fact
sens_norm = (obs_slice[0:6] - sens_deviation) / sens_norm_fact sens_norm = (obs_slice[0:6] - sens_deviation) / sens_norm_fact
obs12 = np.clip(np.hstack([forces_norm, sens_norm]), -1.0, 1.0).astype(np.float32) obs12 = np.clip(np.hstack([forces_norm, sens_norm]), -1.0, 1.0).astype(np.float32)
if s_dim == 14: if s_dim == 14:
# Need target values -- for inference we zero-pad # Reconstruct target_cd, target_cl from harmonics (matching legacy)
target_vals = gen_target_states_at(step, target_harmonics)
target_cd = target_vals[0] / force_norm_fact
target_cl = target_vals[1] / force_norm_fact
target_forces_list.append(target_vals.copy()) # raw lattice-unit target forces
obs = np.zeros(14, dtype=np.float32) obs = np.zeros(14, dtype=np.float32)
obs[:12] = obs12 obs[:12] = obs12
obs[12] = np.clip(target_cd, -1.0, 1.0)
obs[13] = np.clip(target_cl, -1.0, 1.0)
else: else:
obs = obs12 obs = obs12
np.savez(os.path.join(output_root, "controlled.npz"), np.savez(os.path.join(output_root, "controlled.npz"),
sensors=np.array(sens_list, dtype=np.float32), sensors=np.array(sens_list, dtype=np.float32),
forces=np.array(forc_list, dtype=np.float32), forces=np.array(forc_list, dtype=np.float32),
actions=np.array(action_list, dtype=np.float32)) actions=np.array(action_list, dtype=np.float32),
target_forces=np.array(target_forces_list, dtype=np.float32))
# Compute similarity (use the target cylinder's sensor-only signals) # Compute similarity: compare controlled sensors vs target sensors
# For comparison, compute similarity between controlled sensors and target # target_sensors = target_states[:, 2:8] (sensor only, 6 channels)
target_sensors = target_states[:, 0:6] sens_arr = np.array(sens_list, dtype=np.float32)
sim = compute_similarity(target_sensors, sim = compute_similarity(target_sensors, sens_arr, conv_len)
np.array(sens_list, dtype=np.float32), CONV_LEN)
print(f" similarity (vs target cylinder) = {sim:.4f}") print(f" similarity (vs target cylinder) = {sim:.4f}")
result["controlled"] = True result["controlled"] = True
@ -265,7 +340,7 @@ def run_single_illusion(
def main(): def main():
ap = argparse.ArgumentParser(description="Illusion inference") ap = argparse.ArgumentParser(description="Illusion inference (corrected)")
ap.add_argument("--diameter", type=str, default="1.0", ap.add_argument("--diameter", type=str, default="1.0",
help='Diameter: 0.75, 1.0, 1.5, or "all"') help='Diameter: 0.75, 1.0, 1.5, or "all"')
ap.add_argument("--device", type=int, default=0, help="GPU device ID") ap.add_argument("--device", type=int, default=0, help="GPU device ID")
@ -278,7 +353,6 @@ def main():
scene_names = get_scene_list("illusion") scene_names = get_scene_list("illusion")
else: else:
d = float(args.diameter) d = float(args.diameter)
# Match by target_diameter field
scene_names = [] scene_names = []
for sn in get_scene_list("illusion"): for sn in get_scene_list("illusion"):
cfg = get_scene(sn) cfg = get_scene(sn)

View File

@ -0,0 +1,766 @@
#!/usr/bin/env python3
"""Unified SINDy fitting v2 for all scenes.
Uses V2 feature builder: no sin_ua/cos_ua, optional mu, weighted STLSQ.
Generates separate per-scene results and cross-scene comparisons.
Usage:
conda run -n pycuda_3_10 python sindy/run_all_v2.py \\
--scenes karman_re50,karman_re100,steady
# Karman cross-Re only
conda run -n pycuda_3_10 python sindy/run_all_v2.py --family karman
# Cloak family (Karman + steady + vortex)
conda run -n pycuda_3_10 python sindy/run_all_v2.py --family cloak
# All scenes
conda run -n pycuda_3_10 python sindy/run_all_v2.py --family all
"""
from __future__ import annotations
import argparse
import json
import os
import sys
from typing import List, Optional
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
if _REPO not in sys.path:
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from SR_analysis.utils.sindy_fitter import (
fit_sindy_weighted, get_feature_matrix_v2, get_active_support,
get_feature_matrix_deriv, compute_action_deriv,
)
from SR_analysis.utils.feature_builder import CORE_FEAT_KEYS_V2, ALL_FEAT_KEYS_V2, PHYSICS_FEAT_KEYS
from SR_analysis.configs import get_scene, get_scene_list, SCENES
# Base output directory
SINDY_DIR = os.path.join(os.path.dirname(__file__))
# Threshold grid (more granular at low end for better Pareto selection)
THRESHOLDS = [0.0, 0.001, 0.002, 0.003, 0.005, 0.007, 0.01, 0.015, 0.02, 0.03, 0.05, 0.08, 0.1]
FAMILIES = {
"cloak": ["karman", "steady", "vortex"],
"karman": ["karman"],
"illusion": ["illusion"],
"vortex": ["vortex"],
"all": None, # all scenes
}
def load_data(scene_name: str) -> tuple:
"""Load sensors/forces/actions from scene's controlled.npz.
Returns actions in PHYSICAL omega (lattice units).
Also returns target_forces if available (for Illusion scenes).
"""
cfg = get_scene(scene_name)
data_dir = os.path.join(
os.path.dirname(__file__), "..", "data",
cfg["scene_id"], scene_name,
)
npz = np.load(os.path.join(data_dir, "controlled.npz"))
sensors = npz["sensors"].astype(np.float64)
forces = npz["forces"].astype(np.float64)
# actions in .npz are normalized [-1,1]; convert to physical omega
actions_norm = npz["actions"].astype(np.float64)
scale = cfg["action_scale"]
bias = np.array(cfg["action_bias"], dtype=np.float64)
u0 = cfg["u0"]
actions_phys = (actions_norm * scale + bias) * u0
# Load target_forces if available (Illusion scenes)
target_forces = None
if "target_forces" in npz:
target_forces = npz["target_forces"].astype(np.float64)
return sensors, forces, actions_phys, cfg, target_forces
def compute_scene_weight(scene_name: str, cfg: dict) -> float:
"""Compute scene quality weight from PPO similarity."""
result_path = os.path.join(
os.path.dirname(__file__), "..", "data",
cfg["scene_id"], scene_name, "result.json",
)
if os.path.isfile(result_path):
with open(result_path) as f:
r = json.load(f)
sim = r.get("similarity", r.get("avg_reward_last100", 0.5))
return float(sim) ** 2
return 1.0
def run_single_scene(
scene_name: str,
thresholds: Optional[List[float]] = None,
verbose: bool = True,
) -> dict:
"""Run SINDy v2 on a single scene. Returns result dict."""
if thresholds is None:
thresholds = THRESHOLDS
sensors, forces, actions_phys, cfg, target_forces = load_data(scene_name)
mu = cfg["mu"]
u0 = cfg["u0"]
scene_id = cfg["scene_id"]
T = sensors.shape[0]
if verbose:
print(f"\n{'='*60}")
print(f"Scene: {scene_name} ({scene_id}) T={T} mu={mu:.4f} u0={u0}")
print(f"{'='*60}")
# Determine if this is a single-scene (no mu) or cross-Re (with mu)
# For cross-Re, we'll handle separately; for single-scene, no mu
# But steady is special: it needs mu=constant, still no mu in features
# Single scene: no mu, no sin/cos
is_single_re = True # will be overridden for joint fitting
# Build feature matrix (V2: no sin/cos, no mu for single-scene)
# Use dt_c=sample_interval/T0 in T0 units
t0_steps = 2000 # T0 = D/U0 = 20/0.01 = 2000 LBM steps
dt_c = cfg.get("sample_interval", 800) / t0_steps
Theta_f, Theta_r, Y, fn_f, fn_r = get_feature_matrix_v2(
sensors, forces, actions_phys, mu,
u0=u0, include_mu=False, include_cos_sin=False,
use_time_norm=False, dt_c=dt_c, n_warmup=2,
target_forces=target_forces,
)
if verbose:
print(f" Features: front {Theta_f.shape[1]-1} + bias?0, rear {Theta_r.shape[1]-1} + bias?1")
print(f" dt_c={dt_c:.4f} (T0 units)")
# Scene quality weight for STLSQ
w_scene = compute_scene_weight(scene_name, cfg)
if verbose:
print(f" Scene quality weight: {w_scene:.4f}")
# ---------- Front channel (no bias) ----------
if verbose:
print(f"\n --- Front (no bias) ---")
front_results = fit_sindy_weighted(Theta_f, Y[:, 0], thresholds,
n_robust_passes=2)
# Find best by R2 (include th=0.0 for dense coefficients)
front_best = max(front_results, key=lambda r: r["r2"])
if verbose:
_print_active(fn_f, front_best["coef"])
# ---------- Top channel (rear shared-head, with bias) ----------
if verbose:
print(f"\n --- Top (rear shared-head, with bias) ---")
top_results = fit_sindy_weighted(Theta_r, Y[:, 2], thresholds,
n_robust_passes=2)
top_best = max(top_results, key=lambda r: r["r2"])
if verbose:
_print_active(fn_r, top_best["coef"])
# ---------- Bottom (independent, for comparison) ----------
if verbose:
print(f"\n --- Bottom (independent, with bias) ---")
bot_results = fit_sindy_weighted(Theta_r, Y[:, 1], thresholds,
n_robust_passes=2)
bot_best = max(bot_results, key=lambda r: r["r2"])
if verbose:
_print_active(fn_r, bot_best["coef"])
# Compute additional metrics for the best thresholds
n_warmup = 2
n_samples = Theta_f.shape[0]
result = {
"scene": scene_name,
"scene_id": scene_id,
"re_code": cfg.get("re_code", 100),
"mu": mu,
"u0": u0,
"dt_c": dt_c,
"n_samples": n_samples,
"thresholds": thresholds,
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"front": {
"results": _clean_results(front_results),
"best": _clean_single(front_best),
"best_coef": front_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in front_results],
},
"top": {
"results": _clean_results(top_results),
"best": _clean_single(top_best),
"best_coef": top_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in top_results],
},
"bottom": {
"results": _clean_results(bot_results),
"best": _clean_single(bot_best),
"best_coef": bot_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in bot_results],
},
}
if verbose:
print(f"\n Summary:")
print(f" Front: th={front_best['threshold']:.4f} nz={front_best['nz']:2d} R2={front_best['r2']:.4f}")
print(f" Top: th={top_best['threshold']:.4f} nz={top_best['nz']:2d} R2={top_best['r2']:.4f}")
print(f" Bottom: th={bot_best['threshold']:.4f} nz={bot_best['nz']:2d} R2={bot_best['r2']:.4f}")
return result
def _print_active(names: List[str], coef: list, rtol: float = 1e-4):
"""Print active features from coefficient vector."""
for i, c in enumerate(coef):
if abs(c) > rtol and i < len(names):
print(f" {names[i]:20s} = {c:+.6f}")
def _clean_results(results: List[dict]) -> List[dict]:
"""Keep all keys including coef."""
return [dict(r) for r in results]
def _clean_single(entry: dict) -> dict:
"""Remove coef from single entry."""
return {k: v for k, v in entry.items() if k != "coef"}
def run_joint_karman(
scene_names: List[str],
thresholds: Optional[List[float]] = None,
) -> dict:
"""Joint Karman cross-Re SINDy with mu terms and scene weighting."""
if thresholds is None:
thresholds = THRESHOLDS
print(f"\n{'='*60}")
print(f"Joint Karman cross-Re: {scene_names}")
print(f"{'='*60}")
# Load and concatenate all Karman scenes
all_sensors, all_forces, all_actions = [], [], []
all_weights = []
t0_steps = 2000
for sn in scene_names:
sensors, forces, actions, cfg, _ = load_data(sn)
mu = cfg["mu"]
u0 = cfg["u0"]
dt_c = cfg.get("sample_interval", 800) / t0_steps
# Build features WITH mu terms
Theta_f, Theta_r, Y, fn_f, fn_r = get_feature_matrix_v2(
sensors, forces, actions, mu,
u0=u0, include_mu=True, include_cos_sin=False,
use_time_norm=False, dt_c=dt_c, n_warmup=2,
)
w = compute_scene_weight(sn, cfg)
n = Theta_f.shape[0]
print(f" {sn}: {n} samples, weight={w:.4f}")
all_weights.append(np.full(n, w))
all_sensors.append(Theta_f)
all_forces.append(Theta_r) # reuse variable names for concat
all_actions.append(Y)
# Concatenate
Theta_f_j = np.vstack(all_sensors)
Theta_r_j = np.vstack(all_forces)
Y_j = np.vstack(all_actions)
W_j = np.concatenate(all_weights)
print(f" Joint: front {Theta_f_j.shape}, rear {Theta_r_j.shape}")
print(f" Weights: min={W_j.min():.4f} max={W_j.max():.4f}")
# ---------- Front (no bias, all features including mu) ----------
print(f"\n --- Joint Front (no bias, with mu) ---")
front_results = fit_sindy_weighted(Theta_f_j, Y_j[:, 0], thresholds,
sample_weights=W_j, n_robust_passes=2)
front_best = max(front_results, key=lambda r: r["r2"])
_print_active(fn_f, front_best["coef"])
# ---------- Top (with bias, all features including mu) ----------
print(f"\n --- Joint Top (with bias, with mu) ---")
top_results = fit_sindy_weighted(Theta_r_j, Y_j[:, 2], thresholds,
sample_weights=W_j, n_robust_passes=2)
top_best = max(top_results, key=lambda r: r["r2"])
_print_active(fn_r, top_best["coef"])
# ---------- Bottom (independent) ----------
print(f"\n --- Joint Bottom (independent) ---")
bot_results = fit_sindy_weighted(Theta_r_j, Y_j[:, 1], thresholds,
sample_weights=W_j, n_robust_passes=2)
bot_best = max(bot_results, key=lambda r: r["r2"])
_print_active(fn_r, bot_best["coef"])
return {
"joint_scenes": scene_names,
"thresholds": thresholds,
"n_samples_total": Theta_f_j.shape[0],
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"front": {
"results": _clean_results(front_results),
"best": _clean_single(front_best),
"best_coef": front_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in front_results],
},
"top": {
"results": _clean_results(top_results),
"best": _clean_single(top_best),
"best_coef": top_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in top_results],
},
"bottom": {
"results": _clean_results(bot_results),
"best": _clean_single(bot_best),
"best_coef": bot_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in bot_results],
},
}
def run_joint_illusion(
scene_names: List[str],
thresholds: Optional[List[float]] = None,
) -> dict:
"""Joint Illusion cross-diameter SINDy with target features and scene weighting."""
if thresholds is None:
thresholds = THRESHOLDS
print(f"\n{'='*60}")
print(f"Joint Illusion cross-diameter: {scene_names}")
print(f"{'='*60}")
all_sensors, all_forces, all_actions = [], [], []
all_weights = []
t0_steps = 2000
for sn in scene_names:
sensors, forces, actions, cfg, target_forces = load_data(sn)
mu = cfg["mu"]
u0 = cfg["u0"]
dt_c = cfg.get("sample_interval", 600) / t0_steps
# Build features WITH target features, no mu
Theta_f, Theta_r, Y, fn_f, fn_r = get_feature_matrix_v2(
sensors, forces, actions, mu,
u0=u0, include_mu=False, include_cos_sin=False,
use_time_norm=False, dt_c=dt_c, n_warmup=2,
target_forces=target_forces,
)
w = compute_scene_weight(sn, cfg)
n = Theta_f.shape[0]
print(f" {sn}: {n} samples, dt_c={dt_c:.4f}, weight={w:.4f}")
all_weights.append(np.full(n, w))
all_sensors.append(Theta_f)
all_forces.append(Theta_r)
all_actions.append(Y)
Theta_f_j = np.vstack(all_sensors)
Theta_r_j = np.vstack(all_forces)
Y_j = np.vstack(all_actions)
W_j = np.concatenate(all_weights)
print(f" Joint: front {Theta_f_j.shape}, rear {Theta_r_j.shape}")
print(f" Weights: min={W_j.min():.4f} max={W_j.max():.4f}")
# Front (no bias)
print(f"\n --- Joint Front (no bias) ---")
front_results = fit_sindy_weighted(Theta_f_j, Y_j[:, 0], thresholds,
sample_weights=W_j, n_robust_passes=2)
front_best = max(front_results, key=lambda r: r["r2"])
_print_active(fn_f, front_best["coef"])
# Top (with bias)
print(f"\n --- Joint Top (with bias) ---")
top_results = fit_sindy_weighted(Theta_r_j, Y_j[:, 2], thresholds,
sample_weights=W_j, n_robust_passes=2)
top_best = max(top_results, key=lambda r: r["r2"])
_print_active(fn_r, top_best["coef"])
# Bottom (independent)
print(f"\n --- Joint Bottom (independent) ---")
bot_results = fit_sindy_weighted(Theta_r_j, Y_j[:, 1], thresholds,
sample_weights=W_j, n_robust_passes=2)
bot_best = max(bot_results, key=lambda r: r["r2"])
_print_active(fn_r, bot_best["coef"])
return {
"joint_scenes": scene_names,
"thresholds": thresholds,
"n_samples_total": Theta_f_j.shape[0],
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"front": {
"results": _clean_results(front_results),
"best": _clean_single(front_best),
"best_coef": front_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in front_results],
},
"top": {
"results": _clean_results(top_results),
"best": _clean_single(top_best),
"best_coef": top_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in top_results],
},
"bottom": {
"results": _clean_results(bot_results),
"best": _clean_single(bot_best),
"best_coef": bot_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in bot_results],
},
}
def run_single_scene_deriv(
scene_name: str,
thresholds: Optional[List[float]] = None,
verbose: bool = True,
center_diff: bool = False,
augment_level: int = 0,
feat_keys: Optional[List[str]] = None,
output_mode: str = "deriv",
) -> dict:
"""Run SINDy in phase-state mode.
output_mode:
"deriv": Y = d(alpha)/dt (time-normalized). Closed-loop needs integration.
"absolute": Y = alpha (non-dimensional absolute action). No integration.
"""
if thresholds is None:
thresholds = THRESHOLDS
sensors, forces, actions_phys, cfg, target_forces = load_data(scene_name)
mu = cfg["mu"]
u0 = cfg["u0"]
scene_id = cfg["scene_id"]
T = sensors.shape[0]
t0_steps = 2000
dt_c = cfg.get("sample_interval", 800) / t0_steps
if verbose:
print(f"\n{'='*60}")
print(f"{'ABSOLUTE' if output_mode == 'absolute' else 'DERIV'} Scene: {scene_name} dt_c={dt_c:.4f} u0={u0}")
print(f"{'='*60}")
Theta_f, Theta_r, Y, fn_f, fn_r = get_feature_matrix_deriv(
sensors, forces, actions_phys, mu,
u0=u0, dt_c=dt_c,
target_forces=target_forces,
center_diff=center_diff,
augment_level=augment_level,
feat_keys=feat_keys,
output_mode=output_mode,
include_mu=feat_keys is not None and any("mu_" in k for k in feat_keys),
)
y_label = "alpha (abs)" if output_mode == "absolute" else "d(alpha)/dt"
if verbose:
print(f" Features: front {Theta_f.shape[1]} (no bias), rear {Theta_r.shape[1]} (with bias)")
print(f" Fitting {y_label}: {Y.shape}")
w_scene = compute_scene_weight(scene_name, cfg)
# Front (no bias)
if verbose:
print(f"\n --- Front (no bias, target = {y_label}_F) ---")
front_results = fit_sindy_weighted(Theta_f, Y[:, 0], thresholds, n_robust_passes=2)
front_best = max(front_results, key=lambda r: r["r2"])
if verbose:
_print_active(fn_f, front_best["coef"])
print(f" R2={front_best['r2']:.4f} nz={front_best['nz']}")
# Top (with bias, rear shared-head)
if verbose:
print(f"\n --- Top (with bias, target = {y_label}_T) ---")
top_results = fit_sindy_weighted(Theta_r, Y[:, 2], thresholds, n_robust_passes=2)
top_best = max(top_results, key=lambda r: r["r2"])
if verbose:
_print_active(fn_r, top_best["coef"])
print(f" R2={top_best['r2']:.4f} nz={top_best['nz']}")
# Bottom (independent)
if verbose:
print(f"\n --- Bottom (independent, target = {y_label}_B) ---")
bot_results = fit_sindy_weighted(Theta_r, Y[:, 1], thresholds, n_robust_passes=2)
bot_best = max(bot_results, key=lambda r: r["r2"])
if verbose:
_print_active(fn_r, bot_best["coef"])
print(f" R2={bot_best['r2']:.4f} nz={bot_best['nz']}")
result = {
"scene": scene_name,
"scene_id": scene_id,
"mode": output_mode,
"dt_c": dt_c,
"center_diff": center_diff,
"re_code": cfg.get("re_code", 100),
"mu": mu,
"u0": u0,
"n_samples": Theta_f.shape[0],
"thresholds": thresholds,
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"front": {
"results": _clean_results(front_results),
"best": _clean_single(front_best),
"best_coef": front_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in front_results],
},
"top": {
"results": _clean_results(top_results),
"best": _clean_single(top_best),
"best_coef": top_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in top_results],
},
"bottom": {
"results": _clean_results(bot_results),
"best": _clean_single(bot_best),
"best_coef": bot_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in bot_results],
},
}
return result
def run_joint_karman_deriv(
scene_names: List[str],
thresholds: Optional[List[float]] = None,
augment_level: int = 0,
) -> dict:
"""Joint Karman cross-Re derivative-mode SINDy."""
if thresholds is None:
thresholds = THRESHOLDS
print(f"\n{'='*60}")
print(f"DERIV Joint Karman: {scene_names}")
print(f"{'='*60}")
t0_steps = 2000
all_Theta_f, all_Theta_r, all_Y = [], [], []
all_weights = []
for sn in scene_names:
sensors, forces, actions, cfg, _ = load_data(sn)
mu = cfg["mu"]
u0 = cfg["u0"]
dt_c = cfg.get("sample_interval", 800) / t0_steps
Theta_f, Theta_r, Y, fn_f, fn_r = get_feature_matrix_deriv(
sensors, forces, actions, mu,
u0=u0, dt_c=dt_c, include_mu=True,
target_forces=None, augment_level=augment_level,
)
w = compute_scene_weight(sn, cfg)
n = Theta_f.shape[0]
print(f" {sn}: {n} samples, dt_c={dt_c:.4f}, weight={w:.4f}")
all_weights.append(np.full(n, w))
all_Theta_f.append(Theta_f)
all_Theta_r.append(Theta_r)
all_Y.append(Y)
Theta_f_j = np.vstack(all_Theta_f)
Theta_r_j = np.vstack(all_Theta_r)
Y_j = np.vstack(all_Y)
W_j = np.concatenate(all_weights)
print(f" Joint: front {Theta_f_j.shape}, rear {Theta_r_j.shape}")
# Front
print(f"\n --- Joint Front ---")
front_results = fit_sindy_weighted(Theta_f_j, Y_j[:, 0], thresholds, sample_weights=W_j, n_robust_passes=2)
front_best = max(front_results, key=lambda r: r["r2"])
_print_active(fn_f, front_best["coef"])
print(f" R2={front_best['r2']:.4f}")
# Top
print(f"\n --- Joint Top ---")
top_results = fit_sindy_weighted(Theta_r_j, Y_j[:, 2], thresholds, sample_weights=W_j, n_robust_passes=2)
top_best = max(top_results, key=lambda r: r["r2"])
_print_active(fn_r, top_best["coef"])
print(f" R2={top_best['r2']:.4f}")
# Bottom
print(f"\n --- Joint Bottom ---")
bot_results = fit_sindy_weighted(Theta_r_j, Y_j[:, 1], thresholds, sample_weights=W_j, n_robust_passes=2)
bot_best = max(bot_results, key=lambda r: r["r2"])
_print_active(fn_r, bot_best["coef"])
print(f" R2={bot_best['r2']:.4f}")
return {
"joint_scenes": scene_names,
"mode": "deriv",
"thresholds": thresholds,
"n_samples_total": Theta_f_j.shape[0],
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"front": {
"results": _clean_results(front_results),
"best": _clean_single(front_best),
"best_coef": front_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in front_results],
},
"top": {
"results": _clean_results(top_results),
"best": _clean_single(top_best),
"best_coef": top_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in top_results],
},
"bottom": {
"results": _clean_results(bot_results),
"best": _clean_single(bot_best),
"best_coef": bot_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in bot_results],
},
}
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--scenes", type=str, default=None,
help="Comma-separated scene names")
ap.add_argument("--family", type=str, default=None,
choices=list(FAMILIES.keys()) + [None],
help="Scene family")
ap.add_argument("--out", type=str, default=None,
help="Output directory (default: sindy/<scene_id>/")
ap.add_argument("--joint", action="store_true",
help="Run joint Karman cross-Re fitting")
ap.add_argument("--deriv", action="store_true",
help="Use derivative mode: fit d(alpha)/dt = g(physics_state)")
ap.add_argument("--center-diff", action="store_true",
help="Use centered difference for derivative (vs forward diff)")
ap.add_argument("--augment-level", type=int, default=0,
help="Observation augmentation: 0=static, 1=+lags, 2=+derivs, 3=both, 4=+action_lag")
ap.add_argument("--phase", action="store_true",
help="Use phase-state features: [u_a, du_a/dt, Cl_tot, dCl_tot/dt, Cd_tot, Cd_rear]")
ap.add_argument("--karman-expand", action="store_true",
help="Karman expanded: phase-state + u_m/u_c/v_a/Cl_diff (10 dim)")
ap.add_argument("--karman-mu", action="store_true",
help="Karman phase-state + mu modulation: 6-dim + mu*Cl_tot/mu*Cd_tot/mu*u_a (9 dim)")
ap.add_argument("--output-mode", type=str, default="deriv",
choices=["deriv", "absolute"],
help="'deriv': predict d(alpha)/dt; 'absolute': predict alpha directly")
args = ap.parse_args()
# Resolve scene list
if args.scenes:
scene_names = [s.strip() for s in args.scenes.split(",")]
elif args.family and args.family in FAMILIES:
family_ids = FAMILIES[args.family]
if family_ids is None:
scene_names = get_scene_list() # all
else:
scene_names = []
for fid in family_ids:
scene_names.extend(get_scene_list(fid))
else:
scene_names = get_scene_list("karman")
# Determine fitting function based on mode
if args.deriv:
om = args.output_mode
if args.phase:
from SR_analysis.utils.feature_builder import PHASE_STATE_KEYS, ILLUSION_PHASE_KEYS
def _phase_func(sn):
from SR_analysis.configs import get_scene
cfg = get_scene(sn)
feat_keys = ILLUSION_PHASE_KEYS if cfg["scene_id"] == "illusion" else PHASE_STATE_KEYS
return run_single_scene_deriv(sn, center_diff=args.center_diff,
augment_level=args.augment_level,
feat_keys=feat_keys, output_mode=om)
fit_func = _phase_func
elif args.karman_expand:
from SR_analysis.utils.feature_builder import KARMAN_EXPANDED_KEYS
fit_func = lambda sn: run_single_scene_deriv(sn, center_diff=args.center_diff,
augment_level=args.augment_level,
feat_keys=list(KARMAN_EXPANDED_KEYS),
output_mode=om)
elif args.karman_mu:
from SR_analysis.utils.feature_builder import PHASE_STATE_KEYS
# Phase-state + 3 selected mu modulation features
mu_feat = ["mu_Cl_tot", "mu_Cd_tot", "mu_u_a"]
feat_keys = list(PHASE_STATE_KEYS) + mu_feat
fit_func = lambda sn: run_single_scene_deriv(sn, center_diff=args.center_diff,
augment_level=args.augment_level,
feat_keys=feat_keys,
output_mode=om)
else:
fit_func = lambda sn: run_single_scene_deriv(sn, center_diff=args.center_diff,
augment_level=args.augment_level,
output_mode=om)
joint_func_karman = lambda names: run_joint_karman_deriv(names, augment_level=args.augment_level)
else:
fit_func = run_single_scene
joint_func_karman = run_joint_karman
print(f"Running SINDy {'deriv' if args.deriv else 'v2'} for scenes: {scene_names}")
# Run per-scene
for sn in scene_names:
result = fit_func(sn)
# Save per-scene to its scene_id directory
scene_id = result["scene_id"]
out_dir = os.path.join(SINDY_DIR, scene_id)
os.makedirs(out_dir, exist_ok=True)
suffix = "_deriv" if args.deriv else "_v2"
out_path = os.path.join(out_dir, f"sindy_results{suffix}.json")
# Load existing if any, update per_scene
if os.path.isfile(out_path):
with open(out_path) as f:
existing = json.load(f)
else:
existing = {
"thresholds": THRESHOLDS,
"all_feature_names_front": result["feature_names_front"],
"all_feature_names_rear": result["feature_names_rear"],
"per_scene": {},
}
existing["per_scene"][sn] = result
with open(out_path, "w") as f:
json.dump(existing, f, indent=2)
print(f" Saved: {out_path}")
# Joint Karman cross-Re (training Re only: 50, 100, 200, 400)
if args.joint:
karman_train = ["karman_re50", "karman_re100", "karman_re200", "karman_re400"]
joint_result = joint_func_karman(karman_train)
joint_suffix = "_joint_deriv" if args.deriv else "_joint_v2"
out_path = os.path.join(SINDY_DIR, "karman", f"sindy{joint_suffix}.json")
with open(out_path, "w") as f:
json.dump(joint_result, f, indent=2)
print(f"\nSaved joint: {out_path}")
# Joint Illusion (only in non-deriv mode for now)
if args.joint and not args.deriv:
illusion_names = get_scene_list("illusion")
joint_illusion = run_joint_illusion(illusion_names)
out_path = os.path.join(SINDY_DIR, "illusion", "sindy_joint_v2.json")
os.makedirs(os.path.dirname(out_path), exist_ok=True)
with open(out_path, "w") as f:
json.dump(joint_illusion, f, indent=2)
print(f"\nSaved joint illusion: {out_path}")
if __name__ == "__main__":
main()

View File

@ -0,0 +1,197 @@
#!/usr/bin/env python3
"""Restricted PySR on SINDy whitelist features.
Uses the frozen whitelists (whitelist_karman.json, whitelist_illusion.json)
and the controlled CFD data to search for compact closed-form control laws.
Environment: env_sr (conda).
Usage:
conda run -n env_sr python src/SR_analysis/sindy/run_pysr.py --scene karman_re100
conda run -n env_sr python src/SR_analysis/sindy/run_pysr.py --scene illusion_1L
conda run -n env_sr python src/SR_analysis/sindy/run_pysr.py --scene all
"""
from __future__ import annotations
import argparse
import json
import os
import sys
from typing import Dict, List, Optional
import numpy as np
# PySR
from pysr import PySRRegressor
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
sys.path.insert(0, _SRC)
from SR_analysis.utils.feature_builder import (
compute_dimensionless, compute_features, build_feature_matrix,
CORE_FEAT_KEYS_V2,
)
from SR_analysis.configs import get_scene
SINDY_DIR = os.path.join(os.path.dirname(__file__))
def load_controlled_data(scene_name: str):
"""Load controlled data for a scene, returning sensors, forces, actions_phys, target_forces."""
cfg = get_scene(scene_name)
scene_id = cfg["scene_id"]
data_dir = os.path.join(SINDY_DIR, "..", "data", scene_id, scene_name)
npz = np.load(os.path.join(data_dir, "controlled.npz"))
sensors = npz["sensors"].astype(np.float64)
forces = npz["forces"].astype(np.float64)
actions_norm = npz["actions"].astype(np.float64)
actions_phys = (actions_norm * cfg["action_scale"] + np.array(cfg["action_bias"])) * cfg["u0"]
target_forces = npz["target_forces"].astype(np.float64) if "target_forces" in npz else None
return sensors, forces, actions_phys, target_forces, cfg
def run_pysr_scene(scene_name: str, out_dir: str):
"""Run PySR on a single scene."""
# Select whitelist based on scene family
cfg = get_scene(scene_name)
scene_id = cfg["scene_id"]
if scene_id == "karman":
whitelist_path = os.path.join(SINDY_DIR, "whitelist_karman.json")
elif scene_id == "illusion":
whitelist_path = os.path.join(SINDY_DIR, "whitelist_illusion.json")
else:
print(f" SKIP: no whitelist for scene_id={scene_id}")
return
with open(whitelist_path) as f:
wl = json.load(f)
sensors, forces, actions_phys, target_forces, cfg = load_controlled_data(scene_name)
mu = cfg["mu"]
u0 = cfg["u0"]
sample_interval = cfg.get("sample_interval", 800)
t0_steps = 2000 # T0 = D/U0 = 20/0.01 = 2000 LBM steps
dt_c = sample_interval / t0_steps
front_keys = wl["front_active"]
rear_keys = [k for k in wl["rear_active"] if k != "bias"]
# Build features with PROPER lags (matching sindy_fitter.py)
T = sensors.shape[0]
a_prev = np.zeros((T, 3), dtype=np.float64)
a_prev2 = np.zeros((T, 3), dtype=np.float64)
a_prev[1:] = actions_phys[:-1]
a_prev2[2:] = actions_phys[:-2]
dim = compute_dimensionless(sensors, forces, u0=u0, d=20.0)
has_mu = any("mu" in k for k in front_keys)
sym = compute_features(
dim, a_prev, a_prev2, mu,
alpha_mode=False, include_mu=has_mu,
include_cos_sin=False, u0=u0,
target_forces=target_forces,
dt_c=dt_c,
)
n_warmup = 2
X_front = build_feature_matrix(sym, front_keys, add_bias=False)[n_warmup:]
X_rear = build_feature_matrix(sym, rear_keys, add_bias=True)[n_warmup:]
Y = actions_phys[n_warmup:] # target = omega(t) for t >= n_warmup
print(f"\n{'='*60}")
print(f"PySR: {scene_name}")
print(f" Front features ({len(front_keys)}): {front_keys}")
print(f" Rear features ({len(rear_keys)}): {rear_keys}")
print(f" Samples: {X_front.shape[0]}")
print(f"{'='*60}")
# --- Front channel ---
print(f"\n=== PySR: {scene_name} Front ===")
model_f = PySRRegressor(
binary_operators=["+", "-", "*", "/"],
unary_operators=["square"],
niterations=40,
populations=30,
maxsize=15,
complexity_of_constants=2,
parsimony=0.01,
extra_sympy_mappings={},
batching=False,
)
model_f.fit(X_front, Y[:, 0], variable_names=front_keys)
results = {
"scene": scene_name,
"channel": "front",
"feature_names": front_keys,
"equations": model_f.equations_.to_dict(orient="records"),
"best_sympy": str(model_f.sympy()),
"best_score": float(model_f.score(X_front, Y[:, 0])),
}
out_path = os.path.join(out_dir, f"pysr_{scene_name}_front.json")
with open(out_path, "w") as f:
json.dump(results, f, indent=2, default=str)
print(f" Best: {model_f.sympy()}")
print(f" Score: {results['best_score']:.4f}")
print(f" Saved: {out_path}")
# --- Rear/Top channel ---
print(f"\n=== PySR: {scene_name} Top ===")
model_t = PySRRegressor(
binary_operators=["+", "-", "*", "/"],
unary_operators=["square"],
niterations=40,
populations=30,
maxsize=15,
complexity_of_constants=2,
parsimony=0.01,
extra_sympy_mappings={},
batching=False,
)
model_t.fit(X_rear, Y[:, 2], variable_names=["bias"] + rear_keys)
results_t = {
"scene": scene_name,
"channel": "top",
"feature_names": ["bias"] + rear_keys,
"equations": model_t.equations_.to_dict(orient="records"),
"best_sympy": str(model_t.sympy()),
"best_score": float(model_t.score(X_rear, Y[:, 2])),
}
out_path = os.path.join(out_dir, f"pysr_{scene_name}_top.json")
with open(out_path, "w") as f:
json.dump(results_t, f, indent=2, default=str)
print(f" Best: {model_t.sympy()}")
print(f" Score: {results_t['best_score']:.4f}")
print(f" Saved: {out_path}")
def main():
ap = argparse.ArgumentParser(description="Restricted PySR on SINDy whitelists")
ap.add_argument("--scene", type=str, default="all",
help='Scene name like "karman_re100" or "illusion_1L", or "all"')
ap.add_argument("--out", type=str, default=None)
args = ap.parse_args()
if args.out is None:
args.out = os.path.join(SINDY_DIR, "..", "validate", "results")
os.makedirs(args.out, exist_ok=True)
if args.scene == "all":
scenes = ["karman_re100", "illusion_0.75L", "illusion_1L", "illusion_1.5L"]
else:
scenes = [args.scene]
for sn in scenes:
try:
run_pysr_scene(sn, args.out)
except Exception as e:
print(f" ERROR on {sn}: {e}")
import traceback; traceback.print_exc()
print("\nDone.")
if __name__ == "__main__":
main()

View File

@ -0,0 +1,30 @@
{
"scene": "illusion_joint",
"threshold": 0.001,
"selection": "Cleaned features: removed aF_lag1/daF etc, replaced with daF_dt (time-normalized). Features from union across 0.75L/1L/1.5L dense fits with target_Cd/target_Cl.",
"front_active": [
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"daF_dt",
"daB_dt",
"daT_dt",
"target_Cd",
"target_Cl"
],
"rear_active": [
"bias",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"daF_dt",
"daB_dt",
"daT_dt",
"target_Cd",
"target_Cl"
],
"feature_count": 16,
"note": "Cleaned Illusion whitelist: no discrete lags, only time-derivatives + physical features. Separate SINDy closed-loop sims: 0.75L=0.908, 1L=0.962, 1.5L=0.926."
}

View File

@ -0,0 +1,18 @@
{
"scene": "karman_re100",
"threshold": 0.003,
"front_active": [
"daF_dt"
],
"rear_active": [
"bias",
"u_a",
"Cd_rear",
"Cl_diff",
"daF_dt",
"daB_dt",
"daT_dt"
],
"feature_count": 14,
"note": "Cleaned Karman whitelist: removed aF_lag1/daF, replaced with daF_dt (time-normalized). Based on th=0.003 support from joint fit. Separate SINDy closed-loop sim=0.901 (94.4% of PPO)."
}

View File

@ -0,0 +1,66 @@
#!/usr/bin/env python3
"""Wrap joint model into per_scene format.
Usage: python3 src/SR_analysis/sindy/wrap_joint.py [--scene illusion]
"""
import json, sys, os
sys.path.insert(0, os.path.abspath("."))
sys.path.insert(0, os.path.abspath("src"))
import argparse
from SR_analysis.configs import get_scene_list
def wrap_joint(scene_id):
joint_path = f"src/SR_analysis/sindy/{scene_id}/sindy_joint_v2.json"
out_path = f"src/SR_analysis/sindy/{scene_id}/sindy_joint_wrapped.json"
with open(joint_path) as f:
joint = json.load(f)
fn_f = joint["feature_names_front"]
fn_r = joint["feature_names_rear"]
wrapped = {
"thresholds": joint["thresholds"],
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"per_scene": {},
}
# Include all scenes in the family (training + generalization)
all_scene_names = get_scene_list(scene_id)
for scene_name in all_scene_names:
wrapped["per_scene"][scene_name] = {
"scene": scene_name, "re_code": "0", "mu": 0.0,
"feature_names_front": fn_f, "feature_names_rear": fn_r,
"front": {
"results": joint["front"]["results"],
"best": joint["front"]["best"],
"best_coef": joint["front"]["best_coef"],
"sparsity_curve": joint["front"]["sparsity_curve"],
},
"top": {
"results": joint["top"]["results"],
"best": joint["top"]["best"],
"best_coef": joint["top"]["best_coef"],
"sparsity_curve": joint["top"]["sparsity_curve"],
},
"bottom": {
"results": joint["bottom"]["results"],
"best": joint["bottom"]["best"],
"best_coef": joint["bottom"]["best_coef"],
"sparsity_curve": joint["bottom"]["sparsity_curve"],
},
}
with open(out_path, "w") as f:
json.dump(wrapped, f, indent=2)
print(f"Saved {out_path}")
print(f" Scenes ({len(all_scene_names)}): {all_scene_names}")
if __name__ == "__main__":
ap = argparse.ArgumentParser()
ap.add_argument("--scene", type=str, default="karman",
choices=["karman", "illusion"])
args = ap.parse_args()
wrap_joint(args.scene)

View File

@ -0,0 +1,146 @@
# SINDy 与 SR 背景知识
## 文档作用
这份文档只负责一件事:**给正在工作的 coder 提供背景知识、已经确认的经验、已踩过的坑和当前结论强度。**
它不是任务清单,不直接安排"下一步做什么"。凡是执行顺序、阶段划分、最小交付物,统一写在 `sindy_sr_notes`。这份 knowledge 只保留:
- 已确认的技术事实
- 历史错误与纠正
- 结果该如何理解
- 哪些话可以说,哪些话现在还不能说
- 代码和实验上最容易踩的坑
---
## 一、这条线在项目里的位置
SINDy 与 SR 不是独立课题,而是 pinball 后处理主线中的一段工具链。项目真正要解释的是:
\[
\text{obs} \rightarrow \text{act} \rightarrow \text{flow structure} \rightarrow \text{signature}
\]
SINDy 与 SR 当前只直接处理其中的 `obs -> act` 白箱化,但它们的价值在于:
- 检验控制是否真的依赖少数物理量
- 识别不同 cloak 场景中是否复用了同一类反馈结构
- 为后续把控制律与 force、mean wake、observable-related structure 接起来提供接口
因此,任何 SINDy/SR 结果都不应脱离项目总体物理主线单独解读。
---
## 二、当前已经确认的技术事实
### 1. Kármán cloak 的跨 Re_D 统一骨架存在
这是目前最硬的一批证据之一。跨 Re 的 leave-one-Re-out尤其 holdout_200已经显示
- 用 Re50 + Re100 拟合,可高精度预测 Re200
- 这说明统一骨架不是偶然的特征工程产物,而是真实存在于 PPO 策略中的共享结构
### 2. 对称性问题已经纠偏
最重要的 bug 是镜像变换 G 对动作的写法错误。
**错误版本**
\[
[a_F,a_T,a_B] \mapsto [-a_F,a_B,a_T]
\]
**正确版本**
\[
[a_F,a_T,a_B] \mapsto [-a_F,-a_B,-a_T]
\]
修正后rear equivariance 误差从约 100% 降到约 10%,原来"PPO 不尊重交换对称性"的结论应正式撤回。
### 3. front no-bias 被数据支持rear shared-head 是有效结构
### 4. one-step R² 与闭环是两回事
\[
\text{one-step R}^2 \gg 0.95 \not\Rightarrow \text{闭环好}
\]
**关键证据**full-lag 16-dim 模型 R²=0.939 但闭环仅 0.619static 8-dim 模型 R²=0.321 但闭环 0.745。
### 5. 时间特征加剧 rollout mismatch2026-06-15 确认)
训练时使用 `x(t-1)``dx/dt` 特征让 one-step R² 跃升,但闭环验证下降。原因是**分布偏移**:训练时特征来自 PPO 真实轨迹,部署时来自 SINDy 控制轨迹。带时间记忆的模型更容易在自由滚动时暴露分布偏移。
### 6. phase-state + absolute action 路线已被验证2026-06-15 核心结论)
Illusion 0.75L 和 1L 已证明**不需要动作历史特征**也能达到 PPO 的 96%+ 闭环性能:
- 输入:相位状态 `u_a, du_a/dt, Cl_tot, dCl_tot/dt, Cd_tot, Cd_rear` + error-state
- 输出:直接预测 absolute alpha不做导数积分
- 结构v23front no-bias, rear shared-head
### 7. Illusion 1.5L 是独立机制2026-06-15 确认)
1.5L 动作饱和到 [-8,8] 范围,自相关 r=0.07,线性 SINDy 不适用。
---
## 三、哪些结论现在还不能说得太满
### 1. "所有 cloak 已经共享同一骨架"
还不能这么说。all-cloak 统一骨架仍是当前主问题。
### 2. "Karman 新路线已经成功"
不能。Karman 无动作历史最佳仅 0.699(旧 v23 为 0.901),说明状态还不充分。
### 3. "高 Re 退化已经证明是采样率问题"
这是一个强工作假设,需要在时间尺度显式化后重新检验。
### 4. "SR 已经做完一轮"
如果实际做的只是 threshold 网格 + Pareto 分析,不能写成"完整 SR 已完成"。
---
## 四、代码与工程层面的已知经验
### 1. 环境分工
- `pycuda_3_10`CFD、DRL 模型加载、数据采集、SINDypysindy
- `sr_env`PySR 与 SR 相关工作
### 2. 常见坑
- **actions.npz 是归一化动作 [-1,1],不是物理 omega**。所有 SINDy 拟合代码需通过 `(norm * scale + bias) * u0` 转换。
- **Illusion "2U" 误解**"2U" 在模型名中表示 S_DIM=14多了2维目标力观测不是两倍来流速度。u0 始终为 0.01。
- **FIFO bias ≠ DRL action bias**FIFO bias 用 1U `[0, -U0, U0]` 填充历史DRL action decoder 用 2U `action*8+[0,-2,2]`
- **SAMPLE_INTERVAL 因场景而异**0.75L=400, 1L=600, 1.5L=800, Karman=800。
- **验证器 `n_steps` 不能过短**需保证控制传播到传感器S=400→320步, S=600→214步, S=800→160步。
- **保存 `controlled.npz` 时必须包含 `target_forces` 字段**Illusion 场景的 s_dim=14 推理数据)。
- **单步 validatorpredict_v23_deriv中需要传入 `sensors_raw`/`forces_raw`** 以计算相位特征中的导数项。
### 3. 当前可信结果2026-06-15
| 场景 | 路线 | 特征 | 输出 | 是否有动作历史 | 闭环 sim | % of PPO |
|------|------|------|:----:|:-------------:|:--------:|:--------:|
| illusion_0.75L | phase+error | ILLUSION_PHASE (10dim) | alpha | **否** | **0.974** | 100.2% |
| illusion_1L | phase+error | ILLUSION_PHASE (10dim) | alpha | **否** | **0.958** | 98.5% |
| illusion_1.5L | phase+error | ILLUSION_PHASE (10dim) | alpha | **否** | **N/A** | bang-bang |
| karman_re100 | phase-state | PHASE_STATE (6dim) | alpha | **否** | **0.699** | 73.3% |
| karman_re100 | old v23 | CORE_FEAT_KEYS_V2 + a_lag | alpha | **是** | **0.901** | 94.4% |
---
## 五、Bug Audit 补充2026-06-14~15
### 新发现的 Bug
| # | Bug | 文件 | 严重度 | 修复 |
|---|-----|------|--------|------|
| 13 | `run_pysr.py` 滞后构造错误:`actions_phys` 直接当 `actions_prev` 传入 `compute_features`,导致 `aF_lag1 = alpha(t)` 而非 `alpha(t-1)` — 拟合恒等式 | sindy/run_pysr.py | CRITICAL | 改为正确滞后:`a_prev[1:]=actions_phys[:-1]` |
| 14 | `predict_v23_deriv``needs_aug` 只检测 `"lag1"` 关键字,忽略 `"_dt"` 结尾的导数特征 → 闭环中 phase-state 的 `du_a_dt` 等始终为零 | validate/run_closed_loop.py | CRITICAL | `needs_aug = any(k.endswith("_dt") or k.endswith("_lag1") for k in ...)` |
| 15 | `compute_features``target_forces` 未处理 1D 输入(单步验证时 shape 为 (2,) 而非 (1,2) | utils/feature_builder.py | 中 | `if tf.ndim == 1: tf = tf.reshape(1, -1)` |
### 关键方法教训2026-06-15 更新)
1. **PPO 验证优先于 SINDy**。在确认 PPO 推理复现正确之前SINDy 结果无意义。
2. **控制时长必须 >= NX/U0**。50 步验证结果不可靠。
3. **DDF 保存位置决定验证器起始状态**。正确模式:`stabilize → save_ddf(1) → norm → apply_ddf → bias FIFO → save_ddf(2) → apply_ddf`。
4. **状态消融比输出消融更重要**`x_n` 不够 → 加 `x_{n-1}` 大幅提升 → 加 `a_{n-1}` 边际仅 3% → 说明缺的是相位状态,不是动作历史。
5. **离线 rollout 好 ≠ CFD 闭环好**。phase-state 模型 offline 50 步漂移仅 11%,但 CFD 闭环 -7% vs 静态。最终判据是 CFD 短闭环。
6. **如遇数值错误,先检查 obs 切片索引、添加顺序、target_forces 维度**

View File

@ -1,4 +1,4 @@
# SINDy 与 SR 执行计划 # SINDy 与 SR 执行计划2026-06-15 修订版 v4
## 文档作用 ## 文档作用
@ -6,342 +6,183 @@
- 只写执行路线、阶段目标、优先级、输出要求 - 只写执行路线、阶段目标、优先级、输出要求
- 不长篇复述历史争论 - 不长篇复述历史争论
- 不把背景知识和任务顺序混在一起 - 凡是背景知识、已知结论、踩坑记录,统一放到 `sindy_sr_knowledge.md`
- 凡是历史 bug、经验教训、哪些结论已经成立、哪些还不能说统一放到 `sindy_sr_knowladge`
--- ---
## 当前总目标 ## 当前总目标
**所有 cloak 场景** 纳入同一受约束分析框架,使用 **SINDy + SR 并行** 探索控制骨架,先分别拟合,再横向比较,从而判断: 用 SINDy + SR 在所有 cloak / illusion 场景上产出有物理结论的控制律公式,并验证其跨场景泛化能力。
- 哪些项是 shared core 回答三个问题:
- 哪些项是 scene-specific activation 1. **哪些物理量是控制的核心变量?**(力、速度、相位、误差之间的取舍)
- 哪些差异来自时间尺度写法,而不是物理骨架本身 2. **不同场景是否共享同一公式骨架?**family 内 -> family 间)
3. **白箱控制律能否在未训练工况工作?**(泛化性 + 失效边界)
当前默认路线不是“先强行做统一总公式”,而是:
\[
\text{separate fit} \rightarrow \text{compare} \rightarrow \text{shared-backbone test}
\]
--- ---
## 当前工作原则 ## 已完成的工作2026-06-14~15
后续默认遵守: ### 核心成果
1. **所有场景共用同一套 primitive variables 与同一套 \(G\) 规则** 1. **Phase-state 特征体系**`PHASE_STATE_KEYS` (6维) + `ILLUSION_PHASE_KEYS` (10维) — 无需动作历史
2. **front 默认 no-bias + odd structure** 2. **绝对动作输出模式**`get_feature_matrix_deriv(output_mode="absolute")` — 无积分累积
3. **rear 默认 shared-head** 3. **闭合验证器支持**`predict_v23_deriv(output_mode="absolute")` + 闭环 `mode="abs"`
4. **SINDy 与 SR 并行推进**SR 不是替换 SINDy而是并行工具 4. **离线 rollout 评估**`validate/eval_rollout.py` — 多步滚动误差分析
5. **先分场景拟合,再做横向比较** 5. **消融实验**5 种输入配置 × 2 种输出模式 × CFD 闭环验证
6. **闭环验证不可省略**,但不要求每个场景第一轮都做全套闭环
7. **时间尺度问题必须显式进入特征定义** ### Illusion 新路线成功
8. **不再只围绕跨 Re 的 Kármán 单线深挖**;跨 Re 现在是已站住的第一证据,不是全部主线
9. **当前不讨论功率与能量分析**,这不是这条线的优先任务 | 场景 | 闭环 sim | 动作历史 | 结论 |
10. **当前不把时延作为主要矛盾**;对稳定周期控制,当前第一矛盾是变量骨架与时间尺度写法 |------|:--------:|:--------:|------|
| 0.75L | **0.974** | 无 | 新路线成立 |
| 1L | **0.958** | 无 | 新路线成立 |
| 1.5L | N/A | 无 | bang-bang/不同机制 |
### Karman 需要补状态
| 配置 | 闭环 sim | 对比旧 v23 |
|------|:--------:|:----------:|
| 旧 v23 (a_lag) | 0.901 | 基线 |
| phase-state → abs (最佳新) | **0.699** | -22% |
| phase-state → deriv | 0.656 | -27% |
--- ---
## 当前场景优先级 ## 当前最值得优先做的实验
### 第一层:本轮重点场景 ### 第一优先级Illusion 0.75L / 1L → separate PySR → 公式比较
- Kármán cloak 输入特征10 维):
- steady cloak ```python
ILLUSION_PHASE_KEYS = [
"u_a", "du_a_dt", # 振荡相位
"Cl_tot", "dCl_tot_dt", # 升力动力学
"Cd_tot", "Cd_rear", # 阻力反馈
"Cd_err", "Cl_err", # 力误差
"dCd_err_dt", "dCl_err_dt", # 误差动力学
]
```
这两个场景本轮必须做完整输出,因为它们最适合先建立场景间比较模板。 输出absolute alpha非维动作不积分
环境:`conda run -n sr_env`
代码:`sindy/run_pysr.py`
### 第二层:下一轮扩展场景 三个产出要求:
- 最短可接受公式
- 闭环可运行公式
- 0.75L 与 1L 的公式形态比较 / 同形骨架判断
- 单涡 cloakmonopole / taylor / lamb 等已有单涡场景) ### 第二优先级Karman 状态补强CCD/OID 进场)
- erase
- 其他已有 cloak 场景
这批先做轻量版 separate fit再决定哪些值得补完整闭环与深挖。 当前 Karman 新路线最佳 0.699,需要在 phase-state 基础上补充缺失状态量:
- 候选回流区长度、尾迹中心线偏移、POD 模态系数
- 方法CCDCorrelation-based Decomposition或 OIDObserver-based Identification
- 目标:在无动作历史前提下将 Karman 闭环提升到 0.85+
### 第三优先级Illusion 跨直径联合
条件0.75L 和 1L 的 PySR 公式确认同形骨架后
方法error-state 特征 + 目标 signature 描述符target St、amplitude
--- ---
## 阶段 0 ## 当前暂不建议做的事
## 统一接口 - 不回退到动作历史主导Illusion 已经证明不需要)
- 不先做跨 Re Karman 联合拟合(单 Re 还没站稳)
### 0.1 统一 primitive variables - 不先做 vortex 扩展(不是当前瓶颈)
- 不先上 PySR 到 Karman状态还有问题公式不会更物理
所有场景统一输出:
- 无量纲 sensor\(\hat u, \hat v\)
- 力系数:\(C_D, C_L\)
- 无量纲动作:\(\alpha\)
- lagged \(\alpha\)
- \(\Delta \alpha\) 或显式含 \(\Delta t_c\) 的版本
- \(\mu = 1/Re_D\)
- scene metadatascene id、\(Re_D\)、control interval \(\Delta t_c\)、target type、采样设置
### 0.2 固定 \(G\) 算子
统一使用同一套 \(G\) 规则,不允许不同场景临时改写。动作必须满足:
\[
(\alpha_F,\alpha_T,\alpha_B) \mapsto (-\alpha_F,-\alpha_B,-\alpha_T)
\]
### 0.3 统一 feature builder
统一生成三层特征:
| 层级 | 内容 | 用途 |
|---|---|---|
| core | \(\hat u, \hat v, C_D, C_L, \alpha^-, \Delta\alpha^-, \mu\) | 所有场景共用 |
| derived | 对称/反对称组合、总量/差量 | SINDy 主库 |
| time-scale | 显式含 \(\Delta t_c\) 的版本 | 时间尺度分析 |
### 0.4 必做测试
任何场景进入拟合前,先过:
- \(G(G(x)) = x\)
- feature names 与矩阵列顺序一致
- 闭环预测器输入维度一致
- front / rear 的结构约束在数据接口层正确落地
### 阶段 0 输出
- 统一 `feature_builder`
- 统一 `symmetry`
- 统一 `time_scale`
- 统一 `validators`
- 场景级最小数据摘要样本数、变量范围、control interval
--- ---
## 阶段 1 ## 常用命令
## Kármán 与 steady 的第一轮 separate SINDy ### SINDy 拟合
```bash
# Illusion phase-state + absolute action
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes illusion_0.75L,illusion_1L --deriv --phase --output-mode absolute
### 目标 # Karman phase-state + absolute action
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re100 --deriv --phase --output-mode absolute
在统一变量与统一约束下,先得到两个重点场景各自最可信的稀疏 support。 # Karman expanded features
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re100 --deriv --karman-expand --output-mode absolute
### 默认约束 # Karman phase-state + mu modulation
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re100 --deriv --karman-mu --output-mode absolute
- front no-bias # Illusion separate with error-state (old style for comparison)
- front odd structure conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
- rear shared-head --scenes illusion_0.75L,illusion_1L,illusion_1.5L
- correct \(G\) consistency
### 每个场景必须输出 # Karman cross-Re joint
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes "karman_re50,karman_re100,karman_re200,karman_re400" --joint
```
| 输出 | 说明 | ### CFD 闭环验证
|---|---| ```bash
| best support | 最优 support 列表 | # Illusion phase-state + absolute action
| sparsity curve | 稀疏度-误差曲线 | conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop_illusion.py \
| front / rear 主项表 | 主导项与系数 | --scene illusion_0.75L --device 0 --steps 320 \
| one-step metrics | R²、RMSE | --sindy-results src/SR_analysis/sindy/illusion/sindy_results_deriv.json
| key closed-loop result | 至少一个关键闭环指标 |
| support stability | threshold / window / bootstrap 稳定性 |
| contribution table | 主要项贡献度,不只看是否出现 |
### 本阶段比较项 # Karman phase-state + absolute action
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \
--scene karman_re100 --device 1 --steps 200 --mode abs \
--sindy-results src/SR_analysis/sindy/karman/sindy_results_deriv.json
Kármán 与 steady 做第一轮 support overlap # Karman old v23 for comparison
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \
--scene karman_re100 --device 1 --steps 200 --mode v23 \
--sindy-results src/SR_analysis/sindy/karman/sindy_joint_wrapped.json
```
- 哪些项共同出现 ### 离线 rollout 评估
- 哪些项只在 Kármán 激活 ```bash
- steady 是否表现为明显简化版 python3 src/SR_analysis/validate/eval_rollout.py \
- overlap 不能只看布尔出现,还要看贡献量级与稳定性 --sindy-results src/SR_analysis/sindy/karman/sindy_results_deriv.json \
--scene karman_re100
```
### PySR 符号回归
```bash
conda run -n sr_env python src/SR_analysis/sindy/run_pysr.py \
--scene illusion_1L
```
--- ---
## 阶段 2 ## 关键文件索引
## Kármán 与 steady 的第一轮受限 SR | 文件 | 用途 |
|------|------|
### 目标 | `configs.py` | 统一场景元数据10+场景) |
| `utils/feature_builder.py` | 特征工程PHYSICS_FEAT_KEYS, PHASE_STATE_KEYS, ILLUSION_PHASE_KEYS, 导数/滞后特征 |
不是追求最终公式,而是看: | `utils/sindy_fitter.py` | STLSQ拟合get_feature_matrix_v2, get_feature_matrix_deriv, compute_action_deriv |
| `utils/cfd_interface.py` | LegacyCelerisLab封装GPU, 推理, norm |
- 在受限物理库上,是否能压出更短闭式 | `utils/g_operator.py` | 等变性诊断 |
- SINDy 中看似不同的项,是否会被更统一的表达吸收 | `scripts/infer_*.py` | Karman/Illusion/Vortex 推理管线 |
- Kármán 与 steady 是否出现同形公式 | `sindy/run_all_v2.py` | 统一 SINDy 拟合入口(支持 --deriv, --phase, --output-mode, --karman-expand 等) |
| `sindy/run_pysr.py` | 受限 PySR 符号回归 |
### SR 输入规则 | `sindy/wrap_joint.py` | 联合模型 → wrapped 格式 |
| `validate/run_closed_loop.py` | Karman 闭环验证器v23, deriv, abs 模式) |
SR 只能使用: | `validate/run_closed_loop_illusion.py` | Illusion 闭环验证器 |
| `validate/eval_rollout.py` | **新**:离线多步 rollout 评估 |
- 阶段 0 的统一变量 | `compare/support_overlap.py` | 跨场景 support 比较 |
- 阶段 1 的 SINDy 已筛出主项及其邻近项 | `compare/shared_core.py` | 多场景 shared core 检测 |
- 受限运算集合 | `docs/SR_analysis_results.md` | **新**:完整分析报告 |
| `docs/figures/SR_analysis/fig*.png` | **新**结果图表6张 |
### 当前允许的运算
## 环境
- 加减乘
- protected divide | 环境 | conda 名 | 用途 |
- 少量 square |------|---------|------|
- 必要时有限放开 tanh | CFD + DRL + SINDy | `pycuda_3_10` | infer, run_all_v2, closed-loop |
| PySR 符号回归 | `sr_env` | run_pysr.py |
### 当前不允许的运算 | GPU | device 0, device 1 | |
- raw trig 乱搜
- 高次幂
- 指数
- 深层嵌套
### 每个场景必须输出
| 输出 | 说明 |
|---|---|
| shortest acceptable formula | 最短可接受公式 |
| complexity-error pareto | 复杂度 vs 误差 |
| 和 SINDy 的关系 | 压缩了什么、保留了什么、吸收了什么 |
| key closed-loop result | 最佳 SR 公式至少一版关键闭环 |
| formula family notes | 同一场景内是否存在多种同等可接受闭式 |
---
## 阶段 3
## 横向比较
当 Kármán 与 steady 的 SINDy + SR 都出来后,立即做横向比较。
### 3.1 support 比较
不要只看“是否出现”,必须同时看:
- 是否出现
- 系数或贡献量级
- 稳定性
- SR 是否把它吸收到更高层表达里
### 3.2 公式形态比较
至少检查:
- 是否都包含同类 force feedback 核心
- 是否都包含同类 memory 核心
- steady 是否只是删掉了周期相关项
- 是否存在同形结构 + 少量场景激活项
### 3.3 第一轮 shared-backbone 判断
这一轮只回答:
1. 是否存在 Kármán 与 steady 的 shared core
2. steady 是否可以视为 Kármán 的明显简化版
3. 哪些项更像 scene-specific activation而不是 backbone 本身
本阶段不急着拟合 all-cloak 联合总公式。
---
## 阶段 4
## 扩展到单涡与其他 cloak
在 Kármán 与 steady 的第一轮比较完成后,再把单涡 cloak 与其他场景纳入。
### 轻量版输出要求
- separate SINDy
- separate SR
- one-step metrics
- support 形态
- 必要时补关键闭环
### 比较目标
重点看:
- 单涡是否保留 shared core
- 单涡是否主要新增 history / transient 项
- 是否开始出现子家族结构
---
## 阶段 5
## 时间尺度显式化
这条线并行推进,但先不抢在全部场景前面。
### 当前目标
- 让 \(\Delta t_c\) 显式进入特征
- 不再把“1 个采样步”默认当物理可比量
- 比较显式化前后support 是否更收敛
- 为后面严肃讨论采样率影响扫清接口问题
### 第一批测试场景
- Kármán cloak
- steady cloak
### 第一批对比内容
- 旧的 discrete lag / \(\Delta a\)
- 显式含 \(\Delta t_c\) 的版本
### 关注结果
- support 是否变化
- SR 公式是否更统一
- 跨场景比较是否更干净
注意:当前阶段的目标是**时间尺度显式化**,不是立刻给出高采样率优于 PPO 的最终结论。
---
## 本轮必须完成的最小任务
1. 统一 Kármán 与 steady 的 feature builder
2. 统一 Kármán 与 steady 的 \(G\) / time-scale / validators
3. 跑 Kármán 与 steady 的第一轮 separate SINDy
4. 跑 Kármán 与 steady 的第一轮受限 SRPySR
5. 输出 Kármán vs steady 的:
- support overlap
- 公式形态比较
- shared core / scene-specific 的初步分类
---
## 本轮结束时应交付的结果包
每个重点场景至少有一张 summary 表:
| method | sparsity | one-step | closed-loop | key terms | notes |
|---|---:|---:|---:|---|---|
| SINDy | | | | | |
| SR | | | | | |
以及一个跨场景比较表:
| comparison | shared core | scene-enhanced | scene-specific | notes |
|---|---|---|---|---|
| Kármán vs steady | | | | |
---
## 当前不该做的事
- 不继续围绕单一 Kármán across Re 版本号升级
- 不把 threshold 网格或简单 Pareto 扫描直接当成完整 SR
- 不在 raw feature 上做自由 SR
- 不在不同采样间隔下直接复用旧系数并据此下正式结论
- 不在场景比较证据还弱时,提前宣布 all-cloak 统一总公式
- 不把功率、能量、时延这些非当前主矛盾问题拉入本轮 SINDy/SR 主线
---
## 当前这份 notes 的直接收束
接下来核心任务不是“继续优化某个跨 Re 模型”,而是:
\[
\boxed{
\text{把 Kármán 与 steady 先做成可比较的 separate SINDy + separate SR 结果包,然后以它们为模板扩到单涡与其他 cloak。}
}
\]
这一步完成后,才进入更严肃的 all-cloak shared-backbone 判断。

View File

@ -1,9 +1,19 @@
from .feature_builder import ( from .feature_builder import (
compute_dimensionless, compute_features, build_feature_matrix, compute_dimensionless, compute_features, build_feature_matrix,
get_feature_names, apply_G_alpha, apply_G_x, get_feature_names, apply_G_alpha, apply_G_x,
CORE_FEAT_KEYS, MU_FEAT_KEYS, ALL_FEAT_KEYS, CORE_FEAT_KEYS, CORE_FEAT_KEYS_V2, TIME_FEAT_KEYS,
MU_FEAT_KEYS, ALL_FEAT_KEYS, ALL_FEAT_KEYS_V2,
ILLUSION_TARGET_KEYS, ILLUSION_FEAT_KEYS_V2, ILLUSION_ALL_FEAT_KEYS_V2,
PHYSICS_FEAT_KEYS, ILLUSION_ERR_KEYS,
LAG_FEAT_KEYS, DERIV_FEAT_KEYS, ACTION_LAG_KEYS,
AUG_LEVEL_1_KEYS, AUG_LEVEL_2_KEYS, AUG_LEVEL_3_KEYS, AUG_LEVEL_4_KEYS,
ILLUSION_AUG_LEVEL_1_KEYS, ILLUSION_AUG_LEVEL_2_KEYS,
ILLUSION_AUG_LEVEL_3_KEYS, ILLUSION_AUG_LEVEL_4_KEYS,
PHASE_STATE_KEYS, ILLUSION_PHASE_KEYS, KARMAN_EXPANDED_KEYS,
) )
from .sindy_fitter import ( from .sindy_fitter import (
fit_channel, fit_sindy, print_control_law, fit_channel, fit_sindy, print_control_law,
get_active_support, get_feature_matrix_from_data, get_active_support, get_feature_matrix_from_data,
fit_sindy_weighted, get_feature_matrix_v2,
compute_action_deriv, get_feature_matrix_deriv,
) )

View File

@ -208,6 +208,12 @@ def add_pinball(
fifo.append(flow_field.obs.copy()[obs_slice_start:obs_slice_end]) fifo.append(flow_field.obs.copy()[obs_slice_start:obs_slice_end])
save_states = np.array(list(fifo), dtype=data_type) save_states = np.array(list(fifo), dtype=data_type)
# CRITICAL: save DDF again AFTER bias FIFO, so restore_ddf() goes
# to the post-bias state (consistent with saved FIFO).
# Without this, reset() restores to a bare stabilized state with
# no bias history, invalidating the FIFO.
flow_field.get_ddf()
flow_field.save_ddf()
flow_field.apply_ddf() flow_field.apply_ddf()
return { return {
@ -283,7 +289,7 @@ def vorticity_from_ddf(flow_field: FlowField, u0: float) -> np.ndarray:
- ddf[:, :, 3] - ddf[:, :, 6] - ddf[:, :, 7]) / u0 - ddf[:, :, 3] - ddf[:, :, 6] - ddf[:, :, 7]) / u0
uy = (ddf[:, :, 2] + ddf[:, :, 5] + ddf[:, :, 6] uy = (ddf[:, :, 2] + ddf[:, :, 5] + ddf[:, :, 6]
- ddf[:, :, 4] - ddf[:, :, 7] - ddf[:, :, 8]) / u0 - ddf[:, :, 4] - ddf[:, :, 7] - ddf[:, :, 8]) / u0
omega = np.gradient(uy, axis=1) - np.gradient(ux, axis=0) omega = np.gradient(uy, axis=0) - np.gradient(ux, axis=1)
return omega.astype(np.float64) return omega.astype(np.float64)

View File

@ -7,7 +7,7 @@ Copy of analysis_cloak/common/feature_builder.py -- kept as canonical source.
""" """
from __future__ import annotations from __future__ import annotations
from typing import Dict, List, Tuple from typing import Dict, List, Optional, Tuple
import numpy as np import numpy as np
@ -86,6 +86,7 @@ def apply_G_x(dim: Dict[str, np.ndarray],
# -- Feature key definitions ------------------------------------------------- # -- Feature key definitions -------------------------------------------------
# Original feature set (includes sin_ua/cos_ua)
CORE_FEAT_KEYS = [ CORE_FEAT_KEYS = [
"u_m", "u_a", "u_c", "u_m", "u_a", "u_c",
"v_a", "v_a",
@ -96,9 +97,99 @@ CORE_FEAT_KEYS = [
"daF", "daB", "daT", "daF", "daB", "daT",
] ]
# V2 core features: no sin_ua/cos_ua, no mu (for single-scene fitting)
CORE_FEAT_KEYS_V2 = [
"u_m", "u_a", "u_c",
"v_a",
"Cd_tot", "Cd_rear",
"Cl_tot", "Cl_diff",
"aF_lag1", "aB_lag1", "aT_lag1",
"daF", "daB", "daT",
]
# Time-explicit features: da/dt_c (for cross-scene coefficient comparison)
TIME_FEAT_KEYS = [
"daF_dt", "daB_dt", "daT_dt",
]
MU_FEAT_KEYS = ["mu", "mu_u_a", "mu_v_a", "mu_Cd_tot", "mu_Cl_diff"] MU_FEAT_KEYS = ["mu", "mu_u_a", "mu_v_a", "mu_Cd_tot", "mu_Cl_diff"]
# Illusion target force features (for scenes where we know the target Cd/Cl)
ILLUSION_TARGET_KEYS = ["target_Cd", "target_Cl"]
ILLUSION_FEAT_KEYS_V2 = CORE_FEAT_KEYS_V2 + ILLUSION_TARGET_KEYS
ILLUSION_ALL_FEAT_KEYS_V2 = ILLUSION_FEAT_KEYS_V2 + MU_FEAT_KEYS # with mu for joint
# Physics-only feature keys: NO action history terms (no aF_lag1, no daF, etc.)
# These are the clean inputs for learning d(alpha)/dt as a function of physics state.
PHYSICS_FEAT_KEYS = [
"u_m", "u_a", "u_c", "v_a",
"Cd_tot", "Cd_rear", "Cl_tot", "Cl_diff",
]
# Illusion error-state: physics features + force error (not raw target forces)
ILLUSION_ERR_KEYS = PHYSICS_FEAT_KEYS + ["Cd_err", "Cl_err"]
# Observation lag (1-step delayed) and derivative (time-normalized) features
# These add temporal/phasing information to the otherwise static physics features.
LAG_FEAT_KEYS = [
"u_m_lag1", "u_a_lag1", "u_c_lag1", "v_a_lag1",
"Cd_tot_lag1", "Cd_rear_lag1", "Cl_tot_lag1", "Cl_diff_lag1",
]
DERIV_FEAT_KEYS = [
"du_m_dt", "du_a_dt", "du_c_dt", "dv_a_dt",
"dCd_tot_dt", "dCd_rear_dt", "dCl_tot_dt", "dCl_diff_dt",
]
# Action lag (for ablation level 4 — comparing against old approach)
ACTION_LAG_KEYS = ["aF_lag1", "aB_lag1", "aT_lag1"]
# Augmented levels for ablation study:
# Level 0 = x_n (PHYSICS_FEAT_KEYS, no memory at all)
# Level 1 = x_n + x_{n-1} (current + 1-step lag)
# Level 2 = x_n + dx/dt (current + derivative)
# Level 3 = x_n + x_{n-1} + dx/dt (full temporal context)
# Level 4 = Level 3 + a_{n-1} (add action history for comparison)
AUG_LEVEL_1_KEYS = PHYSICS_FEAT_KEYS + LAG_FEAT_KEYS
AUG_LEVEL_2_KEYS = PHYSICS_FEAT_KEYS + DERIV_FEAT_KEYS
AUG_LEVEL_3_KEYS = PHYSICS_FEAT_KEYS + LAG_FEAT_KEYS + DERIV_FEAT_KEYS
AUG_LEVEL_4_KEYS = AUG_LEVEL_3_KEYS + ACTION_LAG_KEYS
# Illusion equivalents with error-state
ILLUSION_LAG_KEYS = [
"Cd_err_lag1", "Cl_err_lag1",
]
ILLUSION_DERIV_KEYS = [
"dCd_err_dt", "dCl_err_dt",
]
ILLUSION_AUG_LEVEL_1_KEYS = ILLUSION_ERR_KEYS + LAG_FEAT_KEYS + ILLUSION_LAG_KEYS
ILLUSION_AUG_LEVEL_2_KEYS = ILLUSION_ERR_KEYS + DERIV_FEAT_KEYS + ILLUSION_DERIV_KEYS
ILLUSION_AUG_LEVEL_3_KEYS = ILLUSION_AUG_LEVEL_1_KEYS + DERIV_FEAT_KEYS + ILLUSION_DERIV_KEYS
ILLUSION_AUG_LEVEL_4_KEYS = ILLUSION_AUG_LEVEL_3_KEYS + ACTION_LAG_KEYS
# Phase-state features: low-dimensional dynamic state z = [phase_obs, phase_deriv, static_force]
# This replaces the full x_n + x_{n-1} with a compact representation.
# The idea: u_a + du_a/dt encodes oscillation phase, Cl_tot + dCl_tot/dt encodes lift dynamics,
# and Cd_tot/Cd_rear provide static force feedback.
PHASE_STATE_KEYS = [
"u_a", "du_a_dt", # oscillation phase (cross-stream asymmetry + rate)
"Cl_tot", "dCl_tot_dt", # lift dynamics
"Cd_tot", "Cd_rear", # static drag feedback
]
# Karman expanded: phase-state + supplementary static quantities
# For Karman, 6-dim phase-state may not capture enough information.
# Adding u_m (mean streamwise), u_c (center sensor), v_a (cross asymmetry), Cl_diff (lift distribution)
KARMAN_EXPANDED_KEYS = PHASE_STATE_KEYS + [
"u_m", "u_c", "v_a", "Cl_diff",
]
# Illusion phase-state: error-based + dynamics
ILLUSION_PHASE_KEYS = PHASE_STATE_KEYS + [
"Cd_err", "Cl_err", "dCd_err_dt", "dCl_err_dt",
]
ALL_FEAT_KEYS = CORE_FEAT_KEYS + MU_FEAT_KEYS ALL_FEAT_KEYS = CORE_FEAT_KEYS + MU_FEAT_KEYS
# V2 all features (for cross-Re joint fitting with mu)
ALL_FEAT_KEYS_V2 = CORE_FEAT_KEYS_V2 + MU_FEAT_KEYS
# -- Feature computation ----------------------------------------------------- # -- Feature computation -----------------------------------------------------
@ -110,7 +201,13 @@ def compute_features(
mu: float, mu: float,
alpha_mode: bool = False, # if True, actions_prev are already nondim alpha alpha_mode: bool = False, # if True, actions_prev are already nondim alpha
include_mu: bool = True, include_mu: bool = True,
include_cos_sin: bool = True, # if True, include sin_ua/cos_ua features
dt_c: float = 1.0, # control interval in T0 units (for time normalization)
u0: float = U0, # inlet velocity for omega->alpha conversion u0: float = U0, # inlet velocity for omega->alpha conversion
target_forces: Optional[np.ndarray] = None, # (T, 2) raw lattice target forces [fx,fy]
rho: float = 1.0, # fluid density for Cd/Cl conversion
sensors_raw: Optional[np.ndarray] = None, # (T, 6) raw lattice sensors, for obs dynamics
forces_raw: Optional[np.ndarray] = None, # (T, 6) raw lattice forces, for obs dynamics
) -> Dict[str, np.ndarray]: ) -> Dict[str, np.ndarray]:
"""Compute unified feature dictionary from dimensionless primitives. """Compute unified feature dictionary from dimensionless primitives.
@ -121,6 +218,8 @@ def compute_features(
mu: 1/Re_D mu: 1/Re_D
alpha_mode: if True, actions are already nondim; else convert alpha_mode: if True, actions are already nondim; else convert
include_mu: include mu modulation terms include_mu: include mu modulation terms
include_cos_sin: include sin_ua/cos_ua phase encoding
dt_c: control interval (in T0 = D/U0 units), for time-normalized deltas
u0: inlet velocity (lattice), used only when alpha_mode=False u0: inlet velocity (lattice), used only when alpha_mode=False
Returns dict with all features as (T,) or (T,3) arrays. Returns dict with all features as (T,) or (T,3) arrays.
@ -153,11 +252,12 @@ def compute_features(
sym["Cl_tot"] = Cl_F + Cl_T + Cl_B sym["Cl_tot"] = Cl_F + Cl_T + Cl_B
sym["Cl_diff"] = Cl_T - Cl_B sym["Cl_diff"] = Cl_T - Cl_B
# Phase # Phase (optional, may obscure linear structure)
sym["sin_ua"] = np.sin(np.pi * sym["u_a"]) if include_cos_sin:
sym["cos_ua"] = np.cos(np.pi * sym["u_a"]) sym["sin_ua"] = np.sin(np.pi * sym["u_a"])
sym["cos_ua"] = np.cos(np.pi * sym["u_a"])
# Memory (nondim alpha) # Memory (nondim alpha) -- discrete version
sym["aF_lag1"] = a[:, 0] sym["aF_lag1"] = a[:, 0]
sym["aB_lag1"] = a[:, 1] sym["aB_lag1"] = a[:, 1]
sym["aT_lag1"] = a[:, 2] sym["aT_lag1"] = a[:, 2]
@ -165,6 +265,22 @@ def compute_features(
sym["daB"] = a[:, 1] - a2[:, 1] sym["daB"] = a[:, 1] - a2[:, 1]
sym["daT"] = a[:, 2] - a2[:, 2] sym["daT"] = a[:, 2] - a2[:, 2]
# Time-normalized deltas (for cross-scene coefficient comparison)
sym["daF_dt"] = sym["daF"] / dt_c
sym["daB_dt"] = sym["daB"] / dt_c
sym["daT_dt"] = sym["daT"] / dt_c
# Target forces (for illusion scenes) — convert lattice forces to Cd/Cl
if target_forces is not None:
tf = np.asarray(target_forces, dtype=np.float64)
if tf.ndim == 1:
tf = tf.reshape(1, -1)
sym["target_Cd"] = 2.0 * tf[:, 0] / (rho * u0**2 * D_CYL)
sym["target_Cl"] = 2.0 * tf[:, 1] / (rho * u0**2 * D_CYL)
# Error-state: deviation between actual and target force
sym["Cd_err"] = sym["Cd_tot"] - sym["target_Cd"]
sym["Cl_err"] = sym["Cl_tot"] - sym["target_Cl"]
# Mu modulation # Mu modulation
if include_mu: if include_mu:
sym["mu"] = np.full(T, mu, dtype=np.float64) sym["mu"] = np.full(T, mu, dtype=np.float64)
@ -172,6 +288,72 @@ def compute_features(
sym["mu_v_a"] = sym["v_a"] * mu sym["mu_v_a"] = sym["v_a"] * mu
sym["mu_Cd_tot"] = sym["Cd_tot"] * mu sym["mu_Cd_tot"] = sym["Cd_tot"] * mu
sym["mu_Cl_diff"] = sym["Cl_diff"] * mu sym["mu_Cl_diff"] = sym["Cl_diff"] * mu
sym["mu_Cl_tot"] = sym["Cl_tot"] * mu # additional for phase-state
# Observation dynamics: 1-step lag + time-normalized derivative
# These add temporal/phase info that static features lack.
if sensors_raw is not None and forces_raw is not None:
sr = np.asarray(sensors_raw, dtype=np.float64)
fr = np.asarray(forces_raw, dtype=np.float64)
# Lag-1 observations (shift by 1, pad first with current)
s_lag1 = np.zeros_like(sr)
f_lag1 = np.zeros_like(fr)
s_lag1[1:] = sr[:-1]
f_lag1[1:] = fr[:-1]
# Compute dim for lag-1
dim_lag1 = compute_dimensionless(s_lag1, f_lag1, u0=u0, d=D_CYL, rho=rho)
# Compute combined features for lag-1
uB_1, uC_1, uT_1 = dim_lag1["u_hat_B"], dim_lag1["u_hat_C"], dim_lag1["u_hat_T"]
vB_1, vC_1, vT_1 = dim_lag1["v_hat_B"], dim_lag1["v_hat_C"], dim_lag1["v_hat_T"]
CdF_1, CdT_1, CdB_1 = dim_lag1["Cd_F"], dim_lag1["Cd_T"], dim_lag1["Cd_B"]
ClF_1, ClT_1, ClB_1 = dim_lag1["Cl_F"], dim_lag1["Cl_T"], dim_lag1["Cl_B"]
u_m_1 = (uB_1 + uC_1 + uT_1) / 3.0
u_a_1 = (uT_1 - uB_1) / 2.0
u_c_1 = uC_1.copy()
v_a_1 = (vT_1 - vB_1) / 2.0
Cd_tot_1 = CdF_1 + CdT_1 + CdB_1
Cd_rear_1 = CdT_1 + CdB_1
Cl_tot_1 = ClF_1 + ClT_1 + ClB_1
Cl_diff_1 = ClT_1 - ClB_1
# Trim to T (in case raw arrays are longer — e.g. validator passes 2 rows)
def _trim(x):
return x[-T:] if x.shape[0] > T else x
# Store lag-1 features (trimmed to T)
sym["u_m_lag1"] = _trim(u_m_1)
sym["u_a_lag1"] = _trim(u_a_1)
sym["u_c_lag1"] = _trim(u_c_1)
sym["v_a_lag1"] = _trim(v_a_1)
sym["Cd_tot_lag1"] = _trim(Cd_tot_1)
sym["Cd_rear_lag1"] = _trim(Cd_rear_1)
sym["Cl_tot_lag1"] = _trim(Cl_tot_1)
sym["Cl_diff_lag1"] = _trim(Cl_diff_1)
# Time-normalized derivatives: (current - lag1) / dt_c
eps = 1e-12
sym["du_m_dt"] = _trim((sym["u_m"] - u_m_1) / (dt_c + eps))
sym["du_a_dt"] = _trim((sym["u_a"] - u_a_1) / (dt_c + eps))
sym["du_c_dt"] = _trim((sym["u_c"] - u_c_1) / (dt_c + eps))
sym["dv_a_dt"] = _trim((sym["v_a"] - v_a_1) / (dt_c + eps))
sym["dCd_tot_dt"] = _trim((sym["Cd_tot"] - Cd_tot_1) / (dt_c + eps))
sym["dCd_rear_dt"] = _trim((sym["Cd_rear"] - Cd_rear_1) / (dt_c + eps))
sym["dCl_tot_dt"] = _trim((sym["Cl_tot"] - Cl_tot_1) / (dt_c + eps))
sym["dCl_diff_dt"] = _trim((sym["Cl_diff"] - Cl_diff_1) / (dt_c + eps))
# Illusion error-state dynamics
if target_forces is not None:
tf = np.asarray(target_forces, dtype=np.float64)
tf_lag1 = np.zeros_like(tf)
tf_lag1[1:] = tf[:-1]
tCd_1 = 2.0 * tf_lag1[:, 0] / (rho * u0**2 * D_CYL)
tCl_1 = 2.0 * tf_lag1[:, 1] / (rho * u0**2 * D_CYL)
sym["Cd_err_lag1"] = _trim(Cd_tot_1 - tCd_1)
sym["Cl_err_lag1"] = _trim(Cl_tot_1 - tCl_1)
sym["dCd_err_dt"] = _trim((sym["Cd_err"] - sym["Cd_err_lag1"]) / (dt_c + eps))
sym["dCl_err_dt"] = _trim((sym["Cl_err"] - sym["Cl_err_lag1"]) / (dt_c + eps))
return sym return sym

View File

@ -179,3 +179,325 @@ def get_feature_matrix_from_data(
return (Theta_f[n_warmup:], Theta_r[n_warmup:], return (Theta_f[n_warmup:], Theta_r[n_warmup:],
actions_phys[n_warmup:], actions_phys[n_warmup:],
feat_names_front, feat_names_rear) feat_names_front, feat_names_rear)
# ---------------------------------------------------------------------------
# Weighted STLSQ with Huber-like robust regression
# ---------------------------------------------------------------------------
def _huber_weights(residuals: np.ndarray, c: float = 1.345) -> np.ndarray:
"""Compute Huber-like weights from residuals.
Args:
residuals: (T,) residual array.
c: tuning constant (default 1.345 gives 95% efficiency for Normal errors).
Returns:
weights: (T,) weight array in [0, 1].
"""
s = np.median(np.abs(residuals)) * 1.4826 # robust scale estimate (MAD)
if s < 1e-12:
s = 1.0
r = np.abs(residuals) / s
w = np.where(r <= c, 1.0, c / r)
return np.asarray(w, dtype=np.float64)
def fit_sindy_weighted(
Theta: np.ndarray,
y: np.ndarray,
thresholds: Optional[List[float]] = None,
alpha: float = 1e-4,
max_iter: int = 25,
sample_weights: Optional[np.ndarray] = None,
n_robust_passes: int = 2,
) -> List[dict]:
"""Run SINDy with threshold grid and optional robust weighting.
Two-stage robust fitting:
1. First pass: OLS, compute residuals, compute Huber weights
2. Second pass: weighted STLSQ with Huber weights
Args:
Theta: (T, N) feature matrix.
y: (T,) target.
thresholds: list of threshold values.
alpha: ridge regularization.
max_iter: max STLSQ iterations.
sample_weights: optional (T,) pre-defined sample weights.
n_robust_passes: number of robust re-weighting passes (1 = skip).
Returns:
results: list of dict per threshold with keys:
threshold, nz, r2, mae, coef, weights_used
"""
import pysindy as ps
if thresholds is None:
thresholds = DEFAULT_THRESHOLDS
# Normalise features for thresholding stability
std = np.std(Theta, axis=0)
std = np.where(std < 1e-8, 1.0, std)
Theta_s = Theta / std
# Initialize weights
if sample_weights is not None:
w = np.asarray(sample_weights, dtype=np.float64).flatten()
w = w / np.mean(w) # normalize to mean=1
else:
w = np.ones(Theta.shape[0], dtype=np.float64)
# Robust re-weighting passes
for _ in range(n_robust_passes - 1):
# OLS on weighted data
Theta_w = Theta_s * np.sqrt(w)[:, None]
y_w = y * np.sqrt(w)
coef_ols, _, _, _ = np.linalg.lstsq(Theta_w, y_w, rcond=None)
coef_ols = coef_ols.flatten() / std
resid = y - Theta @ coef_ols
w_new = _huber_weights(resid)
if sample_weights is not None:
w = w * w_new
else:
w = w_new
w = w / np.mean(w)
results = []
for th in thresholds:
if np.max(w) > 1e-8:
# Weighted STLSQ: apply weights via sample_weight
opt = ps.STLSQ(threshold=th, alpha=alpha, max_iter=max_iter)
opt.fit(Theta_s, y, sample_weight=w)
coef = np.asarray(opt.coef_, dtype=np.float64).flatten() / std
else:
coef = np.zeros(Theta.shape[1], dtype=np.float64)
y_pred = Theta @ coef
# Weighted R2
ssr = float(np.sum(w * (y - y_pred) ** 2))
sst = float(np.sum(w * (y - np.average(y, weights=w)) ** 2) + 1e-12)
r2 = 1.0 - ssr / sst
mae = float(np.mean(np.abs(y - y_pred)))
nz = int(np.sum(np.abs(coef) > 1e-8))
results.append({
"threshold": float(th),
"nz": nz,
"r2": r2,
"mae": mae,
"coef": [float(c) for c in coef],
"weights_min": float(np.min(w)),
"weights_max": float(np.max(w)),
})
return results
# ---------------------------------------------------------------------------
# V2 feature matrix builder (configurable feature sets, time normalization)
# ---------------------------------------------------------------------------
def get_feature_matrix_v2(
sensors: np.ndarray, # (T, 6)
forces: np.ndarray, # (T, 6)
actions_phys: np.ndarray, # (T, 3) physical omega
mu: float,
u0: float = U0,
alpha_mode: bool = False,
include_mu: bool = False, # default False for single-scene
include_cos_sin: bool = False, # default False (avoid masking linear structure)
use_time_norm: bool = False, # if True, use time-normalized deltas (da/dt_c)
feat_keys: Optional[List[str]] = None, # custom feat keys (default CORE_FEAT_KEYS_V2)
dt_c: float = 1.0, # control interval in T0 units
n_warmup: int = 2,
target_forces: Optional[np.ndarray] = None, # (T, 2) raw lattice target forces
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, List[str], List[str]]:
"""Build feature matrices with configurable feature sets.
Uses CORE_FEAT_KEYS_V2 by default (no sin_ua/cos_ua, no mu).
For Illusion scenes with target forces, provide target_forces and
feat_keys=ILLUSION_FEAT_KEYS_V2 (or pass None to auto-detect).
For cross-Re joint fitting, set include_mu=True.
For time-normalized deltas, set use_time_norm=True and provide dt_c.
Returns
-------
Theta_front : (T-warmup, N_front) NO bias column
Theta_rear : (T-warmup, N_rear) WITH bias column
Y : (T-warmup, 3) target actions (physical omega)
feat_names_front, feat_names_rear
"""
from .feature_builder import (
CORE_FEAT_KEYS_V2, TIME_FEAT_KEYS, MU_FEAT_KEYS, ALL_FEAT_KEYS_V2,
ILLUSION_FEAT_KEYS_V2, ILLUSION_ALL_FEAT_KEYS_V2,
)
if feat_keys is None:
if target_forces is not None:
# Illusion scenes: include target force features
feat_keys = ILLUSION_FEAT_KEYS_V2
elif use_time_norm:
# Replace discrete da with time-normalized da/dt
feat_keys = [k for k in CORE_FEAT_KEYS_V2 if not k.startswith("da")]
feat_keys += TIME_FEAT_KEYS
else:
feat_keys = CORE_FEAT_KEYS_V2
if include_mu:
# Add mu features if not already in feat_keys
for mk in MU_FEAT_KEYS:
if mk not in feat_keys:
feat_keys = feat_keys + [mk]
T = sensors.shape[0]
a_prev = np.zeros((T, 3), dtype=np.float64)
a_prev2 = np.zeros((T, 3), dtype=np.float64)
a_prev[1:] = actions_phys[:-1]
a_prev2[2:] = actions_phys[:-2]
dim = compute_dimensionless(sensors, forces, u0=u0, d=20.0)
sym = compute_features(
dim, a_prev, a_prev2, mu,
alpha_mode=alpha_mode,
include_mu=include_mu,
include_cos_sin=include_cos_sin,
dt_c=dt_c,
u0=u0,
target_forces=target_forces,
)
Theta_f = build_feature_matrix(sym, feat_keys, add_bias=False)
Theta_r = build_feature_matrix(sym, feat_keys, add_bias=True)
feat_names_front = get_feature_names(feat_keys, add_bias=False)
feat_names_rear = get_feature_names(feat_keys, add_bias=True)
return (Theta_f[n_warmup:], Theta_r[n_warmup:],
actions_phys[n_warmup:],
feat_names_front, feat_names_rear)
# ---------------------------------------------------------------------------
# Derivative-mode feature builders: fit d(alpha)/dt = g(physics_state)
# No action history in input features; Y is time-normalized action derivative.
# ---------------------------------------------------------------------------
def compute_action_deriv(
actions_phys: np.ndarray, # (T, 3) physical omega
dt_c: float, # control interval in T0 units
u0: float = U0,
center_diff: bool = False, # if True, use (alpha(t) - alpha(t-2))/(2*dt_c)
) -> np.ndarray:
"""Compute time-normalized action derivative d(alpha)/dt.
forward_diff: d(alpha)/dt ~ (alpha(t) - alpha(t-1)) / dt_c
center_diff: d(alpha)/dt ~ (alpha(t) - alpha(t-2)) / (2*dt_c)
Returns (T, 3) array, with first row(s) zero-padded.
"""
alpha = np.asarray(actions_phys, dtype=np.float64) / u0 # non-dim
T = alpha.shape[0]
deriv = np.zeros_like(alpha)
if center_diff and T >= 3:
deriv[2:] = (alpha[2:] - alpha[:-2]) / (2.0 * dt_c)
elif T >= 2:
deriv[1:] = (alpha[1:] - alpha[:-1]) / dt_c
return deriv
def get_feature_matrix_deriv(
sensors: np.ndarray, # (T, 6)
forces: np.ndarray, # (T, 6)
actions_phys: np.ndarray, # (T, 3) physical omega
mu: float,
u0: float = U0,
dt_c: float = 1.0, # control interval in T0 units
feat_keys: Optional[List[str]] = None, # default PHYSICS_FEAT_KEYS
include_mu: bool = False, # add mu modulation features (for cross-Re joint)
target_forces: Optional[np.ndarray] = None, # (T, 2) for Illusion
n_warmup: int = 2,
center_diff: bool = False,
augment_level: int = 0, # 0=static, 1=+lags, 2=+derivs, 3=both, 4=+action_lag
output_mode: str = "deriv", # "deriv": predict d(alpha)/dt; "absolute": predict alpha directly
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, List[str], List[str]]:
"""Build feature matrices for phase-state SINDy.
Input features: physics state with optional temporal context.
Parameters
----------
output_mode : str
"deriv": Y = d(alpha)/dt (time-normalized). Closed-loop needs integration.
"absolute": Y = alpha (non-dimensional action). No integration needed.
augment_level : int
0: PHYSICS_FEAT_KEYS only (static, no memory)
1: + lag-1 obs features (adds temporal context)
2: + obs derivative features (adds rate information)
3: + lag-1 + derivative (full temporal context)
4: + action lag (aF_lag1 etc., for ablation comparison against old approach)
"""
from .feature_builder import (
PHYSICS_FEAT_KEYS, ILLUSION_ERR_KEYS, MU_FEAT_KEYS,
AUG_LEVEL_1_KEYS, AUG_LEVEL_2_KEYS,
AUG_LEVEL_3_KEYS, AUG_LEVEL_4_KEYS,
ILLUSION_AUG_LEVEL_1_KEYS, ILLUSION_AUG_LEVEL_2_KEYS,
ILLUSION_AUG_LEVEL_3_KEYS, ILLUSION_AUG_LEVEL_4_KEYS,
)
# Select feature key set based on augment_level
_AUG_MAP = {0: PHYSICS_FEAT_KEYS, 1: AUG_LEVEL_1_KEYS,
2: AUG_LEVEL_2_KEYS, 3: AUG_LEVEL_3_KEYS, 4: AUG_LEVEL_4_KEYS}
_AUG_ILLUSION_MAP = {0: ILLUSION_ERR_KEYS, 1: ILLUSION_AUG_LEVEL_1_KEYS,
2: ILLUSION_AUG_LEVEL_2_KEYS, 3: ILLUSION_AUG_LEVEL_3_KEYS,
4: ILLUSION_AUG_LEVEL_4_KEYS}
if feat_keys is None:
if target_forces is not None:
feat_keys = _AUG_ILLUSION_MAP.get(augment_level, ILLUSION_AUG_LEVEL_3_KEYS)
else:
feat_keys = _AUG_MAP.get(augment_level, AUG_LEVEL_3_KEYS)
if include_mu:
for mk in MU_FEAT_KEYS:
if mk not in feat_keys:
feat_keys = feat_keys + [mk]
if include_mu:
for mk in MU_FEAT_KEYS:
if mk not in feat_keys:
feat_keys = feat_keys + [mk]
T = sensors.shape[0]
# a_prev/a_prev2 only used for computing the DERIVATIVE target Y, not in Theta
a_prev = np.zeros((T, 3), dtype=np.float64)
a_prev2 = np.zeros((T, 3), dtype=np.float64)
a_prev[1:] = actions_phys[:-1]
a_prev2[2:] = actions_phys[:-2]
dim = compute_dimensionless(sensors, forces, u0=u0, d=20.0)
sym = compute_features(
dim, a_prev, a_prev2, mu,
alpha_mode=False, include_mu=include_mu,
include_cos_sin=False, dt_c=dt_c, u0=u0,
target_forces=target_forces,
sensors_raw=sensors, forces_raw=forces,
)
Theta_f = build_feature_matrix(sym, feat_keys, add_bias=False)
Theta_r = build_feature_matrix(sym, feat_keys, add_bias=True)
feat_names_front = get_feature_names(feat_keys, add_bias=False)
feat_names_rear = get_feature_names(feat_keys, add_bias=True)
# Y: depends on output_mode
if output_mode == "absolute":
Y = np.asarray(actions_phys, dtype=np.float64) / u0 # non-dim alpha (absolute action)
else:
Y = compute_action_deriv(actions_phys, dt_c, u0=u0, center_diff=center_diff)
return (Theta_f[n_warmup:], Theta_r[n_warmup:],
Y[n_warmup:],
feat_names_front, feat_names_rear)

View File

@ -0,0 +1,258 @@
#!/usr/bin/env python3
"""Offline multi-step rollout evaluation for SINDy models.
Evaluates a SINDy model by doing recursive multi-step prediction on
offline PPO data, WITHOUT running CFD. This isolates the model's
rollout behavior from CFD initialization stochasticity.
Usage:
python3 src/SR_analysis/validate/eval_rollout.py \\
--sindy-results src/SR_analysis/sindy/karman/sindy_results_deriv.json \\
--scene karman_re100
# Compare multiple feature presets
python3 src/SR_analysis/validate/eval_rollout.py \\
--scene karman_re100 --preset static --preset phase --preset full_lag
"""
from __future__ import annotations
import argparse
import json
import os
import sys
from typing import List, Optional
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
if _REPO not in sys.path:
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from SR_analysis.utils.feature_builder import (
compute_dimensionless, compute_features, build_feature_matrix,
PHYSICS_FEAT_KEYS, LAG_FEAT_KEYS, DERIV_FEAT_KEYS,
AUG_LEVEL_1_KEYS, PHASE_STATE_KEYS, ILLUSION_ERR_KEYS, ILLUSION_PHASE_KEYS,
)
from SR_analysis.utils.sindy_fitter import compute_action_deriv
from SR_analysis.configs import get_scene
def load_model_data(scene_name: str):
"""Load controlled data and scene config."""
cfg = get_scene(scene_name)
scene_id = cfg["scene_id"]
data_dir = os.path.join(
os.path.dirname(__file__), "..", "data", scene_id, scene_name,
)
npz = np.load(os.path.join(data_dir, "controlled.npz"))
sensors = npz["sensors"].astype(np.float64)
forces = npz["forces"].astype(np.float64)
actions_norm = npz["actions"].astype(np.float64)
scale = cfg["action_scale"]
bias = np.array(cfg["action_bias"], dtype=np.float64)
u0 = cfg["u0"]
actions_phys = (actions_norm * scale + bias) * u0
target_forces = npz["target_forces"].astype(np.float64) if "target_forces" in npz else None
return sensors, forces, actions_phys, target_forces, cfg
def build_sindy_features(sensors, forces, actions_phys, mu, u0, dt_c,
target_forces, feat_keys):
"""Build Theta_f, Theta_r, Y from data using given feat_keys.
Returns (Theta_f, Theta_r, Y_deriv, fn_f, fn_r, a_prev_phys, a_prev2_phys).
"""
T = sensors.shape[0]
a_prev = np.zeros((T, 3), dtype=np.float64)
a_prev2 = np.zeros((T, 3), dtype=np.float64)
a_prev[1:] = actions_phys[:-1]
a_prev2[2:] = actions_phys[:-2]
dim = compute_dimensionless(sensors, forces, u0=u0, d=20.0)
sym = compute_features(
dim, a_prev, a_prev2, mu,
alpha_mode=False, include_mu=False,
include_cos_sin=False, dt_c=dt_c, u0=u0,
target_forces=target_forces,
sensors_raw=sensors, forces_raw=forces,
)
Theta_f = build_feature_matrix(sym, feat_keys, add_bias=False)
Theta_r = build_feature_matrix(sym, feat_keys, add_bias=True)
Y_deriv = compute_action_deriv(actions_phys, dt_c, u0=u0)
return (Theta_f, Theta_r, Y_deriv,
feat_keys, ["bias"] + feat_keys)
def compute_rollout(
sensors, forces, actions_phys, target_forces, cfg,
front_coef, top_coef, bottom_coef,
feat_keys_front, feat_keys_rear,
max_steps: int = 50,
center_diff: bool = False,
) -> dict:
"""Compute offline multi-step rollout error."""
mu = cfg["mu"]
u0 = cfg["u0"]
dt_c = cfg.get("sample_interval", 800) / 2000.0
T = sensors.shape[0]
# True PPO alpha(t)
alpha_true = actions_phys / u0
# SINDy predicted alpha(t) via recursive rollout
alpha_pred = np.zeros_like(alpha_true)
# Step 0, 1: initialize from PPO data (first 2 steps use true values)
alpha_pred[0] = alpha_true[0]
alpha_pred[1] = alpha_true[1]
# Build feature history for augmented features
fbs = {} # feature build state: sensors_raw, forces_raw
for t in range(2, min(T, max_steps + 2)):
# Build current obs: need sensor/force from the actual data
obs_slice = np.hstack([sensors[t], forces[t]])
# For augmented features, also need obs[t-1]
obs_prev = np.hstack([sensors[t-1], forces[t-1]])
# Check if augmented features needed
needs_aug = any(k.endswith("_dt") or k.endswith("_lag1") for k in feat_keys_front)
if needs_aug:
sensors_raw = np.vstack([sensors[t-1:t+1]])
forces_raw = np.vstack([forces[t-1:t+1]])
sensors_curr = sensors[t:t+1]
forces_curr = forces[t:t+1]
else:
sensors_raw = sensors[t:t+1]
forces_raw = forces[t:t+1]
sensors_curr = sensors[t:t+1]
forces_curr = forces[t:t+1]
# Build features for current step
a_prev = alpha_pred[t-1:t] * u0 # use PREDICTED a(t-1)
a_prev2 = alpha_pred[t-2:t-1] * u0 # use PREDICTED a(t-2)
dim = compute_dimensionless(sensors_curr, forces_curr, u0=u0, d=20.0)
tf = target_forces[t:t+1] if target_forces is not None else None
sym = compute_features(dim, a_prev, a_prev2, mu,
alpha_mode=False, include_mu=False,
include_cos_sin=False, dt_c=dt_c, u0=u0,
target_forces=tf,
sensors_raw=sensors_raw, forces_raw=forces_raw)
# Front: no bias
fv_f = build_feature_matrix(sym, feat_keys_front, add_bias=False)[0]
d_alphaF = float(np.dot(fv_f, front_coef))
# Top: with bias
fv_r = build_feature_matrix(sym, feat_keys_rear, add_bias=True)[0]
d_alphaT = float(np.dot(fv_r, top_coef))
# Bottom: from fit (independent since we're rolling out)
d_alphaB = float(np.dot(fv_r, bottom_coef))
# Integrate
alpha_pred[t] = alpha_pred[t-1] + dt_c * np.array([d_alphaF, d_alphaB, d_alphaT])
# Compute errors at various rollout horizons
n = min(T, max_steps + 2)
# Use full-range alpha range for all horizons (stable denominator)
full_alpha_range = alpha_true[2:n].max(axis=0) - alpha_true[2:n].min(axis=0)
results = {}
for horizon in [1, 5, 20, 50]:
if horizon + 2 > n:
horizon = n - 2
if horizon <= 0:
continue
err = alpha_pred[2:2+horizon] - alpha_true[2:2+horizon]
rmse_per_channel = np.sqrt(np.mean(err**2, axis=0))
rel_err = rmse_per_channel / (full_alpha_range + 1e-12)
results[f"{horizon}step_rmse"] = rmse_per_channel.tolist()
results[f"{horizon}step_rel"] = rel_err.tolist()
results["alpha_range"] = full_alpha_range.tolist()
results["n_samples"] = n - 2
return results
def main():
ap = argparse.ArgumentParser(description="Offline SINDy rollout evaluation")
ap.add_argument("--sindy-results", type=str, required=True,
help="Path to sindy_results JSON")
ap.add_argument("--scene", type=str, required=True,
help="Scene name (e.g. karman_re100)")
ap.add_argument("--threshold", type=float, default=None,
help="SINDy threshold (default: best R2)")
args = ap.parse_args()
# Load SINDy model
with open(args.sindy_results) as f:
data = json.load(f)
per = data.get("per_scene", data).get(args.scene)
if per is None:
# Try top-level (joint model format)
per = data
fn_f = per.get("feature_names_front", [])
fn_r = per.get("feature_names_rear", [])
def _coef(ch_name):
ch = per.get(ch_name, {})
if args.threshold is not None:
for r in ch.get("results", []):
if abs(r["threshold"] - args.threshold) < 1e-6:
return np.array(r.get("coef", ch["best_coef"]), dtype=np.float64)
return np.array(ch["best_coef"], dtype=np.float64)
front_coef = _coef("front")[:len(fn_f)]
top_coef = _coef("top")[:len(fn_r)]
bottom_coef = _coef("bottom")[:len(fn_r)]
feat_keys_front = [k for k in fn_f if k != "bias"]
feat_keys_rear = [k for k in fn_r if k != "bias"]
sensors, forces, actions_phys, target_forces, cfg = load_model_data(args.scene)
print(f"Rollout eval: {args.scene}")
print(f" Features: front {len(feat_keys_front)}, rear {len(feat_keys_rear)}")
print(f" Front keys: {feat_keys_front}")
print()
results = compute_rollout(
sensors, forces, actions_phys, target_forces, cfg,
front_coef, top_coef, bottom_coef,
feat_keys_front, feat_keys_rear,
)
print(f" alpha_range = [{results['alpha_range'][0]:.3f}, {results['alpha_range'][1]:.3f}, {results['alpha_range'][2]:.3f}]")
print()
print(f" {'Horizon':>8s} {'Front RMSE':>12s} {'Bottom RMSE':>12s} {'Top RMSE':>12s} {'Front rel':>10s} {'Bottom rel':>10s} {'Top rel':>10s}")
print(f" {'-'*8} {'-'*12} {'-'*12} {'-'*12} {'-'*10} {'-'*10} {'-'*10}")
for h in [1, 5, 20, 50]:
rmse_key = f"{h}step_rmse"
rel_key = f"{h}step_rel"
if rmse_key in results:
rmse = results[rmse_key]
rel = results[rel_key]
print(f" {h:>8d} {rmse[0]:>12.4f} {rmse[1]:>12.4f} {rmse[2]:>12.4f} {rel[0]:>10.3f} {rel[1]:>10.3f} {rel[2]:>10.3f}")
print()
print("Interpretation:")
print(" 1-step = pure one-step prediction quality")
print(" 5-step = short rollout (control propagation ~1/40 of domain)")
print(" 20-step = medium rollout (~1/10 of domain)")
print(" 50-step = long rollout (~1/4 of domain)")
print(f" Relative error = RMSE / alpha_range (lower is better)")
print(f" If 5-step error >> 1-step error: model has rollout instability")
if __name__ == "__main__":
main()

View File

@ -40,7 +40,7 @@ from SR_analysis.utils.cfd_interface import (
from SR_analysis.utils.sindy_fitter import get_feature_matrix_from_data from SR_analysis.utils.sindy_fitter import get_feature_matrix_from_data
from SR_analysis.utils.feature_builder import ( from SR_analysis.utils.feature_builder import (
compute_dimensionless, compute_features, build_feature_matrix, compute_dimensionless, compute_features, build_feature_matrix,
apply_G_x, ALL_FEAT_KEYS, U0, apply_G_x, ALL_FEAT_KEYS, ALL_FEAT_KEYS_V2, CORE_FEAT_KEYS_V2, U0,
) )
from SR_analysis.utils.g_operator import apply_G_raw from SR_analysis.utils.g_operator import apply_G_raw
from SR_analysis.configs import ( from SR_analysis.configs import (
@ -57,10 +57,18 @@ def build_feature_vector(
mu: float, mu: float,
u0: float, u0: float,
add_bias: bool, add_bias: bool,
feat_keys: Optional[List[str]] = None,
use_v2: bool = False,
target_forces: Optional[np.ndarray] = None,
) -> np.ndarray: ) -> np.ndarray:
"""Build a single-row feature vector from raw obs and action state. """Build a single-row feature vector from raw obs and action state.
Matches the feature_builder logic but for a single time step. Matches the feature_builder logic but for a single time step.
Args:
feat_keys: custom feature keys. If None, uses ALL_FEAT_KEYS_V2 or ALL_FEAT_KEYS.
use_v2: if True, use CORE_FEAT_KEYS_V2 (no sin/cos). Ignored if feat_keys given.
target_forces: (2,) raw lattice target forces [fx,fy], for Illusion scenes.
""" """
sensors = obs_slice[0:6].astype(np.float64).reshape(1, 6) sensors = obs_slice[0:6].astype(np.float64).reshape(1, 6)
forces = obs_slice[6:12].astype(np.float64).reshape(1, 6) forces = obs_slice[6:12].astype(np.float64).reshape(1, 6)
@ -68,9 +76,22 @@ def build_feature_vector(
ap2 = a_prev2.astype(np.float64).reshape(1, 3) ap2 = a_prev2.astype(np.float64).reshape(1, 3)
dim = compute_dimensionless(sensors, forces, u0=u0, d=20.0) dim = compute_dimensionless(sensors, forces, u0=u0, d=20.0)
sym = compute_features(dim, ap, ap2, mu, alpha_mode=False, include_mu=True, u0=u0) # Detect whether to include mu terms based on feat_keys
has_mu = feat_keys is not None and "mu" in feat_keys
has_cos_sin = feat_keys is not None and "sin_ua" in feat_keys
tf = target_forces.reshape(1, 2) if target_forces is not None else None
sym = compute_features(dim, ap, ap2, mu, alpha_mode=False,
include_mu=has_mu,
include_cos_sin=has_cos_sin,
u0=u0, target_forces=tf)
feat = build_feature_matrix(sym, ALL_FEAT_KEYS, add_bias=add_bias) if feat_keys is None:
if use_v2:
feat_keys = CORE_FEAT_KEYS_V2
else:
feat_keys = ALL_FEAT_KEYS
feat = build_feature_matrix(sym, feat_keys, add_bias=add_bias)
return feat[0] # single row return feat[0] # single row
@ -84,26 +105,42 @@ def predict_v23(
top_coef: np.ndarray, top_coef: np.ndarray,
feat_names_front: List[str], feat_names_front: List[str],
feat_names_rear: List[str], feat_names_rear: List[str],
feat_keys_front: Optional[List[str]] = None,
feat_keys_rear: Optional[List[str]] = None,
target_forces: Optional[np.ndarray] = None,
) -> np.ndarray: ) -> np.ndarray:
"""Predict actions using v23: front no-bias + rear shared-head. """Predict actions using v23: front no-bias + rear shared-head.
Returns (3,) physical omega array: [front, bottom, top]. Returns (3,) physical omega array: [front, bottom, top].
""" """
# Front channel: no bias # Front channel: no bias
front = float(np.dot( fv_front = build_feature_vector(
build_feature_vector(obs_slice, a_prev, a_prev2, mu, u0, add_bias=False), obs_slice, a_prev, a_prev2, mu, u0, add_bias=False,
front_coef)) feat_keys=feat_keys_front,
target_forces=target_forces,
)
front = float(np.dot(fv_front, front_coef))
# Top channel: with bias # Top channel: with bias
top = float(np.dot( fv_top = build_feature_vector(
build_feature_vector(obs_slice, a_prev, a_prev2, mu, u0, add_bias=True), obs_slice, a_prev, a_prev2, mu, u0, add_bias=True,
top_coef)) feat_keys=feat_keys_rear,
target_forces=target_forces,
)
top = float(np.dot(fv_top, top_coef))
# Bottom = -top(Gx) using shared-head # Bottom = -top(Gx) using shared-head
G_obs, G_a_prev, G_a_prev2 = apply_G_raw(obs_slice, a_prev, a_prev2) G_obs, G_a_prev, G_a_prev2 = apply_G_raw(obs_slice, a_prev, a_prev2)
bottom = -float(np.dot( # G transforms forces: swap bottom<->top, negate fy. Target force G transform:
build_feature_vector(G_obs, G_a_prev, G_a_prev2, mu, u0, add_bias=True), # target_Cd is a drag-like quantity that is even under G (stays),
top_coef)) # target_Cl is a lift-like quantity that is odd under G (negates).
G_target = np.array([target_forces[0], -target_forces[1]]) if target_forces is not None else None
fv_bot = build_feature_vector(
G_obs, G_a_prev, G_a_prev2, mu, u0, add_bias=True,
feat_keys=feat_keys_rear,
target_forces=G_target,
)
bottom = -float(np.dot(fv_bot, top_coef))
return np.array([front, bottom, top], dtype=np.float64) return np.array([front, bottom, top], dtype=np.float64)
@ -118,20 +155,125 @@ def predict_unstructured(
bottom_coef: np.ndarray, bottom_coef: np.ndarray,
top_coef: np.ndarray, top_coef: np.ndarray,
feat_names: List[str], feat_names: List[str],
feat_keys: Optional[List[str]] = None,
) -> np.ndarray: ) -> np.ndarray:
"""Predict actions using unstructured: each channel independent with bias.""" """Predict actions using unstructured: each channel independent with bias."""
feat = build_feature_vector(obs_slice, a_prev, a_prev2, mu, u0, add_bias=True) feat = build_feature_vector(
obs_slice, a_prev, a_prev2, mu, u0, add_bias=True,
feat_keys=feat_keys,
)
front = float(np.dot(feat, front_coef)) front = float(np.dot(feat, front_coef))
bottom = float(np.dot(feat, bottom_coef)) bottom = float(np.dot(feat, bottom_coef))
bottom = float(np.dot(feat, bottom_coef))
top = float(np.dot(feat, top_coef)) top = float(np.dot(feat, top_coef))
return np.array([front, bottom, top], dtype=np.float64) return np.array([front, bottom, top], dtype=np.float64)
def load_sindy_coefs(sindy_path: str, scene_name: str) -> Dict[str, Any]: def predict_v23_deriv(
obs_slice: np.ndarray,
obs_prev: np.ndarray,
a_prev_phys: np.ndarray,
a_prev2_phys: np.ndarray,
mu: float,
u0: float,
dt_c: float,
front_coef: np.ndarray,
top_coef: np.ndarray,
feat_keys_front: Optional[List[str]] = None,
feat_keys_rear: Optional[List[str]] = None,
target_forces: Optional[np.ndarray] = None,
output_mode: str = "deriv",
) -> np.ndarray:
"""Predict action using phase-state features, v23 structure.
Two output modes:
"deriv": Predict d(alpha)/dt, then integrate:
omega_new = (alpha_prev + dt_c * d_alpha/dt_pred) * u0
"absolute": Predict alpha directly:
omega_new = alpha_pred * u0
Returns (3,) physical omega array: [front, bottom, top].
"""
from SR_analysis.utils.feature_builder import PHYSICS_FEAT_KEYS
if feat_keys_front is None:
feat_keys_front = PHYSICS_FEAT_KEYS
if feat_keys_rear is None:
feat_keys_rear = PHYSICS_FEAT_KEYS
# Check if augmented features (derivatives or lags) are needed
needs_aug = any(k.endswith("_dt") or k.endswith("_lag1") for k in feat_keys_front)
# Build sensor/force arrays: if augmented features needed, stack prev+current
if needs_aug:
sensors_raw = np.vstack([obs_prev[0:6], obs_slice[0:6]]).astype(np.float64)
forces_raw = np.vstack([obs_prev[6:12], obs_slice[6:12]]).astype(np.float64)
# Current dim uses only the last row
sensors_curr = obs_slice[0:6].astype(np.float64).reshape(1, 6)
forces_curr = obs_slice[6:12].astype(np.float64).reshape(1, 6)
else:
sensors_raw = sensors_curr = obs_slice[0:6].astype(np.float64).reshape(1, 6)
forces_raw = forces_curr = obs_slice[6:12].astype(np.float64).reshape(1, 6)
ap = a_prev_phys.astype(np.float64).reshape(1, 3)
ap2 = a_prev2_phys.astype(np.float64).reshape(1, 3)
dim = compute_dimensionless(sensors_curr, forces_curr, u0=u0, d=20.0)
has_mu = any("mu" in k for k in feat_keys_front)
tf = target_forces.reshape(1, 2) if target_forces is not None else None
sym = compute_features(dim, ap, ap2, mu, alpha_mode=False,
include_mu=has_mu,
include_cos_sin=False,
u0=u0, target_forces=tf,
sensors_raw=sensors_raw, forces_raw=forces_raw)
# Front: no bias, predict d(alpha_F)/dt
fv_front = build_feature_matrix(sym, feat_keys_front, add_bias=False)[0]
pred_F = float(np.dot(fv_front, front_coef))
# Top: with bias, predict d(alpha_T)/dt
fv_top = build_feature_matrix(sym, feat_keys_rear, add_bias=True)[0]
pred_T = float(np.dot(fv_top, top_coef))
# Bottom = -top(Gx) using shared-head
G_obs, G_a_prev, G_a_prev2 = apply_G_raw(obs_slice, a_prev_phys, a_prev2_phys)
G_target = np.array([target_forces[0], -target_forces[1]]) if target_forces is not None else None
G_sensors = G_obs[0:6].astype(np.float64).reshape(1, 6)
G_forces = G_obs[6:12].astype(np.float64).reshape(1, 6)
G_dim = compute_dimensionless(G_sensors, G_forces, u0=u0, d=20.0)
G_sym = compute_features(G_dim, G_a_prev.reshape(1, 3), G_a_prev2.reshape(1, 3),
mu, alpha_mode=False, include_mu=has_mu,
include_cos_sin=False, u0=u0, target_forces=G_target)
fv_bot = build_feature_matrix(G_sym, feat_keys_rear, add_bias=True)[0]
pred_B = -float(np.dot(fv_bot, top_coef))
# Convert to omega based on output_mode
if output_mode == "absolute":
# Direct prediction: pred = alpha, omega = alpha * u0
alpha_new = np.array([pred_F, pred_B, pred_T])
else:
# Derivative + integrate: pred = d(alpha)/dt
# omega = (alpha_prev + dt_c * d_alpha/dt) * u0
alpha_prev = a_prev_phys / u0
alpha_new = alpha_prev + dt_c * np.array([pred_F, pred_B, pred_T])
omega_new = alpha_new * u0
return omega_new.astype(np.float64)
def load_sindy_coefs(sindy_path: str, scene_name: str,
threshold: Optional[float] = None) -> Dict[str, Any]:
"""Load SINDy coefficients for a scene from results JSON. """Load SINDy coefficients for a scene from results JSON.
If threshold is given (e.g. 0.007), uses coefficients at that threshold
instead of the best_R2 ones.
Auto-detects v1 (21 features, sin_ua/cos_ua, mu) vs v2 (14 features, no sin/cos)
feature sets based on the length of feature_names_front.
Returns dict with keys: front_coef, top_coef, bottom_coef, Returns dict with keys: front_coef, top_coef, bottom_coef,
feat_names_front, feat_names_rear, front_bias_mode. feat_names_front, feat_names_rear, front_bias_mode,
feat_keys_front, feat_keys_rear, is_v2.
""" """
with open(sindy_path) as f: with open(sindy_path) as f:
data = json.load(f) data = json.load(f)
@ -143,20 +285,55 @@ def load_sindy_coefs(sindy_path: str, scene_name: str) -> Dict[str, Any]:
fn_f = per["feature_names_front"] fn_f = per["feature_names_front"]
fn_r = per["feature_names_rear"] fn_r = per["feature_names_rear"]
front_coef = np.array(per["front"]["best_coef"], dtype=np.float64) def _coef_at_threshold(ch_data, th):
top_coef = np.array(per["top"]["best_coef"], dtype=np.float64) """Get coef at given threshold, or best if th=None."""
bottom_coef = np.array(per["bottom"]["best_coef"], dtype=np.float64) if th is None:
return np.array(ch_data["best_coef"], dtype=np.float64)
for r in ch_data.get("results", []):
if abs(r["threshold"] - th) < 1e-6:
return np.array(r.get("coef", ch_data["best_coef"]), dtype=np.float64)
# Fallback: best
return np.array(ch_data["best_coef"], dtype=np.float64)
# Detect if front was fitted with bias (fn_f has "bias") or without front_coef = _coef_at_threshold(per["front"], threshold)
top_coef = _coef_at_threshold(per["top"], threshold)
bottom_coef = _coef_at_threshold(per["bottom"], threshold)
# Trim coefficients to match feature names
front_coef = front_coef[:len(fn_f)]
top_coef = top_coef[:len(fn_r)]
bottom_coef = bottom_coef[:len(fn_r)]
# Detect v1 vs v2: v2 has 14 features (no sin_ua/cos_ua), v1 has 16
# v2 joint has 19 (14 core + 5 mu: mu, mu_u_a, mu_v_a, mu_Cd_tot, mu_Cl_diff)
# v1 joint has 21 (16 core + 5 mu)
# Illusion v2 has 16 (14 core + 2 target: target_Cd, target_Cl)
# Illusion v2 + mu has 21 (14 core + 2 target + 5 mu)
# Deriv mode: PHYSICS_FEAT_KEYS has 8, ILLUSION_ERR_KEYS has 10
is_v2 = len(fn_f) in [8, 10, 14, 16, 19, 21]
# Strip "bias" from feature keys (build_feature_matrix adds it via add_bias flag)
feat_keys_front = [k for k in fn_f if k != "bias"]
feat_keys_rear = [k for k in fn_r if k != "bias"]
# Detect if front was fitted with bias or without
front_has_bias = fn_f[0] == "bias" if len(fn_f) > 0 else False front_has_bias = fn_f[0] == "bias" if len(fn_f) > 0 else False
# Extract mode from per_scene metadata (default: "v23" for old-style results)
sindy_mode = per.get("mode", "v23")
return { return {
"front_coef": front_coef, "front_coef": front_coef,
"top_coef": top_coef, "top_coef": top_coef,
"bottom_coef": bottom_coef, "bottom_coef": bottom_coef,
"feat_names_front": fn_f, "feat_names_front": fn_f,
"feat_names_rear": fn_r, "feat_names_rear": fn_r,
"feat_keys_front": feat_keys_front,
"feat_keys_rear": feat_keys_rear,
"front_has_bias": front_has_bias, "front_has_bias": front_has_bias,
"is_v2": is_v2,
"mode": sindy_mode,
"threshold_used": threshold if threshold is not None else "best",
} }
@ -164,20 +341,14 @@ def run_validation(
scene_name: str, scene_name: str,
coefs: Dict[str, Any], coefs: Dict[str, Any],
device_id: int, device_id: int,
n_steps: int = 100, n_steps: int = 0,
mode: str = "v23", mode: str = "v23",
) -> dict: ) -> dict:
"""Run closed-loop validation using a SINDy control law. """Run closed-loop validation using a SINDy control law.
Parameters Parameters
---------- ----------
scene_name : e.g. "karman_re70" n_steps : int, default=0 (auto: >= NX/U0 / sample_interval)
coefs : dict from load_sindy_coefs()
device_id : GPU device
n_steps : number of closed-loop steps
mode : "v23" (rear shared-head) or "unstructured"
Returns dict with similarity, actions range, etc.
""" """
cfg = get_scene(scene_name) cfg = get_scene(scene_name)
re_code = cfg["re_code"] re_code = cfg["re_code"]
@ -189,8 +360,17 @@ def run_validation(
action_scale = cfg["action_scale"] action_scale = cfg["action_scale"]
action_bias = cfg["action_bias"] action_bias = cfg["action_bias"]
n_obj_total = cfg["n_objects_env"] n_obj_total = cfg["n_objects_env"]
t0_steps = 2000
dt_c = sample_interval / t0_steps
print(f"\n=== Validating {scene_name} (mode={mode}, device={device_id}) ===") # Auto-set steps: control must propagate at least 1*NX/U0 LBM steps
NX = 1280
min_steps = int(1 * NX / u0 / sample_interval)
if n_steps == 0 or n_steps < min_steps:
n_steps = max(min_steps, 200)
print(f" auto-set steps={n_steps} (min_steps={min_steps})")
print(f"\n=== Validating {scene_name} (mode={mode}, device={device_id}, steps={n_steps}) ===")
# Build environment # Build environment
cuda_cfg, field_cfg = load_legacy_configs(LEGACY_CFG_DIR) cuda_cfg, field_cfg = load_legacy_configs(LEGACY_CFG_DIR)
@ -211,58 +391,76 @@ def run_validation(
obs_slice_start=cfg["obs_slice"][0], obs_slice_end=cfg["obs_slice"][1], obs_slice_start=cfg["obs_slice"][0], obs_slice_end=cfg["obs_slice"][1],
) )
# Reset to checkpoint # Now ff is at save_ddf(2) = post-bias state (add_pinball saves after bias rollout).
ff.restore_ddf() # Norm dict contains save_states from the bias FIFO.
ff.apply_ddf()
# Bias FIFO init
fifo = deque(maxlen=FIFO_LEN) fifo = deque(maxlen=FIFO_LEN)
bias_arr = scale_action(np.zeros(3, dtype=np.float32), scale=action_scale, for row in norm["save_states"]:
bias=action_bias, u0=u0, n_total_bodies=n_obj_total) fifo.append(row)
for _ in range(FIFO_LEN):
ff.run(sample_interval, bias_arr)
fifo.append(ff.obs.copy()[2:14])
# Closed-loop with SINDy law # Closed-loop with SINDy law
sens_list = [] sens_list = []
actions_list = [] actions_list = []
a_prev = action_to_physical(np.zeros((1, 3), dtype=np.float32), # For Karman: bias action = DRL action_bias = [0, -4, 4].
scale=action_scale, bias=action_bias, u0=u0).flatten() # Normalized zero action [0,0,0] decodes to the same bias.
a_prev = np.zeros(3, dtype=np.float64)
a_prev2 = a_prev.copy() a_prev2 = a_prev.copy()
for _ in range(n_steps): for _ in range(n_steps):
obs = fifo[-1] if fifo else np.zeros(12, dtype=np.float32) obs = fifo[-1] if fifo else np.zeros(12, dtype=np.float32)
obs_prev = list(fifo)[-2] if len(fifo) >= 2 else obs # for augmented features
# Predict using SINDy law # Convert a_prev from normalized to physical omega for feature builder
# (SINDy was trained on physical omega Y = actions_phys)
a_prev_phys = (a_prev * action_scale + np.array(action_bias, dtype=np.float64)) * u0
a_prev2_phys = (a_prev2 * action_scale + np.array(action_bias, dtype=np.float64)) * u0
# Predict using SINDy law -- returns physical omega
if mode == "v23": if mode == "v23":
omega = predict_v23( omega_phys = predict_v23(
obs, a_prev, a_prev2, mu, u0, obs, a_prev_phys, a_prev2_phys, mu, u0,
coefs["front_coef"], coefs["top_coef"], coefs["front_coef"], coefs["top_coef"],
coefs["feat_names_front"], coefs["feat_names_rear"]) coefs["feat_names_front"], coefs["feat_names_rear"],
feat_keys_front=coefs["feat_keys_front"],
feat_keys_rear=coefs["feat_keys_rear"])
elif mode == "deriv":
omega_phys = predict_v23_deriv(
obs, obs_prev, a_prev_phys, a_prev2_phys, mu, u0, dt_c,
coefs["front_coef"], coefs["top_coef"],
feat_keys_front=coefs["feat_keys_front"],
feat_keys_rear=coefs["feat_keys_rear"],
output_mode="deriv")
elif mode == "abs":
omega_phys = predict_v23_deriv(
obs, obs_prev, a_prev_phys, a_prev2_phys, mu, u0, dt_c,
coefs["front_coef"], coefs["top_coef"],
feat_keys_front=coefs["feat_keys_front"],
feat_keys_rear=coefs["feat_keys_rear"],
output_mode="absolute")
elif mode == "unstructured": elif mode == "unstructured":
omega = predict_unstructured( omega_phys = predict_unstructured(
obs, a_prev, a_prev2, mu, u0, obs, a_prev_phys, a_prev2_phys, mu, u0,
coefs["front_coef"], coefs["bottom_coef"], coefs["top_coef"], coefs["front_coef"], coefs["bottom_coef"], coefs["top_coef"],
coefs["feat_names_rear"]) coefs["feat_names_rear"],
feat_keys=coefs["feat_keys_rear"])
else: else:
raise ValueError(f"Unknown mode: {mode}") raise ValueError(f"Unknown mode: {mode}")
# Clip to valid action range # Convert physical omega -> normalized action -> legacy array
norm_a = (omega / u0 - np.array(action_bias, dtype=np.float64)) / action_scale norm_a = (omega_phys / u0 - np.array(action_bias, dtype=np.float64)) / action_scale
norm_a = np.clip(norm_a, -1.0, 1.0).astype(np.float32) norm_a = np.clip(norm_a, -1.0, 1.0).astype(np.float32)
# Apply to CFD
action_arr = scale_action(norm_a, scale=action_scale, bias=action_bias, action_arr = scale_action(norm_a, scale=action_scale, bias=action_bias,
u0=u0, n_total_bodies=n_obj_total) u0=u0, n_total_bodies=n_obj_total)
ff.run(sample_interval, action_arr) ff.run(sample_interval, action_arr)
obs_new = ff.obs.copy()[2:14] obs_new = ff.obs.copy()[2:14]
fifo.append(obs_new) fifo.append(obs_new)
sens_list.append(obs_new[0:6]) sens_list.append(obs_new[0:6])
actions_list.append(omega.copy()) actions_list.append(omega_phys.copy())
# Keep a_prev in normalized form for next iteration's memory features
a_prev2 = a_prev.copy() a_prev2 = a_prev.copy()
a_prev = omega.copy() a_prev = norm_a.copy().astype(np.float64)
# Evaluate # Evaluate
sens_arr = np.array(sens_list, dtype=np.float32) sens_arr = np.array(sens_list, dtype=np.float32)
@ -289,28 +487,32 @@ def main():
ap.add_argument("--device", type=int, default=2, help="GPU device") ap.add_argument("--device", type=int, default=2, help="GPU device")
ap.add_argument("--steps", type=int, default=100) ap.add_argument("--steps", type=int, default=100)
ap.add_argument("--mode", type=str, default="v23", ap.add_argument("--mode", type=str, default="v23",
choices=["v23", "unstructured"]) choices=["v23", "unstructured", "deriv", "abs"],
help="Control law mode: v23 (default), unstructured, deriv, or abs")
ap.add_argument("--sindy-results", type=str, default=None, ap.add_argument("--sindy-results", type=str, default=None,
help="Path to sindy_results.json") help="Path to sindy_results.json")
ap.add_argument("--threshold", type=float, default=None,
help="SINDy threshold for sparsity (default: best_R2)")
ap.add_argument("--out", type=str, default=None, ap.add_argument("--out", type=str, default=None,
help="Output directory for result JSON") help="Output directory for result JSON")
args = ap.parse_args() args = ap.parse_args()
if args.sindy_results is None: if args.sindy_results is None:
args.sindy_results = os.path.join( args.sindy_results = os.path.join(
os.path.dirname(__file__), "..", "sindy", "karman", "sindy_results.json") os.path.dirname(__file__), "..", "sindy", "karman", "sindy_results_v2.json")
coefs = load_sindy_coefs(args.sindy_results, args.scene) coefs = load_sindy_coefs(args.sindy_results, args.scene, threshold=args.threshold)
result = run_validation(args.scene, coefs, args.device, result = run_validation(args.scene, coefs, args.device,
n_steps=args.steps, mode=args.mode) n_steps=args.steps, mode=args.mode)
th_str = f"_th{args.threshold}" if args.threshold is not None else ""
if args.out is None: if args.out is None:
out_dir = os.path.join(os.path.dirname(__file__), "..", out_dir = os.path.join(os.path.dirname(__file__), "..",
"validate", "results") "validate", "results")
else: else:
out_dir = args.out out_dir = args.out
os.makedirs(out_dir, exist_ok=True) os.makedirs(out_dir, exist_ok=True)
out_path = os.path.join(out_dir, f"{args.scene}_{args.mode}.json") out_path = os.path.join(out_dir, f"{args.scene}_{args.mode}{th_str}.json")
with open(out_path, "w") as f: with open(out_path, "w") as f:
json.dump(result, f, indent=2) json.dump(result, f, indent=2)
print(f"Saved: {out_path}") print(f"Saved: {out_path}")

View File

@ -0,0 +1,329 @@
"""Closed-loop validator for Illusion scenes.
Builds an illusion environment (no disturbance cylinder, 6 objects, obs[0:12]),
loads target data from existing target.npz, and validates SINDy control laws
in closed-loop CFD with v23 mode (front no-bias + rear shared-head).
Usage (from repo root):
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop_illusion.py \\
--scene illusion_1L --device 0 --steps 200
Default steps = max(320, 2*NX/U0/sample_interval) to ensure control propagates.
"""
from __future__ import annotations
import argparse
import json
import os
import sys
from collections import deque
from typing import Any, Dict, Optional
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
if _REPO not in sys.path:
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from LegacyCelerisLab import FlowField # noqa: E402
from SR_analysis.utils.cfd_interface import (
load_legacy_configs, compute_similarity, scale_action,
vorticity_from_ddf, save_vorticity_png,
)
from SR_analysis.configs import get_scene, LEGACY_CFG_DIR, FIFO_LEN
from SR_analysis.validate.run_closed_loop import (
predict_v23, predict_v23_deriv, load_sindy_coefs,
)
# Import target force harmonic reconstruction from inference pipeline
from SR_analysis.scripts.infer_illusion import gen_target_states_at
DATA_TYPE = np.float32
NX = 1280
def run_validation_illusion(
scene_name: str,
sindy_path: str,
device_id: int,
n_steps: int = 0,
threshold: Optional[float] = None,
out_dir: Optional[str] = None,
) -> dict:
"""Run closed-loop validation for an Illusion scene.
Parameters
----------
n_steps : int, default=0 (auto: >= 2*NX/U0 / sample_interval)
"""
cfg = get_scene(scene_name)
u0 = cfg["u0"]
mu = cfg["mu"]
l0 = 20.0
sample_interval = cfg["sample_interval"]
conv_len = cfg.get("conv_len", 36)
action_scale = cfg["action_scale"]
action_bias = cfg["action_bias"]
n_obj_total = cfg["n_objects_env"]
sensor_x = cfg["sensor_x"]
front_x = cfg["pinball_front_x"]
rear_x = cfg["pinball_rear_x"]
# Auto-set steps: control must propagate at least 1*NX/U0 LBM steps
min_steps = int(1 * NX / u0 / sample_interval)
if n_steps == 0 or n_steps < min_steps:
n_steps = max(min_steps, 200)
print(f" auto-set steps={n_steps} (min_steps={min_steps})")
print(f"\n=== Validating {scene_name} (device={device_id}, steps={n_steps}) ===")
if out_dir is None:
out_dir = os.path.join(os.path.dirname(__file__), "results")
os.makedirs(out_dir, exist_ok=True)
# Load SINDy coefficients
coefs = load_sindy_coefs(sindy_path, scene_name, threshold=threshold)
threshold_str = f"_th{threshold}" if threshold is not None else ""
# Load target data (sensors for similarity, harmonics for target force reconstruction)
data_dir = os.path.join(
os.path.dirname(__file__), "..", "data", "illusion", scene_name,
)
target_npz = np.load(os.path.join(data_dir, "target.npz"))
if "target_sensors" in target_npz:
target_states = target_npz["target_sensors"]
else:
target_states = target_npz["target_states"][:, 2:8]
print(f" target loaded: {target_states.shape}")
# Load target harmonics for real-time target force reconstruction
target_harmonics = None
harm_path = os.path.join(data_dir, "target_harmonics.json")
if os.path.isfile(harm_path):
with open(harm_path) as f:
target_harmonics = json.load(f)
print(f" target harmonics loaded: {len(target_harmonics)} channels")
# Build environment
cuda_cfg, field_cfg = load_legacy_configs(LEGACY_CFG_DIR)
field_cfg = field_cfg._replace(viscosity=float(cfg["nu"]))
ff = FlowField(field_cfg, cuda_cfg, device_id=device_id)
ny = ff.FIELD_SHAPE[1]
# Add 3 sensors, then 3 pinball cylinders (illusion positions)
for y_off in [2.0, 0.0, -2.0]:
sc = (sensor_x * l0, (ny - 1) / 2 + y_off * l0, 0.0)
ff.add_sensor(sc, l0 / 4.0)
ff.add_cylinder((front_x * l0, (ny - 1) / 2, 0.0), l0 / 2.0)
ff.add_cylinder((rear_x * l0, (ny - 1) / 2 + 0.75 * l0, 0.0), l0 / 2.0)
ff.add_cylinder((rear_x * l0, (ny - 1) / 2 - 0.75 * l0, 0.0), l0 / 2.0)
n_obj = ff.obs.size // 2
assert n_obj == 6, f"Expected 6 objects, got {n_obj}"
# Stabilize with zero action (4*NX/U0 = 512000 steps)
stabilize_steps = int(4 * NX / u0)
print(f" stabilising pinball ({stabilize_steps} steps)...")
ff.run(stabilize_steps, np.zeros(n_obj, dtype=DATA_TYPE))
# --- save_ddf(1): post-stabilization checkpoint ---
ff.get_ddf()
ff.save_ddf()
# Norm collection (zero action)
print(f" norm collection ({FIFO_LEN} steps)...")
fifo = deque(maxlen=FIFO_LEN)
for _ in range(FIFO_LEN):
ff.run(sample_interval, np.zeros(n_obj, dtype=DATA_TYPE))
fifo.append(ff.obs.copy()[0:12])
temp_states = np.array(fifo, dtype=DATA_TYPE)
force_norm_fact = 6.0 * float(np.max(np.abs(temp_states[:, 6:12])))
sens_deviation = np.mean(temp_states[:, 0:6], axis=0).astype(DATA_TYPE)
sens_norm_fact = np.zeros(6, dtype=DATA_TYPE)
for i in range(6):
sens_norm_fact[i] = 5.0 * float(np.max(np.abs(temp_states[:, i] - sens_deviation[i])))
print(f" norm: force_norm_fact={force_norm_fact:.6f}")
# --- Bias FIFO rollout ---
# Match legacy_env_imit.py line 142-143: [0, -U0, U0]
ff.apply_ddf() # restore save_ddf(1)
bias_arr = np.zeros(n_obj, dtype=DATA_TYPE)
bias_arr[3] = 0.0
bias_arr[4] = -1.0 * u0
bias_arr[5] = 1.0 * u0
fifo.clear()
for _ in range(FIFO_LEN):
ff.run(sample_interval, bias_arr)
fifo.append(ff.obs.copy()[0:12])
# --- save_ddf(2): post-bias checkpoint — consistent with FIFO ---
ff.get_ddf()
ff.save_ddf()
ff.apply_ddf() # restore to post-bias state for inference
# Closed-loop inference
sens_list, actions_list, force_list = [], [], []
bias_norm = np.array([0.0, 0.125, -0.125], dtype=np.float64)
a_prev = bias_norm.copy()
a_prev2 = a_prev.copy()
# Determine mode from SINDy results metadata (default: v23)
sindy_mode = coefs.get("mode", "v23")
for step in range(n_steps):
obs = fifo[-1] if fifo else np.zeros(12, dtype=np.float32)
obs_prev = list(fifo)[-2] if len(fifo) >= 2 else obs
a_prev_phys = (a_prev * action_scale + np.array(action_bias, dtype=np.float64)) * u0
a_prev2_phys = (a_prev2 * action_scale + np.array(action_bias, dtype=np.float64)) * u0
# Reconstruct target forces from harmonics (if available)
target_forces_step = None
if target_harmonics is not None:
target_forces_step = gen_target_states_at(step, target_harmonics)
if sindy_mode == "absolute":
omega_phys = predict_v23_deriv(
obs, obs_prev, a_prev_phys, a_prev2_phys, mu, u0, sample_interval/2000.0,
coefs["front_coef"], coefs["top_coef"],
feat_keys_front=coefs["feat_keys_front"],
feat_keys_rear=coefs["feat_keys_rear"],
target_forces=target_forces_step,
output_mode="absolute")
else:
omega_phys = predict_v23(
obs, a_prev_phys, a_prev2_phys, mu, u0,
coefs["front_coef"], coefs["top_coef"],
coefs["feat_names_front"], coefs["feat_names_rear"],
feat_keys_front=coefs["feat_keys_front"],
feat_keys_rear=coefs["feat_keys_rear"],
target_forces=target_forces_step)
norm_a = (omega_phys / u0 - np.array(action_bias, dtype=np.float64)) / action_scale
norm_a = np.clip(norm_a, -1.0, 1.0).astype(np.float32)
action_arr = scale_action(norm_a, scale=action_scale, bias=action_bias,
u0=u0, n_total_bodies=n_obj_total)
ff.context.push()
ff.run(sample_interval, action_arr)
ff.context.pop()
obs_new = ff.obs.copy()[0:12]
fifo.append(obs_new)
sens_list.append(obs_new[0:6])
force_list.append(obs_new[6:12])
actions_list.append(omega_phys.copy())
a_prev2 = a_prev.copy()
a_prev = norm_a.copy().astype(np.float64)
# Evaluate similarity over full run
sens_arr = np.array(sens_list, dtype=np.float32)
actions_arr = np.array(actions_list, dtype=np.float64)
# Use last 2*conv_len steps for similarity (after control has propagated)
tail_sens = sens_arr[-min(2*conv_len, n_steps):]
tail_target = target_states[-min(2*conv_len, len(target_states)):]
sim_tail = compute_similarity(tail_target, tail_sens, min(conv_len, len(tail_sens)//2))
sim_full = compute_similarity(target_states, sens_arr, conv_len)
action_range = float(np.max(np.abs(actions_arr)))
print(f" full-run similarity={sim_full:.4f} tail similarity={sim_tail:.4f} action_range={action_range:.4f}")
# Save controlled vorticity
try:
omega_ctrl = vorticity_from_ddf(ff, u0)
img_path = os.path.join(out_dir, f"{scene_name}_vorticity{threshold_str}.png")
save_vorticity_png(img_path, omega_ctrl,
title=f"{scene_name} SINDy closed-loop (v23, sim={sim_tail:.3f})")
print(f" controlled vorticity saved: {img_path}")
# Uncontrolled: restore from same post-bias DDF, run zero-action for same duration
ff.restore_ddf()
ff.apply_ddf()
# Run bias steps then zero-action for same duration as controlled
for _ in range(FIFO_LEN):
ff.run(sample_interval, bias_arr)
# Reset FIFO then run zero-action for n_steps
for _ in range(n_steps):
ff.run(sample_interval, np.zeros(n_obj, dtype=DATA_TYPE))
omega_unc = vorticity_from_ddf(ff, u0)
img_unc_path = os.path.join(out_dir, f"{scene_name}_uncontrolled{threshold_str}.png")
save_vorticity_png(img_unc_path, omega_unc,
title=f"{scene_name} uncontrolled")
print(f" uncontrolled vorticity saved: {img_unc_path}")
except Exception as e:
print(f" WARNING: controlled/uncontrolled vorticity export failed: {e}")
import traceback; traceback.print_exc()
del ff
# Target cylinder reference (separate env, fully stabilized)
try:
ff_tgt = FlowField(field_cfg, cuda_cfg, device_id=device_id)
tgt_diam = cfg["target_diameter"]
ff_tgt.add_cylinder((20.0 * l0, (ny - 1) / 2, 0.0), tgt_diam * l0)
for y_off in [2.0, 0.0, -2.0]:
sc = (sensor_x * l0, (ny - 1) / 2 + y_off * l0, 0.0)
ff_tgt.add_sensor(sc, l0 / 4.0)
n_obj_tgt = ff_tgt.obs.size // 2
ff_tgt.run(int(5 * NX / u0), np.zeros(n_obj_tgt, dtype=DATA_TYPE))
omega_tgt = vorticity_from_ddf(ff_tgt, u0)
img_tgt_path = os.path.join(out_dir, f"{scene_name}_target{threshold_str}.png")
save_vorticity_png(img_tgt_path, omega_tgt,
title=f"{scene_name} target cylinder ({tgt_diam}L)")
print(f" target vorticity saved: {img_tgt_path}")
del ff_tgt
except Exception as e:
print(f" WARNING: target vorticity export failed: {e}")
import traceback; traceback.print_exc()
return {
"scene": scene_name,
"mode": "v23",
"similarity_full": sim_full,
"similarity_tail": sim_tail,
"action_range": action_range,
"n_steps": n_steps,
"threshold": threshold if threshold is not None else "best",
}
def main():
ap = argparse.ArgumentParser(description="Illusion closed-loop SINDy validation")
ap.add_argument("--scene", type=str, required=True)
ap.add_argument("--device", type=int, default=0, help="GPU device")
ap.add_argument("--steps", type=int, default=0,
help="Steps (default: auto-set to cover 1*NX/U0)")
ap.add_argument("--threshold", type=float, default=None)
ap.add_argument("--sindy-results", type=str, default=None)
ap.add_argument("--out", type=str, default=None)
args = ap.parse_args()
if args.sindy_results is None:
args.sindy_results = os.path.join(
os.path.dirname(__file__), "..", "sindy", "illusion", "sindy_results_v2.json")
result = run_validation_illusion(
args.scene, args.sindy_results, args.device,
n_steps=args.steps, threshold=args.threshold, out_dir=args.out,
)
th_str = f"_th{args.threshold}" if args.threshold is not None else ""
out_dir = args.out or os.path.join(os.path.dirname(__file__), "results")
os.makedirs(out_dir, exist_ok=True)
out_path = os.path.join(out_dir, f"{args.scene}_v23{th_str}.json")
with open(out_path, "w") as f:
json.dump(result, f, indent=2)
print(f"Saved: {out_path}")
if __name__ == "__main__":
main()