diff --git a/docs/HANDOVER_SR_ANALYSIS.md b/docs/HANDOVER_SR_ANALYSIS.md new file mode 100644 index 0000000..9149682 --- /dev/null +++ b/docs/HANDOVER_SR_ANALYSIS.md @@ -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]`） diff --git a/docs/SR_analysis_results.md b/docs/SR_analysis_results.md new file mode 100644 index 0000000..8fe4cb1 --- /dev/null +++ b/docs/SR_analysis_results.md @@ -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 diff --git a/docs/figures/SR_analysis/fig1_illusion_comparison.png b/docs/figures/SR_analysis/fig1_illusion_comparison.png new file mode 100644 index 0000000..ca8b283 Binary files /dev/null and b/docs/figures/SR_analysis/fig1_illusion_comparison.png differ diff --git a/docs/figures/SR_analysis/fig2_karman_ablation.png b/docs/figures/SR_analysis/fig2_karman_ablation.png new file mode 100644 index 0000000..15074bf Binary files /dev/null and b/docs/figures/SR_analysis/fig2_karman_ablation.png differ diff --git a/docs/figures/SR_analysis/fig3_karman_generalization.png b/docs/figures/SR_analysis/fig3_karman_generalization.png new file mode 100644 index 0000000..9de1a1f Binary files /dev/null and b/docs/figures/SR_analysis/fig3_karman_generalization.png differ diff --git a/docs/figures/SR_analysis/fig4_r2_vs_closedloop.png b/docs/figures/SR_analysis/fig4_r2_vs_closedloop.png new file mode 100644 index 0000000..b2a6956 Binary files /dev/null and b/docs/figures/SR_analysis/fig4_r2_vs_closedloop.png differ diff --git a/docs/figures/SR_analysis/fig5_illusion_coefficients.png b/docs/figures/SR_analysis/fig5_illusion_coefficients.png new file mode 100644 index 0000000..23e088d Binary files /dev/null and b/docs/figures/SR_analysis/fig5_illusion_coefficients.png differ diff --git a/docs/figures/SR_analysis/fig6_roadmap.png b/docs/figures/SR_analysis/fig6_roadmap.png new file mode 100644 index 0000000..e0fe0e8 Binary files /dev/null and b/docs/figures/SR_analysis/fig6_roadmap.png differ diff --git a/scripts/plot_sr_results.py b/scripts/plot_sr_results.py new file mode 100644 index 0000000..f042ee0 --- /dev/null +++ b/scripts/plot_sr_results.py @@ -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}/") diff --git a/src/SR_analysis/README.md b/src/SR_analysis/README.md index a7a96b5..79bda9f 100644 --- a/src/SR_analysis/README.md +++ b/src/SR_analysis/README.md @@ -12,278 +12,174 @@ cloak and illusion scenes, using dimensionless physical features, G-equivariant structural constraints, and STLSQ threshold grids. For background, see: -- `src/sindy_sr_notes.md` -- execution plan -- `src/sindy_sr_knoeledge.md` -- confirmed facts and known pitfalls +- `sindy_sr_notes.md` -- execution plan and task tracking +- `sindy_sr_knowledge.md` -- confirmed facts and known pitfalls +- `../../docs/SR_analysis_results.md` -- comprehensive results report + +--- ## Directory Structure ``` SR_analysis/ - configs.py # Unified scene metadata (all 10 scenes) + configs.py # Unified scene metadata (all 10+ scenes) configs/ - legacy/ # Legacy CFD configs (config_cuda.json, config_flowfield.json) + legacy/ # Legacy CFD configs utils/ __init__.py # Selective exports (no pycuda dependency) - feature_builder.py # Dimensionless features + G-operator (from analysis_cloak) - sindy_fitter.py # STLSQ threshold grid, feature matrix builder + feature_builder.py # Dimensionless features + G-operator + phase-state features + sindy_fitter.py # STLSQ + feature matrices + derivative/absolute modes cfd_interface.py # LegacyCelerisLab wrapper (requires pycuda_3_10) g_operator.py # Equivariance diagnostics data/ - karman/ # Karman cloak: karman_re50, re100, re200, re400 - steady/ # Steady cloak: steady_data.npz - illusion/ # Illusion: illusion_0.75L, illusion_1L, illusion_1.5L - vortex/ # Vortex cloak: vortex_lamb, vortex_taylor + karman/ # Karman cloak: karman_re50/100/200/400 + steady/ # Steady cloak + illusion/ # Illusion: illusion_0.75L/1L/1.5L + vortex/ # Vortex cloak scripts/ infer_karman.py # Inference: LegacyCFD + PPO -> controlled.npz - infer_illusion.py # Inference: for 0.75L, 1L, 1.5L diameters - infer_vortex.py # Inference: for Lamb dipole + Taylor monopole + infer_illusion.py # Inference for illusion scenes + infer_vortex.py # Inference for vortex scenes sindy/ - run_karman.py # SINDy fitting for Karman scenes - run_illusion.py # SINDy fitting for Illusion scenes - run_vortex.py # SINDy fitting for Vortex scenes - run_pareto.py # Pareto-front analysis from SINDy results - karman/ # Output: sindy_results.json, pareto_*.json - illusion/ # Output: sindy_results.json, pareto_*.json - vortex/ # Output: sindy_results.json, pareto_*.json + run_all_v2.py # Unified SINDy fitting (supports --deriv, --phase, --output-mode etc.) + run_pysr.py # Restricted PySR symbolic regression + wrap_joint.py # Joint model -> wrapped format for validator + compare_v2.py # Cross-scene comparison report + karman/illusion/vortex/ # SINDy output JSONs 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/ - support_overlap.py # Pairwise support set comparison - shared_core.py # Multi-scene shared-core detection + support_overlap.py # Support set comparison + shared_core.py # Shared core detection ``` +--- + ## Key Design Decisions ### 1. Scene Metadata Driven -All scene parameters (Re, action scaling, geometry, model paths) are defined -once in `configs.py`, not hard-coded in scripts. Adding a new scene means -adding one dict to `configs.py`. +All scene parameters defined once in `configs.py`. -### 2. Data / Features / Models Separation +### 2. Feature Levels -- `data/` -- raw sensor/force/action arrays (.npz), one-time generation -- `sindy/` -- SINDy fitting results (JSON), reusable for comparison -- `scripts/` -- inference pipelines that produce `data/` +| Level | Features | Dim | Description | +|-------|----------|:---:|-------------| +| 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 -21 dimensionless features from raw lattice-unit sensor/force data: - -**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` +- **deriv**: predict `d(alpha)/dt`, then `alpha(t) = alpha(t-1) + dt_c * dalpha/dt` +- **absolute**: predict `alpha(t)` directly (no integration drift) ### 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 -Top(t) = f_rear(x(t)) # with bias -Bottom(t) = -f_rear(G[x(t)]) # shared-head: bottom = -top(Gx) +Front(t) = f_front(x(t)) # no bias, odd under G +Top(t) = f_rear(x(t)) # with bias +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]` -Per-channel: front (no bias), top (shared-head), bottom (independent, for comparison) +### Illusion — New Route: Phase-state + Error-state + Absolute Action -## 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_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 | +### Karman re100 — Ablation -Note: "Re_code" uses reference length 2*D (code convention). -Physical Re_D = Re_code / 2. E.g. Re_code=100 -> Re_D=50. +| Config | Feat | Output | R2 | Closed-loop | Note | +|--------|:----:|:-----:|:--:|:----------:|------| +| 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 -# Karman cloak -- all 4 training Re -conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_karman.py --re all --device 0 +# Illusion phase-state + absolute (recommended for 0.75L/1L) +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 -conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_karman.py --re 100 --device 0 --steps 200 +# 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 -# Illusion -- all 3 diameters -conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_illusion.py --diameter all --device 0 +# Karman expanded (10 dim) +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 -conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_vortex.py --type all --device 0 +# Karman with 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 + +# 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) - -```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) +### Closed-loop Validation ```bash +# Karman with absolute action conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \ - --scene karman_re70 --device 2 \ - --sindy-results src/SR_analysis/sindy/karman/sindy_results.json + --scene karman_re100 --device 0 --steps 200 --mode abs \ + --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 \ - --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 -# Pairwise support overlap -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 +conda run -n sr_env python src/SR_analysis/sindy/run_pysr.py --scene illusion_1L ``` -## 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 | -|---|---|---|---| -| karman_re50 | 0.998 | 0.989 | 0.996 | -| karman_re100 | 0.995 | 0.993 | 0.997 | -| karman_re200 | 0.957 | 0.914 | 0.918 | -| karman_re400 | 0.991 | 0.979 | 0.969 | -| illusion_0.75L | 0.991 | 0.989 | 0.990 | -| 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 | +- `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] +- "2U" in model name = S_DIM=14 (not 2x velocity), u0 always 0.01 +- SAMPLE_INTERVAL: 0.75L=400, 1L=600, 1.5L/Karman=800 +- Closed-loop steps auto-set: S=400→320, S=600→214, S=800→160 +- One-step R² high ≠ closed-loop good — always validate +- For phase-state features, always pass `sensors_raw`/`forces_raw` to enable derivative computation diff --git a/src/SR_analysis/configs.py b/src/SR_analysis/configs.py index 93ba369..9974f84 100644 --- a/src/SR_analysis/configs.py +++ b/src/SR_analysis/configs.py @@ -27,7 +27,7 @@ NX = 1280 NY = 512 CENTER_Y = (NY - 1) / 2.0 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: @@ -72,7 +72,30 @@ for rc, mn in [(50, "d1a3o12_re50"), (100, "d1a3o12_re100"), "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"] = { "scene_id": "steady", "re_code": 100, @@ -94,26 +117,29 @@ SCENES["steady"] = { "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: """Generate clean illusion scene key.""" s = f"{diam:.3f}".rstrip("0").rstrip(".") return f"illusion_{s}L" -_ILLUSION_1U = [ - (0.75, "d1a3o12_250525_imit_075L_1U"), - (1.0, "d1a3o12_250525_imit_1L_1U"), +_ILLUSION_BEST = [ + (0.75, "d1a3o14_250525_imit_075L_2U_400S", 400, 14), + (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) SCENES[key] = { "scene_id": "illusion", "target_diameter": diam, "re_code": 100, "mu": 2.0 / 100, - "nu": nu_from_re(100), + "nu": 0.004, # default, no Vis suffix in model name "has_disturbance": False, - "sample_interval": 600, + "sample_interval": si, + "conv_len": 36, # Illusion uses 36, not default 30 "action_scale": 8.0, "action_bias": (0.0, -2.0, 2.0), "source": "PPO_inference", @@ -124,19 +150,20 @@ for diam, mn in _ILLUSION_1U: "pinball_front_x": 19.0, "pinball_rear_x": 20.3, "target_type": "periodic", - "s_dim": 12, - "u0": U0, + "s_dim": sd, + "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)] = { "scene_id": "illusion", "target_diameter": 1.5, "re_code": 100, "mu": 2.0 / 100, - "nu": nu_from_re(100, u0=0.02), + "nu": 0.004, # default, no Vis suffix "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_bias": (0.0, -2.0, 2.0), "source": "PPO_inference", @@ -148,7 +175,7 @@ SCENES[_illusion_key(1.5)] = { "pinball_rear_x": 20.3, "target_type": "periodic", "s_dim": 14, - "u0": 0.02, + "u0": U0, # always 0.01 } # -- 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 --------------------------------------------------------- def get_scene(name: str) -> dict: diff --git a/src/SR_analysis/scripts/infer_illusion.py b/src/SR_analysis/scripts/infer_illusion.py index 72dc90c..3278384 100644 --- a/src/SR_analysis/scripts/infer_illusion.py +++ b/src/SR_analysis/scripts/infer_illusion.py @@ -4,9 +4,9 @@ Generates controlled data for a given target cylinder diameter using LegacyCelerisLab + trained PPO model. 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 - 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 """ from __future__ import annotations @@ -37,12 +37,61 @@ from SR_analysis.utils.cfd_interface import ( ) from SR_analysis.configs import ( 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 +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( scene_name: str, device_id: int, @@ -55,13 +104,15 @@ def run_single_illusion( u0 = cfg["u0"] l0 = 20.0 sample_interval = cfg["sample_interval"] + conv_len = cfg.get("conv_len", 30) # Illusion uses 36 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"] sensor_x = cfg["sensor_x"] # 30.0 for illusion front_x = cfg["pinball_front_x"] # 19.0 rear_x = cfg["pinball_rear_x"] # 20.3 target_diam = cfg["target_diameter"] + s_dim = cfg["s_dim"] # 14 for all best illusion models 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)) 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) 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) ------------ ff = FlowField(field_cfg, cuda_cfg, device_id=device_id) 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})") 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)...") 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) - for _ in range(FIFO_LEN): + for _ in range(fifo_len): 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)) print(f" Target recorded: {target_states.shape}") - # Save target - np.savez(os.path.join(output_root, "target.npz"), target_states=target_states) + # Analyze harmonics of force channels (channel 0,1 = force; 2..7 = sensors) + # 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 del ff @@ -123,7 +186,6 @@ def run_single_illusion( ff.add_sensor(sc, l0 / 4.0) # 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((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.save_ddf() - # Norm collection (zero action) - fifo = deque(maxlen=FIFO_LEN) - for _ in range(FIFO_LEN): + # Norm collection (zero action) — matches legacy_env_imit.py + 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]) @@ -161,21 +223,28 @@ def run_single_illusion( } 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() bias_arr = np.zeros(n_obj, dtype=DATA_TYPE) - bias_arr[3] = float(action_bias[0] * u0) - bias_arr[4] = float(action_bias[1] * u0) - bias_arr[5] = float(action_bias[2] * u0) + bias_arr[3] = 0.0 # front = 0 + bias_arr[4] = -1.0 * u0 # bottom = -0.01 (cf. action_bias bottom = -2*u0 = -0.02) + bias_arr[5] = 1.0 * u0 # top = 0.01 (cf. action_bias top = 2*u0 = 0.02) print(f" bias action: {bias_arr}") fifo.clear() - for _ in range(FIFO_LEN): + for _ in range(fifo_len): ff.run(sample_interval, bias_arr) fifo.append(ff.obs.copy()[0:12]) save_states = np.array(list(fifo), dtype=DATA_TYPE) 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 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) # -- 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) if model_path is not None: - s_dim = cfg.get("s_dim", 12) 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.set_random_seed(0) @@ -196,15 +264,16 @@ def run_single_illusion( ff.restore_ddf() ff.apply_ddf() - # Re-bias FIFO - fifo = deque(maxlen=FIFO_LEN) - for _ in range(FIFO_LEN): + # Re-bias FIFO (using same bias action as Phase 2) + fifo = deque(maxlen=fifo_len) + for _ in range(fifo_len): ff.context.push() ff.run(sample_interval, bias_arr) ff.context.pop() fifo.append(ff.obs.copy()[0:12]) sens_list, forc_list, action_list = [], [], [] + target_forces_list = [] obs = np.zeros(s_dim, dtype=np.float32) for step in range(n_infer_steps): @@ -227,28 +296,34 @@ def run_single_illusion( sens_list.append(obs_slice[0:6]) 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 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) 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[:12] = obs12 + obs[12] = np.clip(target_cd, -1.0, 1.0) + obs[13] = np.clip(target_cl, -1.0, 1.0) else: obs = obs12 np.savez(os.path.join(output_root, "controlled.npz"), sensors=np.array(sens_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) - # For comparison, compute similarity between controlled sensors and target - target_sensors = target_states[:, 0:6] - sim = compute_similarity(target_sensors, - np.array(sens_list, dtype=np.float32), CONV_LEN) + # Compute similarity: compare controlled sensors vs target sensors + # target_sensors = target_states[:, 2:8] (sensor only, 6 channels) + sens_arr = np.array(sens_list, dtype=np.float32) + sim = compute_similarity(target_sensors, sens_arr, conv_len) print(f" similarity (vs target cylinder) = {sim:.4f}") result["controlled"] = True @@ -265,7 +340,7 @@ def run_single_illusion( def main(): - ap = argparse.ArgumentParser(description="Illusion inference") + ap = argparse.ArgumentParser(description="Illusion inference (corrected)") ap.add_argument("--diameter", type=str, default="1.0", help='Diameter: 0.75, 1.0, 1.5, or "all"') ap.add_argument("--device", type=int, default=0, help="GPU device ID") @@ -278,7 +353,6 @@ def main(): scene_names = get_scene_list("illusion") else: d = float(args.diameter) - # Match by target_diameter field scene_names = [] for sn in get_scene_list("illusion"): cfg = get_scene(sn) diff --git a/src/SR_analysis/sindy/run_all_v2.py b/src/SR_analysis/sindy/run_all_v2.py new file mode 100644 index 0000000..7a7dd4e --- /dev/null +++ b/src/SR_analysis/sindy/run_all_v2.py @@ -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//") + 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() diff --git a/src/SR_analysis/sindy/run_pysr.py b/src/SR_analysis/sindy/run_pysr.py new file mode 100644 index 0000000..68eb7bb --- /dev/null +++ b/src/SR_analysis/sindy/run_pysr.py @@ -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() diff --git a/src/SR_analysis/sindy/whitelist_illusion.json b/src/SR_analysis/sindy/whitelist_illusion.json new file mode 100644 index 0000000..8257020 --- /dev/null +++ b/src/SR_analysis/sindy/whitelist_illusion.json @@ -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." +} diff --git a/src/SR_analysis/sindy/whitelist_karman.json b/src/SR_analysis/sindy/whitelist_karman.json new file mode 100644 index 0000000..5a070de --- /dev/null +++ b/src/SR_analysis/sindy/whitelist_karman.json @@ -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)." +} diff --git a/src/SR_analysis/sindy/wrap_joint.py b/src/SR_analysis/sindy/wrap_joint.py new file mode 100644 index 0000000..2d2e088 --- /dev/null +++ b/src/SR_analysis/sindy/wrap_joint.py @@ -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) diff --git a/src/SR_analysis/sindy_sr_knowledge.md b/src/SR_analysis/sindy_sr_knowledge.md new file mode 100644 index 0000000..b774b39 --- /dev/null +++ b/src/SR_analysis/sindy_sr_knowledge.md @@ -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.619；static 8-dim 模型 R²=0.321 但闭环 0.745。 + +### 5. 时间特征加剧 rollout mismatch（2026-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，不做导数积分 +- 结构：v23（front 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 模型加载、数据采集、SINDy（pysindy） +- `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 推理数据）。 +- **单步 validator（predict_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 维度**。 \ No newline at end of file diff --git a/src/SR_analysis/sindy_sr_notes.md b/src/SR_analysis/sindy_sr_notes.md index f64dd37..a4075b1 100644 --- a/src/SR_analysis/sindy_sr_notes.md +++ b/src/SR_analysis/sindy_sr_notes.md @@ -1,4 +1,4 @@ -# SINDy 与 SR 执行计划 +# SINDy 与 SR 执行计划（2026-06-15 修订版 v4） ## 文档作用 @@ -6,342 +6,183 @@ - 只写执行路线、阶段目标、优先级、输出要求 - 不长篇复述历史争论 -- 不把背景知识和任务顺序混在一起 -- 凡是历史 bug、经验教训、哪些结论已经成立、哪些还不能说，统一放到 `sindy_sr_knowladge` +- 凡是背景知识、已知结论、踩坑记录，统一放到 `sindy_sr_knowledge.md` --- ## 当前总目标 -把 **所有 cloak 场景** 纳入同一受约束分析框架，使用 **SINDy + SR 并行** 探索控制骨架，先分别拟合，再横向比较，从而判断： +用 SINDy + SR 在所有 cloak / illusion 场景上产出有物理结论的控制律公式，并验证其跨场景泛化能力。 -- 哪些项是 shared core -- 哪些项是 scene-specific activation -- 哪些差异来自时间尺度写法，而不是物理骨架本身 - -当前默认路线不是“先强行做统一总公式”，而是： - -\[ -\text{separate fit} \rightarrow \text{compare} \rightarrow \text{shared-backbone test} -\] +回答三个问题： +1. **哪些物理量是控制的核心变量？**（力、速度、相位、误差之间的取舍） +2. **不同场景是否共享同一公式骨架？**（family 内 -> family 间） +3. **白箱控制律能否在未训练工况工作？**（泛化性 + 失效边界） --- -## 当前工作原则 +## 已完成的工作（2026-06-14~15） -后续默认遵守： +### 核心成果 -1. **所有场景共用同一套 primitive variables 与同一套 \(G\) 规则** -2. **front 默认 no-bias + odd structure** -3. **rear 默认 shared-head** -4. **SINDy 与 SR 并行推进**，SR 不是替换 SINDy，而是并行工具 -5. **先分场景拟合，再做横向比较** -6. **闭环验证不可省略**，但不要求每个场景第一轮都做全套闭环 -7. **时间尺度问题必须显式进入特征定义** -8. **不再只围绕跨 Re 的 Kármán 单线深挖**；跨 Re 现在是已站住的第一证据，不是全部主线 -9. **当前不讨论功率与能量分析**，这不是这条线的优先任务 -10. **当前不把时延作为主要矛盾**；对稳定周期控制，当前第一矛盾是变量骨架与时间尺度写法 +1. **Phase-state 特征体系**：`PHASE_STATE_KEYS` (6维) + `ILLUSION_PHASE_KEYS` (10维) — 无需动作历史 +2. **绝对动作输出模式**：`get_feature_matrix_deriv(output_mode="absolute")` — 无积分累积 +3. **闭合验证器支持**：`predict_v23_deriv(output_mode="absolute")` + 闭环 `mode="abs"` +4. **离线 rollout 评估**：`validate/eval_rollout.py` — 多步滚动误差分析 +5. **消融实验**：5 种输入配置 × 2 种输出模式 × CFD 闭环验证 + +### Illusion 新路线成功 + +| 场景 | 闭环 sim | 动作历史 | 结论 | +|------|:--------:|:--------:|------| +| 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 -- steady cloak +输入特征（10 维）： +```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 的公式形态比较 / 同形骨架判断 -- 单涡 cloak（monopole / taylor / lamb 等已有单涡场景） -- erase -- 其他已有 cloak 场景 +### 第二优先级：Karman 状态补强（CCD/OID 进场） -这批先做轻量版 separate fit，再决定哪些值得补完整闭环与深挖。 +当前 Karman 新路线最佳 0.699，需要在 phase-state 基础上补充缺失状态量： +- 候选：回流区长度、尾迹中心线偏移、POD 模态系数 +- 方法：CCD（Correlation-based Decomposition）或 OID（Observer-based Identification） +- 目标：在无动作历史前提下将 Karman 闭环提升到 0.85+ + +### 第三优先级：Illusion 跨直径联合 + +条件：0.75L 和 1L 的 PySR 公式确认同形骨架后 +方法：error-state 特征 + 目标 signature 描述符（target St、amplitude） --- -## 阶段 0 +## 当前暂不建议做的事 -## 统一接口 - -### 0.1 统一 primitive variables - -所有场景统一输出： - -- 无量纲 sensor：\(\hat u, \hat v\) -- 力系数：\(C_D, C_L\) -- 无量纲动作：\(\alpha\) -- lagged \(\alpha\) -- \(\Delta \alpha\) 或显式含 \(\Delta t_c\) 的版本 -- \(\mu = 1/Re_D\) -- scene metadata：scene 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） +- 不回退到动作历史主导（Illusion 已经证明不需要） +- 不先做跨 Re Karman 联合拟合（单 Re 还没站稳） +- 不先做 vortex 扩展（不是当前瓶颈） +- 不先上 PySR 到 Karman（状态还有问题，公式不会更物理） --- -## 阶段 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 -- front odd structure -- rear shared-head -- correct \(G\) consistency +# Illusion separate with error-state (old style for comparison) +conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \ + --scenes illusion_0.75L,illusion_1L,illusion_1.5L -### 每个场景必须输出 +# 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 +``` -| 输出 | 说明 | -|---|---| -| best support | 最优 support 列表 | -| sparsity curve | 稀疏度-误差曲线 | -| front / rear 主项表 | 主导项与系数 | -| one-step metrics | R²、RMSE | -| key closed-loop result | 至少一个关键闭环指标 | -| support stability | threshold / window / bootstrap 稳定性 | -| contribution table | 主要项贡献度，不只看是否出现 | +### CFD 闭环验证 +```bash +# Illusion phase-state + absolute action +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 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 +``` -- 哪些项共同出现 -- 哪些项只在 Kármán 激活 -- steady 是否表现为明显简化版 -- overlap 不能只看布尔出现，还要看贡献量级与稳定性 +### 离线 rollout 评估 +```bash +python3 src/SR_analysis/validate/eval_rollout.py \ + --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 - -### 目标 - -不是追求最终公式，而是看： - -- 在受限物理库上，是否能压出更短闭式 -- SINDy 中看似不同的项，是否会被更统一的表达吸收 -- Kármán 与 steady 是否出现同形公式 - -### SR 输入规则 - -SR 只能使用： - -- 阶段 0 的统一变量 -- 阶段 1 的 SINDy 已筛出主项及其邻近项 -- 受限运算集合 - -### 当前允许的运算 - -- 加减乘 -- protected divide -- 少量 square -- 必要时有限放开 tanh - -### 当前不允许的运算 - -- 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 的第一轮受限 SR（PySR） -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 判断。 \ No newline at end of file +## 关键文件索引 + +| 文件 | 用途 | +|------|------| +| `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` | 等变性诊断 | +| `scripts/infer_*.py` | Karman/Illusion/Vortex 推理管线 | +| `sindy/run_all_v2.py` | 统一 SINDy 拟合入口（支持 --deriv, --phase, --output-mode, --karman-expand 等） | +| `sindy/run_pysr.py` | 受限 PySR 符号回归 | +| `sindy/wrap_joint.py` | 联合模型 → wrapped 格式 | +| `validate/run_closed_loop.py` | Karman 闭环验证器（v23, deriv, abs 模式） | +| `validate/run_closed_loop_illusion.py` | Illusion 闭环验证器 | +| `validate/eval_rollout.py` | **新**：离线多步 rollout 评估 | +| `compare/support_overlap.py` | 跨场景 support 比较 | +| `compare/shared_core.py` | 多场景 shared core 检测 | +| `docs/SR_analysis_results.md` | **新**：完整分析报告 | +| `docs/figures/SR_analysis/fig*.png` | **新**：结果图表（6张） | + +## 环境 + +| 环境 | conda 名 | 用途 | +|------|---------|------| +| CFD + DRL + SINDy | `pycuda_3_10` | infer, run_all_v2, closed-loop | +| PySR 符号回归 | `sr_env` | run_pysr.py | +| GPU | device 0, device 1 | | diff --git a/src/SR_analysis/utils/__init__.py b/src/SR_analysis/utils/__init__.py index dc9391e..d8f4d98 100644 --- a/src/SR_analysis/utils/__init__.py +++ b/src/SR_analysis/utils/__init__.py @@ -1,9 +1,19 @@ from .feature_builder import ( compute_dimensionless, compute_features, build_feature_matrix, 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 ( fit_channel, fit_sindy, print_control_law, get_active_support, get_feature_matrix_from_data, + fit_sindy_weighted, get_feature_matrix_v2, + compute_action_deriv, get_feature_matrix_deriv, ) diff --git a/src/SR_analysis/utils/cfd_interface.py b/src/SR_analysis/utils/cfd_interface.py index 5c68c6c..a0b0f5e 100644 --- a/src/SR_analysis/utils/cfd_interface.py +++ b/src/SR_analysis/utils/cfd_interface.py @@ -208,6 +208,12 @@ def add_pinball( fifo.append(flow_field.obs.copy()[obs_slice_start:obs_slice_end]) 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() return { @@ -283,7 +289,7 @@ def vorticity_from_ddf(flow_field: FlowField, u0: float) -> np.ndarray: - ddf[:, :, 3] - ddf[:, :, 6] - ddf[:, :, 7]) / u0 uy = (ddf[:, :, 2] + ddf[:, :, 5] + ddf[:, :, 6] - 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) diff --git a/src/SR_analysis/utils/feature_builder.py b/src/SR_analysis/utils/feature_builder.py index 3aaef70..8713f50 100644 --- a/src/SR_analysis/utils/feature_builder.py +++ b/src/SR_analysis/utils/feature_builder.py @@ -7,7 +7,7 @@ Copy of analysis_cloak/common/feature_builder.py -- kept as canonical source. """ from __future__ import annotations -from typing import Dict, List, Tuple +from typing import Dict, List, Optional, Tuple import numpy as np @@ -86,6 +86,7 @@ def apply_G_x(dim: Dict[str, np.ndarray], # -- Feature key definitions ------------------------------------------------- +# Original feature set (includes sin_ua/cos_ua) CORE_FEAT_KEYS = [ "u_m", "u_a", "u_c", "v_a", @@ -96,9 +97,99 @@ CORE_FEAT_KEYS = [ "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"] +# 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 +# V2 all features (for cross-Re joint fitting with mu) +ALL_FEAT_KEYS_V2 = CORE_FEAT_KEYS_V2 + MU_FEAT_KEYS # -- Feature computation ----------------------------------------------------- @@ -110,7 +201,13 @@ def compute_features( mu: float, alpha_mode: bool = False, # if True, actions_prev are already nondim alpha 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 + 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]: """Compute unified feature dictionary from dimensionless primitives. @@ -121,6 +218,8 @@ def compute_features( mu: 1/Re_D alpha_mode: if True, actions are already nondim; else convert 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 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_diff"] = Cl_T - Cl_B - # Phase - sym["sin_ua"] = np.sin(np.pi * sym["u_a"]) - sym["cos_ua"] = np.cos(np.pi * sym["u_a"]) + # Phase (optional, may obscure linear structure) + if include_cos_sin: + 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["aB_lag1"] = a[:, 1] sym["aT_lag1"] = a[:, 2] @@ -165,6 +265,22 @@ def compute_features( sym["daB"] = a[:, 1] - a2[:, 1] 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 if include_mu: 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_Cd_tot"] = sym["Cd_tot"] * 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 diff --git a/src/SR_analysis/utils/sindy_fitter.py b/src/SR_analysis/utils/sindy_fitter.py index 6f97b39..db5cc50 100644 --- a/src/SR_analysis/utils/sindy_fitter.py +++ b/src/SR_analysis/utils/sindy_fitter.py @@ -179,3 +179,325 @@ def get_feature_matrix_from_data( return (Theta_f[n_warmup:], Theta_r[n_warmup:], actions_phys[n_warmup:], 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) diff --git a/src/SR_analysis/validate/eval_rollout.py b/src/SR_analysis/validate/eval_rollout.py new file mode 100644 index 0000000..b941a9b --- /dev/null +++ b/src/SR_analysis/validate/eval_rollout.py @@ -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() diff --git a/src/SR_analysis/validate/run_closed_loop.py b/src/SR_analysis/validate/run_closed_loop.py index 3b56b9f..6f8534f 100644 --- a/src/SR_analysis/validate/run_closed_loop.py +++ b/src/SR_analysis/validate/run_closed_loop.py @@ -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.feature_builder import ( 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.configs import ( @@ -57,10 +57,18 @@ def build_feature_vector( mu: float, u0: float, add_bias: bool, + feat_keys: Optional[List[str]] = None, + use_v2: bool = False, + target_forces: Optional[np.ndarray] = None, ) -> np.ndarray: """Build a single-row feature vector from raw obs and action state. 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) 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) 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 @@ -84,26 +105,42 @@ def predict_v23( top_coef: np.ndarray, feat_names_front: 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: """Predict actions using v23: front no-bias + rear shared-head. Returns (3,) physical omega array: [front, bottom, top]. """ # Front channel: no bias - front = float(np.dot( - build_feature_vector(obs_slice, a_prev, a_prev2, mu, u0, add_bias=False), - front_coef)) + fv_front = build_feature_vector( + obs_slice, a_prev, a_prev2, mu, u0, add_bias=False, + feat_keys=feat_keys_front, + target_forces=target_forces, + ) + front = float(np.dot(fv_front, front_coef)) # Top channel: with bias - top = float(np.dot( - build_feature_vector(obs_slice, a_prev, a_prev2, mu, u0, add_bias=True), - top_coef)) + fv_top = build_feature_vector( + obs_slice, a_prev, a_prev2, mu, u0, add_bias=True, + feat_keys=feat_keys_rear, + target_forces=target_forces, + ) + top = 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, a_prev2) - bottom = -float(np.dot( - build_feature_vector(G_obs, G_a_prev, G_a_prev2, mu, u0, add_bias=True), - top_coef)) + # G transforms forces: swap bottom<->top, negate fy. Target force G transform: + # target_Cd is a drag-like quantity that is even under G (stays), + # 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) @@ -118,20 +155,125 @@ def predict_unstructured( bottom_coef: np.ndarray, top_coef: np.ndarray, feat_names: List[str], + feat_keys: Optional[List[str]] = None, ) -> np.ndarray: """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)) bottom = float(np.dot(feat, bottom_coef)) + bottom = float(np.dot(feat, bottom_coef)) top = float(np.dot(feat, top_coef)) 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. + 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, - 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: 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_r = per["feature_names_rear"] - front_coef = np.array(per["front"]["best_coef"], dtype=np.float64) - top_coef = np.array(per["top"]["best_coef"], dtype=np.float64) - bottom_coef = np.array(per["bottom"]["best_coef"], dtype=np.float64) + def _coef_at_threshold(ch_data, th): + """Get coef at given threshold, or best if th=None.""" + 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 + # Extract mode from per_scene metadata (default: "v23" for old-style results) + sindy_mode = per.get("mode", "v23") + return { "front_coef": front_coef, "top_coef": top_coef, "bottom_coef": bottom_coef, "feat_names_front": fn_f, "feat_names_rear": fn_r, + "feat_keys_front": feat_keys_front, + "feat_keys_rear": feat_keys_rear, "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, coefs: Dict[str, Any], device_id: int, - n_steps: int = 100, + n_steps: int = 0, mode: str = "v23", ) -> dict: """Run closed-loop validation using a SINDy control law. Parameters ---------- - scene_name : e.g. "karman_re70" - 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. + n_steps : int, default=0 (auto: >= NX/U0 / sample_interval) """ cfg = get_scene(scene_name) re_code = cfg["re_code"] @@ -189,8 +360,17 @@ def run_validation( action_scale = cfg["action_scale"] action_bias = cfg["action_bias"] 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 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], ) - # Reset to checkpoint - ff.restore_ddf() - ff.apply_ddf() - - # Bias FIFO init + # Now ff is at save_ddf(2) = post-bias state (add_pinball saves after bias rollout). + # Norm dict contains save_states from the bias FIFO. fifo = deque(maxlen=FIFO_LEN) - bias_arr = scale_action(np.zeros(3, dtype=np.float32), scale=action_scale, - bias=action_bias, u0=u0, n_total_bodies=n_obj_total) - for _ in range(FIFO_LEN): - ff.run(sample_interval, bias_arr) - fifo.append(ff.obs.copy()[2:14]) + for row in norm["save_states"]: + fifo.append(row) # Closed-loop with SINDy law sens_list = [] actions_list = [] - a_prev = action_to_physical(np.zeros((1, 3), dtype=np.float32), - scale=action_scale, bias=action_bias, u0=u0).flatten() + # For Karman: bias action = DRL action_bias = [0, -4, 4]. + # Normalized zero action [0,0,0] decodes to the same bias. + a_prev = np.zeros(3, dtype=np.float64) a_prev2 = a_prev.copy() for _ 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 # 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": - omega = predict_v23( - obs, a_prev, a_prev2, mu, u0, + 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"]) + 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": - omega = predict_unstructured( - obs, a_prev, a_prev2, mu, u0, + omega_phys = predict_unstructured( + obs, a_prev_phys, a_prev2_phys, mu, u0, coefs["front_coef"], coefs["bottom_coef"], coefs["top_coef"], - coefs["feat_names_rear"]) + coefs["feat_names_rear"], + feat_keys=coefs["feat_keys_rear"]) else: raise ValueError(f"Unknown mode: {mode}") - # Clip to valid action range - norm_a = (omega / u0 - np.array(action_bias, dtype=np.float64)) / action_scale + # Convert physical omega -> normalized action -> legacy array + 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) - - # Apply to CFD action_arr = scale_action(norm_a, scale=action_scale, bias=action_bias, u0=u0, n_total_bodies=n_obj_total) + ff.run(sample_interval, action_arr) obs_new = ff.obs.copy()[2:14] fifo.append(obs_new) 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_prev = omega.copy() + a_prev = norm_a.copy().astype(np.float64) # Evaluate 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("--steps", type=int, default=100) 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, 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, help="Output directory for result JSON") args = ap.parse_args() if args.sindy_results is None: 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, 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: out_dir = os.path.join(os.path.dirname(__file__), "..", "validate", "results") else: out_dir = args.out 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: json.dump(result, f, indent=2) print(f"Saved: {out_path}") diff --git a/src/SR_analysis/validate/run_closed_loop_illusion.py b/src/SR_analysis/validate/run_closed_loop_illusion.py new file mode 100644 index 0000000..e1aa0b5 --- /dev/null +++ b/src/SR_analysis/validate/run_closed_loop_illusion.py @@ -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()