fix(oid): confirm FIFO bias bug has no structural impact

- Fixed bias_arr[4] (front 0), bias_arr[5] (bottom -4U0), bias_arr[6] (top +4U0)
- Re-ran full karman pipeline: force-sig overlap unchanged (-0.034)
- Force-OID still beats POD (0.295 vs 0.068, was 0.750 vs 0.418)
- Absolute R2 shifted because corrected FIFO changed PPO trajectory start
- Structural conclusion (force-sig near-orthogonal) is robust

Co-authored-by: Cursor <cursoragent@cursor.com>
This commit is contained in:
Frank14f 2026-06-30 16:26:54 +08:00
parent 5c55c5bdf7
commit 2ae248421d
282 changed files with 27668 additions and 368 deletions

View File

@ -33,5 +33,5 @@
#define V_TAYLOR 0b00000001
// variables
#define N_OBJS 6
#define N_OBJS 7
// #define N_SENS 2

View File

@ -0,0 +1,50 @@
{
"_doc": "Karman Cloak Re100: uniform inlet, free-slip walls, 2000x600 grid. Pinball centered.",
"grid": {
"lattice_model": "D2Q9",
"nx": 2000,
"ny": 600,
"nz": 1
},
"physics": {
"data_type": "FP32",
"viscosity": 0.004,
"velocity": 0.01,
"rho": 1.0
},
"method": {
"collision": "MRT",
"streaming": "double_buffer",
"store_precision": "FP32",
"ddf_shifting": false,
"les": {
"enabled": false,
"cs": 0.16,
"closed_form": true
},
"trt": {
"magic_param": 0.1875
},
"inlet": {
"profile": "uniform",
"scheme": "regularized",
"trt_neq_damp": 0.5,
"regularized_neq_damp": 0.5
},
"outlet": {
"mode": "neq_extrap",
"backflow_clamp": true,
"blend_alpha": 0.7,
"srt_neq_damp": 0.5
},
"y_wall_bc": "free_slip",
"omega_guard": {
"min": 0.01,
"max": 1.99
}
},
"cuda": {
"threads_per_block": 256,
"compute_capability": "auto"
}
}

1000
src/CCD_analysis/Lyu23.md Normal file

File diff suppressed because it is too large Load Diff

View File

@ -1,30 +1,85 @@
# CCD_analysis
# CCD_analysis: Correction-Field CCD Pipeline
Canonical Correlation Decomposition for fluidic pinball control analysis.
Analyzes DRL-controlled fluidic pinball using **correction-field decomposition** + **Canonical Correlation Decomposition (CCD/Lyu23)**. Core question: does `dq_ctl` (what the controller adds) match `dq_tar` (what the target requires)?
## Quick Start
1. Read `ccd_knowledge.md` — 唯一知识库,包含理论、流程、结果
2. Run comparison: `conda run -n pycuda_3_10 python3 correction_analysis/compare_dqctl_scenes.py`
3. Run diagnostics: `conda run -n pycuda_3_10 python3 correction_analysis/diagnose_corrections.py`
```bash
# 1. Panorama comparison figure (primary output)
conda run -n pycuda_3_10 python3 correction_analysis/compare_dqctl_scenes.py
## Directory
# 2. CCD quantitative decomposition
conda run -n pycuda_3_10 python3 correction_analysis/decompose_corrections.py
# 3. Single-scene diagnostics (zone metrics)
conda run -n pycuda_3_10 python3 correction_analysis/diagnose_corrections.py
```
## Pipeline Architecture
```
correction_analysis/ — 所有当前分析代码
scripts/ — GPU 数据采集 + phase alignment
utils/ — 核心算法POD/CCD/场加载/平移)
ccd/ — Round 5 旧基线(已冻结)
scripts/collect_*.py → scripts/{detect_period,replay_fields}.py
(GPU采集) (phase alignment)
correction_analysis/compute_correction_fields.py
(dq_blk, dq_ctl, dq_tar)
┌────┴────┬────────────┐
↓ ↓ ↓
compare_ decompose_ diagnose_
dqctl_ corrections corrections
scenes.py .py .py
(全景对比) (CCD定量) (zone诊断)
```
## Directory Structure
```
CCD_analysis/
README.md # 本文件
ccd_knowledge.md # 唯一知识库 (理论, 结果, 操作流程, bug经验)
Lyu23.md # CCD 方法文献
configs.py # 场景元数据 (统一几何)
utils/
resampling.py # POD, CCD, 场加载
cfd_interface.py # LegacyCelerisLab 封装 (GPU) + build_observation
field_translate.py # 场平移 (备用)
load_vortex_fields.py # 瞬态 vortex 场加载
scripts/
collect_*.py # GPU 数据采集
detect_period.py # 周期检测 → phase_plan.json
replay_fields.py # 场回放 → fields_aligned.npz
correction_analysis/
compute_correction_fields.py # dq 计算 (核心)
compare_dqctl_scenes.py # 全景对比图 (核心输出)
decompose_corrections.py # CCD 定量 (POD + force/action CCD)
diagnose_corrections.py # 单场景 zone 诊断
run_signature_line.py # Signature CCD
run_zone_ccd.py # Zone-restricted CCD
run_steady_metrics.py # Steady cloak 定量
run_15L_correction.py # 1.5L 专项
visualize_action_ccd.py # Action-CCD mode1 可视化
process_legacy_steady.py # 旧格式加载
ccd/ # Round 5 冻结基线 (勿改)
data/
figures/ — 诊断图(无 colorbar
ccd/ — JSON 结果
old_data/ — 归档废弃数据
figures/ # 诊断图 (仅保留核心对比图)
ccd/ # JSON 结果
old_data/ # 归档 (废弃脚本/旧报告/旧数据)
```
## Environment
## Key Documentation
```
conda run -n pycuda_3_10
```
| File | Content |
|------|---------|
| `ccd_knowledge.md` | **Primary entry** — theory, results, conventions, bug history |
| `Lyu23.md` | CCD method paper (Lyu 2023) |
| `data/old_data/ccd_correction_field_report.md` | Old full report (archived reference) |
| `data/old_data/ccd_handover.md` | Old handover notes (archived) |
详见 `ccd_knowledge.md`
## Key Conventions
- **Main analysis object**: `dq_ctl` (not raw `q_ctl`)
- **Observation order**: `[forces/force_norm, sensors/sens_norm]` — force first (see Bug 3 in ccd_knowledge.md §12)
- **All scenes unified geometry**: pinball center at 613px, sensors at 800px
- **Environment**: `conda run -n pycuda_3_10`
- **GPU**: Device 2 for collection

View File

@ -0,0 +1,270 @@
"""LOCO and blocked-split validation for CCD (Round 5).
Reuses shared data loader from resampling.py.
Target-only POD basis. Q_delay=6 for force/action.
0.75L and 1.0L only. 1.5L excluded from validation.
Usage:
conda run -n pycuda_3_10 python ccd/validate.py
"""
from __future__ import annotations
import json
import os
import sys
import time
import numpy as np
_SRC = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", ".."))
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from CCD_analysis.configs import DATA_DIR
from CCD_analysis.utils.resampling import (
compute_reduced_ccd, cumulative_energy,
load_aligned_fields, make_force_obs,
build_field_matrix, project_into_basis,
)
R_LIST = [6, 8, 10]
CCD_Q = 6
N_CYCLES = 4
N_PTS = 24
DIAMETERS = [0.75, 1.0]
def r2(y_true: np.ndarray, y_pred: np.ndarray) -> float:
"""Coefficient of determination."""
ss_r = np.sum((y_true - y_pred) ** 2)
ss_t = np.sum((y_true - np.mean(y_true)) ** 2)
return float(1.0 - ss_r / (ss_t + 1e-12))
def reconstruct_observable(W, sigma, R, a_test, y_train):
"""Reconstruct observable from CCD modes.
Returns dict with 'mode1' and 'm80' reconstructions.
"""
am = np.mean(a_test, axis=1, keepdims=True)
as_ = np.std(a_test, axis=1, keepdims=True) + 1e-12
a_test_z = (a_test - am) / as_
z_test = W.T @ a_test_z
ym = np.mean(y_train, axis=1, keepdims=True)
ys = np.std(y_train, axis=1, keepdims=True) + 1e-12
half = CCD_Q // 2
m_obs = y_train.shape[0]
en = cumulative_energy(sigma)
m80 = int(np.searchsorted(en, 0.80) + 1) if len(en) > 0 else 1
results = {}
# Mode-1
if R.shape[1] >= 1:
pz_1 = R[:, :1] * sigma[:1] @ z_test[:1, :]
yp_1 = pz_1[half * m_obs:(half + 1) * m_obs, :] * ys + ym
results["mode1"] = yp_1
else:
results["mode1"] = np.zeros_like(y_train[:, :a_test.shape[1]])
# M80
n_rm = min(m80, R.shape[1])
if n_rm >= 1:
pz_m = R[:, :n_rm] * sigma[:n_rm] @ z_test[:n_rm, :]
yp_m = pz_m[half * m_obs:(half + 1) * m_obs, :] * ys + ym
results["m80"] = yp_m
else:
results["m80"] = np.zeros_like(y_train[:, :a_test.shape[1]])
return results
def run_single_diameter(diam: float, scene_data: dict) -> dict:
"""Run validation for one diameter. Returns results dict."""
tgt_key = f"target_cylinder_{diam}L"
ill_key = f"illusion_{diam}L"
tgt_d = scene_data[tgt_key]
ill_d = scene_data[ill_key]
unc_d = scene_data["pinball"]
# Pre-build target-only field matrices
tgt_f = build_field_matrix(tgt_d["ux"], tgt_d["uy"])
ill_f = build_field_matrix(ill_d["ux"], ill_d["uy"])
unc_f = build_field_matrix(unc_d["ux"], unc_d["uy"])
diam_results = {}
pod_cache = {}
# Pre-compute target-only POD for each fold
for fold in range(N_CYCLES):
test_cyc = fold
train_cyc = [c for c in range(N_CYCLES) if c != test_cyc]
train_idx = sorted([c * N_PTS + p for c in train_cyc for p in range(N_PTS)])
for r in R_LIST:
Q_ref = tgt_f[:, train_idx]
mf = np.mean(Q_ref, axis=1)
U, _, _ = np.linalg.svd(Q_ref - mf[:, None], full_matrices=False)
pod_cache[("loco", fold, r)] = (mf, U[:, :r])
# Blocked split
train_idx_full = list(range(0, 2 * N_PTS))
for r in R_LIST:
Q_ref = tgt_f[:, train_idx_full]
mf = np.mean(Q_ref, axis=1)
U, _, _ = np.linalg.svd(Q_ref - mf[:, None], full_matrices=False)
pod_cache[("blocked", 0, r)] = (mf, U[:, :r])
# -- LOCO --
print("\n--- LOCO (4-fold) ---", flush=True)
loco_results = {}
for r in R_LIST:
for obs in ["force_fy", "force_fx", "action"]:
fold_r2_m1, fold_r2_m80 = [], []
for fold in range(N_CYCLES):
test_cyc = fold
train_cyc = [c for c in range(N_CYCLES) if c != test_cyc]
train_idx = sorted([c * N_PTS + p for c in train_cyc for p in range(N_PTS)])
test_idx = sorted([c * N_PTS + p for c in [test_cyc] for p in range(N_PTS)])
mf, modes_r = pod_cache[("loco", fold, r)]
for name, d, fld in [
(tgt_key, tgt_d, tgt_f), (ill_key, ill_d, ill_f), ("pinball", unc_d, unc_f)
]:
if obs == "action" and "illusion" not in name:
continue
if d.get("forces") is None and "force" in obs:
continue
a_train = modes_r.T @ (fld[:, train_idx] - mf[:, None])
a_test = modes_r.T @ (fld[:, test_idx] - mf[:, None])
if "force" in obs:
f_mode = obs.split("_")[1] # "fy" or "fx"
y_train = make_force_obs(d["forces"][train_idx], name, mode=f_mode)
y_test = make_force_obs(d["forces"][test_idx], name, mode=f_mode)
else:
y_train = d["actions"][train_idx, :].T
y_test = d["actions"][test_idx, :].T
W, sigma, Rmat, _, _, _ = compute_reduced_ccd(a_train, y_train, Q_delay=CCD_Q)
recon = reconstruct_observable(W, sigma, Rmat, a_test, y_train)
ch_m1 = [r2(y_test[c], recon["mode1"][c]) for c in range(y_test.shape[0])]
ch_m80 = [r2(y_test[c], recon["m80"][c]) for c in range(y_test.shape[0])]
fold_r2_m1.append(float(np.mean(ch_m1)))
fold_r2_m80.append(float(np.mean(ch_m80)))
if fold_r2_m1:
key = f"LOCO_{obs}_r{r}"
loco_results[key] = {
"mode1": {
"mean": float(np.mean(fold_r2_m1)),
"std": float(np.std(fold_r2_m1)),
},
"m80": {
"mean": float(np.mean(fold_r2_m80)),
"std": float(np.std(fold_r2_m80)),
},
}
print(f" {key}: R2_m1={loco_results[key]['mode1']['mean']:.4f}+-"
f"{loco_results[key]['mode1']['std']:.4f} "
f"R2_m80={loco_results[key]['m80']['mean']:.4f}+-"
f"{loco_results[key]['m80']['std']:.4f}", flush=True)
# -- Blocked split --
print("\n--- Blocked Split (train=0-47, test=48-95) ---", flush=True)
blocked = {}
test_idx_full = list(range(2 * N_PTS, 4 * N_PTS))
for r in R_LIST:
for obs in ["force_fy", "force_fx", "action"]:
mf, modes_r = pod_cache[("blocked", 0, r)]
per_case_m1, per_case_m80 = {}, {}
for name, d, fld in [
(tgt_key, tgt_d, tgt_f), (ill_key, ill_d, ill_f), ("pinball", unc_d, unc_f)
]:
if obs == "action" and "illusion" not in name:
continue
if d.get("forces") is None and "force" in obs:
continue
a_train = modes_r.T @ (fld[:, train_idx_full] - mf[:, None])
a_test = modes_r.T @ (fld[:, test_idx_full] - mf[:, None])
if "force" in obs:
f_mode = obs.split("_")[1]
y_train = make_force_obs(d["forces"][train_idx_full], name, mode=f_mode)
y_test = make_force_obs(d["forces"][test_idx_full], name, mode=f_mode)
else:
y_train = d["actions"][train_idx_full, :].T
y_test = d["actions"][test_idx_full, :].T
W, sigma, Rmat, _, _, _ = compute_reduced_ccd(a_train, y_train, Q_delay=CCD_Q)
recon = reconstruct_observable(W, sigma, Rmat, a_test, y_train)
ch_m1 = [r2(y_test[c], recon["mode1"][c]) for c in range(y_test.shape[0])]
ch_m80 = [r2(y_test[c], recon["m80"][c]) for c in range(y_test.shape[0])]
per_case_m1[name] = float(np.mean(ch_m1))
per_case_m80[name] = float(np.mean(ch_m80))
key = f"blocked_{obs}_r{r}"
blocked[key] = {
"mode1": {
"mean": float(np.mean(list(per_case_m1.values()))),
"per_case": per_case_m1,
},
"m80": {
"mean": float(np.mean(list(per_case_m80.values()))),
"per_case": per_case_m80,
},
}
print(f" {key}: R2_m1={blocked[key]['mode1']['mean']:.4f} "
f"R2_m80={blocked[key]['m80']['mean']:.4f}", flush=True)
diam_results["LOCO"] = loco_results
diam_results["blocked_split"] = blocked
return diam_results
def run():
print("=" * 60, flush=True)
print("CCD Validation (Round 5)", flush=True)
print("=" * 60, flush=True)
t_start = time.time()
all_results = {}
for diam in DIAMETERS:
tgt_key = f"target_cylinder_{diam}L"
ill_key = f"illusion_{diam}L"
print(f"\n{'=' * 60}", flush=True)
print(f"Diameter {diam}L", flush=True)
print(f"{'=' * 60}", flush=True)
t0 = time.time()
scene_data = {
tgt_key: load_aligned_fields(tgt_key),
ill_key: load_aligned_fields(ill_key),
"pinball": load_aligned_fields("pinball"),
}
print(f" Data loaded in {time.time() - t0:.0f}s", flush=True)
t1 = time.time()
all_results[f"{diam}L"] = run_single_diameter(diam, scene_data)
print(f" Analysis done in {time.time() - t1:.0f}s", flush=True)
out_dir = os.path.join(DATA_DIR, "ccd")
with open(os.path.join(out_dir, "validation_results.json"), "w") as f:
json.dump(all_results, f, indent=2)
print(f"\nTotal: {time.time() - t_start:.0f}s", flush=True)
print(f"Saved to {out_dir}/validation_results.json", flush=True)
if __name__ == "__main__":
run()

View File

@ -241,17 +241,24 @@ def compute_reduced_ccd(pod_coeffs, observable, Q_delay=6):
## 5. 结果
### 5.1 Correction-field CCD 主表
### 5.1 Correction-field CCD 主表2026-06-28 更新,统一几何后重跑)
| 指标 | 0.75L | 1.0L | 1.5L |
|------|-------|------|------|
| **O(dqctl, dqtar) mode1** | **0.564** | **0.913** | **0.667** |
| force_fy m80 | 2 | **1** (r=8/10) | 2 |
| action sigma1 | 1.39 | 1.13 | **0.28** |
| O(force, sig) tau=0 | **0.413** | **0.551** | — |
| O(force, sig) tau=tau_c | **0.806** | **0.768** | — |
| **O(dqctl, dqtar) mode1 (r=6)** | **0.383** | **0.926** | **0.922** |
| **O(dqctl, dqtar) mode1 (r=10)** | **0.320** | **0.684** | **0.661** |
| force_fy m80 (r=6) | 2 | 2 | **1** |
| action sigma1 (r=6) | 1.49 | 1.17 | **0.20** |
| Phase drift | low | low | **high** |
| Body-wake/sensor KE ratio | 0.73 | 1.17 | **2.58** |
| 1.5L special: rank sensitivity | — | — | O drops 0.922→0.661 (r=6→10) |
**关键变化**vs 6月15日旧版Illusion 旧几何 pinball x≈393px
- 0.75L O 从 0.564 → **0.383**-32%):旧几何的空间错位虚高了 overlap。统一几何后揭示真实匹配度远低于预期。
- 1.0L O 从 0.913 → **0.926**+1.4%基本不变1.0L 的控制修正与目标一致。
- 1.5L O(r=6)=**0.922**首次获得dominant mode 匹配度高,但更高阶 mode 快速发散r=10 时降至 0.661),反映多尺度控制策略。
- 1.5L action sigma1=**0.20**(远低于 0.75L 的 1.49 和 1.0L 的 1.17):确认高频调制机制下动作与流场结构的映射极其分散。
以下 force-sig overlap 和 zone 数据来自 6 月 15 日旧版(待用统一几何和新 zone 定义重跑):
### 5.2 三区域 Force-Signature Overlap
@ -284,8 +291,8 @@ def compute_reduced_ccd(pod_coeffs, observable, Q_delay=6):
|------|------------------|------|
| steady_cloak | 0.196 | 稳态,开环 |
| karman_re100 | 0.397 | 周期PPO 闭环 |
| vortex_lamb | 0.157 | 瞬态PPO 闭环 |
| vortex_taylor | 0.191 | 瞬态PPO 闭环 |
| vortex_lamb | 0.164 | 瞬态PPO 闭环fade-in/out + swapped norm |
| vortex_taylor | 0.203 | 瞬态PPO 闭环fade-in/out + swapped norm |
**结论**:不论上游条件如何(稳态/周期涡街/瞬态涡),控制策略的基本物理机制一致——"通过后两圆柱旋转补偿 pinball 阻塞引起的速度亏损,前端圆柱调节升力"。这与 SR 分析的结论完全吻合。
@ -299,7 +306,7 @@ def compute_reduced_ccd(pod_coeffs, observable, Q_delay=6):
**结论**:开环恒速旋转几乎无法抑制波动。需要闭环 DRL 控制。
### 5.5 Action-CCD Mode 1 可视化
### 5.5 Action-CCD Mode 1
Action-CCD 找出了控制器直接调制的主要结构。对 Cloak 场景action-CCD mode 1 ≈ dq_ctl mean field确认了"控制调制的结构 = correction-field 的主成分"的直觉。
@ -311,14 +318,13 @@ Action-CCD 找出了控制器直接调制的主要结构。对 Cloak 场景ac
所有位于 `data/figures/` 下(无 colorbar干净布局裁剪到 x=300-1100
| 图 | 内容 | 列数 |
|----|------|------|
| `corr_comparison_all_scenes.png` | 7 场景全景4 cloak + 3 illusion | 7 × 4 行 |
| `corr_illusion_comparison_dqctl.png` | Illusion 三直径对比 | 3 × 4 行 |
| `corr_cloak_comparison_dqctl.png` | Cloak 四场景对比 | 4 × 4 行 |
| `corr_illusion_{diam}L_dq_ctl/blk/tar_*.png` | 各场景单独 dq 图 | 3/1 panels |
| `corr_illusion_{diam}L_ctl_vs_tar.png` | dq_ctl vs dq_tar 对比 | 2×2 |
| `action_ccd_mode1_{scene}.png` | Action-CCD mode1 | 1×3 |
| 图 | 内容 |
|----|------|
| `corr_comparison_all_scenes.png` | 7 场景全景4 cloak + 3 illusion4 行 × 7 列 |
| `corr_cloak_comparison_dqctl.png` | Cloak 四场景 dq_ctl 对比 |
| `corr_illusion_comparison_dqctl.png` | Illusion 三直径 dq_ctl 对比 |
| `corr_{scene}_ctl_vs_tar.png` | 单场景 dq_ctl vs dq_tar 对比2×2 |
| `steady_cloak_cancel_test.png` | Steady cloak 抵消检验 |
### 6.2 全景对比图的读法
@ -389,33 +395,37 @@ conda run -n pycuda_3_10 python3 correction_analysis/compare_dqctl_scenes.py
```
src/CCD_analysis/
ccd_knowledge.md <-- 本文档唯一知识库
configs.py -- 场景元数据
ccd_knowledge.md ← 本文档(唯一知识库)
configs.py ← 场景元数据(统一几何)
README.md ← 快速入口
Lyu23.md ← CCD 方法文献
ccd/ ← Round 5 冻结基线(勿改)
utils/
resampling.py -- POD, CCD, 场加载
field_translate.py -- 场平移对齐
load_vortex_fields.py -- 瞬态 vortex 场加载
cfd_interface.py -- LegacyCelerisLab 封装 (GPU)
resampling.py POD, CCD, 场加载
field_translate.py ← 场平移(备用,不参与默认 pipeline
load_vortex_fields.py 瞬态 vortex 场加载
cfd_interface.py LegacyCelerisLab 封装 (GPU)
scripts/
detect_period.py -- 周期检测 → phase_plan.json
replay_fields.py -- 场回放 → fields_aligned.npz
collect_*.py -- GPU 数据采集
detect_period.py 周期检测 → phase_plan.json
replay_fields.py 场回放 → fields_aligned.npz
collect_*.py GPU 数据采集
correction_analysis/
compute_correction_fields.py -- correction-field 计算
diagnose_corrections.py -- 诊断图生成
compare_dqctl_scenes.py -- 多场景对比图
decompose_corrections.py -- 旧 force/action CCD
run_signature_line.py -- 旧 signature CCD
run_15L_correction.py -- 旧 1.5L 分析
run_zone_ccd.py -- 旧 zone CCD
run_steady_metrics.py -- 旧 steady cloak 度量
process_legacy_steady.py -- 旧格式加载
compute_correction_fields.py ← correction-field 计算
diagnose_corrections.py ← 诊断图生成
compare_dqctl_scenes.py ← 多场景对比图
decompose_corrections.py ← CCD 定量分解POD + force/action CCD
run_signature_line.py ← signature CCD
run_15L_correction.py ← 1.5L 专项分析
run_zone_ccd.py ← zone-restricted CCD
run_steady_metrics.py ← steady cloak 定量度量
visualize_action_ccd.py ← action-CCD mode1 可视化
process_legacy_steady.py ← 旧格式加载
data/
{scene_id}/{scene_name}/ -- 各场景数据
resampled/{scene}/ -- phase_plan.json
ccd/ -- JSON 结果文件
figures/ -- PNG 图(无 colorbar 版本)
old_data/ -- 归档废弃数据
{scene_id}/{scene_name}/ 各场景数据
resampled/ phase_plan.json
ccd/ JSON 结果文件
figures/ ← PNG 诊断图
old_data/ ← 归档(旧报告、旧脚本、旧 resampled 数据)
```
---
@ -449,8 +459,51 @@ src/CCD_analysis/
| 方向 | 状态 | 说明 |
|------|------|------|
| Karman cloak CCD 分析 | 数据已齐,分析延后 | 问题定义不同distortion compensation |
| 1.5L force-sig overlap | 待重跑 | SI=800 数据已有SI=200 高频采集待做 |
| SR-CCD-OID 映射 | 草稿 | 归档在 `data/old_data/sr_ccd_oid_mapping.md`,需根据各方向最终报告校正 |
| 1.5L 高频采样 | 待做 | 子步采集SI=200解析高频控制结构 |
| 统一几何后 CCD 重跑 | 待做 | `decompose_corrections.py` 需加 1.5L 后重跑;`correction_ccd_results.json` 仍为 6 月 15 日旧版 |
| 项目目录清理 | 已完成 (2026-06-28) | 移除根级 `old_data/`、`old_scripts/`、`output/`、`steady/`;废弃脚本归档至 `data/old_data/`;文档更新至统一几何 |
| Vortex 数据采集 Bug 修复 | 已完成 | 修复 cylinder order swapBOTTOM=id4, TOP=id5+ FIFO warmup 导致涡量消失;重采集 Taylor/Lamb |
| SR-CCD-OID 映射 | 草稿 | 归档 `data/old_data/sr_ccd_oid_mapping.md`,需最终校正 |
| 统一几何后 CCD 重跑 | 已完成 | `correction_ccd_results.json` 已更新(含 1.5L |
| Vortex 对比图更新 | 已完成 | `compare_dqctl_scenes.py` 已用修正数据重跑 |
| 项目目录清理 + 文档更新 | 已完成 | 移除废弃目录、SI=200 错误数据、旧诊断图 |
---
## 12. Vortex 采集 Bug 排查经验2026-06-29
`collect_vortex.py` 中发现了四个独立 bug按发现顺序
**Bug 1 — 圆柱顺序对调**
- 训练 env `legacy_env_vortex.py` 添加顺序front(id3) → TOP(+y, id4) → BOTTOM(-y, id5)
- action 映射:`temp[3]=front, temp[4]=TOP(bias=-4), temp[5]=BOTTOM(bias=+4)`
- 脚本先后把 TOP/BOTTOM 添加顺序写反 → 后两圆柱旋转方向全错
- 症状Lamb front 剧烈振荡 (std=0.28)dipole 对称性完全破坏
**Bug 2 — 涡量消失**
- FIFO warmup 在涡加入后跑 150×800=120000 步 → 涡从 x=15 漂出域外 (1280 lu)
- 症状:涡量场只有 pinball 尾流,完全看不到 vortex 结构
**Bug 3核心— 观测值归一化顺序错误**
- 训练 env 产出:`obs = [forces/force_norm, sensors/sens_norm]`force 先)
- 脚本将 channel 和 norm 互换:`[sensors/force_norm, forces/sens_norm]`
- 模型收到不匹配分布的反馈
- 症状Lamb cross-corr 仅 0.73, Taylor 仅 0.50;后圆柱同向转而非反向
**Bug 4 — fade-in/out 缺失**
- `uni_test.ipynb` 有 25 步渐入 + 25 步渐出到 steady-cloak bias
- 脚本直接给全量 PPO action无过渡
**最终修正方案**
1. Cylinder order 匹配训练 env
2. Bias 使用 uni_test 值:[-5, +5]FIFO[-5.1, +5.1]fade target
3. Obs 归一化:`forces_norm = obs[6:12]/force_norm``sens_norm = (obs[0:6]-sens_dev)/sens_norm``hstack([forces_norm, sens_norm])`
4. 25-step fade-in / 25-step fade-out
**修正后验证**
| 场景 | sim | cross-corr | active front mean | dq_ctl RMS |
|------|:---:|:----------:|:-----------------:|:----------:|
| vortex_lamb | 0.946 | 0.974 | 0.007 | 0.146 |
| vortex_taylor | 0.923 | 0.953 | -0.030 | 0.188 |
**其他脚本审计**`collect_karman.py` 和 `collect_illusion.py` 使用 `build_observation()` 函数,该函数正确实现了 force-first 归一化,无需修改。`collect_vortex.py` 是唯一手动构建 obs 的脚本。
**代码注释**`collect_vortex.py` 的文件头 docstring 和关键行号均有 BUG HISTORY 和 BUG-FIX 标记。

View File

@ -0,0 +1,188 @@
"""Generate comparison figures across all cloak & illusion scenarios.
All dq_ctl fields use unified geometry (pinball center at ~613px, sensors at ~800px),
set during GPU collection (configs.py UNIFIED coordinates).
Figures zoom into the region around the pinball/cylinder (x=300-1100) to exclude
boundary artifacts.
1. All-scenes panorama: steady_cloak, karman_re100, vortex_lamb, vortex_taylor,
illusion_0.75L, illusion_1.0L, illusion_1.5L
2. Illusion-only comparison: 0.75L, 1.0L, 1.5L
Usage:
conda run -n pycuda_3_10 python correction_analysis/compare_dqctl_scenes.py
"""
from __future__ import annotations
import os
import sys
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
_SRC = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", ".."))
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from CCD_analysis.configs import DATA_DIR, NX, NY
from CCD_analysis.correction_analysis.compute_correction_fields import (
compute_correction,
)
FIG_DIR = os.path.join(DATA_DIR, "figures")
os.makedirs(FIG_DIR, exist_ok=True)
# Display crop region (pixels) — around pinball at x~613
CROP_X0, CROP_X1 = 300, 1100
# Scene groups
CLOAK_SCENES = ["steady_cloak", "karman_re100", "vortex_lamb", "vortex_taylor"]
ILLUSION_SCENES = ["illusion_0.75L", "illusion_1.0L", "illusion_1.5L"]
ALL_SCENES = CLOAK_SCENES + ILLUSION_SCENES
SCENE_LABELS = {
"steady_cloak": "Steady Cloak",
"karman_re100": "Karman Cloak",
"vortex_lamb": "Vortex Lamb",
"vortex_taylor": "Vortex Taylor",
"illusion_0.75L": "Illusion 0.75L",
"illusion_1.0L": "Illusion 1.0L",
"illusion_1.5L": "Illusion 1.5L",
}
FIELD_METRICS = [
("ux_mean", r"mean $u_x$", "RdBu_r", True),
("uy_mean", r"mean $u_y$", "RdBu_r", True),
("rms", "RMS", "viridis", False),
("vorticity", r"$\omega_z$", "RdBu_r", True),
]
def compute_metrics(st: str) -> dict | None:
"""Load dq_ctl for a scene and compute metrics."""
try:
corr = compute_correction(st)
dq = corr.get("dq_ctl")
if dq is None:
return None
ux, uy = dq["ux"], dq["uy"]
return {
"ux_mean": np.mean(ux, axis=0),
"uy_mean": np.mean(uy, axis=0),
"rms": np.sqrt(np.std(ux, axis=0)**2 + np.std(uy, axis=0)**2),
"vorticity": np.gradient(np.mean(uy, axis=0), axis=1)
- np.gradient(np.mean(ux, axis=0), axis=0),
}
except Exception as e:
print(f" SKIP {st}: {e}")
return None
def crop_field(f: np.ndarray) -> np.ndarray:
"""Crop to display region (NY, NX_cropped)."""
return f[:, CROP_X0:CROP_X1]
def plot_comparison(scene_list: str | list, name: str):
"""Generate a grid of dq_ctl metrics for selected scenes."""
if isinstance(scene_list, str):
scene_list = [scene_list]
# Load all fields
fields = {}
for st in scene_list:
print(f" Loading {st}...", flush=True)
m = compute_metrics(st)
if m is not None:
fields[st] = m
n_scenes = len(fields)
if n_scenes == 0:
print(" No valid fields, skipping")
return
scene_names = list(fields.keys())
n_rows = len(FIELD_METRICS)
# Compute global vmax per metric from CROPPED fields
metric_vmax = {}
for mkey, _, _, _ in FIELD_METRICS:
all_vals = np.concatenate(
[abs(crop_field(fields[s][mkey])).ravel() for s in fields])
vmax = float(np.percentile(all_vals[np.isfinite(all_vals)], 99.5))
metric_vmax[mkey] = max(vmax, 1e-12)
fig, axes = plt.subplots(n_rows, n_scenes,
figsize=(3.5 * n_scenes, 3.0 * n_rows))
if n_rows == 1:
axes = [axes]
if n_scenes == 1:
axes = [[a] for a in axes]
nx_crop = CROP_X1 - CROP_X0
extent = (CROP_X0, CROP_X1, 0, NY - 1)
for row, (mkey, mlabel, cmap, symmetric) in enumerate(FIELD_METRICS):
for col, sn in enumerate(scene_names):
ax = axes[row][col]
f = crop_field(fields[sn][mkey])
vmax = metric_vmax[mkey]
kwargs = {"cmap": cmap, "origin": "lower",
"aspect": "equal", "extent": extent}
if symmetric:
kwargs["vmin"] = -vmax
kwargs["vmax"] = vmax
else:
kwargs["vmin"] = 0
kwargs["vmax"] = vmax
ax.imshow(f, **kwargs)
ax.tick_params(left=False, right=False, labelleft=False,
bottom=False, top=False, labelbottom=False)
if row == 0:
ax.set_title(SCENE_LABELS.get(sn, sn), fontsize=10)
if col == 0:
ax.set_ylabel(mlabel, fontsize=10)
plt.suptitle(f"dq_ctl: [{', '.join(SCENE_LABELS.get(s,s) for s in scene_names)}]",
fontsize=12, y=1.01)
plt.tight_layout()
path = os.path.join(FIG_DIR, name)
fig.savefig(path, dpi=150, bbox_inches="tight")
plt.close(fig)
print(f" Saved: {path}", flush=True)
# Stats
print(f"\n --- RMS (cropped region) ---")
for sn in scene_names:
rms_crop = crop_field(fields[sn]["rms"])
rms_val = float(np.sqrt(np.mean(rms_crop**2)))
print(f" {sn:22s}: RMS={rms_val:.6f}")
def main():
print("=" * 60)
print("Comparison: dq_ctl across all cloak & illusion scenes")
print("=" * 60)
# 1. All 7 scenes panorama
print("\n--- All 7 scenes panorama ---")
plot_comparison(ALL_SCENES, "corr_comparison_all_scenes.png")
# 2. Illusion-only (3 diameters)
print("\n--- Illusion-only comparison ---")
plot_comparison(ILLUSION_SCENES, "corr_illusion_comparison_dqctl.png")
# 3. Cloak-only (4 scenes) for reference
print("\n--- Cloak-only comparison ---")
plot_comparison(CLOAK_SCENES, "corr_cloak_comparison_dqctl.png")
print("\nDone!")
if __name__ == "__main__":
sys.exit(main())

View File

@ -30,7 +30,9 @@ from CCD_analysis.correction_analysis.compute_correction_fields import (
R_CANDIDATES = [6, 8, 10]
CCD_Q = 6
SCENE_TYPES = ["illusion_0.75L", "illusion_1.0L", "steady_cloak"]
SCENE_TYPES = ["illusion_0.75L", "illusion_1.0L", "illusion_1.5L", "steady_cloak"]
DIAMETERS_MAIN = [0.75, 1.0]
DIAMETER_SPECIAL = 1.5 # flagged as special_mechanism (high-freq modulation)
def compute_modal_overlap(W_dict, scene_label, r, obs_label="force_fy"):
@ -101,7 +103,9 @@ def main():
continue
diam = corr.get("diam")
print(f"\n--- {st} (diam={diam}) ---", flush=True)
is_special = (diam is not None and diam >= DIAMETER_SPECIAL)
flag = " [SPECIAL MECHANISM — high-freq modulation]" if is_special else ""
print(f"\n--- {st} (diam={diam}){flag} ---", flush=True)
Q_ctl = dict_to_field_matrix(dq_ctl)
N = Q_ctl.shape[1]
@ -146,6 +150,7 @@ def main():
"scene": st, "diam": diam, "obs": flabel, "r": r,
"m80": m80, "N": sig.size,
"sigma_top3": [float(sig[i]) for i in range(min(3,len(sig)))],
"special_mechanism": is_special,
}
if fmode == "fy":
print(f" {key}: m80={m80} s1={float(sig[0]):.4f}", flush=True)
@ -163,6 +168,7 @@ def main():
"scene": st, "diam": diam, "obs": "action", "r": r,
"m80": m80, "N": sig.size,
"sigma_top3": [float(sig[i]) for i in range(min(3,len(sig)))],
"special_mechanism": is_special,
}
print(f" {key}: m80={m80} s1={float(sig[0]):.4f}", flush=True)
@ -181,6 +187,7 @@ def main():
"m80": int(np.searchsorted(cumulative_energy(sig_t), 0.80)+1) if len(sig_t) > 0 else 0,
"N": sig_t.size,
"sigma_top3": [float(sig_t[i]) for i in range(min(3,len(sig_t)))],
"special_mechanism": is_special,
}
# Overlap: dq_ctl vs dq_tar
ck = f"{st}_dqctl_force_fy_r{r}"

View File

@ -47,46 +47,27 @@ SCENE_TYPES = [
# ---------------------------------------------------------------------------
# Three-zone masks
# Three-zone masks (unified geometry: pinball center at 613 px, sensors at 800 px)
# ---------------------------------------------------------------------------
def define_zones_illusion() -> dict:
"""Define three-zone masks for illusion layout (sensors at x=30*L0)."""
def define_zones() -> dict:
"""Define three-zone masks for all scenes (unified geometry, 2026-06-28).
All scenes now use the same pinball/sensor positions after unified collection.
Zone ranges:
near_body: 580-720 px (around pinball at x613)
body_wake: 720-850 px (near wake downstream)
sensor_zone: 780-850 px (around sensors at x=800)
"""
zones = {}
# Zone 1: near-body — envelope around pinball cylinders
# pinball front at x=380 (19*L0), rear at x=406 (20.3*L0)
# Extend to x=[350, 500], full height
mask = np.zeros((NY, NX), dtype=bool)
mask[:, 350:500] = True
zones["near_body"] = mask
# Zone 2: body-connected near wake — immediate downstream
mask = np.zeros((NY, NX), dtype=bool)
mask[:, 500:700] = True
zones["body_wake"] = mask
# Zone 3: downstream sensor zone — around sensors at x=600
mask = np.zeros((NY, NX), dtype=bool)
mask[:, 580:650] = True
zones["sensor_zone"] = mask
return zones
def define_zones_karman() -> dict:
"""Define three-zone masks for Karman layout (sensors at x=40*L0=800)."""
zones = {}
# Zone 1: near-body — around pinball at x=600
mask = np.zeros((NY, NX), dtype=bool)
mask[:, 580:720] = True
zones["near_body"] = mask
# Zone 2: body-connected near wake
mask = np.zeros((NY, NX), dtype=bool)
mask[:, 720:850] = True
zones["body_wake"] = mask
# Zone 3: downstream sensor zone — around sensors at x=800
mask = np.zeros((NY, NX), dtype=bool)
mask[:, 780:850] = True
zones["sensor_zone"] = mask
@ -272,8 +253,7 @@ def run():
print("Phase 2: Baseline Diagnostics (correction fields)", flush=True)
print("=" * 60, flush=True)
zones_ill = define_zones_illusion()
zones_karman = define_zones_karman()
zones = define_zones()
all_metrics = {}
@ -292,11 +272,7 @@ def run():
print(f" SKIP: no valid data (N=0)", flush=True)
continue
# Determine which zones to use
is_illusion = "illusion" in scene_type
is_karman = "karman" in scene_type or "vortex" in scene_type
zones = zones_ill if is_illusion else (zones_karman if is_karman else zones_ill)
# Unified geometry — same zones for all scenes
for dq_key, dq_label in [
("dq_blk", "dq_blk (pinball blockage)"),
("dq_ctl", "dq_ctl (control correction)"),
@ -313,8 +289,8 @@ def run():
metrics = zone_metrics(dq, zones, dq_label)
all_metrics[f"{scene_type}_{dq_key}"] = metrics
# For illusion/karman/vortex, also plot dq_tar if available
if dq_key == "dq_ctl" and (is_illusion or is_karman):
# For scenes with a target, also plot dq_tar if available
if dq_key == "dq_ctl" and corr.get("dq_tar") is not None:
dq_tar = corr.get("dq_tar")
if dq_tar is not None:
plot_mean_rms(dq_tar, "dq_tar (target correction)", prefix, zones)

View File

@ -0,0 +1,126 @@
"""Action-CCD mode 1 visualization for cloak scenes.
Action-CCD finds correction-field structures most correlated with cylinder
rotation speeds. For cloak scenes (steady/karman/vortex), this should reveal
the structures the controller directly modulates clean velocity deficit
compensation and cylinder dipoles, excluding upstream disturbance structures.
Usage:
conda run -n pycuda_3_10 python correction_analysis/visualize_action_ccd.py
"""
from __future__ import annotations
import os
import sys
import numpy as np
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
_SRC = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", ".."))
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from CCD_analysis.configs import DATA_DIR, NX, NY, L0
from CCD_analysis.utils.resampling import (
compute_pod, cumulative_energy, e95_index, compute_reduced_ccd,
unstack_velocity_modes,
)
from CCD_analysis.correction_analysis.compute_correction_fields import (
compute_correction, dict_to_field_matrix,
)
FIG_DIR = os.path.join(DATA_DIR, "figures")
os.makedirs(FIG_DIR, exist_ok=True)
CLOAK_SCENES = ["steady_cloak", "vortex_lamb", "vortex_taylor"]
# karman_re100 excluded due to 72 vs 96 frame mismatch
R = 10
CCD_Q = 6
CROP_X0, CROP_X1 = 300, 1100
def main():
print("=" * 60)
print("Action-CCD Mode 1: Cloak scenes")
print("=" * 60)
for st in CLOAK_SCENES:
print(f"\n--- {st} ---", flush=True)
# Load correction fields
corr = compute_correction(st)
dq_ctl = corr.get("dq_ctl")
if dq_ctl is None or dq_ctl.get("actions") is None:
print(f" SKIP: no dq_ctl or no actions")
continue
# Build snapshot matrix and compute POD
Q = dict_to_field_matrix(dq_ctl)
N = Q.shape[1]
mf, modes, sv, coeffs = compute_pod(Q)
e95 = e95_index(cumulative_energy(sv))
print(f" POD: E95={e95}, N_modes={len(sv)}")
# Action-CCD: find structures correlated with actions
a_r = coeffs[:R, :]
actions = dq_ctl["actions"][:N].T # (3, N)
W, sigma, _, _, _, _ = compute_reduced_ccd(a_r, actions, Q_delay=CCD_Q)
print(f" Action-CCD: sigma[0]={sigma[0]:.4f}, sigma_top3={sigma[:3]}")
# Reconstruct CCD mode 1 in physical space
# z1 = W[:, 0] @ A_z → CCD temporal coefficient
# CCD mode = sum over POD modes of (CCD direction weights * POD mode)
w1 = W[:, 0] / (np.linalg.norm(W[:, 0]) + 1e-12)
ccd_mode1 = modes[:, :R] @ w1 # (2*NX*NY,)
# Unstack into ux, uy
half = NX * NY
ux_mode = ccd_mode1[:half].reshape(NY, NX)
uy_mode = ccd_mode1[half:].reshape(NY, NX)
# Plot mode 1: ux + uy + vorticity, cropped
vor = np.gradient(uy_mode, axis=1) - np.gradient(ux_mode, axis=0)
fig, axes = plt.subplots(1, 3, figsize=(14, 4))
extent = (CROP_X0, CROP_X1, 0, NY - 1)
# ux
vmax_ux = max(abs(ux_mode).max(), 1e-12)
axes[0].imshow(ux_mode[:, CROP_X0:CROP_X1], cmap="RdBu_r",
vmin=-vmax_ux, vmax=vmax_ux,
origin="lower", aspect="equal", extent=extent)
axes[0].set_title(f"{st}: Action-CCD mode 1 ux")
# uy
vmax_uy = max(abs(uy_mode).max(), 1e-12)
axes[1].imshow(uy_mode[:, CROP_X0:CROP_X1], cmap="RdBu_r",
vmin=-vmax_uy, vmax=vmax_uy,
origin="lower", aspect="equal", extent=extent)
axes[1].set_title(f"{st}: Action-CCD mode 1 uy")
# vorticity
vmax_vor = max(np.percentile(abs(vor), 99), 1e-12)
axes[2].imshow(vor[:, CROP_X0:CROP_X1], cmap="RdBu_r",
vmin=-vmax_vor, vmax=vmax_vor,
origin="lower", aspect="equal", extent=extent)
axes[2].set_title(f"{st}: Action-CCD mode 1 vorticity")
for ax in axes:
ax.tick_params(left=False, right=False, labelleft=False,
bottom=False, top=False, labelbottom=False)
plt.tight_layout()
path = os.path.join(FIG_DIR, f"action_ccd_mode1_{st}.png")
fig.savefig(path, dpi=150, bbox_inches="tight")
plt.close(fig)
print(f" Saved: {path}")
print("\nDone!")
if __name__ == "__main__":
sys.exit(main())

View File

@ -0,0 +1,509 @@
"""Collect vortex cloak data for CCD correction-field analysis.
Collects field snapshots for three scene types per vortex type:
- vortex_target_{type}: vortex only (no pinball), the "ideal" target flow
- vortex_uncontrolled_{type}: vortex + pinball, zero control (q_blk)
- vortex_{type}: vortex + pinball + PPO control (q_ctl)
All use LegacyCelerisLab (matching existing CCD data convention).
Usage:
conda run -n pycuda_3_10 python scripts/collect_vortex.py \\
--type lamb --device 2
conda run -n pycuda_3_10 python scripts/collect_vortex.py \\
--type taylor --device 2
conda run -n pycuda_3_10 python scripts/collect_vortex.py \\
--type all --device 2
--- BUG HISTORY (2026-06-29) ---
Three independent bugs in this script were discovered and fixed after
Lamb dipole showed non-physical front cylinder oscillation:
BUG 1 - Cylinder order swap (lines 166-168, 254-257):
Original training env adds cylinders: front(id3) -> TOP(+y, id4) -> BOTTOM(-y, id5).
action = [aF, aT(+bias), aB(-bias)] -> temp[3]=aF, temp[4]=aT(+bias), temp[5]=aB(-bias).
We had TOP and BOTTOM swapped, so bias -4 went to the wrong cylinder.
Consequence: rear cylinders rotated in opposite directions; front overcompensated.
BUG 2 - Observation normalization swap (lines 368-371):
Training env produces obs = [forces/force_norm, sensors/sens_norm] (force-first).
We incorrectly fed [sensors/force_norm, forces/sens_norm] (channel-swapped + wrong norms).
Consequence: model received garbage feedback; on Lamb, similarity appeared 0.94 from
the similarity computation but cross-correlation was only 0.73.
BUG 3 - Missing fade-in/out transitions (lines 314-357):
uni_test.ipynb uses 25-step fade-in from steady-cloak bias [-5.1U0, +5.1U0] to PPO action,
and 25-step fade-out back to steady-cloak. We applied PPO action immediately at full scale.
Consequence: flow instability from abrupt control changes.
Fix summary:
- Cylinder order: TOP(+y) at id4, BOTTOM(-y) at id5.
- Bias: [-5, +5] U0 for both Lamb and Taylor (matching uni_test).
- Observation: forces_norm=obs[6:12]/force_norm_fact, sens_norm=(obs[0:6]-sens_dev)/sens_norm_fact.
Assembled as hstack([forces_norm, sens_norm]) force first.
- Fade-in/out: 25 steps linear interpolation between steady_bias and PPO action.
Comparison with other collect scripts:
collect_karman.py and collect_illusion.py use build_observation() from cfd_interface.py
which does force-first correctly. Only this script built obs manually.
"""
from __future__ import annotations
import argparse
import json
import os
import sys
import time
from collections import deque
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
from CCD_analysis.configs import (
get_scene, data_dir_for_scene, model_path_for_scene,
LEGACY_CFG_DIR, L0, CENTER_Y, U0,
)
from CCD_analysis.utils.cfd_interface import (
load_legacy_configs, get_velocity_field,
load_ppo_model,
)
DATA_TYPE = np.float32
FIFO_LEN = 150
# ---------------------------------------------------------------------------
# Vortex configuration
# ---------------------------------------------------------------------------
_VORTEX_CFG = {
"lamb": {"vortex_type": "lamb", "vortex_strength": 0.5 * U0},
"taylor": {"vortex_type": "taylor", "vortex_strength": 0.03 * U0},
}
# ---------------------------------------------------------------------------
# Target: vortex only (no pinball), save fields.npz
# ---------------------------------------------------------------------------
def collect_target(vtype: str, device_id: int, out_dir: str,
n_steps: int = 150) -> dict:
"""Collect field snapshots for vortex-only target flow.
Records the vortex evolving through a sensor-only environment
at x=40*L0. Saves fields.npz and sensors.npz.
"""
cfg = _VORTEX_CFG[vtype]
cuda_cfg, field_cfg = load_legacy_configs(LEGACY_CFG_DIR)
field_cfg = field_cfg._replace(viscosity=float(0.004))
ff = FlowField(field_cfg, cuda_cfg, device_id=device_id)
ny = ff.FIELD_SHAPE[1]
n_sensors = 3
# Add 3 sensors at x=40*L0
sensor_positions = [2.0, 0.0, -2.0]
for y_off in sensor_positions:
ff.add_sensor((40.0 * L0, CENTER_Y + y_off * L0, 0.0), L0 / 4.0)
# Short stabilization (1xNX/U0 for vortex, not 4x)
stabilize_steps = int(1 * ff.FIELD_SHAPE[0] / U0)
ff.run(stabilize_steps, np.zeros(n_sensors, dtype=DATA_TYPE))
# Save clean flow DDF
ff.get_ddf()
ff.save_ddf()
# Add vortex at x=10*L0
vc = cfg["vortex_type"]
vs = cfg["vortex_strength"]
ff.add_vortex((10.0 * L0, CENTER_Y, 0.0),
2.0 * L0, vs, 0, vc)
# Record vortex evolution
sens_list, forc_list = [], []
ux_list, uy_list = [], []
for step in range(n_steps):
ff.run(800, np.zeros(n_sensors, dtype=DATA_TYPE))
obs = ff.obs.copy() # (n_sensors*2,) = 6 sensor channels
sens_list.append(obs)
forc_list.append(np.zeros(6, dtype=DATA_TYPE)) # placeholder
ux, uy = get_velocity_field(ff, u0=U0)
ux_list.append(ux)
uy_list.append(uy)
# Save
np.savez_compressed(os.path.join(out_dir, "fields.npz"),
ux=np.stack(ux_list), uy=np.stack(uy_list))
np.savez(os.path.join(out_dir, "sensors.npz"),
sensors=np.array(sens_list, dtype=np.float32),
forces=np.array(forc_list, dtype=np.float32))
del ff
return {"scene": f"vortex_target_{vtype}", "n_steps": n_steps}
# ---------------------------------------------------------------------------
# Uncontrolled: vortex + pinball, zero control (q_blk)
# ---------------------------------------------------------------------------
def collect_uncontrolled(vtype: str, device_id: int, out_dir: str,
n_steps: int = 150) -> dict:
"""Collect field snapshots for vortex + pinball with zero control.
Records the transient interaction of vortex with the pinball.
Target phases are recorded alongside field snapshots.
"""
cfg = _VORTEX_CFG[vtype]
cuda_cfg, field_cfg = load_legacy_configs(LEGACY_CFG_DIR)
field_cfg = field_cfg._replace(viscosity=float(0.004))
ff = FlowField(field_cfg, cuda_cfg, device_id=device_id)
ny = ff.FIELD_SHAPE[1]
n_sensors = 3
# ---- Phase 1: Sensors + target recording ----
for y_off in [2.0, 0.0, -2.0]:
ff.add_sensor((40.0 * L0, CENTER_Y + y_off * L0, 0.0), L0 / 4.0)
stabilize_steps = int(1 * ff.FIELD_SHAPE[0] / U0)
ff.run(stabilize_steps, np.zeros(n_sensors, dtype=DATA_TYPE))
# Save clean DDF (pre-pinball, pre-vortex)
ff.get_ddf()
ff.save_ddf()
# Record target (vortex only, for similarity reference)
ff.add_vortex((10.0 * L0, CENTER_Y, 0.0),
2.0 * L0, cfg["vortex_strength"], 0, cfg["vortex_type"])
target_states = np.empty((0, 6), dtype=DATA_TYPE)
for _ in range(min(FIFO_LEN, n_steps)):
ff.run(800, np.zeros(n_sensors, dtype=DATA_TYPE))
target_states = np.vstack((target_states, ff.obs.copy()))
np.savez(os.path.join(out_dir, "target.npz"), target_states=target_states)
# ---- Phase 2: Add pinball, record uncontrolled flow ----
ff.restore_ddf()
ff.apply_ddf()
# BUG-FIX (2026-06-29): cylinder order MUST match training env.
# Original env adds: front(3) -> TOP(+y,id4) -> BOTTOM(-y,id5).
# Swapping TOP/BOTTOM causes rear cylinders to rotate wrong direction.
ff.add_cylinder((30.0 * L0, CENTER_Y, 0.0), L0 / 2.0) # id 3: front
ff.add_cylinder((31.3 * L0, CENTER_Y + 0.75 * L0, 0.0), L0 / 2.0) # id 4: TOP (+y)
ff.add_cylinder((31.3 * L0, CENTER_Y - 0.75 * L0, 0.0), L0 / 2.0) # id 5: BOTTOM (-y)
n_obj = ff.obs.size // 2
assert n_obj == 6, f"Expected 6, got {n_obj}"
# Bias action stabilization (matching uni_test: [-5, +5] for both Lamb and Taylor)
ff.run(stabilize_steps, np.zeros(n_obj, dtype=DATA_TYPE))
ff.run(stabilize_steps, np.array([0.0, 0.0, 0.0, 0.0, -5.0 * U0, 5.0 * U0], dtype=DATA_TYPE))
# Add vortex at x=15*L0
ff.add_vortex((15.0 * L0, CENTER_Y, 0.0),
2.0 * L0, cfg["vortex_strength"], 0, cfg["vortex_type"])
# Record uncontrolled flow (zero actions on pinball)
sens_list, forc_list = [], []
ux_list, uy_list = [], []
for step in range(n_steps):
# Zero control on pinball
ff.run(800, np.zeros(n_obj, dtype=DATA_TYPE))
obs = ff.obs.copy()
sens_list.append(obs[0:6])
forc_list.append(obs[6:12])
ux, uy = get_velocity_field(ff, u0=U0)
ux_list.append(ux)
uy_list.append(uy)
np.savez_compressed(os.path.join(out_dir, "fields.npz"),
ux=np.stack(ux_list), uy=np.stack(uy_list))
np.savez(os.path.join(out_dir, "sensors.npz"),
sensors=np.array(sens_list, dtype=np.float32),
forces=np.array(forc_list, dtype=np.float32))
del ff
return {"scene": f"vortex_uncontrolled_{vtype}", "n_steps": n_steps}
# ---------------------------------------------------------------------------
# Controlled: vortex + pinball + PPO (q_ctl)
# ---------------------------------------------------------------------------
def collect_controlled(vtype: str, device_id: int, out_dir: str,
n_steps: int = 150) -> dict:
"""Collect field snapshots for vortex + pinball with PPO control.
Loads the trained PPO model, runs inference, saves fields and telemetry.
Also saves norm.json and ddf/fifo checkpoints for possible replay.
"""
scene_name = f"vortex_{vtype}"
cfg_src = get_scene(scene_name)
cfg_v = _VORTEX_CFG[vtype]
cuda_cfg, field_cfg = load_legacy_configs(LEGACY_CFG_DIR)
field_cfg = field_cfg._replace(viscosity=float(0.004))
ff = FlowField(field_cfg, cuda_cfg, device_id=device_id)
ny = ff.FIELD_SHAPE[1]
n_sensors = 3
# Save config
with open(os.path.join(out_dir, "config.json"), "w") as f:
json.dump({k: str(v) if not isinstance(v, (int, float, list, bool)) else v
for k, v in cfg_src.items()}, f, indent=2)
# ---- Phase 1: Sensor-only target recording with vortex ----
for y_off in [2.0, 0.0, -2.0]:
ff.add_sensor((40.0 * L0, CENTER_Y + y_off * L0, 0.0), L0 / 4.0)
stabilize_steps_short = int(1 * ff.FIELD_SHAPE[0] / U0)
ff.run(stabilize_steps_short, np.zeros(n_sensors, dtype=DATA_TYPE))
ff.get_ddf()
ff.save_ddf()
ff.add_vortex((10.0 * L0, CENTER_Y, 0.0),
2.0 * L0, cfg_v["vortex_strength"], 0, cfg_v["vortex_type"])
target_states = np.empty((0, 6), dtype=DATA_TYPE)
for _ in range(min(FIFO_LEN, n_steps)):
ff.run(800, np.zeros(n_sensors, dtype=DATA_TYPE))
target_states = np.vstack((target_states, ff.obs.copy()))
np.savez(os.path.join(out_dir, "target.npz"), target_states=target_states)
# ---- Phase 2: Add pinball, compute norm ----
ff.restore_ddf()
ff.apply_ddf()
# Object order MUST match training env: front(id3), TOP(+y,id4), BOTTOM(-y,id5)
ff.add_cylinder((30.0 * L0, CENTER_Y, 0.0), L0 / 2.0) # id 3: front
ff.add_cylinder((31.3 * L0, CENTER_Y + 0.75 * L0, 0.0), L0 / 2.0) # id 4: TOP (+y)
ff.add_cylinder((31.3 * L0, CENTER_Y - 0.75 * L0, 0.0), L0 / 2.0) # id 5: BOTTOM (-y)
n_obj = ff.obs.size // 2
assert n_obj == 6, f"Expected 6, got {n_obj}"
# Stabilize with zero action, then bias action matching training env
ff.run(stabilize_steps_short, np.zeros(n_obj, dtype=DATA_TYPE))
ff.run(stabilize_steps_short, np.array([0.0, 0.0, 0.0, 0.0, -5.0 * U0, 5.0 * U0], dtype=DATA_TYPE))
# Add vortex at x=15*L0
ff.add_vortex((15.0 * L0, CENTER_Y, 0.0),
2.0 * L0, cfg_v["vortex_strength"], 0, cfg_v["vortex_type"])
# Save DDF checkpoint (vortex at x=15, pinball stabilized with bias)
ff.get_ddf()
ff.save_ddf()
# Norm collection (zero action on vortex+pinball)
fifo = deque(maxlen=FIFO_LEN)
for _ in range(FIFO_LEN):
ff.run(800, np.zeros(n_obj, dtype=DATA_TYPE))
fifo.append(ff.obs.copy())
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])))
norm = {
"force_norm_fact": force_norm_fact,
"sens_deviation": sens_deviation.tolist(),
"sens_norm_fact": sens_norm_fact.tolist(),
"action_bias": [0.0, -4.0, 4.0],
}
with open(os.path.join(out_dir, "norm.json"), "w") as f:
json.dump(norm, f, indent=2)
# Bias FIFO init (restore DDF so vortex starts from x=15)
ff.restore_ddf()
ff.apply_ddf()
fifo.clear()
# BUG-FIX (2026-06-29): use uni_test bias [-5,+5], not training env bias [-4,+4].
# The model was tested with [-5,+5] in uni_test.ipynb.
bias_fifo_arr = np.array([0.0, 0.0, 0.0, 0.0, -5.0 * U0, 5.0 * U0], dtype=DATA_TYPE)
for _ in range(FIFO_LEN):
ff.run(800, bias_fifo_arr)
fifo.append(ff.obs.copy())
save_states = np.array(list(fifo), dtype=DATA_TYPE)
# Restore DDF back to vortex-at-x=15 checkpoint
ff.restore_ddf()
ff.apply_ddf()
# Save checkpoints for replay (vortex at initial position)
np.save(os.path.join(out_dir, "ddf_checkpoint.npy"), ff.ddf)
np.save(os.path.join(out_dir, "fifo_checkpoint.npy"), save_states)
# ---- Phase 3: Controlled PPO inference (fade-in/out matching uni_test) ----
model_path = model_path_for_scene(scene_name)
if model_path is None:
raise FileNotFoundError(f"No model found for scene: {scene_name}")
model = load_ppo_model(model_path, device=f"cuda:{device_id}", s_dim=12)
model.set_random_seed(0)
# Restore DDF to vortex-at-x=15 checkpoint
ff.restore_ddf()
ff.apply_ddf()
# Start with zeros observation (matching uni_test)
obs = np.zeros(12, dtype=np.float32)
sens_list, forc_list, act_list, rew_list = [], [], [], []
ux_list = []
uy_list = []
fifo = deque(maxlen=FIFO_LEN)
# Steady-cloak bias for transition (matching uni_test: [-5.1, +5.1])
# BUG-FIX (2026-06-29): uni_test uses [-5.1,+5.1], not training env [-4,+4].
steady_bias = np.array([0.0, -5.1 * U0, 5.1 * U0], dtype=DATA_TYPE)
for step in range(n_steps):
action, _ = model.predict(obs, deterministic=True)
action = action.astype(np.float32).flatten()
act_list.append(action.copy())
temp_action = np.array(action * 4.0 + np.array([0.0, -4.0, 4.0]), dtype=DATA_TYPE)
# Fade-in (0-24), active (25-44), fade-out (45-69), steady-cloak (70+)
if step < 25:
w = step / 25.0
temp_val = temp_action * w * U0 + steady_bias * (1.0 - w)
elif 45 <= step < 70:
w = (step - 45) / 25.0
temp_val = temp_action * (1.0 - w) * U0 + steady_bias * w
elif step >= 70:
temp_val = steady_bias
else:
temp_val = temp_action * U0
temp = np.zeros(n_obj, dtype=DATA_TYPE)
temp[3:6] = temp_val
ff.context.push()
ff.run(800, temp)
ff.context.pop()
obs_slice = ff.obs.copy()
fifo.append(obs_slice)
sens_list.append(obs_slice[0:6])
forc_list.append(obs_slice[6:12])
# Build observation for next step
# BUG-FIX (2026-06-29): force-first with CORRECT norms.
# OLD (broken): sens_raw = obs[0:6]/force_norm, force_raw = (obs[6:12]-sens_dev)/sens_norm
# obs = hstack([sens_raw, force_raw]) ← completely wrong!
# NEW (correct): forces_norm = obs[6:12]/force_norm, sens_norm = (obs[0:6]-sens_dev)/sens_norm
# obs = hstack([forces_norm, sens_norm]) ← force first
# Training env produces: obs = hstack([forces/force_norm, sensors/sens_norm]).
# collect_karman and collect_illusion use build_observation() which does this correctly.
forces_norm = obs_slice[6:12] / force_norm_fact
sens_norm = (obs_slice[0:6] - sens_deviation) / sens_norm_fact
obs = np.clip(np.hstack([forces_norm, sens_norm]), -1.0, 1.0).astype(np.float32)
# Save field snapshot
ux, uy = get_velocity_field(ff, u0=U0)
ux_list.append(ux)
uy_list.append(uy)
# Save field snapshots
np.savez_compressed(os.path.join(out_dir, "fields.npz"),
ux=np.stack(ux_list), uy=np.stack(uy_list))
# Save telemetry
np.savez(os.path.join(out_dir, "controlled.npz"),
sensors=np.array(sens_list, dtype=np.float32),
forces=np.array(forc_list, dtype=np.float32),
actions=np.array(act_list, dtype=np.float32))
# Compute similarity
from CCD_analysis.utils.cfd_interface import compute_similarity
states_arr = np.array(sens_list, dtype=np.float32)
n_align = min(states_arr.shape[0], target_states.shape[0])
if n_align >= 30:
sim = compute_similarity(target_states, states_arr[:n_align], 30)
else:
sim = 0.0
result = {"scene": scene_name, "similarity": float(sim), "n_steps": n_steps}
with open(os.path.join(out_dir, "result.json"), "w") as f:
json.dump(result, f, indent=2)
del ff, model
return result
# ---------------------------------------------------------------------------
# Main
# ---------------------------------------------------------------------------
def main():
ap = argparse.ArgumentParser(description="Collect vortex cloak fields for CCD")
ap.add_argument("--type", type=str, default="lamb",
help='Vortex type: lamb, taylor, or "all"')
ap.add_argument("--device", type=int, default=2, help="GPU device ID")
ap.add_argument("--steps", type=int, default=150,
help="Number of recording steps (max 150 for transient)")
ap.add_argument("--skip-target", action="store_true",
help="Skip vortex_target collection")
ap.add_argument("--skip-uncontrolled", action="store_true",
help="Skip uncontrolled collection")
ap.add_argument("--skip-controlled", action="store_true",
help="Skip controlled collection")
args = ap.parse_args()
if args.type.lower() == "all":
vtypes = ["lamb", "taylor"]
else:
vtypes = [args.type.lower()]
t_start = time.time()
for vtype in vtypes:
print(f"\n{'=' * 60}")
print(f"Vortex type: {vtype}")
print(f"{'=' * 60}")
# --- Target: vortex only ---
if not args.skip_target:
scene_name = f"vortex_target_{vtype}"
print(f"\n--- Collecting target: {scene_name} ---")
out_dir = data_dir_for_scene(scene_name)
r = collect_target(vtype, args.device, out_dir, args.steps)
print(f" Done: {r['scene']} -> {out_dir}")
# --- Uncontrolled: vortex + pinball, zero control ---
if not args.skip_uncontrolled:
scene_name = f"vortex_uncontrolled_{vtype}"
print(f"\n--- Collecting uncontrolled: {scene_name} ---")
out_dir = data_dir_for_scene(scene_name)
r = collect_uncontrolled(vtype, args.device, out_dir, args.steps)
print(f" Done: {r['scene']} -> {out_dir}")
# --- Controlled: vortex + pinball + PPO ---
if not args.skip_controlled:
scene_name = f"vortex_{vtype}"
print(f"\n--- Collecting controlled: {scene_name} ---")
out_dir = data_dir_for_scene(scene_name)
r = collect_controlled(vtype, args.device, out_dir, args.steps)
print(f" Done: {r['scene']} -> {out_dir} sim={r['similarity']:.4f}")
elapsed = time.time() - t_start
print(f"\nTotal time: {elapsed:.1f}s")
if __name__ == "__main__":
sys.exit(main())

View File

@ -0,0 +1,113 @@
"""Field translation (spatial shifting) utilities for CCD analysis.
When computing correction fields across scenes with different object positions
(e.g., illusion pinball at x=393 vs cloak pinball at x=613), fields must be
translated so that reference points align before subtraction.
All functions operate on field arrays with shape (NY, NX) or (N, NY, NX).
NX = 1280, NY = 512 in the standard grid.
"""
from __future__ import annotations
import numpy as np
def translate_field_x(field: np.ndarray, shift_x: int,
fill_edge: bool = True) -> np.ndarray:
"""Horizontally shift a field by a given number of pixels.
Parameters
----------
field : (NY, NX) or (N, NY, NX) ndarray
Velocity field(s) in legacy (NY, NX) order.
shift_x : int
Positive = shift right, negative = shift left.
fill_edge : bool
If True, fill vacated columns with edge values (smooth).
If False, fill with zeros (creates boundary artifacts).
Returns
-------
shifted : ndarray with same shape as input.
"""
if shift_x == 0:
return field.copy()
if field.ndim == 2:
ny, nx = field.shape
result = np.zeros_like(field)
shift = shift_x
if shift > 0:
result[:, shift:] = field[:, :-shift]
if fill_edge:
# Fill left vacated columns with leftmost column value
result[:, :shift] = field[:, :1]
else:
s = -shift
result[:, :-s] = field[:, s:]
if fill_edge:
# Fill right vacated columns with rightmost column value
result[:, -s:] = field[:, -1:]
return result
elif field.ndim == 3:
n, ny, nx = field.shape
result = np.zeros_like(field)
shift = shift_x
if shift > 0:
result[:, :, shift:] = field[:, :, :-shift]
if fill_edge:
# Broadcast leftmost column across vacated columns
result[:, :, :shift] = field[:, :, :1]
else:
s = -shift
result[:, :, :-s] = field[:, :, s:]
if fill_edge:
result[:, :, -s:] = field[:, :, -1:]
return result
else:
raise ValueError(f"Unsupported field ndim: {field.ndim}")
def translate_fields_dict(q: dict, shift_x: int) -> dict:
"""Translate ux and uy fields in a data dict, leaving telemetry unchanged."""
if shift_x == 0:
return q
return {
"ux": translate_field_x(q["ux"], shift_x),
"uy": translate_field_x(q["uy"], shift_x),
"forces": q.get("forces"),
"sensors": q.get("sensors"),
"actions": q.get("actions"),
"meta": {**q.get("meta", {}), "translated_by": shift_x},
"step_indices": q.get("step_indices"),
}
# ---------------------------------------------------------------------------
# Reference positions (pixel coordinates)
# ---------------------------------------------------------------------------
# All scenes now use UNIFIED geometry:
# pinball at (30, 31.3) x L0 -> center = 613 px
# sensors at 40 x L0
# target cylinder at 30.65 x L0 = 613 px (same as pinball center)
# field_translate kept for optional cross-comparison or future use.
ILLUSION_REF_X = 613 # UNIFIED: was 393
CLOAK_REF_X = 613 # unchanged
TARGET_CYL_REF_X = 613 # UNIFIED: was 400 (now same as pinball center)
def get_scene_ref_x(scene_name: str) -> int | None:
"""Get the reference x-position (pinball/cylinder center) for a scene.
Returns pixel coordinate, or None if no reference (e.g. target_channel).
"""
if "illusion" in scene_name:
return ILLUSION_REF_X
if "target_cylinder" in scene_name:
return TARGET_CYL_REF_X
if scene_name in ("pinball", "steady_cloak", "karman_re100", "karman_q_blk",
"vortex_lamb", "vortex_taylor",
"vortex_uncontrolled_lamb", "vortex_uncontrolled_taylor"):
return CLOAK_REF_X
return None

View File

@ -0,0 +1,114 @@
"""Load vortex scene fields.npz format for transient CCD analysis.
Vortex scenes (vortex_lamb, vortex_taylor, vortex_target_*, vortex_uncontrolled_*)
have field snapshots saved as raw fields.npz with shape (N, NX, NY) and no phase
plan (transient, not periodic).
Converts to the same convention as load_aligned_fields():
- Transposes fields from (N, NX, NY) -> (N, NY, NX)
- Loads telemetry from controlled.npz or sensors.npz
- Returns dict with identical key structure
"""
from __future__ import annotations
import json
import os
from typing import Any
import numpy as np
from CCD_analysis.configs import DATA_DIR, NX, NY, SCENES
def load_vortex_fields(scene_name: str) -> dict:
"""Load vortex scene field data from raw fields.npz.
Parameters
----------
scene_name : str one of the vortex scene names registered in configs.
Returns
-------
dict with same keys as load_aligned_fields():
ux, uy : (N, NY, NX) ndarray field snapshots (transposed)
forces : (N, 6) ndarray or None
sensors : (N, 6) ndarray or None
actions : (N, 3) ndarray or None
meta : dict with scene info
step_indices : list of int (sequential 0..N-1 for transient)
"""
if scene_name not in SCENES:
raise KeyError(f"Unknown scene: {scene_name}")
cfg = SCENES[scene_name]
scene_id = cfg["scene_id"]
data_dir = os.path.join(DATA_DIR, scene_id, scene_name)
if not os.path.isdir(data_dir):
raise FileNotFoundError(f"Vortex scene directory not found: {data_dir}")
# -- fields.npz (native simulation order: NX first) --
fields_path = os.path.join(data_dir, "fields.npz")
if not os.path.isfile(fields_path):
raise FileNotFoundError(f"{fields_path} not found")
fd = np.load(fields_path)
ux_raw = fd["ux"] # (N, NX, NY)
uy_raw = fd["uy"]
N = ux_raw.shape[0]
fd.close()
# Transpose (N, NX, NY) -> (N, NY, NX) to match load_aligned_fields convention
ux = np.ascontiguousarray(ux_raw.transpose(0, 2, 1))
uy = np.ascontiguousarray(uy_raw.transpose(0, 2, 1))
# -- Telemetry (controlled.npz or sensors.npz) --
tele_path = None
for p in [
os.path.join(data_dir, "controlled.npz"),
os.path.join(data_dir, "sensors.npz"),
]:
if os.path.isfile(p):
tele_path = p
break
sensors, forces, actions = None, None, None
if tele_path is not None:
td = np.load(tele_path)
if "sensors" in td:
sensors = td["sensors"] # (N, 6)
if "forces" in td:
forces = td["forces"] # (N, 6)
if "actions" in td:
actions = td["actions"] # (N, 3)
td.close()
# Verify N matches
if sensors is not None and sensors.shape[0] != N:
raise ValueError(
f"sensors ({sensors.shape[0]}) != fields ({N})"
)
if forces is not None and forces.shape[0] != N:
raise ValueError(
f"forces ({forces.shape[0]}) != fields ({N})"
)
meta = {
"scene": scene_name,
"scene_id": scene_id,
"source": "vortex_transient",
"target_type": "transient",
"n_frames": N,
}
result: dict[str, Any] = {
"ux": ux,
"uy": uy,
"forces": forces,
"actions": actions,
"sensors": sensors,
"meta": meta,
"step_indices": list(range(N)), # sequential — no phase plan
}
return result

View File

@ -4,10 +4,10 @@
"force-oid_m2": 0.43535277661654437,
"force-oid_m3": 0.43535277661654437,
"force-oid_m5": 0.43535277661654437,
"sig-oid_m1": 0.017340158959890415,
"sig-oid_m2": 0.30180744104636803,
"sig-oid_m3": 0.09790190212866957,
"sig-oid_m5": 0.02065845627823741,
"sig-oid_m1": -0.15117050299050566,
"sig-oid_m2": 0.3261349291477396,
"sig-oid_m3": 0.21184489728386308,
"sig-oid_m5": 0.025177590823196313,
"sig-pcd_m1": -0.03521104267963627,
"sig-pcd_m2": 0.20732705633252407,
"sig-pcd_m3": 0.11988520616706469,
@ -18,21 +18,21 @@
"pod_m5": -3.0459602065033202
},
"future_sig": {
"force-oid_m1": 0.013511471239959334,
"force-oid_m2": 0.07098337174417249,
"force-oid_m3": 0.07098337174417249,
"force-oid_m5": 0.07098337174417249,
"sig-oid_m1": 0.3740715508827751,
"sig-oid_m2": 0.6608883811088201,
"sig-oid_m3": 0.5592259563419594,
"sig-oid_m5": 0.533343435056657,
"sig-pcd_m1": 0.20205641028110888,
"sig-pcd_m2": 0.4672590946761527,
"sig-pcd_m3": 0.4468990482184305,
"sig-pcd_m5": 0.41968732641466205,
"pod_m1": -0.2540973280602146,
"pod_m2": -0.0339567217960513,
"pod_m3": 0.054785729407538376,
"pod_m5": 0.3000378545113639
"force-oid_m1": 0.02808159096495584,
"force-oid_m2": 0.16866913136184078,
"force-oid_m3": 0.16866913136184078,
"force-oid_m5": 0.16866913136184078,
"sig-oid_m1": 0.33101756431711893,
"sig-oid_m2": 0.5494672616793316,
"sig-oid_m3": 0.530079016273497,
"sig-oid_m5": 0.32401702949087063,
"sig-pcd_m1": 0.16309737465172808,
"sig-pcd_m2": 0.5183971881220766,
"sig-pcd_m3": 0.5204063038036019,
"sig-pcd_m5": 0.4265691501067544,
"pod_m1": -0.052450576630182974,
"pod_m2": 0.04555274170153097,
"pod_m3": 0.0926186828625095,
"pod_m5": 0.33872788235759366
}
}

View File

@ -4,10 +4,10 @@
"force-oid_m2": 0.6705941647225692,
"force-oid_m3": 0.6705941647225692,
"force-oid_m5": 0.6705941647225692,
"sig-oid_m1": -2.7646669027021566,
"sig-oid_m2": -2.539151608216361,
"sig-oid_m3": -1.47692206321327,
"sig-oid_m5": -1.5110272636915942,
"sig-oid_m1": -2.590735492534273,
"sig-oid_m2": -1.746340572676032,
"sig-oid_m3": -0.505449315098804,
"sig-oid_m5": 0.20254103672707288,
"sig-pcd_m1": -1.6681874811094162,
"sig-pcd_m2": -1.342076853642359,
"sig-pcd_m3": 0.04249559758899473,
@ -18,21 +18,21 @@
"pod_m5": -0.09978939220221528
},
"future_sig": {
"force-oid_m1": -0.688583875677999,
"force-oid_m2": 0.0977498946249901,
"force-oid_m3": 0.0977498946249901,
"force-oid_m5": 0.0977498946249901,
"sig-oid_m1": 0.3400013732837159,
"sig-oid_m2": 0.5855599713349928,
"sig-oid_m3": 0.6757301882995801,
"sig-oid_m5": 0.6051731015609549,
"sig-pcd_m1": -0.045583209960261946,
"sig-pcd_m2": -0.07349660070560707,
"sig-pcd_m3": 0.534579564047348,
"sig-pcd_m5": 0.6365887267870092,
"pod_m1": -0.3737860307570295,
"pod_m2": -0.1596051084511593,
"pod_m3": 0.08266261398865987,
"pod_m5": -0.33155756162258626
"force-oid_m1": 0.045801892403349205,
"force-oid_m2": 0.2629802484043238,
"force-oid_m3": 0.2629802484043238,
"force-oid_m5": 0.2629802484043238,
"sig-oid_m1": 0.49884985568656176,
"sig-oid_m2": 0.7228144172245711,
"sig-oid_m3": 0.800731534358897,
"sig-oid_m5": 0.7557779433603652,
"sig-pcd_m1": 0.15914778050377293,
"sig-pcd_m2": 0.41216776425952434,
"sig-pcd_m3": 0.7363453012643849,
"sig-pcd_m5": 0.8203723624009522,
"pod_m1": -0.21134397189909435,
"pod_m2": 0.3003051304179969,
"pod_m3": 0.46937683250351936,
"pod_m5": 0.38547893664057714
}
}

View File

@ -4,10 +4,10 @@
"force-oid_m2": 0.6397818250190341,
"force-oid_m3": 0.6397818250190341,
"force-oid_m5": 0.6397818250190341,
"sig-oid_m1": 0.5371119596459986,
"sig-oid_m2": 0.5689626851549741,
"sig-oid_m3": 0.5480702090166246,
"sig-oid_m5": 0.49764490473273426,
"sig-oid_m1": -0.052064205105397096,
"sig-oid_m2": 0.5476475353492796,
"sig-oid_m3": 0.5104783631016963,
"sig-oid_m5": 0.5006514705153909,
"sig-pcd_m1": 0.02922230950403174,
"sig-pcd_m2": 0.4747650262191032,
"sig-pcd_m3": 0.5480885671190363,
@ -18,21 +18,21 @@
"pod_m5": 0.5163077019241664
},
"future_sig": {
"force-oid_m1": 0.25720592794565883,
"force-oid_m2": 0.07069504954059229,
"force-oid_m3": 0.07069504954059229,
"force-oid_m5": 0.07069504954059229,
"sig-oid_m1": 0.3378310203787158,
"sig-oid_m2": 0.3147990569733715,
"sig-oid_m3": 0.34429262568108926,
"sig-oid_m5": 0.33509730308714486,
"sig-pcd_m1": -0.002980254539846315,
"sig-pcd_m2": 0.35229352094431515,
"sig-pcd_m3": 0.3046193933085676,
"sig-pcd_m5": 0.3332866518761705,
"pod_m1": -0.01505551568170228,
"pod_m2": 0.05972906227204257,
"pod_m3": 0.050175749864584104,
"pod_m5": 0.2244801105502782
"force-oid_m1": 0.17421268011080493,
"force-oid_m2": 0.14691686742584273,
"force-oid_m3": 0.14691686742584273,
"force-oid_m5": 0.14691686742584273,
"sig-oid_m1": 0.5380880846698864,
"sig-oid_m2": 0.7415813766705678,
"sig-oid_m3": 0.7457945827838065,
"sig-oid_m5": 0.743474115106988,
"sig-pcd_m1": 0.45784877152463377,
"sig-pcd_m2": 0.6555116659772753,
"sig-pcd_m3": 0.6986463965973865,
"sig-pcd_m5": 0.7043777394599675,
"pod_m1": 0.43491763790324245,
"pod_m2": 0.5396007454438271,
"pod_m3": 0.5888838546919968,
"pod_m5": 0.675748834414234
}
}

View File

@ -1,21 +1,21 @@
{
"force": {
"force-oid_m1": 0.3973693481528069,
"force-oid_m2": 0.7503722371594272,
"force-oid_m3": 0.7503722371594272,
"force-oid_m5": 0.7503722371594272,
"sig-oid_m1": 0.047626492192117884,
"sig-oid_m2": -0.0899320785087113,
"sig-oid_m3": -0.06793031697290859,
"sig-oid_m5": 0.050723754778942164,
"sig-pcd_m1": -0.032869874061091105,
"sig-pcd_m2": -0.03470568581697929,
"sig-pcd_m3": -0.0024393867643178763,
"sig-pcd_m5": 0.20808527695100557,
"pod_m1": -0.028581812678008658,
"pod_m2": 0.41796895591108846,
"pod_m3": 0.3922853200314628,
"pod_m5": 0.5941700935980355
"force-oid_m1": 0.11266023422325092,
"force-oid_m2": 0.29502997457263447,
"force-oid_m3": 0.29502997457263447,
"force-oid_m5": 0.29502997457263447,
"sig-oid_m1": -0.04664427253668195,
"sig-oid_m2": -0.1264864272302698,
"sig-oid_m3": -0.5831767238197324,
"sig-oid_m5": -0.468171814213285,
"sig-pcd_m1": -0.01403430119734832,
"sig-pcd_m2": -0.10619108268969935,
"sig-pcd_m3": -0.4769605173029764,
"sig-pcd_m5": -0.5750085233764719,
"pod_m1": -0.2578359217654613,
"pod_m2": 0.06838543401441043,
"pod_m3": -0.04891972521400708,
"pod_m5": 0.14664737200706499
},
"future_sig": {
"force-oid_m1": 0.0,

View File

@ -2,13 +2,11 @@
"scene": "karman_re100",
"n_snapshots": 500,
"dof": 67200,
"ranks_computed": [
6,
8,
10,
12,
16
],
"energy_r10_5modes": 0.999034936679307,
"energy_r10_10modes": 1.0
"roi": {
"x_start": 400,
"x_end": 1000,
"y_start": 100,
"y_end": 400
},
"energy_r10_5modes": 0.9971600541895967
}

View File

@ -1,6 +1,6 @@
{
"obs_act_train": 0.9560973048210144,
"oid_m3_act_train": 0.22518758634108724,
"oid_m5_act_train": 0.22518758634108724,
"oid_force_act_train": 0.23275911247341877
"obs_act_train": 0.7495932579040527,
"oid_m3_act_train": 0.09605772105033679,
"oid_m5_act_train": 0.09605772105033679,
"oid_force_act_train": 0.11671916873232673
}

View File

@ -1,20 +1,20 @@
{
"force_norm_fact": 0.019301448483020067,
"force_norm_fact": 0.0192322488874197,
"sens_deviation": [
0.7973665595054626,
-0.12411565333604813,
0.25158512592315674,
-0.008040801621973515,
0.8032128214836121,
0.11381910741329193
0.8234996795654297,
-0.12610766291618347,
0.2525006830692291,
-0.0017381409415975213,
0.7810041308403015,
0.12196661531925201
],
"sens_norm_fact": [
2.679999828338623,
3.1077914237976074,
1.8671960830688477,
3.329291582107544,
2.9416096210479736,
3.267305374145508
3.3339574337005615,
3.35296630859375,
1.8934848308563232,
3.4226977825164795,
3.0502705574035645,
3.111417293548584
],
"action_bias": [
0.0,

View File

@ -1,6 +1,7 @@
{
"scene_id": "karman",
"re_code": 100,
"nu": 0.004,
"has_disturbance": true,
"sample_interval": 800,
"action_scale": 8.0,
@ -15,6 +16,5 @@
"pinball_rear_x": 31.3,
"target_type": "periodic",
"s_dim": 12,
"u0": 0.01,
"nu": 0.004
"u0": 0.01
}

View File

@ -1,5 +1,5 @@
{
"scene": "karman_re100",
"similarity": 0.9587061844662661,
"avg_reward": 0.6636279821395874
"n_steps": 500,
"similarity": 0.0
}

144
src/SR_analysis/PIPELINE.md Normal file
View File

@ -0,0 +1,144 @@
# SR Analysis Pipeline
> Symbolic regression pipeline for extracting interpretable DRL control laws (obs → act) from the fluidic pinball.
> Four independent stages: inference → fitting → validation → analysis.
## Pipeline Architecture
```
[PPO 推理] [PySR 拟合] [CFD 验证] [分析/画图]
stage_1_infer.py stage_2_fit.py stage_3_validate.py stage_4_analyze.py
↓ ↓ ↓ ↓
controlled.npz formulas/*.json validations/*.json data/figures/*.png
```
## Quick Start
```bash
# 1. Generate PPO data
conda run -n pycuda_3_10 python stage_1_infer.py --scene karman_re100 --device 2
# 2. Fit PySR formula
conda run -n sr_env python stage_2_fit.py --scene karman_re100 --mode per-scene
# 3. Validate in CFD
conda run -n pycuda_3_10 python stage_3_validate.py \
--scene karman_re100 --device 2 --mode pysr \
--formula-front results/formulas/karman_re100_front.json \
--formula-top results/formulas/karman_re100_top.json
# 4. Analyze results
conda run -n pycuda_3_10 python stage_4_analyze.py --scene karman_re100 --mode ppo-viz
```
## Environments
| Env | Used For |
|-----|----------|
| `pycuda_3_10` | Stage 1 (CFD inference), Stage 3 (CFD validation), Stage 4 (analysis) |
| `sr_env` | Stage 2 (PySR symbolic regression) |
GPU: device 2 recommended (device 0 may conflict with PyTorch).
## Key Conventions
### Reynolds Number
- Code Re uses reference length 2D = 40: `Re = U0 * 40 / nu`
- Physical Re_D uses D = 20: `Re_D = Re / 2`
- Default: Re_code=100 → Re_D=50, nu=0.004
### Action
- `controlled.npz` stores actions as **normalized [-1, +1]** (not physical omega)
- Physical omega: `omega = (action * scale + bias) * U0`, then divided by radius for angular velocity
- Fitting target: **non-dimensional alpha = omega / U0** (not omega)
### Action Bias vs FIFO Bias
- **DRL action decoder bias**: `action * scale + bias` → physical omega. Karman: [0,-4,4], Illusion: [0,-2,2], Vortex: [0,-4,4]
- **FIFO initialization bias** (environment warmup): different values! Illusion FIFO uses [0, -U0, U0] (1U scale), not [0, -2U0, 2U0]
### Inlet
- Parabolic velocity profile (not uniform). Top/bottom walls are no-slip bounce-back.
- U0 = 0.01 at centerline (lattice units)
### G-mirror
- Correct: `[aF, aT, aB] → [-aF, -aB, -aT]` (not the old buggy version)
- v23 structure: Front no-bias (α_F = 0 when features = 0), rear shared-head (α_B = -Top∘G)
### Norm
- Each scene computes its own force/sensor normalization during environment initialization
- Must use the **same norm values** during inference and during validation
- Norm values are saved in `norm.json` in each data directory
### Sample Interval & Steps
| Scene | SI | Validation Steps |
|-------|:--:|:----------------:|
| Karman | 800 | 160-200 |
| Illusion 0.75L | 400 | 320 |
| Illusion 1L | 600 | 214 |
| Illusion 1.5L | 800 | 160 |
| Vortex | 800 | 150 (transient) |
Rule: steps ≥ NX/U0/SI = 1280/0.01/SI ≈ 128000/SI
## Results Index
All results indexed in [`scene_registry.json`](scene_registry.json). Canonical formulas in `results/formulas/`, CFD validations in `results/validations/`.
### Canonical Formulas
| Formula File | Scene | Formula |
|-------------|-------|---------|
| `results/formulas/karman_joint_front.json` | Karman cross-Re (joint) | `daF_dt - 14.952*mu*Cl_tot` |
| `results/formulas/karman_joint_top.json` | Karman cross-Re (joint) | `3.414` (constant) |
| `results/formulas/illusion_joint_front.json` | Illusion joint (0.75L+1L) | `Cd_tot - (Cd_err + 5.428) - 0.00978*(du_a_dt + u_a)` |
| `results/formulas/illusion_joint_top.json` | Illusion joint (0.75L+1L) | `(Cd_err - (Cd_rear - Cl_err))*0.535 + 2.782` |
### Key CFD Results
| Scene | Formula | Similarity |
|-------|---------|:----------:|
| Karman cross-Re avg | Joint | 0.847 |
| Illusion 0.75L | Joint | 0.978 |
| Illusion 1L | Joint | 0.970 |
| Vortex lamb | Karman joint | 0.949 (exceeds PPO 0.942) |
| Illusion 0.6L | Joint (generalization) | 0.939 |
| Illusion 0.8L | Joint (generalization) | 0.908 |
| Illusion 1.2L | Joint (generalization) | 0.849 |
| Illusion 2L | Joint (generalization) | 0.676 |
### Feature Sets
| Name | Features | Dim | Used For |
|------|----------|:---:|----------|
| PHASE_STATE_KEYS | u_a, du_a_dt, Cl_tot, dCl_tot_dt, Cd_tot, Cd_rear | 6 | Karman per-Re |
| ILLUSION_PHASE_KEYS | phase-state + Cd_err, Cl_err, dCd_err_dt, dCl_err_dt | 10 | Illusion |
| PHYS_DADT | physics + daF_dt, daB_dt, daT_dt + mu | 17 | Karman joint/deep |
## Key Documentation
| File | Content |
|------|---------|
| `PIPELINE.md` | This file — overview, environment, conventions |
| `sindy_sr_knowledge.md` | Bug history, confirmed facts, known limitations |
| `sindy_sr_notes.md` | Task list, current status |
| `STAGE_1_INFER.md` | Stage 1: PPO data generation |
| `STAGE_2_FIT.md` | Stage 2: PySR fitting |
| `STAGE_3_VALIDATE.md` | Stage 3: CFD validation |
| `STAGE_4_ANALYZE.md` | Stage 4: Analysis & visualization |
| `docs/SR_analysis_report.md` | Full report (465+ lines) |
| `docs/illusion_joint_formula_analysis.md` | Illusion joint formula deep dive |
## Stage 0 Audit (2026-06-28)
All imports verified in both conda environments:
| Script | pycuda_3_10 | sr_env |
|--------|:-----------:|:------:|
| `stage_1_infer.py` (infer_karman / infer_illusion / infer_vortex) | OK | — |
| `stage_2_fit.py` (PySR) | — | OK |
| `stage_3_validate.py` (closed-loop) | OK | — |
| `stage_4_analyze.py` (analysis) | OK | — |
| `core/features.py` (feature_builder) | OK | OK |
| `core/cfd.py` (cfd_interface) | OK | — |
No broken imports. All dependencies available.

View File

@ -4,181 +4,81 @@ Extracts interpretable control laws (`obs -> act`) from DRL-trained policies for
fluidic pinball. Uses **PySR symbolic regression** on dimensionless physical features with
G-equivariant structural constraints (v23: front no-bias, rear shared-head).
## Current Results (2026-06-25)
## Quick Start
### Karman Cloak — Cross-Re Unified Formula
```bash
# 1. Generate PPO data
conda run -n pycuda_3_10 python stage_1_infer.py --scene karman_re100 --device 2
| Scene | Front Formula | Top Formula | CFD Closed-Loop |
|-------|--------------|-------------|:---------------:|
| Joint (Re50-400) | `daF_dt - 14.952*mu*Cl_tot` | `alpha_T = 3.414` (const) | **0.847 avg** |
| Re50 independent | PySR per-Re best | — | **0.895** |
| Re100 independent | PySR per-Re best | — | **0.888** |
| Re200 independent | PySR per-Re best | — | **0.916** |
| Re400 (SI=400 opt) | Joint formula | Joint formula | **0.819** |
# 2. Fit formula
conda run -n sr_env python stage_2_fit.py --scene karman_re100 --mode per-scene
### Illusion
# 3. Validate in CFD
conda run -n pycuda_3_10 python stage_3_validate.py \\
--scene karman_re100 --device 2 --mode pysr \\
--formula-front results/formulas/karman_joint_front.json \\
--formula-top results/formulas/karman_joint_top.json
| Scene | Front Formula | CFD Closed-Loop | % of PPO |
|-------|--------------|:---------------:|:--------:|
| 0.75L | `-0.169*(Cl_tot + dCl_tot_dt) - 1.240` | **0.979** | 100.7% |
| 1L | `(du_a_dt + u_a + 26.5)*0.0123` | **0.957** | 98.4% |
| **Joint (0.75L+1L)** | `target_Cd - 5.428 + 0.0098*(du_a_dt + u_a)` | **0.978 / 0.970** | — |
| 1.5L | High-freq periodic modulation (not SR-amenable) | — | — |
# 4. Analyze
conda run -n pycuda_3_10 python stage_4_analyze.py --scene karman_re100 --mode ppo-viz
```
**Key finding**: 0.75L and 1L formulas have fundamentally different skeletons (Cl_tot vs u_a
dominant). Joint formula still achieves excellent CFD results on both although the underlying
mechanisms differ.
## Pipeline Architecture
### Illusion Generalization (Joint Formula, No PPO)
```
stage_1_infer.py → stage_2_fit.py → stage_3_validate.py → stage_4_analyze.py
(PPO数据) (PySR拟合) (CFD闭环验证) (分析/画图)
```
| Diameter | Similarity | Notes |
|:--------:|:----------:|-------|
| 0.5L | 0.854 | Signal weak, noise-dominated |
| 0.6L | **0.939** | Generalizes well |
| 0.8L | **0.908** | Generalizes well |
| 1.2L | 0.849 | Begins to degrade |
| 1.5L | N/A | High-frequency regime, different mechanism |
| 2.0L | 0.676 | Degraded, near 1.5L regime |
## Directory Structure
Valid range: 0.6L-1.0L (similarity > 0.90).
### Vortex Cloak (Generalization)
Karman joint formula tested on vortex scenes (no retraining):
| Scene | Karman Joint Formula | PPO Baseline |
|-------|:-------------------:|:------------:|
| vortex_lamb | **0.949** | 0.942 |
| vortex_taylor | **0.905** | 0.916 |
```
SR_analysis/
PIPELINE.md # 总览文档 (入口)
README.md # 本文件
sindy_sr_knowledge.md # 知识库 (bugs, 事实, 结果)
sindy_sr_notes.md # 任务清单
scene_registry.json # 所有场景的规范结果索引
configs.py # 场景注册表
core/ # 共享工具库
features.py # 特征构建 (无量纲化+phase-state)
fitting.py # STLSQ拟合+特征矩阵
cfd.py # LegacyCelerisLab接口
g_operator.py # G-mirror变换
data/ # 运行时生成的.npz数据
results/
formulas/ # 规范公式JSON
validations/ # CFD闭环验证结果
archive/ # 归档的中间文件
stage_1_infer.py # Stage 1: 统一PPO推理入口
STAGE_1_INFER.md
stage_2_fit.py # Stage 2: 统一PySR拟合入口
STAGE_2_FIT.md
stage_3_validate.py # Stage 3: 统一CFD闭环验证入口
STAGE_3_VALIDATE.md
stage_4_analyze.py # Stage 4: 统一分析/画图入口
STAGE_4_ANALYZE.md
```
---
## Pipeline Overview
```
controlled.npz (PPO rollout)
|
v
compute_features() --> dimensionless physics features (ILLUSION_PHASE_KEYS, etc.)
|
v
PySR symbolic regression --> sparse interpretable formulas
|
v
CFD closed-loop validation --> final similarity score
```
### Key Design Decisions
## Key Design Decisions
1. **Feature levels**: Static (8-dim) -> Phase-state (6-dim) -> Illusion-phase (10-dim)
2. **Output target**: Non-dimensional alpha, not physical omega
3. **v23 structure**: Front no-bias, rear shared-head (Bottom = -Top(Gx))
4. **Final judge**: CFD closed-loop similarity, not one-step R2
---
## Directory Structure
```
SR_analysis/
configs.py # Scene metadata (Karman, Illusion, Vortex)
configs/legacy/ # Legacy CFD configs (config_cuda.json, config_flowfield.json)
utils/
__init__.py # Exports (no pycuda dependency)
feature_builder.py # Dimensionless features, G-operator, phase-state features
sindy_fitter.py # STLSQ fitting + feature matrices
cfd_interface.py # LegacyCelerisLab wrapper (requires pycuda_3_10)
g_operator.py # Equivariance diagnostics
data/ # Inference output data (controlled.npz, target.npz)
karman/ karman_re50..400/
illusion/ illusion_0.75L,1L,1.5L/
vortex/ vortex_lamb,taylor/
scripts/
infer_karman.py # PPO inference -> controlled.npz
infer_illusion.py # PPO inference -> controlled.npz
infer_vortex.py # PPO inference -> controlled.npz
gen_illusion_target.py # Target data generation for generalization scenes
visualize_ppo_illusion.py# PPO visualization with vorticity
sindy/
run_pysr.py # PySR symbolic regression (niter=40)
run_pysr_deep.py # Karman deep PySR (niter=120, Re independent + joint)
run_pysr_deep_illusion.py# Illusion deep+joint PySR (niter=120)
validate/
run_closed_loop.py # Karman closed-loop validator
run_closed_loop_illusion.py # Illusion closed-loop validator
run_closed_loop_vortex.py # Vortex closed-loop validator
run_closed_loop_re400_si.py # Karman re400 short-SI validator
predict_pysr.py # PySR formula sympy.lambdify wrapper
eval_rollout.py # Offline multi-step rollout evaluation
launch_pysr_validation.py # Batch CFD validation launcher
batch_illusion_generalization.sh# Batch generalization CFD validation
results/ # 136 JSON files — canonical + intermediate
results/README.md # Result file reference table
results/archive/ # Archived intermediate search attempts
```
---
## Usage
### PySR Symbolic Regression (conda: sr_env)
```bash
# Illusion
conda run -n sr_env python src/SR_analysis/sindy/run_pysr_deep_illusion.py --individual
# Karman deep (cross-Re independent + joint)
conda run -n sr_env python src/SR_analysis/sindy/run_pysr_deep.py --both
```
### CFD Closed-Loop Validation (conda: pycuda_3_10, GPU 1 or 2)
```bash
# Illusion PySR formula
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop_illusion.py \
--scene illusion_1L --device 2 --steps 320 --mode pysr \
--pysr-front validate/results/pysr_illusion_1L_front.json \
--pysr-top validate/results/pysr_illusion_1L_top.json
# Karman joint formula
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \
--scene karman_re100 --device 2 --steps 200 --mode pysr \
--pysr-front validate/results/karman_joint_deep_front.json \
--pysr-top validate/results/karman_joint_deep_top.json
# Vortex (generalization test)
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop_vortex.py \
--scene vortex_lamb --device 2 --steps 150 --mode pysr \
--pysr-front validate/results/karman_joint_deep_front.json \
--pysr-top validate/results/karman_joint_deep_top.json
```
### PPO Inference (generate controlled.npz)
```bash
conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_karman.py --re 100 --device 2
conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_illusion.py --diameter 1.0 --device 2
conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_vortex.py --type lamb --device 2
```
---
## Critical Reminders
- **actions.npz are normalized [-1,1]**, not physical omega. Convert: `(action * scale + bias) * u0`
- **PySR needs `sensors_raw`/`forces_raw`** passed to `compute_features()` or derivative features are zero
- **Output target must be alpha** (non-dim): `Y = actions_phys / u0`
- **One-step R2 high != closed-loop good** -- always validate in CFD
- **Controls must propagate**: steps >= NX/u0/SI (S=400->320, S=600->214, S=800->160)
- **FIFO bias != DRL action bias** for Illusion: FIFO=[0,-U0,U0], decode=[0,-2,2]*U0
- **Joint formula must be manually reviewed** for spurious terms (e.g. `daB_dt` is constant=0 at deployment)
---
## Key Documentation
| File | Content |
|------|---------|
| `src/SR_analysis/sindy_sr_knowledge.md` | Background knowledge, bug history, known pitfalls (for coder reference) |
| `src/SR_analysis/sindy_sr_notes.md` | Task list, phase breakdown, current status |
| `docs/SR_analysis_report.md` | **Single consolidated report** — all formulas, results, methodology, structural analysis |
| `docs/illusion_joint_formula_analysis.md` | Illusion joint formula deep dive — physical interpretation, generalization curve |
| `PIPELINE.md` | **Primary entry** — pipeline overview, environment, conventions |
| `sindy_sr_knowledge.md` | Bug history, confirmed facts, known limitations |
| `sindy_sr_notes.md` | Task list, current status |
| `docs/SR_analysis_report.md` | Full report (465+ lines) |
| `docs/illusion_joint_formula_analysis.md` | Illusion joint formula deep dive |
## Core Files (≤20)
`stage_1_infer.py`, `stage_2_fit.py`, `stage_3_validate.py`, `stage_4_analyze.py`, `configs.py`, `scene_registry.json`, `core/features.py`, `core/fitting.py`, `core/cfd.py`, `core/g_operator.py` + 8 docs.

View File

@ -0,0 +1,31 @@
# Stage 1: PPO Inference
Generates `controlled.npz` data by running trained PPO models in LegacyCelerisLab CFD.
## Usage
```bash
conda run -n pycuda_3_10 python stage_1_infer.py --scene karman_re100 --device 2
conda run -n pycuda_3_10 python stage_1_infer.py --group karman_trained --device 2
conda run -n pycuda_3_10 python stage_1_infer.py --scene illusion_0.6L --target-only --device 2
```
## Scene Groups
karman_trained (re50-400), illusion_trained (0.75L/1L/1.5L), illusion_generalization (0.5L-2L), vortex_all.
## Output per Scene
`data/{scene_id}/{scene_name}/`: target.npz, controlled.npz, config.json, norm.json, result.json, target_harmonics.json (Illusion only).
actions in controlled.npz are normalized [-1,+1]. Physical omega = (action*scale+bias)*U0.
## Per-Scene Notes
- **Karman**: 7 objects, obs_slice=(2,14), action_bias=[0,-4,4], s_dim=12
- **Illusion**: 6 objects, obs_slice=(0,12), action_bias=[0,-2,2], s_dim=14, FIFO bias=[0,-U0,U0] differs from DRL bias
- **Vortex**: 6 objects, MAX_STEPS=150, action_scale=4, vortex added after DDF checkpoint
## Expected Similarities
Karman re100 ~0.90, Illusion 1L ~0.97, Illusion 0.75L ~0.98, Illusion 1.5L ~0.95

View File

@ -0,0 +1,45 @@
# Stage 2: PySR Symbolic Regression
Fits interpretable formulas (obs → act) from controlled.npz data using PySR.
## Usage
```bash
# Per-scene fitting
conda run -n sr_env python stage_2_fit.py --scene karman_re100 --mode per-scene
# Joint cross-scene
conda run -n sr_env python stage_2_fit.py --scenes karman_re50,karman_re100,karman_re200,karman_re400 --mode joint
# Deep search
conda run -n sr_env python stage_2_fit.py --scenes illusion_0.75L,illusion_1L --mode joint --deep
```
## Output
`results/formulas/{label}_{front,top}.json`: best_sympy formula, feature_keys, R2 score.
Fitting target is non-dimensional **alpha = omega/U0** (not physical omega).
## Feature Sets
- **PHASE_STATE_KEYS** (6): u_a, du_a_dt, Cl_tot, dCl_tot_dt, Cd_tot, Cd_rear — Karman per-Re
- **ILLUSION_PHASE_KEYS** (10): above + Cd_err, Cl_err, dCd_err_dt, dCl_err_dt — Illusion
- **PHYS_DADT + mu** (17): physics + daF/dt + mu modulation — Karman joint
## Per-Scene vs Joint
- **per-scene**: Fit one formula per scene. Use for individual Re/diameter analysis.
- **joint**: Concatenate multiple scenes' data, fit single formula. Use for cross-scene generalization.
## Formula Review Checklist
After fitting, check:
1. **No spurious terms**: `daB_dt` = 0 at deployment (rear constant), remove if present
2. **Front no-bias**: α_F ≈ 0 when all features ≈ 0
3. **Rear shared-head**: α_B = -α_T when flow is symmetric (G-mirror applied)
4. **One-step R2 ≠ closed-loop**: Always validate in Stage 3
## Environments
`conda run -n sr_env` (PySR installed, separate from pycuda to avoid CUDA conflicts).

View File

@ -0,0 +1,52 @@
# Stage 3: CFD Closed-Loop Validation
Validates PySR formulas or PPO baselines in closed-loop CFD. Final judge is DTW similarity (not one-step R2).
## Usage
```bash
# PySR formula validation
conda run -n pycuda_3_10 python stage_3_validate.py \
--scene karman_re100 --device 2 --mode pysr \
--formula-front results/formulas/karman_joint_front.json \
--formula-top results/formulas/karman_joint_top.json
# PPO baseline
conda run -n pycuda_3_10 python stage_3_validate.py --scene illusion_1L --device 2 --mode ppo
# Batch generalization
conda run -n pycuda_3_10 python stage_3_validate.py \
--group illusion_generalization --device 2 --mode pysr \
--formula-front results/formulas/illusion_joint_front.json \
--formula-top results/formulas/illusion_joint_top.json
```
## Modes
- **pysr**: Deploy PySR formula in closed-loop CFD (v23: front no-bias, rear shared-head)
- **ppo**: Run trained PPO model as baseline
- **uncontrolled**: Zero-action baseline
## Steps
Auto-calculated from sample interval: SI=400→320, SI=600→214, SI=800→160.
## Output
`results/validations/{scene_name}.json`: similarity score, mode, n_steps.
## v23 Structure
- Front: α_F = f_front(x), no bias — should be ~0 when features ~0
- Top: α_T = f_rear(x), with bias
- Bottom: α_B = -f_rear(G[x]) — shared-head via G-mirror
G-mirror: [aF,aT,aB]→[-aF,-aB,-aT], sensors swap top↔bottom with v sign flip.
## Similarity Interpretation
- > 0.90: Excellent (cloak/illusion effective)
- 0.70-0.90: Partial control
- < 0.50: Poor
Tail similarity isolates the far-wake portion.

View File

@ -0,0 +1,38 @@
# Stage 4: Analysis & Visualization
Analyzes PPO policies and SR formulas. Generates FFT spectra, action timeseries, and degradation curves.
## Usage
```bash
# PPO action visualization (timeseries + FFT)
conda run -n pycuda_3_10 python stage_4_analyze.py --scene illusion_1L --mode ppo-viz
# Cross-diameter degradation analysis
conda run -n pycuda_3_10 python stage_4_analyze.py \
--scenes illusion_0.75L,illusion_1L,illusion_1.5L --mode degradation
# Formula comparison (coming soon)
conda run -n pycuda_3_10 python stage_4_analyze.py --scene illusion_1L --mode formula-compare
```
## Modes
- **ppo-viz**: Action timeseries plot + FFT spectrum for each cylinder
- **degradation**: Cross-scene comparison (alpha std, Cd_tot, action range) — finds transition points where control regime changes
- **formula-compare**: (placeholder) Compare PySR formula predictions vs PPO actions
## Output
Figures written to `data/figures/`:
- `ppo_viz_{scene}.png` — timeseries + FFT
- `degradation_metrics.png` — cross-diameter comparison panels
## Regime Detection
The degradation mode detects control regime transitions by:
1. Action amplitude jump (alpha std > 4x)
2. Frequency shift (FFT dominant frequency > 5x)
3. Autocorrelation pattern change (lag-2 ≈ -0.9 = high-frequency switching)
Known transition: Illusion 1.5L shifts from phase-lead compensation to high-frequency periodic modulation.

403
src/SR_analysis/core/cfd.py Normal file
View File

@ -0,0 +1,403 @@
"""CFD interface for LegacyCelerisLab (pycuda_3_10 env).
All functions use the LegacyCelerisLab (old) CFD API via:
from LegacyCelerisLab import FlowField
Must be run inside: conda run -n pycuda_3_10
NOTE: This module should be imported directly, not through SR_analysis.utils
because it requires pycuda. Other utils (sindy_fitter, feature_builder, g_operator)
do NOT require pycuda and can be imported from the __init__.
"""
from __future__ import annotations
import json
import os
import sys
from collections import deque
from typing import Any, Dict, List, Optional, Tuple
import numpy as np
# -- Import legacy CFD -------------------------------------------------------
# LegacyCelerisLab lives at the repo root; SR_analysis is at repo_root/src/SR_analysis.
_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 LegacyCelerisLab import utils as legacy_utils # noqa: E402
# ---------------------------------------------------------------------------
# Action-smoothing constant (legacy run() internal)
# ---------------------------------------------------------------------------
ACTION_SMOOTH_WEIGHT = 0.1 # used by FlowField.run() internally
def nu_from_re(re_code: float, u0: float = 0.01, d_ref: float = 40.0) -> float:
"""Return kinematic viscosity for a given code Reynolds number.
``re_code`` uses reference length *2*D* = 40.0 (matching model file naming).
"""
return u0 * d_ref / re_code
def load_legacy_configs(config_dir: str) -> Tuple[Any, Any]:
"""Load and return legacy (cuda_config, field_config) from *config_dir*."""
cuda_cfg = legacy_utils.load_cuda_config(
os.path.join(config_dir, "config_cuda.json")
)
field_cfg = legacy_utils.load_flow_field_config(
os.path.join(config_dir, "config_flowfield.json")
)
return cuda_cfg, field_cfg
# ---------------------------------------------------------------------------
# Environment helpers -- Karman cloak (disturbance cylinder + pinball)
# ---------------------------------------------------------------------------
def build_karman_cloak_env(
flow_field: FlowField,
*,
u0: float,
l0: float,
sample_interval: int,
fifo_len: int,
data_type: type,
) -> Tuple[np.ndarray, dict]:
"""Phase 0-1: add dist-cylinder & 3 sensors, stabilize, record target.
Steps (mirrors env_karman_cloak_standard.__init__):
1. add dist_cylinder (id=0)
2. add 3 sensors (id=1,2,3)
3. stabilize run(4*NX/U0, zero-action[4])
4. record FIFO_LEN x run(SAMPLE_INTERVAL, zero[4]), collect obs[2:8]
Returns
-------
target_states : ndarray (FIFO_LEN, 6) -- sensor0/1/2 ux,uy
info : dict with n_objects, NX, NY
"""
# dist cylinder
center = (10.0 * l0, (flow_field.FIELD_SHAPE[1] - 1) / 2, 0.0)
flow_field.add_cylinder(center, l0)
# sensors
for y_off in [2.0, 0.0, -2.0]:
sc = (40.0 * l0, (flow_field.FIELD_SHAPE[1] - 1) / 2 + y_off * l0, 0.0)
flow_field.add_sensor(sc, l0 / 4.0)
n_obj = flow_field.obs.size // 2
# stabilize
stabilize_steps = int(4 * flow_field.FIELD_SHAPE[0] / u0)
print(f" stabilising ({stabilize_steps} steps)...")
flow_field.run(stabilize_steps, np.zeros(n_obj, dtype=data_type))
# record target (only sensor signals = obs[2:8])
target_states = np.empty((0, 6), dtype=data_type)
for _ in range(fifo_len):
flow_field.run(sample_interval, np.zeros(n_obj, dtype=data_type))
new_state = flow_field.obs.copy()[2:8]
target_states = np.vstack((target_states, new_state))
print(f" target recorded: {target_states.shape}")
return target_states, {"n_objects": n_obj, "NX": flow_field.FIELD_SHAPE[0],
"NY": flow_field.FIELD_SHAPE[1]}
def add_pinball(
flow_field: FlowField,
*,
l0: float,
u0: float,
sample_interval: int,
fifo_len: int,
data_type: type,
action_bias: Optional[Tuple[float, float, float]] = None,
pinball_front_x: float = 30.0,
pinball_rear_x: float = 31.3,
obs_slice_start: int = 2,
obs_slice_end: int = 14,
n_objects_total: Optional[int] = None,
) -> dict:
"""Add pinball cylinders, stabilize, compute norm, bias rollout.
Steps:
1. add front, bottom, top cylinders
2. stabilize run(4*NX/U0, zero-action)
3. get_ddf() + save_ddf() (checkpoint)
4. FIFO_LEN x run(SAMPLE_INTERVAL, zero) -> compute norm
5. apply_ddf() (restore pre-bias state)
6. FIFO_LEN x run(SAMPLE_INTERVAL, bias-action) -> save_states
7. apply_ddf()
Parameters
----------
pinball_front_x, pinball_rear_x : pinball geometry (L0 units).
Default 30.0/31.3 for Karman; 19.0/20.3 for Illusion.
obs_slice_start, obs_slice_end : slice of obs for norm.
Default [2:14] for Karman (7 objects); [0:12] for Illusion (6 objects).
n_objects_total : if provided, used for bias array length.
Default: inferred from flow_field after adding cylinders.
Returns dict with norm values.
"""
if action_bias is None:
action_bias = (0.0, -4.0, 4.0)
u0_float = float(u0)
# add 3 pinball cylinders
ny = flow_field.FIELD_SHAPE[1]
centers = [
(pinball_front_x * l0, (ny - 1) / 2, 0.0),
(pinball_rear_x * l0, (ny - 1) / 2 + 0.75 * l0, 0.0),
(pinball_rear_x * l0, (ny - 1) / 2 - 0.75 * l0, 0.0),
]
for c in centers:
flow_field.add_cylinder(c, l0 / 2.0)
n_obj = flow_field.obs.size // 2 if n_objects_total is None else n_objects_total
print(f" bodies after pinball: {n_obj}")
# stabilize
stabilize_steps = int(4 * flow_field.FIELD_SHAPE[0] / u0_float)
print(f" stabilising pinball ({stabilize_steps} steps)...")
flow_field.run(stabilize_steps, np.zeros(n_obj, dtype=data_type))
# checkpoint DDF
flow_field.get_ddf()
flow_field.save_ddf()
# --- norm phase (zero-action) ---
fifo = deque(maxlen=fifo_len)
for _ in range(fifo_len):
flow_field.run(sample_interval, np.zeros(n_obj, dtype=data_type))
fifo.append(flow_field.obs.copy()[obs_slice_start:obs_slice_end])
temp_states = np.array(fifo, dtype=data_type)
# forces are at indices [6:12] relative to the slice end
force_start = obs_slice_end - obs_slice_start - 6
force_end = force_start + 6
force_norm_fact = 6.0 * float(np.max(np.abs(temp_states[:, force_start:force_end])))
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}")
print(f" norm: sens_deviation={sens_deviation}")
print(f" norm: sens_norm_fact={sens_norm_fact}")
# --- bias-action rollout ---
flow_field.apply_ddf()
bias = np.zeros(n_obj, dtype=data_type)
bias[n_obj - 3] = float(action_bias[0] * u0_float)
bias[n_obj - 2] = float(action_bias[1] * u0_float)
bias[n_obj - 1] = float(action_bias[2] * u0_float)
print(f" bias action: {bias}")
fifo.clear()
for _ in range(fifo_len):
flow_field.run(sample_interval, bias)
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 {
"force_norm_fact": force_norm_fact,
"sens_deviation": sens_deviation.tolist(),
"sens_norm_fact": sens_norm_fact.tolist(),
"action_bias": list(action_bias),
"save_states": save_states,
}
def build_observation(
obs_slice: np.ndarray,
norm: dict,
) -> np.ndarray:
"""Assemble normalised DRL observation (12-dim) from a single obs slice.
``obs_slice`` is 12-element: sensor[0:6] + force[6:12].
Returns clipped 12-dim array in [-1, 1].
"""
forces = obs_slice[6:12] / norm["force_norm_fact"]
sens = (obs_slice[0:6] - norm["sens_deviation"]) / norm["sens_norm_fact"]
obs = np.clip(np.hstack([forces, sens]), -1.0, 1.0).astype(np.float32)
return obs
def action_to_physical(
action_norm: np.ndarray,
*,
scale: float = 8.0,
bias: Tuple[float, float, float] = (0.0, -4.0, 4.0),
u0: float = 0.01,
) -> np.ndarray:
"""Convert normalized action [-1,1] to physical omega (lattice units).
physical_omega[i] = (action_norm[i] * scale + bias[i]) * u0
"""
action_norm = np.asarray(action_norm, dtype=np.float64).reshape(-1, 3)
bias_arr = np.array(bias, dtype=np.float64)
return (action_norm * scale + bias_arr) * u0
def scale_action(
action_norm: np.ndarray,
*,
scale: float = 8.0,
bias: Tuple[float, float, float] = (0.0, -4.0, 4.0),
u0: float = 0.01,
n_total_bodies: int = 7,
) -> np.ndarray:
"""Convert normalised action ([-1,1]^3) to legacy CFD action array.
Returns array of length *n_total_bodies* with cylinders' omegas at the
last 3 slots.
"""
a = np.zeros(n_total_bodies, dtype=np.float32)
omega = (np.array(action_norm, dtype=np.float32) * scale + np.array(bias, dtype=np.float32)) * u0
a[n_total_bodies - 3:] = omega
return a
# ---------------------------------------------------------------------------
# Vorticity & field export
# ---------------------------------------------------------------------------
def vorticity_from_ddf(flow_field: FlowField, u0: float) -> np.ndarray:
"""Compute z-vorticity from current DDF on host."""
flow_field.get_ddf()
ddf = flow_field.ddf.copy().reshape((9, flow_field.FIELD_SHAPE[1],
flow_field.FIELD_SHAPE[0])).transpose(2, 1, 0)
ux = (ddf[:, :, 1] + ddf[:, :, 5] + ddf[:, :, 8]
- 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=0) - np.gradient(ux, axis=1)
return omega.astype(np.float64)
def save_vorticity_png(path: str, omega: np.ndarray, title: str = ""):
"""Save vorticity field as a PNG with symmetric colour bar."""
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
abs_o = np.abs(omega[np.isfinite(omega)])
vmax = float(np.percentile(abs_o, 99.5)) if abs_o.size > 0 else 1.0
if vmax <= 0:
vmax = 1.0
ny, nx = omega.shape
fig, ax = plt.subplots(figsize=(min(18, max(8, nx / 60)), min(10, max(3, ny / 40))))
im = ax.imshow(omega, origin="lower", aspect="equal", cmap="RdBu_r",
vmin=-vmax, vmax=vmax, extent=(0, nx - 1, 0, ny - 1))
ax.set_xlabel("x (lattice)")
ax.set_ylabel("y (lattice)")
if title:
ax.set_title(title)
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label=r"$\omega_z$")
fig.tight_layout()
fig.savefig(path, dpi=150, bbox_inches="tight")
plt.close(fig)
# ---------------------------------------------------------------------------
# DTW similarity
# ---------------------------------------------------------------------------
def calc_lag(target: np.ndarray, state: np.ndarray) -> int:
"""Find lag that maximises cross-correlation between two 1-D signals."""
t = target - np.mean(target)
s = state - np.mean(state)
corr = np.correlate(t, s, mode="full")
lags = np.arange(-len(target) + 1, len(target))
return int(lags[np.argmax(corr)])
def calc_dtw_sim(target: np.ndarray, state: np.ndarray) -> float:
"""DTW-based similarity: 1 - (DTW distance / len(target))."""
n, m = len(target), len(state)
dtw = np.full((n + 1, m + 1), np.inf)
dtw[0, 0] = 0.0
for i in range(1, n + 1):
for j in range(1, m + 1):
cost = abs(float(target[i - 1]) - float(state[j - 1]))
dtw[i, j] = cost + min(dtw[i - 1, j], dtw[i, j - 1], dtw[i - 1, j - 1])
return float(1.0 - dtw[n, m] / n)
def compute_similarity(
target_states: np.ndarray,
state_series: np.ndarray,
conv_len: int,
) -> float:
"""Compute lag-compensated DTW similarity over *conv_len* window."""
ref = target_states[conv_len:2 * conv_len, 1]
cur = state_series[-conv_len:, 1]
lag = calc_lag(ref, cur)
sim_sum = 0.0
for i in range(6):
target_seq = np.roll(target_states[:, i], -lag)[conv_len:2 * conv_len]
state_seq = state_series[-conv_len:, i]
sim_sum += calc_dtw_sim(target_seq, state_seq) / 6.0
return float(sim_sum)
# ---------------------------------------------------------------------------
# Dummy env for loading SB3 models
# ---------------------------------------------------------------------------
def create_dummy_env(s_dim: int = 12, a_dim: int = 3):
"""Return a gym.Env with correct observation/action spaces for model loading."""
import gymnasium as gym
from gymnasium import spaces
class DummyEnv(gym.Env):
def __init__(self):
super().__init__()
self.observation_space = spaces.Box(low=-1, high=1, shape=(s_dim,), dtype=np.float32)
self.action_space = spaces.Box(low=-1, high=1, shape=(a_dim,), dtype=np.float32)
def reset(self, seed=None):
return np.zeros(s_dim, dtype=np.float32), {}
def step(self, action):
return np.zeros(s_dim, dtype=np.float32), 0.0, False, False, {}
def render(self):
pass
return DummyEnv()
def load_ppo_model(model_path: str, device: str = "cuda:0", s_dim: int = 12, a_dim: int = 3):
"""Load a PPO model with Sin activation."""
import torch
from torch.nn import Module
from stable_baselines3 import PPO
class Sin(Module):
def forward(self, x):
return torch.sin(x)
dummy_env = create_dummy_env(s_dim, a_dim)
model = PPO.load(model_path, env=dummy_env, device=device)
return model

View File

@ -0,0 +1,386 @@
"""Unified feature builder for all cloak scenes.
Produces dimensionless features with consistent G-equivariant structure.
All scenes (Karman, steady, vortex, illusion) use this same builder.
Copy of analysis_cloak/common/feature_builder.py -- kept as canonical source.
"""
from __future__ import annotations
from typing import Dict, List, Optional, Tuple
import numpy as np
# -- Physical constants ------------------------------------------------------
U0 = 0.01 # inlet velocity (lattice units)
D_CYL = 20.0 # cylinder diameter (lattice)
# -- Dimensionless conversion ------------------------------------------------
def compute_dimensionless(
sensors: np.ndarray, # (T, 6) raw lattice [s0_ux,s0_uy, s1_ux,s1_uy, s2_ux,s2_uy]
forces: np.ndarray, # (T, 6) raw lattice [f0_fx,f0_fy, f1_fx,f1_fy, f2_fx,f2_fy]
u0: float = U0,
d: float = D_CYL,
rho: float = 1.0,
) -> Dict[str, np.ndarray]:
"""Convert raw lattice CFD data to dimensionless quantities.
Sensor order: [s0_ux,s0_uy, s1_ux,s1_uy, s2_ux,s2_uy]
where s0=top(y=+2L0), s1=mid(y=0), s2=bottom(y=-2L0)
Force order: [front_fx,front_fy, bottom_fx,bottom_fy, top_fx,top_fy]
Returns:
u_hat_B, u_hat_C, u_hat_T: nondim streamwise velocity (bottom/centre/top)
v_hat_B, v_hat_C, v_hat_T: nondim crosswise velocity
Cd_F, Cd_T, Cd_B: drag coefficient per cylinder
Cl_F, Cl_T, Cl_B: lift coefficient per cylinder
"""
s = np.asarray(sensors, dtype=np.float64)
f = np.asarray(forces, dtype=np.float64)
# Sensor positions: s0=top, s1=centre, s2=bottom
# Convention: B=bottom=s2, C=centre=s1, T=top=s0
return {
"u_hat_T": s[:, 0] / u0,
"v_hat_T": s[:, 1] / u0,
"u_hat_C": s[:, 2] / u0,
"v_hat_C": s[:, 3] / u0,
"u_hat_B": s[:, 4] / u0,
"v_hat_B": s[:, 5] / u0,
"Cd_F": 2.0 * f[:, 0] / (rho * u0**2 * d),
"Cl_F": 2.0 * f[:, 1] / (rho * u0**2 * d),
"Cd_B": 2.0 * f[:, 2] / (rho * u0**2 * d),
"Cl_B": 2.0 * f[:, 3] / (rho * u0**2 * d),
"Cd_T": 2.0 * f[:, 4] / (rho * u0**2 * d),
"Cl_T": 2.0 * f[:, 5] / (rho * u0**2 * d),
}
# -- G operator (corrected) --------------------------------------------------
def apply_G_alpha(alpha: np.ndarray) -> np.ndarray:
"""Apply mirror G to action: [aF, aT, aB] -> [-aF, -aB, -aT]."""
return np.array([-alpha[0], -alpha[2], -alpha[1]], dtype=alpha.dtype)
def apply_G_x(dim: Dict[str, np.ndarray],
a_prev: np.ndarray,
a_prev2: np.ndarray) -> Tuple[Dict, np.ndarray, np.ndarray]:
"""Apply G to dimensionless state.
Returns (G_dim, G_a_prev, G_a_prev2) with corrected sign rules.
"""
G_dim = {
"u_hat_B": dim["u_hat_T"], "u_hat_C": dim["u_hat_C"], "u_hat_T": dim["u_hat_B"],
"v_hat_B": -dim["v_hat_T"], "v_hat_C": -dim["v_hat_C"], "v_hat_T": -dim["v_hat_B"],
"Cd_F": dim["Cd_F"], "Cd_T": dim["Cd_B"], "Cd_B": dim["Cd_T"],
"Cl_F": -dim["Cl_F"], "Cl_T": -dim["Cl_B"], "Cl_B": -dim["Cl_T"],
}
G_a_prev = np.column_stack([-a_prev[:, 0], -a_prev[:, 2], -a_prev[:, 1]])
G_a_prev2 = np.column_stack([-a_prev2[:, 0], -a_prev2[:, 2], -a_prev2[:, 1]])
return G_dim, G_a_prev, G_a_prev2
# -- Feature key definitions -------------------------------------------------
# Original feature set (includes sin_ua/cos_ua)
CORE_FEAT_KEYS = [
"u_m", "u_a", "u_c",
"v_a",
"Cd_tot", "Cd_rear",
"Cl_tot", "Cl_diff",
"sin_ua", "cos_ua",
"aF_lag1", "aB_lag1", "aT_lag1",
"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 -----------------------------------------------------
def compute_features(
dim: Dict[str, np.ndarray],
actions_prev: np.ndarray, # (T, 3) physical omega(t-1) or nondim alpha(t-1)
actions_prev2: np.ndarray, # (T, 3) physical omega(t-2)
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.
Args:
dim: from compute_dimensionless()
actions_prev: lagged actions (physical omega or nondim alpha)
actions_prev2: twice-lagged actions
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.
"""
T = actions_prev.shape[0]
u_B, u_C, u_T = dim["u_hat_B"], dim["u_hat_C"], dim["u_hat_T"]
v_B, v_C, v_T = dim["v_hat_B"], dim["v_hat_C"], dim["v_hat_T"]
Cd_F, Cd_T, Cd_B = dim["Cd_F"], dim["Cd_T"], dim["Cd_B"]
Cl_F, Cl_T, Cl_B = dim["Cl_F"], dim["Cl_T"], dim["Cl_B"]
# If actions are in physical omega, convert to nondim alpha
if alpha_mode:
a = actions_prev.astype(np.float64)
a2 = actions_prev2.astype(np.float64)
else:
a = actions_prev.astype(np.float64) / u0
a2 = actions_prev2.astype(np.float64) / u0
sym = {}
# Sensor combinations (nondim)
sym["u_m"] = (u_B + u_C + u_T) / 3.0
sym["u_a"] = (u_T - u_B) / 2.0
sym["u_c"] = u_C.copy()
sym["v_a"] = (v_T - v_B) / 2.0
# Force combinations (dimensionless Cd/Cl)
sym["Cd_tot"] = Cd_F + Cd_T + Cd_B
sym["Cd_rear"] = Cd_T + Cd_B
sym["Cl_tot"] = Cl_F + Cl_T + Cl_B
sym["Cl_diff"] = Cl_T - Cl_B
# 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) -- discrete version
sym["aF_lag1"] = a[:, 0]
sym["aB_lag1"] = a[:, 1]
sym["aT_lag1"] = a[:, 2]
sym["daF"] = a[:, 0] - a2[:, 0]
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)
sym["mu_u_a"] = sym["u_a"] * mu
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
def build_feature_matrix(
sym: Dict[str, np.ndarray],
feat_keys: List[str],
add_bias: bool = True,
) -> np.ndarray:
"""Build feature matrix (T, N) from symbol dict."""
cols = []
if add_bias:
cols.append(np.ones(sym[feat_keys[0]].shape[0], dtype=np.float64))
for k in feat_keys:
if k in sym:
cols.append(np.asarray(sym[k], dtype=np.float64))
else:
# Missing key -> zero
T = sym.get("u_m", np.ones(1)).shape[0]
cols.append(np.zeros(T, dtype=np.float64))
return np.column_stack(cols)
def get_feature_names(feat_keys: List[str], add_bias: bool = True) -> List[str]:
"""Get feature names matching build_feature_matrix output."""
names = []
if add_bias:
names.append("bias")
names.extend(feat_keys)
return names

View File

@ -0,0 +1,503 @@
"""SINDy fitting utilities: STLSQ threshold grid, feature matrix building.
All features are built using the unified feature_builder module.
"""
from __future__ import annotations
from typing import Dict, List, Optional, Tuple
import numpy as np
from .feature_builder import (
compute_dimensionless, compute_features, build_feature_matrix,
get_feature_names, ALL_FEAT_KEYS, U0,
)
# Default thresholds used across all scenes
DEFAULT_THRESHOLDS = [0.0, 0.001, 0.002, 0.005, 0.01, 0.015, 0.02, 0.03, 0.05, 0.1]
def fit_channel(
Theta: np.ndarray,
y: np.ndarray,
thresholds: Optional[List[float]] = None,
alpha: float = 1e-4,
max_iter: int = 25,
) -> Tuple[List[dict], dict]:
"""Fit a single channel (one cylinder) with STLSQ threshold grid.
Returns
-------
rows : list of dict per threshold
best : dict with best threshold entry (highest R2)
"""
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
best = None
rows = []
for th in thresholds:
opt = ps.STLSQ(threshold=th, alpha=alpha, max_iter=max_iter)
opt.fit(Theta_s, y)
coef = np.asarray(opt.coef_, dtype=np.float64).flatten() / std
y_pred = Theta @ coef
ssr = float(np.sum((y - y_pred) ** 2))
sst = float(np.sum((y - np.mean(y)) ** 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))
entry = {"threshold": float(th), "nz": nz, "r2": r2, "mae": mae, "coef": coef}
rows.append(entry)
if best is None or r2 > best["r2"]:
best = entry
return rows, best
def fit_sindy(
Theta: np.ndarray,
y: np.ndarray,
thresholds: Optional[List[float]] = None,
) -> List[dict]:
"""Run SINDy with threshold grid, return results list.
Each result dict has keys: threshold, nz, r2, mae, coef.
"""
if thresholds is None:
thresholds = DEFAULT_THRESHOLDS
std = np.std(Theta, axis=0)
std = np.where(std < 1e-8, 1.0, std)
Theta_s = Theta / std
results = []
for th in thresholds:
import pysindy as ps
opt = ps.STLSQ(threshold=th, alpha=1e-4, max_iter=25)
opt.fit(Theta_s, y)
coef = np.asarray(opt.coef_, dtype=np.float64).flatten() / std
y_pred = Theta @ coef
ssr = float(np.sum((y - y_pred) ** 2))
sst = float(np.sum((y - np.mean(y)) ** 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],
})
return results
def print_control_law(feature_names: List[str], coef: np.ndarray, channel_label: str = "ch"):
"""Pretty-print a sparse control law."""
terms = []
for i, c in enumerate(coef):
if abs(c) > 1e-8:
terms.append(f"{c:.6f} * {feature_names[i]}")
print(f" {channel_label}: {' + '.join(terms)}")
nz = sum(1 for c in coef if abs(c) > 1e-8)
print(f" non-zero terms: {nz}")
def get_active_support(
coef: np.ndarray,
feat_names: List[str],
relative_threshold: float = 0.02,
) -> Dict[str, float]:
"""Extract active features from coefficient vector.
Features with |coef| / max(|coef|) >= relative_threshold are considered active.
"""
max_c = np.max(np.abs(coef))
if max_c < 1e-12:
return {}
active = {}
for name, c in zip(feat_names, coef):
if abs(c) / max_c >= relative_threshold:
active[name] = float(c)
return active
def get_feature_matrix_from_data(
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 = True,
n_warmup: int = 2,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, List[str], List[str]]:
"""Build feature matrices from raw CFD data.
Constructs dimensionless features via feature_builder, creates front (no bias)
and rear (with bias) feature matrices, and returns them aligned with Y.
Parameters
----------
sensors, forces, actions_phys : raw data arrays.
mu : 1/Re_D.
u0 : inlet velocity (lattice units).
alpha_mode : if True, actions_phys are already nondim alpha.
include_mu : include mu modulation features.
n_warmup : number of warmup steps to discard (default 2 for lag/da).
Returns
-------
Theta_front : (T-warmup, N_front) feature matrix, NO bias column
Theta_rear : (T-warmup, N_rear) feature matrix, WITH bias column
Y : (T-warmup, 3) target action matrix
feat_names_front : list of feature names for front
feat_names_rear : list of feature names for rear
"""
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, u0=u0)
Theta_f = build_feature_matrix(sym, ALL_FEAT_KEYS, add_bias=False)
Theta_r = build_feature_matrix(sym, ALL_FEAT_KEYS, add_bias=True)
feat_names_front = get_feature_names(ALL_FEAT_KEYS, add_bias=False)
feat_names_rear = get_feature_names(ALL_FEAT_KEYS, add_bias=True)
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)

View File

@ -0,0 +1,191 @@
"""G-operator and equivariance tools.
Provides G-operator transformations, dimensionless conversion,
and equivariance diagnostics for PPO control laws.
"""
from __future__ import annotations
from typing import Any, Dict, Optional, Tuple
import numpy as np
from .feature_builder import compute_dimensionless as _compute_dimless
def apply_G_alpha(alpha: np.ndarray) -> np.ndarray:
"""Apply mirror G to action: [aF, aT, aB] -> [-aF, -aB, -aT]."""
return np.array([-alpha[0], -alpha[2], -alpha[1]], dtype=alpha.dtype)
def apply_G_raw(obs_slice: np.ndarray,
a_prev: np.ndarray,
a_prev2: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Apply G to raw obs slice [sensor(6)+force(6)] and action arrays.
Parameters
----------
obs_slice : (12,) raw obs [s0_ux,uy, s1_ux,uy, s2_ux,uy, f0_fx,fy, f1_fx,fy, f2_fx,fy]
a_prev : (3,) physical omega at t-1
a_prev2 : (3,) physical omega at t-2
Returns
-------
G_obs : (12,) transformed obs slice
G_a_prev : (3,) transformed a_prev
G_a_prev2 : (3,) transformed a_prev2
"""
G_obs = np.zeros(12, dtype=np.float64)
# sensors: swap top(0,1) <-> bottom(4,5), negate v
G_obs[0] = obs_slice[4]
G_obs[1] = -obs_slice[5]
G_obs[2] = obs_slice[2]
G_obs[3] = -obs_slice[3]
G_obs[4] = obs_slice[0]
G_obs[5] = -obs_slice[1]
# forces: swap bottom(2,3) <-> top(4,5), negate fy
G_obs[6] = obs_slice[6]
G_obs[7] = -obs_slice[7]
G_obs[8] = obs_slice[10]
G_obs[9] = -obs_slice[11]
G_obs[10] = obs_slice[8]
G_obs[11] = -obs_slice[9]
G_a_prev = np.array([-a_prev[0], -a_prev[2], -a_prev[1]], dtype=np.float64)
G_a_prev2 = np.array([-a_prev2[0], -a_prev2[2], -a_prev2[1]], dtype=np.float64)
return G_obs, G_a_prev, G_a_prev2
def check_equivariance(
model: Any,
obs_slice_series: np.ndarray, # (T, 12) raw obs
actions_phys: np.ndarray, # (T, 3) physical omega
norm: dict,
action_scale: float = 8.0,
action_bias: Tuple[float, float, float] = (0.0, -4.0, 4.0),
u0: float = 0.01,
) -> Dict[str, float]:
"""Check G-equivariance of a PPO model over a time series.
Returns dict with front/rear equivariance errors.
"""
from .cfd_interface import build_observation, action_to_physical
T = min(obs_slice_series.shape[0], actions_phys.shape[0])
ef, eb, et = [], [], []
for t in range(2, T):
# Get current obs
osl = obs_slice_series[t]
a_prev = actions_phys[t - 1] if t > 0 else actions_phys[0]
a_prev2 = actions_phys[t - 2] if t > 1 else actions_phys[0]
# Predict action for current state
obs = build_observation(osl, norm)
act, _ = model.predict(obs, deterministic=True)
act = act.astype(np.float32).flatten()
alpha = action_to_physical(act.reshape(1, 3),
scale=action_scale, bias=action_bias, u0=u0).flatten()
# Apply G to state
G_obs, _, _ = apply_G_raw(osl, a_prev, a_prev2)
obs_G = build_observation(G_obs, norm)
act_G, _ = model.predict(obs_G, deterministic=True)
act_G = act_G.astype(np.float32).flatten()
alpha_G = action_to_physical(act_G.reshape(1, 3),
scale=action_scale, bias=action_bias, u0=u0).flatten()
# Expected: G(alpha) = [-aF, -aB, -aT]
expected = apply_G_alpha(alpha)
ef.append(abs(float(alpha_G[0]) - float(expected[0])))
eb.append(abs(float(alpha_G[1]) - float(expected[1])))
et.append(abs(float(alpha_G[2]) - float(expected[2])))
ef_arr = np.array(ef)
eb_arr = np.array(eb)
et_arr = np.array(et)
alpha_range = float(np.max(np.abs(actions_phys[2:])))
return {
"front_mean_abs_error": float(np.mean(ef_arr)),
"front_rel_error": float(np.mean(ef_arr) / (alpha_range + 1e-12)),
"rear_bottom_rel_error": float(np.mean(eb_arr) / (alpha_range + 1e-12)),
"rear_top_rel_error": float(np.mean(et_arr) / (alpha_range + 1e-12)),
"alpha_range": alpha_range,
}
def diagnose_one_re(model, ff, target_states, norm, config, n_steps=150) -> dict:
"""Run PPO inference and check equivariance.
Parameters
----------
model : loaded PPO model
ff : FlowField instance (must be at saved checkpoint state)
target_states : (FIFO_LEN, 6) target sensor signals
norm : norm dict
config : scene config dict with action_scale, action_bias, u0, etc.
Returns
-------
dict with equivariance metrics.
"""
from collections import deque
from .cfd_interface import (build_observation, scale_action,
action_to_physical, compute_similarity)
action_scale = config.get("action_scale", 8.0)
action_bias = config.get("action_bias", (0.0, -4.0, 4.0))
u0 = config.get("u0", 0.01)
sample_interval = config.get("sample_interval", 800)
fifo_len = config.get("fifo_len", 150)
n_obj_total = config.get("n_objects_total", 7)
ff.restore_ddf()
ff.apply_ddf()
# Bias FIFO init
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])
# Inference
obs_array = []
action_array = []
obs = np.zeros(12, dtype=np.float32)
for _ in range(n_steps):
act, _ = model.predict(obs, deterministic=True)
act = act.astype(np.float32).flatten()
action_array.append(act.copy())
action_arr = scale_action(act, 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_slice = ff.obs.copy()[2:14]
fifo.append(obs_slice)
obs_array.append(obs_slice)
obs = build_observation(obs_slice, norm)
obs_series = np.array(obs_array, dtype=np.float64)
actions_phys = action_to_physical(np.array(action_array),
scale=action_scale, bias=action_bias, u0=u0)
states_arr = np.array(list(fifo), dtype=np.float32)
sim = compute_similarity(target_states, states_arr[:, 0:6],
config.get("conv_len", 30))
# Equivariance check
eq = check_equivariance(model, obs_series, actions_phys, norm,
action_scale, action_bias, u0)
return {
"similarity": sim,
"equivariance": eq,
}

Binary file not shown.

After

Width:  |  Height:  |  Size: 412 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 189 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 234 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 169 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 185 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 167 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 135 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 299 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 151 KiB

View File

@ -0,0 +1,50 @@
[
{
"dc": 0.0029594360229869684,
"amps": [
4.482508629307246e-06,
1.7559728832214555e-06,
1.695399969073208e-06,
4.863359419766791e-07,
4.0744896335315875e-07
],
"freqs": [
0.08666666666666667,
0.18000000000000002,
0.09333333333333334,
0.1,
0.08
],
"phases": [
-2.223160683015736,
-3.090016544658663,
-1.4396296733151448,
-1.3066393251376214,
1.6772756869554482
]
},
{
"dc": 1.3870471496678268e-05,
"amps": [
0.00028510443996025985,
0.0002599505102795033,
9.365825828740738e-05,
8.877059656510341e-05,
5.6681594533049977e-05
],
"freqs": [
0.04666666666666667,
0.04,
0.05333333333333334,
0.03333333333333333,
0.060000000000000005
],
"phases": [
-2.111421674988832,
0.931445743532081,
-2.026139500541211,
0.8174729094507539,
-1.9494987485100541
]
}
]

View File

@ -0,0 +1,50 @@
[
{
"dc": 0.0034837019443511963,
"amps": [
1.1376483917607967e-05,
7.248102177867125e-06,
5.838593601492066e-06,
3.099357973227853e-06,
2.721364398772006e-06
],
"freqs": [
0.08,
0.07333333333333333,
0.09333333333333334,
0.18000000000000002,
0.1
],
"phases": [
-1.0408841490876113,
1.9485158983384934,
-1.3904977741558484,
2.7148922736497454,
-1.2577909435667654
]
},
{
"dc": 2.0249660774425138e-05,
"amps": [
0.0006472307684973089,
0.0001580755533972986,
9.359867367951902e-05,
7.772244469336163e-05,
5.558302481685001e-05
],
"freqs": [
0.04,
0.03333333333333333,
0.04666666666666667,
0.02666666666666667,
0.02
],
"phases": [
3.026701709868741,
-0.09658000009988132,
3.0095201121212978,
-0.07819925121021988,
-0.05948277368734109
]
}
]

View File

@ -0,0 +1,50 @@
[
{
"dc": 0.004293269868940115,
"amps": [
3.8063591256693715e-05,
8.663235430632004e-06,
7.901108470774483e-06,
7.1510384113964075e-06,
6.272449978203118e-06
],
"freqs": [
0.06666666666666667,
0.09333333333333334,
0.08666666666666667,
0.18000000000000002,
0.006666666666666667
],
"phases": [
2.9599042697964997,
-1.7903030042984585,
1.1028484616237284,
1.731865132582606,
-2.5926944354040096
]
},
{
"dc": -1.1502189348296572e-06,
"amps": [
0.0011732713800496743,
8.715048457232978e-06,
8.113829684732725e-06,
2.805951507723943e-06,
2.2956073119401994e-06
],
"freqs": [
0.03333333333333333,
0.02666666666666667,
0.04,
0.02,
0.1
],
"phases": [
-1.225940564415875,
2.5874876850956996,
-1.9715232735875732,
-1.0313721340317885,
-2.2437335289130216
]
}
]

View File

@ -0,0 +1,50 @@
[
{
"dc": 0.004459296974043052,
"amps": [
4.365837863945567e-05,
1.6627461404015874e-05,
1.499327969764015e-05,
1.0882240141500203e-05,
8.648375991472251e-06
],
"freqs": [
0.09333333333333334,
0.1,
0.13333333333333333,
0.08666666666666667,
0.10666666666666667
],
"phases": [
-0.955645707880348,
2.2853500968785534,
0.39850354801733023,
-0.6758960563529793,
2.500473431154046
]
},
{
"dc": -9.51454075887644e-05,
"amps": [
0.0012646521264823782,
0.00021337454745985233,
0.00018204900807997863,
0.0001045493570771146,
9.486869355170616e-05
],
"freqs": [
0.04666666666666667,
0.05333333333333334,
0.04,
0.03333333333333333,
0.060000000000000005
],
"phases": [
3.0652144428661194,
-0.07281059956203381,
3.0763553688568286,
3.0963827460231994,
-0.08166873030294015
]
}
]

View File

@ -0,0 +1,50 @@
[
{
"dc": 0.006754336698601643,
"amps": [
0.0001205644768172709,
0.00010325204771730541,
6.1032211125949833e-05,
4.133092864305172e-05,
3.322629507619568e-05
],
"freqs": [
0.06666666666666667,
0.07333333333333333,
0.13333333333333333,
0.060000000000000005,
0.08
],
"phases": [
0.21557592230485592,
-2.886865293479548,
-1.017309664173338,
0.18116041037477737,
-2.833486698409445
]
},
{
"dc": -0.00011794015166136282,
"amps": [
0.0027813392111168526,
0.0008110970246664206,
0.0005849442269674323,
0.00035855985851303973,
0.00033022754730012617
],
"freqs": [
0.03333333333333333,
0.04,
0.02666666666666667,
0.02,
0.04666666666666667
],
"phases": [
-2.728322858690305,
0.49031498750676245,
-2.808700524748656,
-2.886553982431292,
0.565211864386275
]
}
]

View File

@ -0,0 +1,50 @@
[
{
"dc": 0.008447802712519964,
"amps": [
0.0002532245295392234,
0.00014147298112950578,
0.0001297400791941424,
6.434300108105502e-05,
5.7950849123282916e-05
],
"freqs": [
0.08,
0.07333333333333333,
0.18000000000000002,
0.08666666666666667,
0.36000000000000004
],
"phases": [
0.5998824737031915,
-2.5827112781492545,
-3.0945367379195816,
0.6391142016853608,
0.10751503340383528
]
},
{
"dc": 0.00011803933939518174,
"amps": [
0.004330276189305638,
0.000976254738854663,
0.0006186459734055383,
0.0004722765324672463,
0.000331904307497245
],
"freqs": [
0.04,
0.03333333333333333,
0.04666666666666667,
0.02666666666666667,
0.02
],
"phases": [
-2.551505303958517,
0.5007165195151991,
-2.466226000290465,
0.4069960130268731,
0.3090395499543929
]
}
]

View File

@ -0,0 +1,50 @@
[
{
"dc": 0.005645164093002677,
"amps": [
9.94330093591307e-05,
3.232916972111346e-05,
1.3215355231141013e-05,
1.1981542249045906e-05,
1.1909839518599968e-05
],
"freqs": [
0.08,
0.13333333333333333,
0.26666666666666666,
0.2733333333333334,
0.08666666666666667
],
"phases": [
1.7159609635464452,
-0.2507805590458965,
1.0324605783803387,
-1.6430097798845071,
-1.2507241172270498
]
},
{
"dc": -7.872594342188676e-06,
"amps": [
0.0020951717499817367,
0.00011698728994078175,
9.559502592412231e-05,
5.94451334704438e-05,
4.6248624934531934e-05
],
"freqs": [
0.04,
0.04666666666666667,
0.03333333333333333,
0.05333333333333334,
0.02666666666666667
],
"phases": [
-1.9167863100531028,
1.3060842750219288,
-2.0155661953529633,
1.3740336240558457,
-2.1415896469819797
]
}
]

View File

@ -0,0 +1,50 @@
[
{
"dc": 0.011301154401153327,
"amps": [
0.0005916813028667388,
0.0003228542347087601,
0.00011299686050540704,
0.00010933004657394343,
9.57751063779222e-05
],
"freqs": [
0.060000000000000005,
0.18000000000000002,
0.35333333333333333,
0.06666666666666667,
0.36000000000000004
],
"phases": [
-1.2029089114719096,
0.10667750417539099,
0.8128110625734376,
1.9355161635025424,
-1.7526560628277699
]
},
{
"dc": 0.000209219001474897,
"amps": [
0.00560199847829251,
0.0038799719967085497,
0.0017754993231967535,
0.0012951582037757647,
0.001094199761464694
],
"freqs": [
0.03333333333333333,
0.02666666666666667,
0.04,
0.02,
0.04666666666666667
],
"phases": [
-1.8630748640530403,
1.1749393426220804,
-1.7802277933394777,
1.028375164210408,
-1.7139048570401034
]
}
]

View File

@ -0,0 +1,24 @@
{
"force_norm_fact": 1.0,
"sens_deviation": [
0.0,
0.0,
0.0,
0.0,
0.0,
0.0
],
"sens_norm_fact": [
1.0,
1.0,
1.0,
1.0,
1.0,
1.0
],
"action_bias": [
0.0,
-5.1,
5.1
]
}

View File

@ -0,0 +1,11 @@
{
"scene": "steady",
"controlled": true,
"source": "open_loop",
"action_constant": [
0.0,
-0.051,
0.051
],
"note": "open-loop constant rotation, no PPO model"
}

View File

@ -0,0 +1,118 @@
{
"scene": "karman_re100",
"channels": [
{
"channel": "front",
"best_r2": 0.9950013415086292,
"best_nz": 21,
"pareto": [
{
"nz": 2,
"r2": 0.8681678005649112
},
{
"nz": 3,
"r2": 0.9692414821446799
},
{
"nz": 4,
"r2": 0.989378747581318
},
{
"nz": 6,
"r2": 0.9908017426802015
},
{
"nz": 12,
"r2": 0.9939890287801123
},
{
"nz": 13,
"r2": 0.9947660566975605
},
{
"nz": 18,
"r2": 0.9949963098276752
},
{
"nz": 21,
"r2": 0.9950013415086292
}
]
},
{
"channel": "top",
"best_r2": 0.9928439394510484,
"best_nz": 22,
"pareto": [
{
"nz": 1,
"r2": 0.44188145385324007
},
{
"nz": 2,
"r2": 0.9285296099294669
},
{
"nz": 6,
"r2": 0.9812946866555053
},
{
"nz": 12,
"r2": 0.9909377708950154
},
{
"nz": 17,
"r2": 0.9927843560231949
},
{
"nz": 20,
"r2": 0.992824536423302
},
{
"nz": 22,
"r2": 0.9928439394510484
}
]
},
{
"channel": "bottom",
"best_r2": 0.9965335484991993,
"best_nz": 22,
"pareto": [
{
"nz": 1,
"r2": 0.8456588970543057
},
{
"nz": 2,
"r2": 0.9030386794798415
},
{
"nz": 8,
"r2": 0.9947803148283133
},
{
"nz": 10,
"r2": 0.9962551509064745
},
{
"nz": 13,
"r2": 0.9963877299663628
},
{
"nz": 16,
"r2": 0.9964335809851655
},
{
"nz": 18,
"r2": 0.9965311727162228
},
{
"nz": 22,
"r2": 0.9965335484991993
}
]
}
]
}

View File

@ -0,0 +1,110 @@
{
"scene": "vortex_lamb",
"channels": [
{
"channel": "front",
"best_r2": 0.9035051679446613,
"best_nz": 21,
"pareto": [
{
"nz": 3,
"r2": 0.8536388609680176
},
{
"nz": 8,
"r2": 0.8957625139671737
},
{
"nz": 9,
"r2": 0.8985847902462778
},
{
"nz": 16,
"r2": 0.9017891237504362
},
{
"nz": 20,
"r2": 0.9035028044287374
},
{
"nz": 21,
"r2": 0.9035051679446613
}
]
},
{
"channel": "top",
"best_r2": 0.979810012198325,
"best_nz": 22,
"pareto": [
{
"nz": 0,
"r2": -0.2926478447353347
},
{
"nz": 1,
"r2": 0.9262664550660489
},
{
"nz": 2,
"r2": 0.9602153300731789
},
{
"nz": 5,
"r2": 0.9685659511773673
},
{
"nz": 6,
"r2": 0.9752589669369754
},
{
"nz": 19,
"r2": 0.9796029724451923
},
{
"nz": 20,
"r2": 0.9797408009698567
},
{
"nz": 22,
"r2": 0.979810012198325
}
]
},
{
"channel": "bottom",
"best_r2": 0.9334422885170565,
"best_nz": 22,
"pareto": [
{
"nz": 0,
"r2": -5.4349952892154265
},
{
"nz": 1,
"r2": 0.6935730348615337
},
{
"nz": 5,
"r2": 0.8721441125489685
},
{
"nz": 7,
"r2": 0.8963211646824344
},
{
"nz": 11,
"r2": 0.9294742132351927
},
{
"nz": 15,
"r2": 0.9318911728011019
},
{
"nz": 22,
"r2": 0.9334422885170565
}
]
}
]
}

View File

@ -0,0 +1,58 @@
{
"scene": "vortex_taylor",
"channels": [
{
"channel": "front",
"best_r2": 0.9603622630700551,
"best_nz": 21,
"pareto": [
{
"nz": 0,
"r2": -1762.7026716770156
},
{
"nz": 21,
"r2": 0.9603622630700551
}
]
},
{
"channel": "top",
"best_r2": 0.809824114603052,
"best_nz": 22,
"pareto": [
{
"nz": 0,
"r2": -3813.3981416722418
},
{
"nz": 1,
"r2": 4.909409545561516e-09
},
{
"nz": 22,
"r2": 0.809824114603052
}
]
},
{
"channel": "bottom",
"best_r2": 0.6431303693566448,
"best_nz": 22,
"pareto": [
{
"nz": 0,
"r2": -11389.175969107222
},
{
"nz": 1,
"r2": 2.5346330034814457e-08
},
{
"nz": 22,
"r2": 0.6431303693566448
}
]
}
]
}

Binary file not shown.

After

Width:  |  Height:  |  Size: 144 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 143 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 139 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 142 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 145 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 145 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 142 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 141 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 141 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 140 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 128 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 138 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 147 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 134 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 134 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 138 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 157 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 148 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 144 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 146 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 141 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 147 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 146 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 152 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 148 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 151 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 148 KiB

Binary file not shown.

After

Width:  |  Height:  |  Size: 148 KiB

View File

@ -0,0 +1,210 @@
{
"karman_re100_front": {
"output_label": "karman_re100_front",
"n_samples": 198,
"n_features": 1,
"feature_names": [
"aF_lag1"
],
"niterations": 50,
"best_equation": "(x0 * 0.01106541) / 1.106541",
"best_loss": 1.4959569e-19,
"r2": 1.0,
"equations": [
{
"complexity": 1,
"loss": 0.00018302062,
"score": 0.0,
"equation": "-0.003915151",
"sympy_format": "-0.00391515100000000"
},
{
"complexity": 3,
"loss": 1.7232441999999998e-19,
"score": 17.299498177983697,
"equation": "x0 * 0.01",
"sympy_format": "x0*0.01"
},
{
"complexity": 5,
"loss": 1.4959569e-19,
"score": 0.07072130403388514,
"equation": "(x0 * 0.01106541) / 1.106541",
"sympy_format": "x0*0.01106541/1.106541"
}
]
},
"karman_re100_top": {
"output_label": "karman_re100_top",
"n_samples": 198,
"n_features": 7,
"feature_names": [
"bias",
"u_a",
"Cd_rear",
"Cl_diff",
"aF_lag1",
"aB_lag1",
"aT_lag1"
],
"niterations": 50,
"best_equation": "(x6 * -0.00835853) / -0.835853",
"best_loss": 2.8386383999999997e-18,
"r2": 1.0,
"equations": [
{
"complexity": 1,
"loss": 0.00031669735,
"score": 0.0,
"equation": "0.041159406",
"sympy_format": "0.0411594060000000"
},
{
"complexity": 3,
"loss": 3.732284e-18,
"score": 16.035973661600742,
"equation": "x6 * 0.01",
"sympy_format": "x6*0.01"
},
{
"complexity": 5,
"loss": 2.8386383999999997e-18,
"score": 0.13684793905473897,
"equation": "(x6 * -0.00835853) / -0.835853",
"sympy_format": "x6*(-0.00835853)/(-0.835853)"
}
]
},
"illusion_1L_front": {
"output_label": "illusion_1L_front",
"n_samples": 198,
"n_features": 14,
"feature_names": [
"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"
],
"niterations": 100,
"best_equation": "((0.0039467034 * (x8 + x8)) - ((x8 / (-0.9984451 - (x4 + x8))) * square(2.5222103e-6))) / 0.7893407",
"best_loss": 8.728418e-19,
"r2": 1.0,
"equations": [
{
"complexity": 1,
"loss": 0.004748166,
"score": 0.0,
"equation": "0.0012775909",
"sympy_format": "0.00127759090000000"
},
{
"complexity": 3,
"loss": 2.6305633999999998e-18,
"score": 17.564668394678097,
"equation": "x8 * 0.01",
"sympy_format": "x8*0.01"
},
{
"complexity": 5,
"loss": 1.3032371999999998e-18,
"score": 0.35117336038102565,
"equation": "(0.007208903 * x8) / 0.7208903",
"sympy_format": "0.007208903*x8/0.7208903"
},
{
"complexity": 7,
"loss": 1.1532010999999999e-18,
"score": 0.06115484118442259,
"equation": "((x8 + x8) * 0.0039467034) / 0.7893407",
"sympy_format": "(x8 + x8)*0.0039467034/0.7893407"
},
{
"complexity": 9,
"loss": 1.1488888999999999e-18,
"score": 0.001873169616955362,
"equation": "(((x8 + x8) - 1.3038937e-7) * 0.0039467034) / 0.7893407",
"sympy_format": "(x8 + x8 - 1*1.3038937e-7)*0.0039467034/0.7893407"
},
{
"complexity": 12,
"loss": 1.1395758999999998e-18,
"score": 0.0027130419227974806,
"equation": "(((x8 + x8) * 0.0039467034) - square(x7 * -3.0155393e-6)) / 0.7893407",
"sympy_format": "(-9.09347726984449e-12*x7**2 + (x8 + x8)*0.0039467034)/0.7893407"
},
{
"complexity": 14,
"loss": 1.0790729e-18,
"score": 0.02727696454758814,
"equation": "((x8 + (x8 - square((x13 + x13) * -1.5024679e-5))) * 0.0039467034) / 0.7893407",
"sympy_format": "(x8 + x8 - 2.25740979053041e-10*(x13 + x13)**2)*0.0039467034/0.7893407"
},
{
"complexity": 18,
"loss": 8.728418e-19,
"score": 0.05302580007984551,
"equation": "((0.0039467034 * (x8 + x8)) - ((x8 / (-0.9984451 - (x4 + x8))) * square(2.5222103e-6))) / 0.7893407",
"sympy_format": "(-6.36154479742609e-12*x8/(-(x4 + x8) - 0.9984451) + 0.0039467034*(x8 + x8))/0.7893407"
}
]
},
"illusion_1L_top": {
"output_label": "illusion_1L_top",
"n_samples": 198,
"n_features": 15,
"feature_names": [
"bias",
"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"
],
"niterations": 100,
"best_equation": "(x11 * 0.0121531105) / 1.215311",
"best_loss": 2.9856626999999997e-18,
"r2": 1.0,
"equations": [
{
"complexity": 1,
"loss": 0.0034831127,
"score": 0.0,
"equation": "0.01056136",
"sympy_format": "0.0105613600000000"
},
{
"complexity": 3,
"loss": 3.6416598e-18,
"score": 17.247131588103095,
"equation": "x11 * 0.01",
"sympy_format": "x11*0.01"
},
{
"complexity": 5,
"loss": 2.9856626999999997e-18,
"score": 0.09930891723700508,
"equation": "(x11 * 0.0121531105) / 1.215311",
"sympy_format": "x11*0.0121531105/1.215311"
}
]
}
}

File diff suppressed because it is too large Load Diff

View File

@ -0,0 +1,618 @@
{
"thresholds": [
0.0,
0.001,
0.002,
0.005,
0.01,
0.015,
0.02,
0.03,
0.05,
0.1
],
"feature_names_front": [
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff"
],
"feature_names_rear": [
"bias",
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff"
],
"per_scene": {
"karman_re100": {
"scene": "karman_re100",
"re_code": 100,
"mu": 0.02,
"feature_names_front": [
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff"
],
"feature_names_rear": [
"bias",
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff"
],
"front": {
"results": [
{
"threshold": 0.0,
"nz": 21,
"r2": 0.9963468671522658,
"mae": 0.008495019589563084,
"weights_min": 0.2538872797382758,
"weights_max": 1.0371907088754313
},
{
"threshold": 0.001,
"nz": 18,
"r2": 0.9963454010383349,
"mae": 0.008500938429721638,
"weights_min": 0.2538872797382758,
"weights_max": 1.0371907088754313
},
{
"threshold": 0.002,
"nz": 18,
"r2": 0.9963454010383349,
"mae": 0.008500938429721638,
"weights_min": 0.2538872797382758,
"weights_max": 1.0371907088754313
},
{
"threshold": 0.005,
"nz": 12,
"r2": 0.9947081910813219,
"mae": 0.010856555791902431,
"weights_min": 0.2538872797382758,
"weights_max": 1.0371907088754313
},
{
"threshold": 0.01,
"nz": 9,
"r2": 0.9944015269767014,
"mae": 0.010728158860757616,
"weights_min": 0.2538872797382758,
"weights_max": 1.0371907088754313
},
{
"threshold": 0.015,
"nz": 4,
"r2": 0.9805810563712329,
"mae": 0.02109193834223528,
"weights_min": 0.2538872797382758,
"weights_max": 1.0371907088754313
},
{
"threshold": 0.02,
"nz": 3,
"r2": 0.969677527911751,
"mae": 0.024398378538262957,
"weights_min": 0.2538872797382758,
"weights_max": 1.0371907088754313
},
{
"threshold": 0.03,
"nz": 3,
"r2": 0.969677527911751,
"mae": 0.024398378538262957,
"weights_min": 0.2538872797382758,
"weights_max": 1.0371907088754313
},
{
"threshold": 0.05,
"nz": 3,
"r2": 0.969677527911751,
"mae": 0.024398378538262957,
"weights_min": 0.2538872797382758,
"weights_max": 1.0371907088754313
},
{
"threshold": 0.1,
"nz": 2,
"r2": 0.8671246137520265,
"mae": 0.0536574860231316,
"weights_min": 0.2538872797382758,
"weights_max": 1.0371907088754313
}
],
"best": {
"threshold": 0.002,
"nz": 18,
"r2": 0.9963454010383349,
"mae": 0.008500938429721638,
"weights_min": 0.2538872797382758,
"weights_max": 1.0371907088754313
},
"best_coef": [
0.003380802864171953,
-0.0008809022653997198,
0.0,
0.00208717851683782,
0.009307953320832238,
0.012225442968436223,
-0.007383343207562919,
0.019594876551595624,
0.0,
0.0,
0.0071310792485739126,
0.0014536838339802387,
-0.004129923355485657,
0.0014498851666906254,
-0.0027671395780320186,
0.0021240428248806,
0.9081582144711916,
-0.044045113030994884,
0.10435892480150129,
0.4653976700088291,
0.979743839707262
],
"sparsity_curve": [
[
0.0,
21,
0.9963468671522658
],
[
0.001,
18,
0.9963454010383349
],
[
0.002,
18,
0.9963454010383349
],
[
0.005,
12,
0.9947081910813219
],
[
0.01,
9,
0.9944015269767014
],
[
0.015,
4,
0.9805810563712329
],
[
0.02,
3,
0.969677527911751
],
[
0.03,
3,
0.969677527911751
],
[
0.05,
3,
0.969677527911751
],
[
0.1,
2,
0.8671246137520265
]
]
},
"top": {
"results": [
{
"threshold": 0.0,
"nz": 22,
"r2": 0.9947943202493076,
"mae": 0.013882836661429927,
"weights_min": 0.36423882973791505,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.001,
"nz": 20,
"r2": 0.9947724287185891,
"mae": 0.013998132434573892,
"weights_min": 0.36423882973791505,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.002,
"nz": 18,
"r2": 0.9947437375025526,
"mae": 0.013973729698665602,
"weights_min": 0.36423882973791505,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.005,
"nz": 18,
"r2": 0.9947437375025526,
"mae": 0.013973729698665602,
"weights_min": 0.36423882973791505,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.01,
"nz": 12,
"r2": 0.9930561323691804,
"mae": 0.015350435484475112,
"weights_min": 0.36423882973791505,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.015,
"nz": 12,
"r2": 0.9930561323691804,
"mae": 0.015350435484475112,
"weights_min": 0.36423882973791505,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.02,
"nz": 12,
"r2": 0.9930561323691804,
"mae": 0.015350435484475112,
"weights_min": 0.36423882973791505,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.03,
"nz": 6,
"r2": 0.9829012294927744,
"mae": 0.024013187353778373,
"weights_min": 0.36423882973791505,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.05,
"nz": 2,
"r2": 0.9310887280966665,
"mae": 0.04519335614359589,
"weights_min": 0.36423882973791505,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.1,
"nz": 1,
"r2": 0.42576849886018275,
"mae": 0.12579383759812146,
"weights_min": 0.36423882973791505,
"weights_max": 1.1067365435175973
}
],
"best": {
"threshold": 0.002,
"nz": 18,
"r2": 0.9947437375025526,
"mae": 0.013973729698665602,
"weights_min": 0.36423882973791505,
"weights_max": 1.1067365435175973
},
"best_coef": [
0.9115729610071152,
-0.006805683305581021,
-0.00407714864334097,
0.0028603164071186993,
0.0,
-0.020649798768545318,
0.033591198632074694,
0.09195198469780816,
0.0232186348784368,
0.0,
0.0,
0.0035421651717378612,
0.009890204105631987,
0.0028849152878799955,
-0.000946773567394545,
-0.00391306181767762,
0.005116314559440586,
0.018231459167452172,
-0.20385743262256592,
0.0,
-1.0324899193623733,
1.1609316346427214
],
"sparsity_curve": [
[
0.0,
22,
0.9947943202493076
],
[
0.001,
20,
0.9947724287185891
],
[
0.002,
18,
0.9947437375025526
],
[
0.005,
18,
0.9947437375025526
],
[
0.01,
12,
0.9930561323691804
],
[
0.015,
12,
0.9930561323691804
],
[
0.02,
12,
0.9930561323691804
],
[
0.03,
6,
0.9829012294927744
],
[
0.05,
2,
0.9310887280966665
],
[
0.1,
1,
0.42576849886018275
]
]
},
"bottom": {
"results": [
{
"threshold": 0.0,
"nz": 22,
"r2": 0.9974702578481651,
"mae": 0.011121305785433843,
"weights_min": 0.2863478575619962,
"weights_max": 1.0399071561860243
},
{
"threshold": 0.001,
"nz": 18,
"r2": 0.9974635520604825,
"mae": 0.011155185750820898,
"weights_min": 0.2863478575619962,
"weights_max": 1.0399071561860243
},
{
"threshold": 0.002,
"nz": 16,
"r2": 0.9974006411842374,
"mae": 0.011251386324348665,
"weights_min": 0.2863478575619962,
"weights_max": 1.0399071561860243
},
{
"threshold": 0.005,
"nz": 16,
"r2": 0.9974006411842374,
"mae": 0.011251386324348665,
"weights_min": 0.2863478575619962,
"weights_max": 1.0399071561860243
},
{
"threshold": 0.01,
"nz": 10,
"r2": 0.9971310266820465,
"mae": 0.011717975996893105,
"weights_min": 0.2863478575619962,
"weights_max": 1.0399071561860243
},
{
"threshold": 0.015,
"nz": 8,
"r2": 0.9955140894896422,
"mae": 0.01509793962669375,
"weights_min": 0.2863478575619962,
"weights_max": 1.0399071561860243
},
{
"threshold": 0.02,
"nz": 8,
"r2": 0.9955140894896422,
"mae": 0.01509793962669375,
"weights_min": 0.2863478575619962,
"weights_max": 1.0399071561860243
},
{
"threshold": 0.03,
"nz": 5,
"r2": 0.9920045223882177,
"mae": 0.01968881107737399,
"weights_min": 0.2863478575619962,
"weights_max": 1.0399071561860243
},
{
"threshold": 0.05,
"nz": 2,
"r2": 0.9027764309517692,
"mae": 0.06593190744450401,
"weights_min": 0.2863478575619962,
"weights_max": 1.0399071561860243
},
{
"threshold": 0.1,
"nz": 1,
"r2": 0.8435989700357883,
"mae": 0.08688057754133625,
"weights_min": 0.2863478575619962,
"weights_max": 1.0399071561860243
}
],
"best": {
"threshold": 0.002,
"nz": 16,
"r2": 0.9974006411842374,
"mae": 0.011251386324348665,
"weights_min": 0.2863478575619962,
"weights_max": 1.0399071561860243
},
"best_coef": [
-0.3062392108066253,
0.005083824376570689,
-0.0019100164756928636,
-0.0032417665760502293,
0.0,
0.0,
-0.016413610065012466,
0.062146484641493534,
-0.005222530364610314,
0.0,
0.0,
0.003517044018032268,
0.00988827788555667,
0.0023305542256132016,
-0.0015468136369196132,
0.004824491806330017,
-0.0032570158106502616,
-0.00612478415924937,
-0.09550082351895658,
0.0,
0.0,
-0.26112655898206505
],
"sparsity_curve": [
[
0.0,
22,
0.9974702578481651
],
[
0.001,
18,
0.9974635520604825
],
[
0.002,
16,
0.9974006411842374
],
[
0.005,
16,
0.9974006411842374
],
[
0.01,
10,
0.9971310266820465
],
[
0.015,
8,
0.9955140894896422
],
[
0.02,
8,
0.9955140894896422
],
[
0.03,
5,
0.9920045223882177
],
[
0.05,
2,
0.9027764309517692
],
[
0.1,
1,
0.8435989700357883
]
]
}
}
}
}

View File

@ -0,0 +1,420 @@
{
"thresholds": [
0.0,
0.001,
0.002,
0.005,
0.01,
0.015,
0.02,
0.03,
0.05,
0.1
],
"feature_names_front": [
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT"
],
"feature_names_rear": [
"bias",
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT"
],
"per_scene": {
"karman_re100": {
"scene": "karman_re100",
"re_code": 100,
"mu": 0.02,
"feature_names_front": [
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT"
],
"feature_names_rear": [
"bias",
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT"
],
"front": {
"results": [
{
"threshold": 0.0,
"nz": 16,
"r2": 0.9964173297164212,
"mae": 0.008496979102955536,
"weights_min": 0.2778368062626606,
"weights_max": 1.037970618101812
},
{
"threshold": 0.001,
"nz": 13,
"r2": 0.9964170471154649,
"mae": 0.00849877674299794,
"weights_min": 0.2778368062626606,
"weights_max": 1.037970618101812
},
{
"threshold": 0.002,
"nz": 13,
"r2": 0.9964170471154649,
"mae": 0.00849877674299794,
"weights_min": 0.2778368062626606,
"weights_max": 1.037970618101812
},
{
"threshold": 0.005,
"nz": 3,
"r2": 0.9885903068210383,
"mae": 0.014540007299834038,
"weights_min": 0.2778368062626606,
"weights_max": 1.037970618101812
},
{
"threshold": 0.01,
"nz": 3,
"r2": 0.9801597431201753,
"mae": 0.02107344007147761,
"weights_min": 0.2778368062626606,
"weights_max": 1.037970618101812
},
{
"threshold": 0.015,
"nz": 3,
"r2": 0.9801597431201753,
"mae": 0.02107344007147761,
"weights_min": 0.2778368062626606,
"weights_max": 1.037970618101812
},
{
"threshold": 0.02,
"nz": 2,
"r2": 0.9680925119257522,
"mae": 0.024303837665639722,
"weights_min": 0.2778368062626606,
"weights_max": 1.037970618101812
},
{
"threshold": 0.03,
"nz": 2,
"r2": 0.9680925119257522,
"mae": 0.024303837665639722,
"weights_min": 0.2778368062626606,
"weights_max": 1.037970618101812
},
{
"threshold": 0.05,
"nz": 2,
"r2": 0.9680925119257522,
"mae": 0.024303837665639722,
"weights_min": 0.2778368062626606,
"weights_max": 1.037970618101812
},
{
"threshold": 0.1,
"nz": 1,
"r2": 0.8667444738394041,
"mae": 0.053551706949699907,
"weights_min": 0.2778368062626606,
"weights_max": 1.037970618101812
}
],
"best": {
"threshold": 0.002,
"nz": 13,
"r2": 0.9964170471154649,
"mae": 0.00849877674299794,
"weights_min": 0.2778368062626606,
"weights_max": 1.037970618101812
},
"best_coef": [
0.0035693198413540234,
-0.0018962945486085396,
0.0,
0.004246223224503542,
0.01867909228324871,
0.011809790151579671,
-0.006291860766813562,
0.03938466562589553,
0.0,
0.0,
0.0071131542686216345,
0.0015118413457375216,
-0.004145051172507872,
0.0014265611111992293,
-0.0028705905475650928,
0.0021914873908442335
],
"sparsity_curve": [
[
0.0,
16,
0.9964173297164212
],
[
0.001,
13,
0.9964170471154649
],
[
0.002,
13,
0.9964170471154649
],
[
0.005,
3,
0.9885903068210383
],
[
0.01,
3,
0.9801597431201753
],
[
0.015,
3,
0.9801597431201753
],
[
0.02,
2,
0.9680925119257522
],
[
0.03,
2,
0.9680925119257522
],
[
0.05,
2,
0.9680925119257522
],
[
0.1,
1,
0.8667444738394041
]
]
},
"top": {
"results": [
{
"threshold": 0.0,
"nz": 17,
"r2": 0.9947943200325069,
"mae": 0.01388285452459937,
"weights_min": 0.3642388297379152,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.001,
"nz": 15,
"r2": 0.9947724285032213,
"mae": 0.013998153842830187,
"weights_min": 0.3642388297379152,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.002,
"nz": 15,
"r2": 0.9947724285032213,
"mae": 0.013998153842830187,
"weights_min": 0.3642388297379152,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.005,
"nz": 14,
"r2": 0.9947437372987523,
"mae": 0.013973763318069576,
"weights_min": 0.3642388297379152,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.01,
"nz": 11,
"r2": 0.994527427094902,
"mae": 0.014375559061365345,
"weights_min": 0.3642388297379152,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.015,
"nz": 10,
"r2": 0.9930561323361152,
"mae": 0.015350427162982516,
"weights_min": 0.3642388297379152,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.02,
"nz": 10,
"r2": 0.9930561323361152,
"mae": 0.015350427162982516,
"weights_min": 0.3642388297379152,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.03,
"nz": 9,
"r2": 0.9923352206869693,
"mae": 0.016117400302495512,
"weights_min": 0.3642388297379152,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.05,
"nz": 5,
"r2": 0.9832113360018212,
"mae": 0.02344478050006772,
"weights_min": 0.3642388297379152,
"weights_max": 1.1067365435175973
},
{
"threshold": 0.1,
"nz": 2,
"r2": 0.5190549772066559,
"mae": 0.1224813760307,
"weights_min": 0.3642388297379152,
"weights_max": 1.1067365435175973
}
],
"best": {
"threshold": 0.002,
"nz": 15,
"r2": 0.9947724285032213,
"mae": 0.013998153842830187,
"weights_min": 0.3642388297379152,
"weights_max": 1.1067365435175973
},
"best_coef": [
1.0021483455333764,
-0.008212536312529608,
-0.008095982403317927,
0.0033883289805039874,
-0.0014027415625294862,
-0.043634070706600914,
0.034508992736236595,
0.09250796489630796,
0.046370682051450104,
0.0,
0.0,
0.0035461180448679397,
0.009932878542230005,
0.002841799103744343,
-0.0009874050179528864,
-0.0039603133924009494,
0.0051583627341490485
],
"sparsity_curve": [
[
0.0,
17,
0.9947943200325069
],
[
0.001,
15,
0.9947724285032213
],
[
0.002,
15,
0.9947724285032213
],
[
0.005,
14,
0.9947437372987523
],
[
0.01,
11,
0.994527427094902
],
[
0.015,
10,
0.9930561323361152
],
[
0.02,
10,
0.9930561323361152
],
[
0.03,
9,
0.9923352206869693
],
[
0.05,
5,
0.9832113360018212
],
[
0.1,
2,
0.5190549772066559
]
]
}
}
}
}

View File

@ -0,0 +1,996 @@
{
"thresholds": [
0.0,
0.001,
0.002,
0.005,
0.01,
0.015,
0.02,
0.03,
0.05,
0.1
],
"per_scene": {
"vortex_lamb": {
"scene": "vortex_lamb",
"re_code": 100,
"mu": 0.02,
"n_samples": 148,
"feature_names_front": [
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff"
],
"feature_names_rear": [
"bias",
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff"
],
"front": {
"results": [
{
"threshold": 0.0,
"nz": 21,
"r2": 0.9035051679446613,
"mae": 0.05502827225008196
},
{
"threshold": 0.001,
"nz": 20,
"r2": 0.9035028044287374,
"mae": 0.054966606007031876
},
{
"threshold": 0.002,
"nz": 20,
"r2": 0.9035028044287374,
"mae": 0.054966606007031876
},
{
"threshold": 0.005,
"nz": 16,
"r2": 0.9017891237504362,
"mae": 0.053710513734883655
},
{
"threshold": 0.01,
"nz": 9,
"r2": 0.8985847902462778,
"mae": 0.054491226319598914
},
{
"threshold": 0.015,
"nz": 8,
"r2": 0.8957625139671737,
"mae": 0.05306021520916354
},
{
"threshold": 0.02,
"nz": 8,
"r2": 0.8957625139671737,
"mae": 0.05306021520916354
},
{
"threshold": 0.03,
"nz": 8,
"r2": 0.8957625139671737,
"mae": 0.05306021520916354
},
{
"threshold": 0.05,
"nz": 3,
"r2": 0.8536388609680176,
"mae": 0.048637995761707624
},
{
"threshold": 0.1,
"nz": 3,
"r2": 0.8536388609680176,
"mae": 0.048637995761707624
}
],
"best": {
"threshold": 0.0,
"nz": 21,
"r2": 0.9035051679446613,
"mae": 0.05502827225008196
},
"best_coef": [
0.0072499249196507284,
-0.013406647625693383,
-0.00038101504054493135,
-0.01577043426582744,
-0.021357793291808133,
0.09538149016633536,
-0.17597430070013861,
0.01105460767691196,
0.013467531810937875,
-0.008447877114831191,
-0.009842767133851753,
0.019209968175028087,
-0.036801729600496,
0.0036358387699529284,
0.006079481319894673,
0.003342294121453103,
-0.2971343137386559,
-0.6703323813926846,
-0.7885217073102587,
-1.0678896759007668,
0.5527301613931733
],
"sparsity_curve": [
[
0.0,
21,
0.9035051679446613
],
[
0.001,
20,
0.9035028044287374
],
[
0.002,
20,
0.9035028044287374
],
[
0.005,
16,
0.9017891237504362
],
[
0.01,
9,
0.8985847902462778
],
[
0.015,
8,
0.8957625139671737
],
[
0.02,
8,
0.8957625139671737
],
[
0.03,
8,
0.8957625139671737
],
[
0.05,
3,
0.8536388609680176
],
[
0.1,
3,
0.8536388609680176
]
]
},
"top": {
"results": [
{
"threshold": 0.0,
"nz": 22,
"r2": 0.979810012198325,
"mae": 0.009405479490894075
},
{
"threshold": 0.001,
"nz": 20,
"r2": 0.9797408009698567,
"mae": 0.009464924709826631
},
{
"threshold": 0.002,
"nz": 19,
"r2": 0.9796029724451923,
"mae": 0.009503048001896178
},
{
"threshold": 0.005,
"nz": 6,
"r2": 0.9752589669369754,
"mae": 0.009083428201966606
},
{
"threshold": 0.01,
"nz": 5,
"r2": 0.9685659511773673,
"mae": 0.01027038394304665
},
{
"threshold": 0.015,
"nz": 2,
"r2": 0.9602153300731789,
"mae": 0.012152564325363832
},
{
"threshold": 0.02,
"nz": 2,
"r2": 0.9602153300731789,
"mae": 0.012152564325363832
},
{
"threshold": 0.03,
"nz": 2,
"r2": 0.9602153300731789,
"mae": 0.012152564325363832
},
{
"threshold": 0.05,
"nz": 1,
"r2": 0.9262664550660489,
"mae": 0.013270022046144844
},
{
"threshold": 0.1,
"nz": 0,
"r2": -0.2926478447353347,
"mae": 0.10377883511859723
}
],
"best": {
"threshold": 0.0,
"nz": 22,
"r2": 0.979810012198325,
"mae": 0.009405479490894075
},
"best_coef": [
1.7326565320927485,
-0.021510750075876994,
0.0007162636287450304,
0.003690011764588065,
-0.001255780396435901,
0.008069876064011276,
-0.027165643005235728,
0.019732103606457083,
-0.007023240411992188,
-0.00413143280625655,
0.0016423995129504244,
0.0008079534031552331,
-0.005050487340309005,
0.011582526033145867,
0.0010455718346845003,
0.0012590379370475096,
0.0019097817325384912,
0.0346531306245988,
0.03581318435812541,
-0.06278853147020641,
0.40349380438892807,
-0.3511614171629877
],
"sparsity_curve": [
[
0.0,
22,
0.979810012198325
],
[
0.001,
20,
0.9797408009698567
],
[
0.002,
19,
0.9796029724451923
],
[
0.005,
6,
0.9752589669369754
],
[
0.01,
5,
0.9685659511773673
],
[
0.015,
2,
0.9602153300731789
],
[
0.02,
2,
0.9602153300731789
],
[
0.03,
2,
0.9602153300731789
],
[
0.05,
1,
0.9262664550660489
],
[
0.1,
0,
-0.2926478447353347
]
]
},
"bottom": {
"results": [
{
"threshold": 0.0,
"nz": 22,
"r2": 0.9334422885170565,
"mae": 0.0062053698726041735
},
{
"threshold": 0.001,
"nz": 15,
"r2": 0.9318911728011019,
"mae": 0.00612106433470266
},
{
"threshold": 0.002,
"nz": 11,
"r2": 0.9294742132351927,
"mae": 0.006263040564511156
},
{
"threshold": 0.005,
"nz": 7,
"r2": 0.8963211646824344,
"mae": 0.006779203944819282
},
{
"threshold": 0.01,
"nz": 5,
"r2": 0.8721441125489685,
"mae": 0.007291293800356046
},
{
"threshold": 0.015,
"nz": 1,
"r2": 0.6935730348615337,
"mae": 0.009771975563999018
},
{
"threshold": 0.02,
"nz": 1,
"r2": 0.6935730348615337,
"mae": 0.009771975563999018
},
{
"threshold": 0.03,
"nz": 1,
"r2": 3.300804074513053e-12,
"mae": 0.029445380091027484
},
{
"threshold": 0.05,
"nz": 1,
"r2": 3.300804074513053e-12,
"mae": 0.029445380091027484
},
{
"threshold": 0.1,
"nz": 0,
"r2": -5.4349952892154265,
"mae": 0.0796928089615461
}
],
"best": {
"threshold": 0.0,
"nz": 22,
"r2": 0.9334422885170565,
"mae": 0.0062053698726041735
},
"best_coef": [
0.6365299396095444,
-0.00680014903595817,
-0.0014763845224131282,
0.0008331648138976364,
-0.0013963333722837737,
0.0026865900945637236,
-0.014425589645304323,
-0.008472370593784707,
-6.004332296000475e-05,
-0.0014926889561804521,
-0.0004258139447154254,
-0.0018644416734531306,
0.006364757842372448,
-0.003199874412701356,
0.001370932025212518,
0.002115447057480394,
0.0019992859663057654,
0.012730598783987558,
-0.07381922521108229,
-0.069816526655473,
0.13432950141437106,
-0.0030019765287261236
],
"sparsity_curve": [
[
0.0,
22,
0.9334422885170565
],
[
0.001,
15,
0.9318911728011019
],
[
0.002,
11,
0.9294742132351927
],
[
0.005,
7,
0.8963211646824344
],
[
0.01,
5,
0.8721441125489685
],
[
0.015,
1,
0.6935730348615337
],
[
0.02,
1,
0.6935730348615337
],
[
0.03,
1,
3.300804074513053e-12
],
[
0.05,
1,
3.300804074513053e-12
],
[
0.1,
0,
-5.4349952892154265
]
]
}
},
"vortex_taylor": {
"scene": "vortex_taylor",
"re_code": 100,
"mu": 0.02,
"n_samples": 148,
"feature_names_front": [
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff"
],
"feature_names_rear": [
"bias",
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff"
],
"front": {
"results": [
{
"threshold": 0.0,
"nz": 21,
"r2": 0.9603622630700551,
"mae": 5.5722956333111334e-05
},
{
"threshold": 0.001,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.002,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.005,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.01,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.015,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.02,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.03,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.05,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.1,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
}
],
"best": {
"threshold": 0.0,
"nz": 21,
"r2": 0.9603622630700551,
"mae": 5.5722956333111334e-05
},
"best_coef": [
0.00013185241334928363,
0.0013869486093058166,
-2.0689454061431898e-05,
0.0003650495654675631,
-7.646260510377518e-05,
0.0016834710794748663,
0.002138870811614627,
0.0009672545588015969,
-0.0011958663539923654,
-0.00015017143903531285,
0.007403929744140764,
-0.0015484862269750206,
-0.0003240821985992223,
-0.0010978795962499714,
-0.0002774795287867582,
-0.00041159658606479107,
0.00011190538695352297,
0.06934742736821846,
0.018252511507361284,
-0.0038231311262574906,
0.04836265570630274
],
"sparsity_curve": [
[
0.0,
21,
0.9603622630700551
],
[
0.001,
0,
-1762.7026716770156
],
[
0.002,
0,
-1762.7026716770156
],
[
0.005,
0,
-1762.7026716770156
],
[
0.01,
0,
-1762.7026716770156
],
[
0.015,
0,
-1762.7026716770156
],
[
0.02,
0,
-1762.7026716770156
],
[
0.03,
0,
-1762.7026716770156
],
[
0.05,
0,
-1762.7026716770156
],
[
0.1,
0,
-1762.7026716770156
]
]
},
"top": {
"results": [
{
"threshold": 0.0,
"nz": 22,
"r2": 0.809824114603052,
"mae": 0.00019280734797529785
},
{
"threshold": 0.001,
"nz": 1,
"r2": 4.909409545561516e-09,
"mae": 0.0004627442799109909
},
{
"threshold": 0.002,
"nz": 1,
"r2": 4.909409545561516e-09,
"mae": 0.0004627442799109909
},
{
"threshold": 0.005,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
},
{
"threshold": 0.01,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
},
{
"threshold": 0.015,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
},
{
"threshold": 0.02,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
},
{
"threshold": 0.03,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
},
{
"threshold": 0.05,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
},
{
"threshold": 0.1,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
}
],
"best": {
"threshold": 0.0,
"nz": 22,
"r2": 0.809824114603052,
"mae": 0.00019280734797529785
},
"best_coef": [
-0.003855648809544881,
-0.00016402917649747364,
-0.00930469925984968,
0.0003725286141146794,
-0.0018143741612594137,
-0.0006228672933056743,
-0.004123657830294623,
0.009617379353839065,
-0.0020693580870050163,
0.0064472830467429834,
-0.0009363904689749538,
0.002848747579355727,
0.008718629599353982,
0.006611618826853361,
0.0026442625498493562,
-0.000822158875902874,
-0.00019526846815097482,
-7.71129751901901e-05,
-0.4652349584088788,
-0.09071890000133181,
-0.031143453311218216,
-0.10346632929440941
],
"sparsity_curve": [
[
0.0,
22,
0.809824114603052
],
[
0.001,
1,
4.909409545561516e-09
],
[
0.002,
1,
4.909409545561516e-09
],
[
0.005,
0,
-3813.3981416722418
],
[
0.01,
0,
-3813.3981416722418
],
[
0.015,
0,
-3813.3981416722418
],
[
0.02,
0,
-3813.3981416722418
],
[
0.03,
0,
-3813.3981416722418
],
[
0.05,
0,
-3813.3981416722418
],
[
0.1,
0,
-3813.3981416722418
]
]
},
"bottom": {
"results": [
{
"threshold": 0.0,
"nz": 22,
"r2": 0.6431303693566448,
"mae": 0.00012115305583366485
},
{
"threshold": 0.001,
"nz": 1,
"r2": 2.5346330034814457e-08,
"mae": 0.00029202565622359403
},
{
"threshold": 0.002,
"nz": 1,
"r2": 2.5346330034814457e-08,
"mae": 0.00029202565622359403
},
{
"threshold": 0.005,
"nz": 1,
"r2": 2.5346330034814457e-08,
"mae": 0.00029202565622359403
},
{
"threshold": 0.01,
"nz": 1,
"r2": 2.5346330034814457e-08,
"mae": 0.00029202565622359403
},
{
"threshold": 0.015,
"nz": 1,
"r2": 2.5346330034814457e-08,
"mae": 0.00029202565622359403
},
{
"threshold": 0.02,
"nz": 1,
"r2": 2.5346330034814457e-08,
"mae": 0.00029202565622359403
},
{
"threshold": 0.03,
"nz": 0,
"r2": -11389.175969107222,
"mae": 0.05019249195686063
},
{
"threshold": 0.05,
"nz": 0,
"r2": -11389.175969107222,
"mae": 0.05019249195686063
},
{
"threshold": 0.1,
"nz": 0,
"r2": -11389.175969107222,
"mae": 0.05019249195686063
}
],
"best": {
"threshold": 0.0,
"nz": 22,
"r2": 0.6431303693566448,
"mae": 0.00012115305583366485
},
"best_coef": [
0.024327554221764875,
-0.0005249278989068122,
-0.0012993577286068104,
2.2680966036791235e-07,
-0.00010725143675978186,
-6.83864821520398e-05,
0.0010096871646746155,
0.010716542986185212,
0.0009908500661735154,
0.0011374615793514392,
0.0001758268322788275,
-0.0008997304730805852,
-0.0002053415257431248,
0.0003868365622255378,
0.005272105474899691,
0.0033010584355117824,
0.0027208004849191268,
0.0004865510823514911,
-0.06496786884052097,
-0.005364174292120229,
-0.003419346453559459,
0.04954205911186436
],
"sparsity_curve": [
[
0.0,
22,
0.6431303693566448
],
[
0.001,
1,
2.5346330034814457e-08
],
[
0.002,
1,
2.5346330034814457e-08
],
[
0.005,
1,
2.5346330034814457e-08
],
[
0.01,
1,
2.5346330034814457e-08
],
[
0.015,
1,
2.5346330034814457e-08
],
[
0.02,
1,
2.5346330034814457e-08
],
[
0.03,
0,
-11389.175969107222
],
[
0.05,
0,
-11389.175969107222
],
[
0.1,
0,
-11389.175969107222
]
]
}
}
}
}

View File

@ -0,0 +1,8 @@
{
"scene": "illusion_0.5L",
"mode": "pysr",
"similarity_full": 0.8539735529506018,
"similarity_tail": 0.8411474326454078,
"action_range": 0.06714382640888465,
"n_steps": 320
}

View File

@ -0,0 +1,8 @@
{
"scene": "illusion_0.6L",
"mode": "pysr",
"similarity_full": 0.9388972369369757,
"similarity_tail": 0.8608654259522963,
"action_range": 0.044579733956375364,
"n_steps": 320
}

View File

@ -0,0 +1,8 @@
{
"scene": "illusion_0.75L",
"mode": "pysr",
"similarity_full": 0.9821949471763758,
"similarity_tail": 0.8071620435178003,
"action_range": 0.029841950903943048,
"n_steps": 320
}

View File

@ -0,0 +1,9 @@
{
"scene": "illusion_0.75L",
"mode": "v23",
"similarity_full": 0.9741372518705987,
"similarity_tail": 0.9726050667188786,
"action_range": 0.03237645857182351,
"n_steps": 320,
"threshold": "best"
}

View File

@ -0,0 +1,8 @@
{
"scene": "illusion_0.8L",
"mode": "pysr",
"similarity_full": 0.9076366235277858,
"similarity_tail": 0.9003715798241534,
"action_range": 0.02705657196270742,
"n_steps": 214
}

View File

@ -0,0 +1,8 @@
{
"scene": "illusion_1.2L",
"mode": "pysr",
"similarity_full": 0.8485137980740682,
"similarity_tail": 0.7692440834631539,
"action_range": 0.026140261986299245,
"n_steps": 214
}

View File

@ -0,0 +1,9 @@
{
"scene": "illusion_1.5L",
"mode": "v23",
"similarity_full": 0.9296044171184163,
"similarity_tail": 0.925839316423258,
"action_range": 0.035437546689072105,
"n_steps": 200,
"threshold": "best"
}

View File

@ -0,0 +1,8 @@
{
"scene": "illusion_1L",
"mode": "pysr",
"similarity_full": 0.9580609580562278,
"similarity_tail": 0.8857815231916633,
"action_range": 0.011634223590664046,
"n_steps": 214
}

View File

@ -0,0 +1,8 @@
{
"scene": "illusion_1L",
"mode": "pysr",
"similarity_full": 0.9580609580562278,
"similarity_tail": 0.8857815231916633,
"action_range": 0.011634223590664046,
"n_steps": 214
}

View File

@ -0,0 +1,9 @@
{
"scene": "illusion_1L",
"mode": "v23",
"similarity_full": 0.9581664533499054,
"similarity_tail": 0.9575937863460072,
"action_range": 0.014234293670783398,
"n_steps": 320,
"threshold": "best"
}

View File

@ -0,0 +1,8 @@
{
"scene": "illusion_2L",
"mode": "pysr",
"similarity_full": 0.6754607445898861,
"similarity_tail": 0.6292225479264744,
"action_range": 0.3686004659586289,
"n_steps": 160
}

View File

@ -0,0 +1,27 @@
{
"scene": "karman_joint",
"scene_id": "karman",
"channel": "front",
"output": "alpha",
"feature_names": [
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"daF_dt",
"daB_dt",
"daT_dt",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff",
"mu_Cl_tot"
],
"best_sympy": "daB_dt*0.38504666 + daF_dt + mu_Cl_tot*(-14.951645)",
"best_score": 1.0
}

View File

@ -0,0 +1,28 @@
{
"scene": "karman_joint",
"scene_id": "karman",
"channel": "top",
"output": "alpha",
"feature_names": [
"bias",
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"daF_dt",
"daB_dt",
"daT_dt",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff",
"mu_Cl_tot"
],
"best_sympy": "3.4138296",
"best_score": 1.0
}

Some files were not shown because too many files have changed in this diff Show More