Frank14f/DynamisLab
1
0
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases Wiki Activity

chore: project-wide cleanup — consolidate docs, remove obsolete code, update .gitignore

Browse Source 
- Remove obsolete docs (OID_handover, SR_analysis_results, ccd_* handover)
- Remove CCD legacy output_redux and old scripts
- Remove SR old sindy scripts and compare modules
- Update .gitignore to cover all analysis-generated outputs
- Retain all active code in OID/SR/CCD analysis directories

Co-authored-by: Cursor <cursoragent@cursor.com>
This commit is contained in:
Frank14f
2026-06-28 16:52:48 +08:00
co-authored by Cursor
parent c918ac0de4
commit 56e3c78a83
73 changed files with 1365 additions and 11072 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/SR_analysis/README.md
+145 -146
View File
@@ -1,20 +1,78 @@
# SR_analysis: Unified SINDy-SR Analysis Pipeline
# SR_analysis: Symbolic Regression Analysis Pipeline
## Overview
Extracts interpretable control laws (`obs -> act`) from DRL-trained policies for the
fluidic pinball. Uses **PySR symbolic regression** on dimensionless physical features with
G-equivariant structural constraints (v23: front no-bias, rear shared-head).
This directory consolidates the SINDy-and-symbolic-regression analysis pipeline
for the DynamisLab fluidic pinball project. It replaces the old
`src/analysis_crossre/` and `src/analysis_cloak/` directories with a unified
structure.
## Current Results (2026-06-25)
The pipeline fits **sparse interpretable control laws** (`obs -> act`) for all
cloak and illusion scenes, using dimensionless physical features,
G-equivariant structural constraints, and STLSQ threshold grids.
### Karman Cloak — Cross-Re Unified Formula
For background, see:
- `sindy_sr_notes.md` -- execution plan and task tracking
- `sindy_sr_knowledge.md` -- confirmed facts and known pitfalls
- `../../docs/SR_analysis_results.md` -- comprehensive results report
| 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** |
### Illusion
| 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) | — | — |
**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.
### Illusion Generalization (Joint Formula, No PPO)
| 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 |
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 |
---
## 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
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
---
@@ -22,164 +80,105 @@ For background, see:
```
SR_analysis/
configs.py # Unified scene metadata (all 10+ scenes)
configs/
legacy/ # Legacy CFD configs
configs.py # Scene metadata (Karman, Illusion, Vortex)
configs/legacy/ # Legacy CFD configs (config_cuda.json, config_flowfield.json)
utils/
__init__.py # Selective exports (no pycuda dependency)
feature_builder.py # Dimensionless features + G-operator + phase-state features
sindy_fitter.py # STLSQ + feature matrices + derivative/absolute modes
__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/
karman/ # Karman cloak: karman_re50/100/200/400
steady/ # Steady cloak
illusion/ # Illusion: illusion_0.75L/1L/1.5L
vortex/ # Vortex cloak
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 # Inference: LegacyCFD + PPO -> controlled.npz
infer_illusion.py # Inference for illusion scenes
infer_vortex.py # Inference for vortex scenes
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_all_v2.py # Unified SINDy fitting (supports --deriv, --phase, --output-mode etc.)
run_pysr.py # Restricted PySR symbolic regression
wrap_joint.py # Joint model -> wrapped format for validator
compare_v2.py # Cross-scene comparison report
karman/illusion/vortex/ # SINDy output JSONs
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 (v23/deriv/abs modes)
run_closed_loop_illusion.py # Illusion closed-loop validator
eval_rollout.py # Offline multi-step rollout evaluation
results/ # Validation result JSONs
compare/
support_overlap.py # Support set comparison
shared_core.py # Shared core detection
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
```
---
## Key Design Decisions
## Usage
### 1. Scene Metadata Driven
All scene parameters defined once in `configs.py`.
### 2. Feature Levels
| Level | Features | Dim | Description |
|-------|----------|:---:|-------------|
| Static | u_m, u_a, u_c, v_a, Cd_tot, Cd_rear, Cl_tot, Cl_diff | 8 | Current-step only |
| **Phase-state** | u_a, du_a/dt, Cl_tot, dCl_tot/dt, Cd_tot, Cd_rear | **6** | Oscillation phase + rate |
| Illusion-phase | Phase-state + Cd_err, Cl_err, dCd_err/dt, dCl_err/dt | **10** | Phase + error-state |
| Karman-expanded | Phase-state + u_m, u_c, v_a, Cl_diff | **10** | Phase + supplementary |
| Full-lag | Static + lag-1 | 16 | Full temporal context |
### 3. Output Modes
- **deriv**: predict `d(alpha)/dt`, then `alpha(t) = alpha(t-1) + dt_c * dalpha/dt`
- **absolute**: predict `alpha(t)` directly (no integration drift)
### 4. G-Equivariant Structure (v23)
```
Front(t) = f_front(x(t)) # no bias, odd under G
Top(t) = f_rear(x(t)) # with bias
Bottom(t) = -f_rear(G[x(t)]) # shared-head
```
---
## Current Best Results (2026-06-15)
### Illusion — New Route: Phase-state + Error-state + Absolute Action
| Scene | Closed-loop | % of PPO | Action history? | Features |
|-------|:----------:|:--------:|:---------------:|----------|
| 0.75L | **0.974** | 100.2% | **No** | ILLUSION_PHASE (10dim) |
| 1L | **0.958** | 98.5% | **No** | ILLUSION_PHASE (10dim) |
| 1.5L | N/A | — | **No** | Bang-bang regime |
### Karman re100 — Ablation
| Config | Feat | Output | R2 | Closed-loop | Note |
|--------|:----:|:-----:|:--:|:----------:|------|
| old v23 (a_lag) | 14+3 | alpha | 0.996 | **0.901** | Baseline |
| **Phase->abs** | **6** | **alpha** | **0.965** | **0.699** | Best new route |
| Phase->deriv | 6 | dalpha/dt | 0.837 | 0.656 | |
| Phase+mu->abs | 9 | alpha | 0.979 | 0.700 | mu helps cross-Re |
| Expanded->abs | 10 | alpha | 0.980 | 0.580 | Overfitting |
---
## Commands
All from repo root (`/home/frank14f/DynamisLab`).
### SINDy Fitting
### PySR Symbolic Regression (conda: sr_env)
```bash
# Illusion phase-state + absolute (recommended for 0.75L/1L)
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes illusion_0.75L,illusion_1L --deriv --phase --output-mode absolute
# Illusion
conda run -n sr_env python src/SR_analysis/sindy/run_pysr_deep_illusion.py --individual
# Karman phase-state + absolute
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re100 --deriv --phase --output-mode absolute
# Karman expanded (10 dim)
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re100 --deriv --karman-expand --output-mode absolute
# Karman with mu modulation
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re100 --deriv --karman-mu --output-mode absolute
# Old-style (v2, with action history)
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re50,karman_re100 --joint
# Karman deep (cross-Re independent + joint)
conda run -n sr_env python src/SR_analysis/sindy/run_pysr_deep.py --both
```
### Closed-loop Validation
### CFD Closed-Loop Validation (conda: pycuda_3_10, GPU 1 or 2)
```bash
# Karman with absolute action
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \
--scene karman_re100 --device 0 --steps 200 --mode abs \
--sindy-results src/SR_analysis/sindy/karman/sindy_results_deriv.json
# Karman old v23
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \
--scene karman_re100 --device 0 --steps 200 --mode v23 \
--sindy-results src/SR_analysis/sindy/karman/sindy_joint_wrapped.json
# Illusion with absolute action
# Illusion PySR formula
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop_illusion.py \
--scene illusion_1L --device 0 --steps 320 \
--sindy-results src/SR_analysis/sindy/illusion/sindy_results_deriv.json
--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
```
### PySR Symbolic Regression
### PPO Inference (generate controlled.npz)
```bash
conda run -n sr_env python src/SR_analysis/sindy/run_pysr.py --scene illusion_1L
```
### Offline Rollout Evaluation
```bash
python3 src/SR_analysis/validate/eval_rollout.py \
--sindy-results src/SR_analysis/sindy/karman/sindy_results_deriv.json \
--scene karman_re100
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
```
---
## Important Reminders
## Critical Reminders
- `controlled.npz` actions are **normalized [-1,1]** — must convert via `(norm * scale + bias) * u0`
- **FIFO bias ≠ DRL action bias** for Illusion: FIFO=[0, -0.01, 0.01], decode=[0, -0.02, 0.02]
- "2U" in model name = S_DIM=14 (not 2x velocity), u0 always 0.01
- SAMPLE_INTERVAL: 0.75L=400, 1L=600, 1.5L/Karman=800
- Closed-loop steps auto-set: S=400→320, S=600→214, S=800→160
- One-step R² high ≠ closed-loop good — always validate
- For phase-state features, always pass `sensors_raw`/`forces_raw` to enable derivative computation
- **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 |
src/SR_analysis/compare/shared_core.py
-152
View File
@@ -1,152 +0,0 @@
"""Shared core detection across scenes.
Finds features that are active across ALL scenes in a group (e.g. all Karman Re,
all Illusion diameters) and identifies the cross-scene shared core.
Usage:
python compare/shared_core.py --sindy-results sindy/karman/sindy_results.json \\
--scenes karman_re50 karman_re100 karman_re200 karman_re400
python compare/shared_core.py \\
--sindy-results sindy/results.json \\
--scenes karman_re100 illusion_1L vortex_lamb steady
"""
from __future__ import annotations
import argparse
import json
import os
import sys
from typing import Dict, List, Tuple
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
if _REPO not in sys.path:
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from SR_analysis.utils.sindy_fitter import get_active_support
RELATIVE_THRESHOLD = 0.02
def feat_group(name: str) -> str:
if name == "bias":
return "bias"
if name in ("u_m", "u_a", "u_c", "v_a", "sin_ua", "cos_ua"):
return "sensor"
if name.startswith("Cd") or name.startswith("Cl"):
return "force"
if "lag1" in name:
return "memory_lag"
if name.startswith("da"):
return "memory_delta"
if name == "mu" or name.startswith("mu_"):
return "mu_mod"
return "other"
def detect_core(scene_data: Dict[str, dict], channels: List[str],
threshold: float) -> dict:
"""Find features active in ALL scenes for each channel."""
scene_names = list(scene_data.keys())
results = {}
for ch_name in channels:
fn_key = f"feature_names_{'front' if ch_name == 'front' else 'rear'}"
# Collect active sets per scene
active_per_scene = {}
for sn in scene_names:
fn = scene_data[sn][fn_key]
coef = scene_data[sn][ch_name]["best_coef"]
active = get_active_support(np.array(coef, dtype=np.float64)[:len(fn)],
fn, threshold)
active_per_scene[sn] = set(active.keys())
# Intersection = shared core
core_keys = set.intersection(*active_per_scene.values()) if active_per_scene else set()
# Union for scene-specific detection
all_keys = set.union(*active_per_scene.values()) if active_per_scene else set()
scene_specific = {}
for sn in scene_names:
others = set.union(*[v for k, v in active_per_scene.items() if k != sn])
diff = active_per_scene[sn] - others
if diff:
scene_specific[sn] = sorted(diff)
# Coef means for core features
core_coefs = {}
for k in sorted(core_keys):
vals = []
for sn in scene_names:
fn = scene_data[sn][fn_key]
coef = scene_data[sn][ch_name]["best_coef"]
if k in fn:
idx = fn.index(k)
vals.append(float(coef[idx]))
core_coefs[k] = {
"mean": float(np.mean(vals)),
"std": float(np.std(vals)),
"per_scene": {sn: vals[i] for i, sn in enumerate(scene_names)},
}
results[ch_name] = {
"n_scenes": len(scene_names),
"n_core": len(core_keys),
"core_features": {k: {"group": feat_group(k), "coef": v}
for k, v in core_coefs.items()},
"scene_specific": {sn: sorted(v) for sn, v in scene_specific.items()},
}
return results
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--sindy-results", type=str, required=True)
ap.add_argument("--scenes", type=str, nargs="+", required=True)
ap.add_argument("--threshold", type=float, default=RELATIVE_THRESHOLD)
ap.add_argument("--out", type=str, default=None)
args = ap.parse_args()
with open(args.sindy_results) as f:
all_data = json.load(f)
per = all_data.get("per_scene", {})
scene_data = {sn: per[sn] for sn in args.scenes if sn in per}
if len(scene_data) < 2:
print(f"Need >=2 scenes. Found: {list(scene_data.keys())}")
return 1
print(f"Shared Core Detection: {len(scene_data)} scenes")
for sn in scene_data:
print(f" {sn}")
print(f" threshold={args.threshold}")
results = detect_core(scene_data, ["front", "top", "bottom"], args.threshold)
for ch_name, ch_data in results.items():
print(f"\n--- {ch_name} ---")
print(f" Core features: {ch_data['n_core']}")
for k, v in ch_data["core_features"].items():
c = v["coef"]
print(f" {k:20s} mean={c['mean']:+.6f} std={c['std']:.6f} [{v['group']}]")
for sn, keys in ch_data["scene_specific"].items():
if keys:
print(f" {sn} specific: {', '.join(keys)}")
if args.out:
output = {"scenes": args.scenes, "threshold": args.threshold,
"channels": results}
os.makedirs(os.path.dirname(args.out), exist_ok=True)
with open(args.out, "w") as f:
json.dump(output, f, indent=2)
print(f"\nSaved: {args.out}")
if __name__ == "__main__":
main()
src/SR_analysis/compare/support_overlap.py
-158
View File
@@ -1,158 +0,0 @@
"""Cross-scene support overlap analysis.
Compares SINDy support sets across scenes at a given relative threshold.
Usage:
python compare/support_overlap.py --sindy-results sindy/karman/sindy_results.json \\
--scenes karman_re100 karman_re200 illusion_1L vortex_lamb
"""
from __future__ import annotations
import argparse
import json
import os
import sys
from typing import Dict, List, Tuple
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
if _REPO not in sys.path:
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from SR_analysis.utils.sindy_fitter import get_active_support
RELATIVE_THRESHOLD = 0.02 # default: 2% of max coefficient
def load_sindy_scenes(sindy_path: str, scenes: List[str]) -> dict:
"""Load sindy results for the given scene names."""
with open(sindy_path) as f:
data = json.load(f)
result = {}
for sn in scenes:
per = data["per_scene"].get(sn)
if per is None:
print(f"WARNING: {sn} not found in {sindy_path}")
continue
result[sn] = per
return result
def feat_group(name: str) -> str:
"""Classify a feature into a group."""
if name == "bias":
return "bias"
if name in ("u_m", "u_a", "u_c", "v_a", "sin_ua", "cos_ua"):
return "sensor"
if name.startswith("Cd") or name.startswith("Cl"):
return "force"
if "lag1" in name:
return "memory_lag"
if name.startswith("da"):
return "memory_delta"
if name == "mu" or name.startswith("mu_"):
return "mu_mod"
return "other"
def classify(
a_active: Dict[str, float],
b_active: Dict[str, float],
) -> Tuple[List[Tuple[str, float, float]], List[Tuple[str, float]], List[Tuple[str, float]]]:
"""Classify features as shared, A-only, B-only.
Returns (shared, A_only, B_only) where shared has feature name + both coeffs.
"""
a_keys = set(a_active.keys())
b_keys = set(b_active.keys())
shared = sorted(a_keys & b_keys)
a_only = sorted(a_keys - b_keys)
b_only = sorted(b_keys - a_keys)
shared_out = [(k, a_active[k], b_active[k]) for k in shared]
a_out = [(k, a_active[k]) for k in a_only]
b_out = [(k, b_active[k]) for k in b_only]
return shared_out, a_out, b_out
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--sindy-results", type=str, required=True)
ap.add_argument("--scenes", type=str, nargs="+", required=True,
help="Scene names to compare")
ap.add_argument("--threshold", type=float, default=RELATIVE_THRESHOLD)
ap.add_argument("--channels", type=str, nargs="+",
default=["front", "top", "bottom"],
help="Which channels to compare")
ap.add_argument("--out", type=str, default=None)
args = ap.parse_args()
data = load_sindy_scenes(args.sindy_results, args.scenes)
if len(data) < 2:
print("Need at least 2 scenes to compare")
return 1
scene_names = list(data.keys())
scene_a, scene_b = scene_names[0], scene_names[1]
print(f"Support Overlap: {scene_a} vs {scene_b} (th={args.threshold})")
print("=" * 60)
all_results = {}
for ch_name in args.channels:
# Map "front" -> "feature_names_front", etc
fn_key = f"feature_names_{'front' if ch_name == 'front' else 'rear'}"
fn_a = data[scene_a][fn_key]
fn_b = data[scene_b][fn_key]
fn_min = min(len(fn_a), len(fn_b))
fn_a_trim = fn_a[:fn_min]
fn_b_trim = fn_b[:fn_min]
ch_a = get_active_support(np.array(data[scene_a][ch_name]["best_coef"])[:fn_min],
fn_a_trim, args.threshold)
ch_b = get_active_support(np.array(data[scene_b][ch_name]["best_coef"])[:fn_min],
fn_b_trim, args.threshold)
shared, a_only, b_only = classify(ch_a, ch_b)
print(f"\n--- {ch_name} ---")
print(f" {scene_a} nz={len(ch_a)} {scene_b} nz={len(ch_b)} Shared={len(shared)}")
for name, ca, cb in shared:
print(f" {name:20s} A={ca:+9.6f} B={cb:+9.6f} [{feat_group(name)}]")
for name, ca in a_only:
print(f" {scene_a[:10]:>10s} {name:20s} A={ca:+9.6f} [{feat_group(name)}]")
for name, cb in b_only:
print(f" {scene_b[:10]:>10s} {name:20s} B={cb:+9.6f} [{feat_group(name)}]")
all_results[ch_name] = {
"scene_a_nz": len(ch_a),
"scene_b_nz": len(ch_b),
"shared_nz": len(shared),
"shared": [{"name": n, "coef_a": ca, "coef_b": cb} for n, ca, cb in shared],
f"{scene_a}_only": [{"name": n, "coef": ca} for n, ca in a_only],
f"{scene_b}_only": [{"name": n, "coef": cb} for n, cb in b_only],
}
if args.out:
output = {"scene_a": scene_a, "scene_b": scene_b,
"threshold": args.threshold,
"channels": all_results}
os.makedirs(os.path.dirname(args.out), exist_ok=True)
with open(args.out, "w") as f:
json.dump(output, f, indent=2)
print(f"\nSaved: {args.out}")
if __name__ == "__main__":
main()
src/SR_analysis/configs.py
+2
View File
@@ -213,6 +213,8 @@ for vtype, mn, strength in _SCENES_VORTEX:
# -- Illusion generalization (unseen diameters, no PPO model) ----------------
for diam, mn, si in [
(0.5, None, 400),
(0.6, None, 400), # added 2026-06-25 for generalization test
(0.8, None, 600), # added 2026-06-25 for generalization test
(1.2, None, 600),
(2.0, None, 800),
]:
src/SR_analysis/data/illusion/illusion_0.75L/config.json
+4 -3
View File
@@ -5,17 +5,18 @@
"mu": 0.02,
"nu": 0.004,
"has_disturbance": false,
"sample_interval": 600,
"sample_interval": 400,
"conv_len": 36,
"action_scale": 8.0,
"action_bias": "(0.0, -2.0, 2.0)",
"source": "PPO_inference",
"model_name": "d1a3o12_250525_imit_075L_1U",
"model_name": "d1a3o14_250525_imit_075L_2U_400S",
"n_objects_env": 6,
"obs_slice": "(0, 12)",
"sensor_x": 30.0,
"pinball_front_x": 19.0,
"pinball_rear_x": 20.3,
"target_type": "periodic",
"s_dim": 12,
"s_dim": 14,
"u0": 0.01
}
src/SR_analysis/data/illusion/illusion_0.75L/norm.json
+13 -13
View File
@@ -1,20 +1,20 @@
{
"force_norm_fact": 0.013476977124810219,
"force_norm_fact": 0.013459767680615187,
"sens_deviation": [
0.962617814540863,
-0.12039308249950409,
0.6415857672691345,
0.011103342287242413,
0.9339056611061096,
0.11935960501432419
0.9466423392295837,
-0.1479361653327942,
0.6496055126190186,
-0.06913633644580841,
0.9421071410179138,
0.08593743294477463
],
"sens_norm_fact": [
2.0483264923095703,
2.5809006690979004,
0.7443606853485107,
3.4969263076782227,
2.1811583042144775,
2.5745153427124023
2.1083030700683594,
2.7172951698303223,
0.7220026850700378,
3.7818100452423096,
2.1336052417755127,
2.4038033485412598
],
"action_bias": [
0.0,
src/SR_analysis/data/illusion/illusion_0.75L/result.json
+1 -1
View File
@@ -1,5 +1,5 @@
{
"scene": "illusion_0.75L",
"controlled": true,
"similarity": 0.18393706196948187
"similarity": 0.980430435808052
}
src/SR_analysis/data/illusion/illusion_1.5L/config.json
+4 -3
View File
@@ -3,9 +3,10 @@
"target_diameter": 1.5,
"re_code": 100,
"mu": 0.02,
"nu": 0.008,
"nu": 0.004,
"has_disturbance": false,
"sample_interval": 600,
"sample_interval": 800,
"conv_len": 36,
"action_scale": 8.0,
"action_bias": "(0.0, -2.0, 2.0)",
"source": "PPO_inference",
@@ -17,5 +18,5 @@
"pinball_rear_x": 20.3,
"target_type": "periodic",
"s_dim": 14,
"u0": 0.02
"u0": 0.01
}
src/SR_analysis/data/illusion/illusion_1.5L/norm.json
+13 -13
View File
@@ -1,20 +1,20 @@
{
"force_norm_fact": 0.054184265434741974,
"force_norm_fact": 0.013502256944775581,
"sens_deviation": [
1.9146838188171387,
-0.23843951523303986,
1.305143117904663,
0.0009254463366232812,
1.8709181547164917,
0.2255263477563858
0.9501011371612549,
-0.1168946698307991,
0.6514688730239868,
0.0026063816621899605,
0.9466588497161865,
0.1172625869512558
],
"sens_norm_fact": [
4.165492057800293,
5.171088695526123,
1.5217405557632446,
6.904999732971191,
4.384937763214111,
5.106513023376465
2.11934757232666,
2.5621178150177,
0.7280007004737854,
3.4596621990203857,
2.1354496479034424,
2.5682904720306396
],
"action_bias": [
0.0,
src/SR_analysis/data/illusion/illusion_1.5L/result.json
+2 -2
View File
@@ -1,5 +1,5 @@
{
"scene": "illusion_15L",
"scene": "illusion_1.5L",
"controlled": true,
"similarity": 0.30179914651024675
"similarity": 0.9453180853373244
}
src/SR_analysis/data/illusion/illusion_1L/config.json
+3 -2
View File
@@ -6,16 +6,17 @@
"nu": 0.004,
"has_disturbance": false,
"sample_interval": 600,
"conv_len": 36,
"action_scale": 8.0,
"action_bias": "(0.0, -2.0, 2.0)",
"source": "PPO_inference",
"model_name": "d1a3o12_250525_imit_1L_1U",
"model_name": "d1a3o14_250525_imit_1L_2U_600S",
"n_objects_env": 6,
"obs_slice": "(0, 12)",
"sensor_x": 30.0,
"pinball_front_x": 19.0,
"pinball_rear_x": 20.3,
"target_type": "periodic",
"s_dim": 12,
"s_dim": 14,
"u0": 0.01
}
src/SR_analysis/data/illusion/illusion_1L/norm.json
+13 -13
View File
@@ -1,20 +1,20 @@
{
"force_norm_fact": 0.013487594202160835,
"force_norm_fact": 0.01349810091778636,
"sens_deviation": [
0.9391862154006958,
-0.09573575109243393,
0.6525679230690002,
0.03420928493142128,
0.949753999710083,
0.13210755586624146
0.9381087422370911,
-0.12299147248268127,
0.648796796798706,
-0.022100985050201416,
0.9596271514892578,
0.11032649874687195
],
"sens_norm_fact": [
2.1654911041259766,
2.459106683731079,
0.7431825995445251,
3.613541603088379,
2.1102607250213623,
2.6434059143066406
2.1753337383270264,
2.5954699516296387,
0.707850456237793,
3.554702043533325,
2.0688724517822266,
2.5344114303588867
],
"action_bias": [
0.0,
src/SR_analysis/data/illusion/illusion_1L/result.json
+2 -2
View File
@@ -1,5 +1,5 @@
{
"scene": "illusion_1.0L",
"scene": "illusion_1L",
"controlled": true,
"similarity": 0.5543174409436081
"similarity": 0.9753732477509874
}
src/SR_analysis/pysr.md
-842
View File
@@ -1,842 +0,0 @@
# Toy Examples with Code
## Preamble
```python
import numpy as np
from pysr import *
```
## 1. Simple search
Here's a simple example where we
find the expression `2 cos(x3) + x0^2 - 2`.
```python
X = 2 * np.random.randn(100, 5)
y = 2 * np.cos(X[:, 3]) + X[:, 0] ** 2 - 2
model = PySRRegressor(binary_operators=["+", "-", "*", "/"])
model.fit(X, y)
print(model)
```
## 2. Custom operator
Here, we define a custom operator and use it to find an expression:
```python
X = 2 * np.random.randn(100, 5)
y = 1 / X[:, 0]
model = PySRRegressor(
binary_operators=["+", "*"],
unary_operators=["inv(x) = 1/x"],
extra_sympy_mappings={"inv": lambda x: 1/x},
)
model.fit(X, y)
print(model)
```
## 3. Multiple outputs
Here, we do the same thing, but with multiple expressions at once,
each requiring a different feature.
```python
X = 2 * np.random.randn(100, 5)
y = 1 / X[:, [0, 1, 2]]
model = PySRRegressor(
binary_operators=["+", "*"],
unary_operators=["inv(x) = 1/x"],
extra_sympy_mappings={"inv": lambda x: 1/x},
)
model.fit(X, y)
```
## 4. Plotting an expression
For now, let's consider the expressions for output 0.
We can see the LaTeX version of this with:
```python
model.latex()[0]
```
or output 1 with `model.latex()[1]`.
Let's plot the prediction against the truth:
```python
from matplotlib import pyplot as plt
plt.scatter(y[:, 0], model.predict(X)[:, 0])
plt.xlabel('Truth')
plt.ylabel('Prediction')
plt.show()
```
Which gives us:
![Truth vs Prediction](/images/example_plot.png)
We may also plot the output of a particular expression
by passing the index of the expression to `predict` (or
`sympy` or `latex` as well)
## 5. Feature selection
PySR and evolution-based symbolic regression in general performs
very poorly when the number of features is large.
Even, say, 10 features might be too much for a typical equation search.
If you are dealing with high-dimensional data with a particular type of structure,
you might consider using deep learning to break the problem into
smaller "chunks" which can then be solved by PySR, as explained in the paper
[2006.11287](https://arxiv.org/abs/2006.11287).
For tabular datasets, this is a bit trickier. Luckily, PySR has a built-in feature
selection mechanism. Simply declare the parameter `select_k_features=5`, for selecting
the most important 5 features.
Here is an example. Let's say we have 30 input features and 300 data points, but only 2
of those features are actually used:
```python
X = np.random.randn(300, 30)
y = X[:, 3]**2 - X[:, 19]**2 + 1.5
```
Let's create a model with the feature selection argument set up:
```python
model = PySRRegressor(
binary_operators=["+", "-", "*", "/"],
unary_operators=["exp"],
select_k_features=5,
)
```
Now let's fit this:
```python
model.fit(X, y)
```
Before the Julia backend is launched, you can see the string:
```text
Using features ['x3', 'x5', 'x7', 'x19', 'x21']
```
which indicates that the feature selection (powered by a gradient-boosting tree)
has successfully selected the relevant two features.
This fit should find the solution quickly, whereas with the huge number of features,
it would have struggled.
This simple preprocessing step is enough to simplify our tabular dataset,
but again, for more structured datasets, you should try the deep learning
approach mentioned above.
## 6. Denoising
Many datasets, especially in the observational sciences,
contain intrinsic noise. PySR is noise robust itself, as it is simply optimizing a loss function,
but there are still some additional steps you can take to reduce the effect of noise.
One thing you could do, which we won't detail here, is to create a custom log-likelihood
given some assumed noise model. By passing weights to the fit function, and
defining a custom loss function such as `elementwise_loss="myloss(x, y, w) = w * (x - y)^2"`,
you can define any sort of log-likelihood you wish. (However, note that it must be bounded at zero)
However, the simplest thing to do is preprocessing, just like for feature selection. To do this,
set the parameter `denoise=True`. This will fit a Gaussian process (containing a white noise kernel)
to the input dataset, and predict new targets (which are assumed to be denoised) from that Gaussian process.
For example:
```python
X = np.random.randn(100, 5)
noise = np.random.randn(100) * 0.1
y = np.exp(X[:, 0]) + X[:, 1] + X[:, 2] + noise
```
Let's create and fit a model with the denoising argument set up:
```python
model = PySRRegressor(
binary_operators=["+", "-", "*", "/"],
unary_operators=["exp"],
denoise=True,
)
model.fit(X, y)
print(model)
```
If all goes well, you should find that it predicts the correct input equation, without the noise term!
## 7. Julia packages and types
PySR uses [SymbolicRegression.jl](https://github.com/MilesCranmer/SymbolicRegression.jl)
as its search backend. This is a pure Julia package, and so can interface easily with any other
Julia package.
For some tasks, it may be necessary to load such a package.
For example, let's say we wish to discovery the following relationship:
$$ y = p_{3x + 1} - 5, $$
where $p_i$ is the $i$th prime number, and $x$ is the input feature.
Let's see if we can discover this using
the [Primes.jl](https://github.com/JuliaMath/Primes.jl) package.
First, let's get the Julia backend:
```python
from pysr import jl
```
`jl` stores the Julia runtime.
Now, let's run some Julia code to add the Primes.jl
package to the PySR environment:
```python
jl.seval("""
import Pkg
Pkg.add("Primes")
""")
```
This imports the Julia package manager, and uses it to install
`Primes.jl`. Now let's import `Primes.jl`:
```python
jl.seval("import Primes")
```
Now, we define a custom operator:
```python
jl.seval("""
function p(i::T) where T
if (0.5 < i < 1000)
return T(Primes.prime(round(Int, i)))
else
return T(NaN)
end
end
""")
```
We have created a a function `p`, which takes an arbitrary number as input.
`p` first checks whether the input is between 0.5 and 1000.
If out-of-bounds, it returns `NaN`.
If in-bounds, it rounds it to the nearest integer, compures the corresponding prime number, and then
converts it to the same type as input.
Next, let's generate a list of primes for our test dataset.
Since we are using juliacall, we can just call `p` directly to do this:
```python
primes = {i: jl.p(i*1.0) for i in range(1, 999)}
```
Next, let's use this list of primes to create a dataset of $x, y$ pairs:
```python
import numpy as np
X = np.random.randint(0, 100, 100)[:, None]
y = [primes[3*X[i, 0] + 1] - 5 + np.random.randn()*0.001 for i in range(100)]
```
Note that we have also added a tiny bit of noise to the dataset.
Finally, let's create a PySR model, and pass the custom operator. We also need to define the sympy equivalent, which we can leave as a placeholder for now:
```python
from pysr import PySRRegressor
import sympy
class sympy_p(sympy.Function):
pass
model = PySRRegressor(
binary_operators=["+", "-", "*", "/"],
unary_operators=["p"],
niterations=100,
extra_sympy_mappings={"p": sympy_p}
)
```
We are all set to go! Let's see if we can find the true relation:
```python
model.fit(X, y)
```
if all works out, you should be able to see the true relation (note that the constant offset might not be exactly 1, since it is allowed to round to the nearest integer).
You can get the sympy version of the best equation with:
```python
model.sympy()
```
## 8. Complex numbers
PySR can also search for complex-valued expressions. Simply pass
data with a complex datatype (e.g., `np.complex128`),
and PySR will automatically search for complex-valued expressions:
```python
import numpy as np
X = np.random.randn(100, 1) + 1j * np.random.randn(100, 1)
y = (1 + 2j) * np.cos(X[:, 0] * (0.5 - 0.2j))
model = PySRRegressor(
binary_operators=["+", "-", "*"], unary_operators=["cos"], niterations=100,
)
model.fit(X, y)
```
You can see that all of the learned constants are now complex numbers.
We can get the sympy version of the best equation with:
```python
model.sympy()
```
We can also make predictions normally, by passing complex data:
```python
model.predict(X, -1)
```
to make predictions with the most accurate expression.
## 9. Custom objectives
You can also pass a custom objectives as a snippet of Julia code,
which might include symbolic manipulations or custom functional forms.
These do not even need to be differentiable! First, let's look at the
default objective used (a simplified version, without weights
and with mean square error), so that you can see how to write your own:
```julia
function default_objective(tree, dataset::Dataset{T,L}, options)::L where {T,L}
(prediction, completion) = eval_tree_array(tree, dataset.X, options)
if !completion
return L(Inf)
end
diffs = prediction .- dataset.y
return sum(diffs .^ 2) / length(diffs)
end
```
Here, the `where {T,L}` syntax defines the function for arbitrary types `T` and `L`.
If you have `precision=32` (default) and pass in regular floating point data,
then both `T` and `L` will be equal to `Float32`. If you pass in complex data,
then `T` will be `ComplexF32` and `L` will be `Float32` (since we need to return
a real number from the loss function). But, you don't need to worry about this, just
make sure to return a scalar number of type `L`.
The `tree` argument is the current expression being evaluated. You can read
about the `tree` fields [here](https://ai.damtp.cam.ac.uk/symbolicregression/stable/types/).
For example, let's fix a symbolic form of an expression,
as a rational function. i.e., $P(X)/Q(X)$ for polynomials $P$ and $Q$.
```python
objective = """
function my_custom_objective(tree, dataset::Dataset{T,L}, options) where {T,L}
# Require root node to be binary, so we can split it,
# otherwise return a large loss:
tree.degree != 2 && return L(Inf)
P = tree.l
Q = tree.r
# Evaluate numerator:
P_prediction, flag = eval_tree_array(P, dataset.X, options)
!flag && return L(Inf)
# Evaluate denominator:
Q_prediction, flag = eval_tree_array(Q, dataset.X, options)
!flag && return L(Inf)
# Impose functional form:
prediction = P_prediction ./ Q_prediction
diffs = prediction .- dataset.y
return sum(diffs .^ 2) / length(diffs)
end
"""
model = PySRRegressor(
niterations=100,
binary_operators=["*", "+", "-"],
loss_function=objective,
)
```
> **Warning**: When using a custom objective like this that performs symbolic
> manipulations, many functionalities of PySR will not work, such as `.sympy()`,
> `.predict()`, etc. This is because the SymPy parsing does not know about
> how you are manipulating the expression, so you will need to do this yourself.
Note how we did not pass `/` as a binary operator; it will just be implicit
in the functional form.
Let's generate an equation of the form $\frac{x_0^2 x_1 - 2}{x_2^2 + 1}$:
```python
X = np.random.randn(1000, 3)
y = (X[:, 0]**2 * X[:, 1] - 2) / (X[:, 2]**2 + 1)
```
Finally, let's fit:
```python
model.fit(X, y)
```
> Note that the printed equation is not the same as the evaluated equation,
> because the printing functionality does not know about the functional form.
We can get the string format with:
```python
model.get_best().equation
```
(or, you could use `model.equations_.iloc[-1].equation`)
For me, this equation was:
```text
(((2.3554819 + -0.3554746) - (x1 * (x0 * x0))) - (-1.0000019 - (x2 * x2)))
```
looking at the bracket structure of the equation, we can see that the outermost
bracket is split at the `-` operator (note that we ignore the root operator in
the evaluation, as we simply evaluated each argument and divided the result) into
`((2.3554819 + -0.3554746) - (x1 * (x0 * x0)))` and
`(-1.0000019 - (x2 * x2))`, meaning that our discovered equation is
equal to:
$\frac{x_0^2 x_1 - 2.0000073}{x_2^2 + 1.0000019}$, which
is nearly the same as the true equation!
## 10. Dimensional constraints
One other feature we can exploit is dimensional analysis.
Say that we know the physical units of each feature and output,
and we want to find an expression that is dimensionally consistent.
We can do this as follows, using `DynamicQuantities.jl` to assign units,
passing a string specifying the units for each variable.
First, let's make some data on Newton's law of gravitation, using
astropy for units:
```python
import numpy as np
from astropy import units as u, constants as const
M = (np.random.rand(100) + 0.1) * const.M_sun
m = 100 * (np.random.rand(100) + 0.1) * u.kg
r = (np.random.rand(100) + 0.1) * const.R_earth
G = const.G
F = G * M * m / r**2
```
We can see the units of `F` with `F.unit`.
Now, let's create our model.
Since this data has such a large dynamic range,
let's also create a custom loss function
that looks at the error in log-space:
```python
elementwise_loss = """function loss_fnc(prediction, target)
scatter_loss = abs(log((abs(prediction)+1e-20) / (abs(target)+1e-20)))
sign_loss = 10 * (sign(prediction) - sign(target))^2
return scatter_loss + sign_loss
end
"""
```
Now let's define our model:
```python
model = PySRRegressor(
binary_operators=["+", "-", "*", "/"],
unary_operators=["square"],
elementwise_loss=elementwise_loss,
complexity_of_constants=2,
maxsize=25,
niterations=100,
populations=50,
# Amount to penalize dimensional violations:
dimensional_constraint_penalty=10**5,
)
```
and fit it, passing the unit information.
To do this, we need to use the format of [DynamicQuantities.jl](https://symbolicml.org/DynamicQuantities.jl/dev/#Usage).
```python
# Get numerical arrays to fit:
X = pd.DataFrame(dict(
M=M.to("M_sun").value,
m=m.to("kg").value,
r=r.to("R_earth").value,
))
y = F.value
model.fit(
X,
y,
X_units=["Constants.M_sun", "kg", "Constants.R_earth"],
y_units="kg * m / s^2"
)
```
You can observe that all expressions with a loss under
our penalty are dimensionally consistent!
(The `"[⋅]"` indicates free units in a constant, which can cancel out other units in the expression.)
For example,
```julia
"y[m s⁻² kg] = (M[kg] * 2.6353e-22[⋅])"
```
would indicate that the expression is dimensionally consistent, with
a constant `"2.6353e-22[m s⁻²]"`.
Note that this expression has a large dynamic range so may be difficult to find. Consider searching with a larger `niterations` if needed.
Note that you can also search for exclusively dimensionless constants by settings
`dimensionless_constants_only` to `true`.
## 11. Expression Specifications
PySR 1.0 introduces powerful expression specifications that allow you to define structured equations. Here are two examples:
### Template Expressions
`TemplateExpressionSpec` allows you to define a specific structure for the equation.
For example, let's say we want to learn an equation of the form:
$$ y = \sin(f(x_1, x_2)) + g(x_3) $$
We can do this as follows:
```python
import numpy as np
from pysr import PySRRegressor, TemplateExpressionSpec
# Create data
X = np.random.randn(1000, 3)
y = np.sin(X[:, 0] + X[:, 1]) + X[:, 2]**2
# Define template: we want sin(f(x1, x2)) + g(x3)
template = TemplateExpressionSpec(
expressions=["f", "g"],
variable_names=["x1", "x2", "x3"],
combine="sin(f(x1, x2)) + g(x3)",
)
model = PySRRegressor(
expression_spec=template,
binary_operators=["+", "*", "-", "/"],
unary_operators=["sin"],
maxsize=10,
)
model.fit(X, y)
```
### Parametric Expressions
When your data has categories with shared equation structure but different parameters,
you can use the `parameters` argument of `TemplateExpressionSpec` to specify learned category-specific parameters.
For example, let's say we want to learn an equation of the form:
$$ y = \alpha \sin(x_1) + \beta $$
where $\alpha$ and $\beta$ are different for each category.
Further, let's say we have 3 categories,
with $\alpha \in \{0.1, 1.5, -0.5\}$ and $\beta \in \{1.0, 2.0, 0.5\}$.
```python
import numpy as np
from pysr import PySRRegressor, TemplateExpressionSpec
# Create data with 2 features and 3 categories
X = np.random.uniform(-3, 3, (1000, 2))
category = np.random.randint(0, 3, 1000)
# Parameters for each category
offsets = [0.1, 1.5, -0.5]
scales = [1.0, 2.0, 0.5]
# y = scale[category] * sin(x1) + offset[category]
y = np.array([
scales[c] * np.sin(x1) + offsets[c]
for x1, c in zip(X[:, 0], category)
])
```
Now, let's define our parametric expression:
```python
template = TemplateExpressionSpec(
expressions=["f"],
variable_names=["x1", "x2", "category"],
parameters={"p1": 3, "p2": 3}, # One parameter per category
combine="f(x1, x2, p1[category], p2[category])"
)
```
Next, we pass the category as a _column_ in `X`
corresponding to the index we defined in `variable_names`.
**Note that because Julia is 1-indexed, we need to add 1 to the category index.**
```python
category_p_one = category + 1
X_with_category = np.column_stack([X, category])
```
Now, we can fit our model:
```python
model = PySRRegressor(
expression_spec=template,
binary_operators=["+", "*", "-", "/"],
unary_operators=["sin"],
maxsize=10,
)
model.fit(X_with_category, y)
# Predicting on new data
# model.predict(X_test_with_category)
```
See [Expression Specifications](/api/#expression-specifications) for more details.
You can use this approach for more complex cases,
where you have multiple expressions in the template and parameters that vary by category.
## 12. Using TensorBoard for Logging
You can use TensorBoard to visualize the search progress, as well as
record hyperparameters and final metrics (like `min_loss` and `pareto_volume` - the latter of which
is a performance measure of the entire Pareto front).
```python
import numpy as np
from pysr import PySRRegressor, TensorBoardLoggerSpec
rstate = np.random.RandomState(42)
# Uniform dist between -3 and 3:
X = rstate.uniform(-3, 3, (1000, 2))
y = np.exp(X[:, 0]) + X[:, 1]
# Create a logger that writes to "logs/run*":
logger_spec = TensorBoardLoggerSpec(
log_dir="logs/run",
log_interval=10, # Log every 10 iterations
)
model = PySRRegressor(
binary_operators=["+", "*", "-", "/"],
logger_spec=logger_spec,
)
model.fit(X, y)
```
You can then view the logs with:
```bash
tensorboard --logdir logs/
```
## 13. Vector-valued expressions
You can use `TemplateExpressionSpec` to find expressions for vector-valued data,
where each component might share a common structure.
The trick is to put each vector element into your feature matrix `X`,
and then use a template expression to define the relationships.
For example, say we have 3-dimensional vectors where each component
follows a pattern with a shared term. Say the true model is:
$$\begin{align*}
y_1 &= \exp(x_1) + x_2^2 \\
y_2 &= \exp(x_1) + \sin(x_3) \\
y_3 &= \exp(x_1) + x_1 \cdot x_2
\end{align*}$$
Let's set this up:
```python
import numpy as np
from pysr import PySRRegressor, TemplateExpressionSpec
n = 200
rstate = np.random.RandomState(0)
x1 = rstate.uniform(-2, 2, n)
x2 = rstate.uniform(-2, 2, n)
x3 = rstate.uniform(-2, 2, n)
# True model with shared component exp(x1):
y1 = np.exp(x1) + x2**2
y2 = np.exp(x1) + np.sin(x3)
y3 = np.exp(x1) + x1 * x2
# Add some noise
y1 += 0.05 * rstate.randn(n)
y2 += 0.05 * rstate.randn(n)
y3 += 0.05 * rstate.randn(n)
```
Now, we put everything in `X`; BOTH features and targets:
```python
X = np.column_stack([x1, x2, x3, y1, y2, y3])
```
Now, we can define our template expression:
```python
spec = TemplateExpressionSpec(
expressions=["f1", "f2", "f3", "shared"],
variable_names=["x1", "x2", "x3", "y1", "y2", "y3"],
combine="""
v = shared(x1, x2, x3)
y1_predicted = v + f1(x1, x2, x3)
y2_predicted = v + f2(x1, x2, x3)
y3_predicted = v + f3(x1, x2, x3)
residuals = (
abs2(y1 - y1_predicted) +
abs2(y2 - y2_predicted) +
abs2(y3 - y3_predicted)
)
residuals
"""
)
```
Now, we can fit our model using this template. Since
we already computed the per-row squared error inside the template,
we can pass a dummy `y` to the `fit` method, and also define
an `elementwise_loss` that simply returns the residuals (which get
summed over the data):
```python
model = PySRRegressor(
expression_spec=spec,
binary_operators=["+", "-", "*", "/"],
unary_operators=["exp", "sin"],
maxsize=20,
niterations=50,
elementwise_loss="(pred, target) -> pred",
)
dummy_y = np.zeros(n)
model.fit(X, dummy_y)
```
After running, PySR should find both the shared component (`exp(x1)`) as well as individual components (`square(x2)`, `sin(x3)`, and `x1 * x2`).
You can access the individual expressions through the Julia objects:
```python
# Simply get the expression with the highest score:
idx = model.equations_.score.idxmax()
# Extract the Julia object:
julia_expr = model.equations_.loc[idx, 'julia_expression']
# Access individual subexpressions:
for name in ['f1', 'f2', 'f3', 'shared']:
tree = getattr(julia_expr.trees, name)
print(f"{name}: {tree}")
```
We can also evaluate individual expressions:
```python
from pysr import jl
from pysr.julia_helpers import jl_array
SR = jl.SymbolicRegression
# Get individual trees
f1_tree = julia_expr.trees.f1
shared_tree = julia_expr.trees.shared
# Evaluate at specific points (x1=1, x2=2, x3=3)
test_inputs = jl_array(np.array([[1.0], [2.0], [3.0]]))
f1_result, _ = SR.eval_tree_array(f1_tree, test_inputs, model.julia_options_)
shared_result, _ = SR.eval_tree_array(shared_tree, test_inputs, model.julia_options_)
print(f"f1 at (1,2,3): {f1_result[0]}") # Should be ~4.0 for x2^2
print(f"shared at (1,2,3): {shared_result[0]}") # Should be ~2.718 for exp(1)
```
## 14. Using differential operators
As part of the `TemplateExpressionSpec` described above,
you can also use differential operators within the template.
The operator for this is `D` which takes an expression as the first argument,
and the argument _index_ we are differentiating as the second argument.
This lets you compute integrals via evolution.
For example, let's say we wish to find the integral of $\frac{1}{x^2 \sqrt{x^2 - 1}}$
in the range $x > 1$.
We can compute the derivative of a function $f(x)$, and compare that
to numerical samples of $\frac{1}{x^2\sqrt{x^2-1}}$. Then, by extension,
$f(x)$ represents the indefinite integral of it with some constant offset!
```python
import numpy as np
from pysr import PySRRegressor, TemplateExpressionSpec
x = np.random.uniform(1, 10, (1000,)) # Integrand sampling points
y = 1 / (x**2 * np.sqrt(x**2 - 1)) # Evaluation of the integrand
expression_spec = TemplateExpressionSpec(
expressions=["f"],
variable_names=["x"],
combine="df = D(f, 1); df(x)",
)
model = PySRRegressor(
binary_operators=["+", "-", "*", "/"],
unary_operators=["sqrt"],
expression_spec=expression_spec,
maxsize=20,
)
model.fit(x[:, np.newaxis], y)
```
If everything works, you should find something that simplifies to $\frac{\sqrt{x^2 - 1}}{x}$.
Here, we write out a full function in Julia.
## 15. Additional features
For the many other features available in PySR, please
read the [Options section](options.md).
src/SR_analysis/sindy/illusion/pareto_illusion_1.5L.json
-102
View File
@@ -1,102 +0,0 @@
{
"scene": "illusion_1.5L",
"channels": [
{
"channel": "front",
"best_r2": 0.9594009004329207,
"best_nz": 21,
"pareto": [
{
"nz": 5,
"r2": 0.9393791282052957
},
{
"nz": 6,
"r2": 0.9467043546594699
},
{
"nz": 7,
"r2": 0.9468038158796819
},
{
"nz": 12,
"r2": 0.958013912668809
},
{
"nz": 17,
"r2": 0.9593925359990215
},
{
"nz": 18,
"r2": 0.9593997378572243
},
{
"nz": 21,
"r2": 0.9594009004329207
}
]
},
{
"channel": "top",
"best_r2": 0.9283646632651096,
"best_nz": 22,
"pareto": [
{
"nz": 2,
"r2": 0.8636623835462103
},
{
"nz": 5,
"r2": 0.9176885699701556
},
{
"nz": 6,
"r2": 0.9196961922610852
},
{
"nz": 11,
"r2": 0.9227131504032893
},
{
"nz": 13,
"r2": 0.926100716890473
},
{
"nz": 22,
"r2": 0.9283646632651096
}
]
},
{
"channel": "bottom",
"best_r2": 0.9318647363961834,
"best_nz": 21,
"pareto": [
{
"nz": 3,
"r2": 0.7872211556048209
},
{
"nz": 6,
"r2": 0.9223640246817683
},
{
"nz": 7,
"r2": 0.9258419638948643
},
{
"nz": 11,
"r2": 0.929107795801212
},
{
"nz": 16,
"r2": 0.9315566670173666
},
{
"nz": 21,
"r2": 0.9318647363961834
}
]
}
]
}
src/SR_analysis/sindy/illusion/pareto_illusion_1L.json
-110
View File
@@ -1,110 +0,0 @@
{
"scene": "illusion_1L",
"channels": [
{
"channel": "front",
"best_r2": 0.9793036424523165,
"best_nz": 21,
"pareto": [
{
"nz": 4,
"r2": 0.9718953011650306
},
{
"nz": 6,
"r2": 0.9752101462860064
},
{
"nz": 11,
"r2": 0.9785538576893853
},
{
"nz": 15,
"r2": 0.9791278336272012
},
{
"nz": 18,
"r2": 0.9792726941557117
},
{
"nz": 21,
"r2": 0.9793036424523165
}
]
},
{
"channel": "top",
"best_r2": 0.9838580289472151,
"best_nz": 22,
"pareto": [
{
"nz": 7,
"r2": 0.9786186299815743
},
{
"nz": 10,
"r2": 0.9828568398768021
},
{
"nz": 11,
"r2": 0.9833618417657272
},
{
"nz": 12,
"r2": 0.9835181072357578
},
{
"nz": 16,
"r2": 0.9837742843659423
},
{
"nz": 22,
"r2": 0.9838580289472151
}
]
},
{
"channel": "bottom",
"best_r2": 0.983658612471297,
"best_nz": 22,
"pareto": [
{
"nz": 4,
"r2": 0.9698703453821834
},
{
"nz": 6,
"r2": 0.9816905531339832
},
{
"nz": 7,
"r2": 0.9817463485828739
},
{
"nz": 8,
"r2": 0.9822685326747084
},
{
"nz": 11,
"r2": 0.9831202415012674
},
{
"nz": 12,
"r2": 0.983291741316277
},
{
"nz": 15,
"r2": 0.9836126332791877
},
{
"nz": 18,
"r2": 0.9836582680775273
},
{
"nz": 22,
"r2": 0.983658612471297
}
]
}
]
}
src/SR_analysis/sindy/illusion/sindy_results.json
-1486
View File
File diff suppressed because it is too large Load Diff
src/SR_analysis/sindy/karman/pareto_karman_re100.json
-118
View File
@@ -1,118 +0,0 @@
{
"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
}
]
}
]
}
src/SR_analysis/sindy/karman/sindy_results.json
-2031
View File
File diff suppressed because it is too large Load Diff
src/SR_analysis/sindy/run_all_v2.py
-766
View File
@@ -1,766 +0,0 @@
#!/usr/bin/env python3
"""Unified SINDy fitting v2 for all scenes.
Uses V2 feature builder: no sin_ua/cos_ua, optional mu, weighted STLSQ.
Generates separate per-scene results and cross-scene comparisons.
Usage:
conda run -n pycuda_3_10 python sindy/run_all_v2.py \\
--scenes karman_re50,karman_re100,steady
# Karman cross-Re only
conda run -n pycuda_3_10 python sindy/run_all_v2.py --family karman
# Cloak family (Karman + steady + vortex)
conda run -n pycuda_3_10 python sindy/run_all_v2.py --family cloak
# All scenes
conda run -n pycuda_3_10 python sindy/run_all_v2.py --family all
"""
from __future__ import annotations
import argparse
import json
import os
import sys
from typing import List, Optional
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
if _REPO not in sys.path:
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from SR_analysis.utils.sindy_fitter import (
fit_sindy_weighted, get_feature_matrix_v2, get_active_support,
get_feature_matrix_deriv, compute_action_deriv,
)
from SR_analysis.utils.feature_builder import CORE_FEAT_KEYS_V2, ALL_FEAT_KEYS_V2, PHYSICS_FEAT_KEYS
from SR_analysis.configs import get_scene, get_scene_list, SCENES
# Base output directory
SINDY_DIR = os.path.join(os.path.dirname(__file__))
# Threshold grid (more granular at low end for better Pareto selection)
THRESHOLDS = [0.0, 0.001, 0.002, 0.003, 0.005, 0.007, 0.01, 0.015, 0.02, 0.03, 0.05, 0.08, 0.1]
FAMILIES = {
"cloak": ["karman", "steady", "vortex"],
"karman": ["karman"],
"illusion": ["illusion"],
"vortex": ["vortex"],
"all": None, # all scenes
}
def load_data(scene_name: str) -> tuple:
"""Load sensors/forces/actions from scene's controlled.npz.
Returns actions in PHYSICAL omega (lattice units).
Also returns target_forces if available (for Illusion scenes).
"""
cfg = get_scene(scene_name)
data_dir = os.path.join(
os.path.dirname(__file__), "..", "data",
cfg["scene_id"], scene_name,
)
npz = np.load(os.path.join(data_dir, "controlled.npz"))
sensors = npz["sensors"].astype(np.float64)
forces = npz["forces"].astype(np.float64)
# actions in .npz are normalized [-1,1]; convert to physical omega
actions_norm = npz["actions"].astype(np.float64)
scale = cfg["action_scale"]
bias = np.array(cfg["action_bias"], dtype=np.float64)
u0 = cfg["u0"]
actions_phys = (actions_norm * scale + bias) * u0
# Load target_forces if available (Illusion scenes)
target_forces = None
if "target_forces" in npz:
target_forces = npz["target_forces"].astype(np.float64)
return sensors, forces, actions_phys, cfg, target_forces
def compute_scene_weight(scene_name: str, cfg: dict) -> float:
"""Compute scene quality weight from PPO similarity."""
result_path = os.path.join(
os.path.dirname(__file__), "..", "data",
cfg["scene_id"], scene_name, "result.json",
)
if os.path.isfile(result_path):
with open(result_path) as f:
r = json.load(f)
sim = r.get("similarity", r.get("avg_reward_last100", 0.5))
return float(sim) ** 2
return 1.0
def run_single_scene(
scene_name: str,
thresholds: Optional[List[float]] = None,
verbose: bool = True,
) -> dict:
"""Run SINDy v2 on a single scene. Returns result dict."""
if thresholds is None:
thresholds = THRESHOLDS
sensors, forces, actions_phys, cfg, target_forces = load_data(scene_name)
mu = cfg["mu"]
u0 = cfg["u0"]
scene_id = cfg["scene_id"]
T = sensors.shape[0]
if verbose:
print(f"\n{'='*60}")
print(f"Scene: {scene_name} ({scene_id}) T={T} mu={mu:.4f} u0={u0}")
print(f"{'='*60}")
# Determine if this is a single-scene (no mu) or cross-Re (with mu)
# For cross-Re, we'll handle separately; for single-scene, no mu
# But steady is special: it needs mu=constant, still no mu in features
# Single scene: no mu, no sin/cos
is_single_re = True # will be overridden for joint fitting
# Build feature matrix (V2: no sin/cos, no mu for single-scene)
# Use dt_c=sample_interval/T0 in T0 units
t0_steps = 2000 # T0 = D/U0 = 20/0.01 = 2000 LBM steps
dt_c = cfg.get("sample_interval", 800) / t0_steps
Theta_f, Theta_r, Y, fn_f, fn_r = get_feature_matrix_v2(
sensors, forces, actions_phys, mu,
u0=u0, include_mu=False, include_cos_sin=False,
use_time_norm=False, dt_c=dt_c, n_warmup=2,
target_forces=target_forces,
)
if verbose:
print(f" Features: front {Theta_f.shape[1]-1} + bias?0, rear {Theta_r.shape[1]-1} + bias?1")
print(f" dt_c={dt_c:.4f} (T0 units)")
# Scene quality weight for STLSQ
w_scene = compute_scene_weight(scene_name, cfg)
if verbose:
print(f" Scene quality weight: {w_scene:.4f}")
# ---------- Front channel (no bias) ----------
if verbose:
print(f"\n --- Front (no bias) ---")
front_results = fit_sindy_weighted(Theta_f, Y[:, 0], thresholds,
n_robust_passes=2)
# Find best by R2 (include th=0.0 for dense coefficients)
front_best = max(front_results, key=lambda r: r["r2"])
if verbose:
_print_active(fn_f, front_best["coef"])
# ---------- Top channel (rear shared-head, with bias) ----------
if verbose:
print(f"\n --- Top (rear shared-head, with bias) ---")
top_results = fit_sindy_weighted(Theta_r, Y[:, 2], thresholds,
n_robust_passes=2)
top_best = max(top_results, key=lambda r: r["r2"])
if verbose:
_print_active(fn_r, top_best["coef"])
# ---------- Bottom (independent, for comparison) ----------
if verbose:
print(f"\n --- Bottom (independent, with bias) ---")
bot_results = fit_sindy_weighted(Theta_r, Y[:, 1], thresholds,
n_robust_passes=2)
bot_best = max(bot_results, key=lambda r: r["r2"])
if verbose:
_print_active(fn_r, bot_best["coef"])
# Compute additional metrics for the best thresholds
n_warmup = 2
n_samples = Theta_f.shape[0]
result = {
"scene": scene_name,
"scene_id": scene_id,
"re_code": cfg.get("re_code", 100),
"mu": mu,
"u0": u0,
"dt_c": dt_c,
"n_samples": n_samples,
"thresholds": thresholds,
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"front": {
"results": _clean_results(front_results),
"best": _clean_single(front_best),
"best_coef": front_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in front_results],
},
"top": {
"results": _clean_results(top_results),
"best": _clean_single(top_best),
"best_coef": top_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in top_results],
},
"bottom": {
"results": _clean_results(bot_results),
"best": _clean_single(bot_best),
"best_coef": bot_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in bot_results],
},
}
if verbose:
print(f"\n Summary:")
print(f" Front: th={front_best['threshold']:.4f} nz={front_best['nz']:2d} R2={front_best['r2']:.4f}")
print(f" Top: th={top_best['threshold']:.4f} nz={top_best['nz']:2d} R2={top_best['r2']:.4f}")
print(f" Bottom: th={bot_best['threshold']:.4f} nz={bot_best['nz']:2d} R2={bot_best['r2']:.4f}")
return result
def _print_active(names: List[str], coef: list, rtol: float = 1e-4):
"""Print active features from coefficient vector."""
for i, c in enumerate(coef):
if abs(c) > rtol and i < len(names):
print(f" {names[i]:20s} = {c:+.6f}")
def _clean_results(results: List[dict]) -> List[dict]:
"""Keep all keys including coef."""
return [dict(r) for r in results]
def _clean_single(entry: dict) -> dict:
"""Remove coef from single entry."""
return {k: v for k, v in entry.items() if k != "coef"}
def run_joint_karman(
scene_names: List[str],
thresholds: Optional[List[float]] = None,
) -> dict:
"""Joint Karman cross-Re SINDy with mu terms and scene weighting."""
if thresholds is None:
thresholds = THRESHOLDS
print(f"\n{'='*60}")
print(f"Joint Karman cross-Re: {scene_names}")
print(f"{'='*60}")
# Load and concatenate all Karman scenes
all_sensors, all_forces, all_actions = [], [], []
all_weights = []
t0_steps = 2000
for sn in scene_names:
sensors, forces, actions, cfg, _ = load_data(sn)
mu = cfg["mu"]
u0 = cfg["u0"]
dt_c = cfg.get("sample_interval", 800) / t0_steps
# Build features WITH mu terms
Theta_f, Theta_r, Y, fn_f, fn_r = get_feature_matrix_v2(
sensors, forces, actions, mu,
u0=u0, include_mu=True, include_cos_sin=False,
use_time_norm=False, dt_c=dt_c, n_warmup=2,
)
w = compute_scene_weight(sn, cfg)
n = Theta_f.shape[0]
print(f" {sn}: {n} samples, weight={w:.4f}")
all_weights.append(np.full(n, w))
all_sensors.append(Theta_f)
all_forces.append(Theta_r) # reuse variable names for concat
all_actions.append(Y)
# Concatenate
Theta_f_j = np.vstack(all_sensors)
Theta_r_j = np.vstack(all_forces)
Y_j = np.vstack(all_actions)
W_j = np.concatenate(all_weights)
print(f" Joint: front {Theta_f_j.shape}, rear {Theta_r_j.shape}")
print(f" Weights: min={W_j.min():.4f} max={W_j.max():.4f}")
# ---------- Front (no bias, all features including mu) ----------
print(f"\n --- Joint Front (no bias, with mu) ---")
front_results = fit_sindy_weighted(Theta_f_j, Y_j[:, 0], thresholds,
sample_weights=W_j, n_robust_passes=2)
front_best = max(front_results, key=lambda r: r["r2"])
_print_active(fn_f, front_best["coef"])
# ---------- Top (with bias, all features including mu) ----------
print(f"\n --- Joint Top (with bias, with mu) ---")
top_results = fit_sindy_weighted(Theta_r_j, Y_j[:, 2], thresholds,
sample_weights=W_j, n_robust_passes=2)
top_best = max(top_results, key=lambda r: r["r2"])
_print_active(fn_r, top_best["coef"])
# ---------- Bottom (independent) ----------
print(f"\n --- Joint Bottom (independent) ---")
bot_results = fit_sindy_weighted(Theta_r_j, Y_j[:, 1], thresholds,
sample_weights=W_j, n_robust_passes=2)
bot_best = max(bot_results, key=lambda r: r["r2"])
_print_active(fn_r, bot_best["coef"])
return {
"joint_scenes": scene_names,
"thresholds": thresholds,
"n_samples_total": Theta_f_j.shape[0],
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"front": {
"results": _clean_results(front_results),
"best": _clean_single(front_best),
"best_coef": front_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in front_results],
},
"top": {
"results": _clean_results(top_results),
"best": _clean_single(top_best),
"best_coef": top_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in top_results],
},
"bottom": {
"results": _clean_results(bot_results),
"best": _clean_single(bot_best),
"best_coef": bot_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in bot_results],
},
}
def run_joint_illusion(
scene_names: List[str],
thresholds: Optional[List[float]] = None,
) -> dict:
"""Joint Illusion cross-diameter SINDy with target features and scene weighting."""
if thresholds is None:
thresholds = THRESHOLDS
print(f"\n{'='*60}")
print(f"Joint Illusion cross-diameter: {scene_names}")
print(f"{'='*60}")
all_sensors, all_forces, all_actions = [], [], []
all_weights = []
t0_steps = 2000
for sn in scene_names:
sensors, forces, actions, cfg, target_forces = load_data(sn)
mu = cfg["mu"]
u0 = cfg["u0"]
dt_c = cfg.get("sample_interval", 600) / t0_steps
# Build features WITH target features, no mu
Theta_f, Theta_r, Y, fn_f, fn_r = get_feature_matrix_v2(
sensors, forces, actions, mu,
u0=u0, include_mu=False, include_cos_sin=False,
use_time_norm=False, dt_c=dt_c, n_warmup=2,
target_forces=target_forces,
)
w = compute_scene_weight(sn, cfg)
n = Theta_f.shape[0]
print(f" {sn}: {n} samples, dt_c={dt_c:.4f}, weight={w:.4f}")
all_weights.append(np.full(n, w))
all_sensors.append(Theta_f)
all_forces.append(Theta_r)
all_actions.append(Y)
Theta_f_j = np.vstack(all_sensors)
Theta_r_j = np.vstack(all_forces)
Y_j = np.vstack(all_actions)
W_j = np.concatenate(all_weights)
print(f" Joint: front {Theta_f_j.shape}, rear {Theta_r_j.shape}")
print(f" Weights: min={W_j.min():.4f} max={W_j.max():.4f}")
# Front (no bias)
print(f"\n --- Joint Front (no bias) ---")
front_results = fit_sindy_weighted(Theta_f_j, Y_j[:, 0], thresholds,
sample_weights=W_j, n_robust_passes=2)
front_best = max(front_results, key=lambda r: r["r2"])
_print_active(fn_f, front_best["coef"])
# Top (with bias)
print(f"\n --- Joint Top (with bias) ---")
top_results = fit_sindy_weighted(Theta_r_j, Y_j[:, 2], thresholds,
sample_weights=W_j, n_robust_passes=2)
top_best = max(top_results, key=lambda r: r["r2"])
_print_active(fn_r, top_best["coef"])
# Bottom (independent)
print(f"\n --- Joint Bottom (independent) ---")
bot_results = fit_sindy_weighted(Theta_r_j, Y_j[:, 1], thresholds,
sample_weights=W_j, n_robust_passes=2)
bot_best = max(bot_results, key=lambda r: r["r2"])
_print_active(fn_r, bot_best["coef"])
return {
"joint_scenes": scene_names,
"thresholds": thresholds,
"n_samples_total": Theta_f_j.shape[0],
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"front": {
"results": _clean_results(front_results),
"best": _clean_single(front_best),
"best_coef": front_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in front_results],
},
"top": {
"results": _clean_results(top_results),
"best": _clean_single(top_best),
"best_coef": top_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in top_results],
},
"bottom": {
"results": _clean_results(bot_results),
"best": _clean_single(bot_best),
"best_coef": bot_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in bot_results],
},
}
def run_single_scene_deriv(
scene_name: str,
thresholds: Optional[List[float]] = None,
verbose: bool = True,
center_diff: bool = False,
augment_level: int = 0,
feat_keys: Optional[List[str]] = None,
output_mode: str = "deriv",
) -> dict:
"""Run SINDy in phase-state mode.
output_mode:
"deriv": Y = d(alpha)/dt (time-normalized). Closed-loop needs integration.
"absolute": Y = alpha (non-dimensional absolute action). No integration.
"""
if thresholds is None:
thresholds = THRESHOLDS
sensors, forces, actions_phys, cfg, target_forces = load_data(scene_name)
mu = cfg["mu"]
u0 = cfg["u0"]
scene_id = cfg["scene_id"]
T = sensors.shape[0]
t0_steps = 2000
dt_c = cfg.get("sample_interval", 800) / t0_steps
if verbose:
print(f"\n{'='*60}")
print(f"{'ABSOLUTE' if output_mode == 'absolute' else 'DERIV'} Scene: {scene_name} dt_c={dt_c:.4f} u0={u0}")
print(f"{'='*60}")
Theta_f, Theta_r, Y, fn_f, fn_r = get_feature_matrix_deriv(
sensors, forces, actions_phys, mu,
u0=u0, dt_c=dt_c,
target_forces=target_forces,
center_diff=center_diff,
augment_level=augment_level,
feat_keys=feat_keys,
output_mode=output_mode,
include_mu=feat_keys is not None and any("mu_" in k for k in feat_keys),
)
y_label = "alpha (abs)" if output_mode == "absolute" else "d(alpha)/dt"
if verbose:
print(f" Features: front {Theta_f.shape[1]} (no bias), rear {Theta_r.shape[1]} (with bias)")
print(f" Fitting {y_label}: {Y.shape}")
w_scene = compute_scene_weight(scene_name, cfg)
# Front (no bias)
if verbose:
print(f"\n --- Front (no bias, target = {y_label}_F) ---")
front_results = fit_sindy_weighted(Theta_f, Y[:, 0], thresholds, n_robust_passes=2)
front_best = max(front_results, key=lambda r: r["r2"])
if verbose:
_print_active(fn_f, front_best["coef"])
print(f" R2={front_best['r2']:.4f} nz={front_best['nz']}")
# Top (with bias, rear shared-head)
if verbose:
print(f"\n --- Top (with bias, target = {y_label}_T) ---")
top_results = fit_sindy_weighted(Theta_r, Y[:, 2], thresholds, n_robust_passes=2)
top_best = max(top_results, key=lambda r: r["r2"])
if verbose:
_print_active(fn_r, top_best["coef"])
print(f" R2={top_best['r2']:.4f} nz={top_best['nz']}")
# Bottom (independent)
if verbose:
print(f"\n --- Bottom (independent, target = {y_label}_B) ---")
bot_results = fit_sindy_weighted(Theta_r, Y[:, 1], thresholds, n_robust_passes=2)
bot_best = max(bot_results, key=lambda r: r["r2"])
if verbose:
_print_active(fn_r, bot_best["coef"])
print(f" R2={bot_best['r2']:.4f} nz={bot_best['nz']}")
result = {
"scene": scene_name,
"scene_id": scene_id,
"mode": output_mode,
"dt_c": dt_c,
"center_diff": center_diff,
"re_code": cfg.get("re_code", 100),
"mu": mu,
"u0": u0,
"n_samples": Theta_f.shape[0],
"thresholds": thresholds,
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"front": {
"results": _clean_results(front_results),
"best": _clean_single(front_best),
"best_coef": front_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in front_results],
},
"top": {
"results": _clean_results(top_results),
"best": _clean_single(top_best),
"best_coef": top_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in top_results],
},
"bottom": {
"results": _clean_results(bot_results),
"best": _clean_single(bot_best),
"best_coef": bot_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in bot_results],
},
}
return result
def run_joint_karman_deriv(
scene_names: List[str],
thresholds: Optional[List[float]] = None,
augment_level: int = 0,
) -> dict:
"""Joint Karman cross-Re derivative-mode SINDy."""
if thresholds is None:
thresholds = THRESHOLDS
print(f"\n{'='*60}")
print(f"DERIV Joint Karman: {scene_names}")
print(f"{'='*60}")
t0_steps = 2000
all_Theta_f, all_Theta_r, all_Y = [], [], []
all_weights = []
for sn in scene_names:
sensors, forces, actions, cfg, _ = load_data(sn)
mu = cfg["mu"]
u0 = cfg["u0"]
dt_c = cfg.get("sample_interval", 800) / t0_steps
Theta_f, Theta_r, Y, fn_f, fn_r = get_feature_matrix_deriv(
sensors, forces, actions, mu,
u0=u0, dt_c=dt_c, include_mu=True,
target_forces=None, augment_level=augment_level,
)
w = compute_scene_weight(sn, cfg)
n = Theta_f.shape[0]
print(f" {sn}: {n} samples, dt_c={dt_c:.4f}, weight={w:.4f}")
all_weights.append(np.full(n, w))
all_Theta_f.append(Theta_f)
all_Theta_r.append(Theta_r)
all_Y.append(Y)
Theta_f_j = np.vstack(all_Theta_f)
Theta_r_j = np.vstack(all_Theta_r)
Y_j = np.vstack(all_Y)
W_j = np.concatenate(all_weights)
print(f" Joint: front {Theta_f_j.shape}, rear {Theta_r_j.shape}")
# Front
print(f"\n --- Joint Front ---")
front_results = fit_sindy_weighted(Theta_f_j, Y_j[:, 0], thresholds, sample_weights=W_j, n_robust_passes=2)
front_best = max(front_results, key=lambda r: r["r2"])
_print_active(fn_f, front_best["coef"])
print(f" R2={front_best['r2']:.4f}")
# Top
print(f"\n --- Joint Top ---")
top_results = fit_sindy_weighted(Theta_r_j, Y_j[:, 2], thresholds, sample_weights=W_j, n_robust_passes=2)
top_best = max(top_results, key=lambda r: r["r2"])
_print_active(fn_r, top_best["coef"])
print(f" R2={top_best['r2']:.4f}")
# Bottom
print(f"\n --- Joint Bottom ---")
bot_results = fit_sindy_weighted(Theta_r_j, Y_j[:, 1], thresholds, sample_weights=W_j, n_robust_passes=2)
bot_best = max(bot_results, key=lambda r: r["r2"])
_print_active(fn_r, bot_best["coef"])
print(f" R2={bot_best['r2']:.4f}")
return {
"joint_scenes": scene_names,
"mode": "deriv",
"thresholds": thresholds,
"n_samples_total": Theta_f_j.shape[0],
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"front": {
"results": _clean_results(front_results),
"best": _clean_single(front_best),
"best_coef": front_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in front_results],
},
"top": {
"results": _clean_results(top_results),
"best": _clean_single(top_best),
"best_coef": top_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in top_results],
},
"bottom": {
"results": _clean_results(bot_results),
"best": _clean_single(bot_best),
"best_coef": bot_best["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in bot_results],
},
}
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--scenes", type=str, default=None,
help="Comma-separated scene names")
ap.add_argument("--family", type=str, default=None,
choices=list(FAMILIES.keys()) + [None],
help="Scene family")
ap.add_argument("--out", type=str, default=None,
help="Output directory (default: sindy/<scene_id>/")
ap.add_argument("--joint", action="store_true",
help="Run joint Karman cross-Re fitting")
ap.add_argument("--deriv", action="store_true",
help="Use derivative mode: fit d(alpha)/dt = g(physics_state)")
ap.add_argument("--center-diff", action="store_true",
help="Use centered difference for derivative (vs forward diff)")
ap.add_argument("--augment-level", type=int, default=0,
help="Observation augmentation: 0=static, 1=+lags, 2=+derivs, 3=both, 4=+action_lag")
ap.add_argument("--phase", action="store_true",
help="Use phase-state features: [u_a, du_a/dt, Cl_tot, dCl_tot/dt, Cd_tot, Cd_rear]")
ap.add_argument("--karman-expand", action="store_true",
help="Karman expanded: phase-state + u_m/u_c/v_a/Cl_diff (10 dim)")
ap.add_argument("--karman-mu", action="store_true",
help="Karman phase-state + mu modulation: 6-dim + mu*Cl_tot/mu*Cd_tot/mu*u_a (9 dim)")
ap.add_argument("--output-mode", type=str, default="deriv",
choices=["deriv", "absolute"],
help="'deriv': predict d(alpha)/dt; 'absolute': predict alpha directly")
args = ap.parse_args()
# Resolve scene list
if args.scenes:
scene_names = [s.strip() for s in args.scenes.split(",")]
elif args.family and args.family in FAMILIES:
family_ids = FAMILIES[args.family]
if family_ids is None:
scene_names = get_scene_list() # all
else:
scene_names = []
for fid in family_ids:
scene_names.extend(get_scene_list(fid))
else:
scene_names = get_scene_list("karman")
# Determine fitting function based on mode
if args.deriv:
om = args.output_mode
if args.phase:
from SR_analysis.utils.feature_builder import PHASE_STATE_KEYS, ILLUSION_PHASE_KEYS
def _phase_func(sn):
from SR_analysis.configs import get_scene
cfg = get_scene(sn)
feat_keys = ILLUSION_PHASE_KEYS if cfg["scene_id"] == "illusion" else PHASE_STATE_KEYS
return run_single_scene_deriv(sn, center_diff=args.center_diff,
augment_level=args.augment_level,
feat_keys=feat_keys, output_mode=om)
fit_func = _phase_func
elif args.karman_expand:
from SR_analysis.utils.feature_builder import KARMAN_EXPANDED_KEYS
fit_func = lambda sn: run_single_scene_deriv(sn, center_diff=args.center_diff,
augment_level=args.augment_level,
feat_keys=list(KARMAN_EXPANDED_KEYS),
output_mode=om)
elif args.karman_mu:
from SR_analysis.utils.feature_builder import PHASE_STATE_KEYS
# Phase-state + 3 selected mu modulation features
mu_feat = ["mu_Cl_tot", "mu_Cd_tot", "mu_u_a"]
feat_keys = list(PHASE_STATE_KEYS) + mu_feat
fit_func = lambda sn: run_single_scene_deriv(sn, center_diff=args.center_diff,
augment_level=args.augment_level,
feat_keys=feat_keys,
output_mode=om)
else:
fit_func = lambda sn: run_single_scene_deriv(sn, center_diff=args.center_diff,
augment_level=args.augment_level,
output_mode=om)
joint_func_karman = lambda names: run_joint_karman_deriv(names, augment_level=args.augment_level)
else:
fit_func = run_single_scene
joint_func_karman = run_joint_karman
print(f"Running SINDy {'deriv' if args.deriv else 'v2'} for scenes: {scene_names}")
# Run per-scene
for sn in scene_names:
result = fit_func(sn)
# Save per-scene to its scene_id directory
scene_id = result["scene_id"]
out_dir = os.path.join(SINDY_DIR, scene_id)
os.makedirs(out_dir, exist_ok=True)
suffix = "_deriv" if args.deriv else "_v2"
out_path = os.path.join(out_dir, f"sindy_results{suffix}.json")
# Load existing if any, update per_scene
if os.path.isfile(out_path):
with open(out_path) as f:
existing = json.load(f)
else:
existing = {
"thresholds": THRESHOLDS,
"all_feature_names_front": result["feature_names_front"],
"all_feature_names_rear": result["feature_names_rear"],
"per_scene": {},
}
existing["per_scene"][sn] = result
with open(out_path, "w") as f:
json.dump(existing, f, indent=2)
print(f" Saved: {out_path}")
# Joint Karman cross-Re (training Re only: 50, 100, 200, 400)
if args.joint:
karman_train = ["karman_re50", "karman_re100", "karman_re200", "karman_re400"]
joint_result = joint_func_karman(karman_train)
joint_suffix = "_joint_deriv" if args.deriv else "_joint_v2"
out_path = os.path.join(SINDY_DIR, "karman", f"sindy{joint_suffix}.json")
with open(out_path, "w") as f:
json.dump(joint_result, f, indent=2)
print(f"\nSaved joint: {out_path}")
# Joint Illusion (only in non-deriv mode for now)
if args.joint and not args.deriv:
illusion_names = get_scene_list("illusion")
joint_illusion = run_joint_illusion(illusion_names)
out_path = os.path.join(SINDY_DIR, "illusion", "sindy_joint_v2.json")
os.makedirs(os.path.dirname(out_path), exist_ok=True)
with open(out_path, "w") as f:
json.dump(joint_illusion, f, indent=2)
print(f"\nSaved joint illusion: {out_path}")
if __name__ == "__main__":
main()
src/SR_analysis/sindy/run_illusion.py
-133
View File
@@ -1,133 +0,0 @@
"""SINDy fitting for Illusion scenes.
Usage:
conda run -n pycuda_3_10 python sindy/run_illusion.py
conda run -n pycuda_3_10 python sindy/run_illusion.py --diameters 0.75,1.0
"""
from __future__ import annotations
import argparse
import json
import os
import sys
from typing import List, Optional
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
if _REPO not in sys.path:
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from SR_analysis.utils.sindy_fitter import fit_sindy, get_feature_matrix_from_data
from SR_analysis.configs import get_scene, get_scene_list
SINDY_DIR = os.path.join(os.path.dirname(__file__), "..", "sindy", "illusion")
THRESHOLDS = [0.0, 0.001, 0.002, 0.005, 0.01, 0.015, 0.02, 0.03, 0.05, 0.1]
def load_data(scene_name: str) -> tuple:
data_dir = os.path.join(os.path.dirname(__file__), "..", "data", "illusion", scene_name)
npz = np.load(os.path.join(data_dir, "controlled.npz"))
sensors = npz["sensors"].astype(np.float64)
forces = npz["forces"].astype(np.float64)
actions = npz["actions"].astype(np.float64)
return sensors, forces, actions
def run(scene_names: Optional[List[str]] = None):
if scene_names is None:
scene_names = get_scene_list("illusion")
per_scene = {}
for sn in scene_names:
print(f"\n{'='*60}")
print(f"Scene: {sn}")
print(f"{'='*60}")
cfg = get_scene(sn)
sensors, forces, actions_phys = load_data(sn)
mu = cfg["mu"]
print(f" T={sensors.shape[0]}, mu={mu:.6f}")
Theta_f, Theta_r, Y, fn_f, fn_r = get_feature_matrix_from_data(
sensors, forces, actions_phys, mu, u0=cfg["u0"],
alpha_mode=False, include_mu=True, n_warmup=2,
)
print(f" Front: {Theta_f.shape}, Rear: {Theta_r.shape}")
# Front channel
print(f"\n --- Front (no bias) ---")
front_results = fit_sindy(Theta_f, Y[:, 0], THRESHOLDS)
best_f = max(front_results, key=lambda r: r["r2"])
print(f" Best: th={best_f['threshold']:.4f} nz={best_f['nz']:2d} R2={best_f['r2']:.6f}")
# Top channel (rear shared-head)
print(f"\n --- Top (rear shared-head) ---")
top_results = fit_sindy(Theta_r, Y[:, 2], THRESHOLDS)
best_t = max(top_results, key=lambda r: r["r2"])
print(f" Best: th={best_t['threshold']:.4f} nz={best_t['nz']:2d} R2={best_t['r2']:.6f}")
# Bottom (independent)
print(f"\n --- Bottom (independent) ---")
bot_results = fit_sindy(Theta_r, Y[:, 1], THRESHOLDS)
best_b = max(bot_results, key=lambda r: r["r2"])
print(f" Best: th={best_b['threshold']:.4f} nz={best_b['nz']:2d} R2={best_b['r2']:.6f}")
per_scene[sn] = {
"scene": sn,
"re_code": cfg["re_code"],
"mu": mu,
"n_samples": Theta_f.shape[0],
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"front": {
"results": [{k: v for k, v in r.items() if k != "coef"} for r in front_results],
"best": {k: v for k, v in best_f.items() if k != "coef"},
"best_coef": best_f["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in front_results],
},
"top": {
"results": [{k: v for k, v in r.items() if k != "coef"} for r in top_results],
"best": {k: v for k, v in best_t.items() if k != "coef"},
"best_coef": best_t["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in top_results],
},
"bottom": {
"results": [{k: v for k, v in r.items() if k != "coef"} for r in bot_results],
"best": {k: v for k, v in best_b.items() if k != "coef"},
"best_coef": best_b["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in bot_results],
},
}
os.makedirs(SINDY_DIR, exist_ok=True)
out_path = os.path.join(SINDY_DIR, "sindy_results.json")
result = {"thresholds": THRESHOLDS, "per_scene": per_scene}
with open(out_path, "w") as f:
json.dump(result, f, indent=2)
print(f"\nSaved: {out_path}")
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--diameters", type=str, default=None,
help="Comma-separated diameters (e.g. 0.75,1.0,1.5)")
ap.add_argument("--scene-names", type=str, default=None)
args = ap.parse_args()
if args.scene_names:
names = [s.strip() for s in args.scene_names.split(",")]
elif args.diameters:
names = [f"illusion_{d.strip()}L" for d in args.diameters.split(",")]
else:
names = None
run(names)
if __name__ == "__main__":
main()
src/SR_analysis/sindy/run_karman.py
-159
View File
@@ -1,159 +0,0 @@
"""SINDy fitting for Karman cloak scenes.
Runs STLSQ threshold grid for Karman scenes (training Re or all Re).
Usage:
conda run -n pycuda_3_10 python sindy/run_karman.py --re-codes 50,100,200
conda run -n pycuda_3_10 python sindy/run_karman.py
"""
from __future__ import annotations
import argparse
import json
import os
import sys
from typing import List, Optional
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
if _REPO not in sys.path:
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from SR_analysis.utils.sindy_fitter import fit_sindy, get_feature_matrix_from_data
from SR_analysis.utils.feature_builder import ALL_FEAT_KEYS, U0
from SR_analysis.configs import get_scene, get_scene_list, SCENES
SINDY_DIR = os.path.join(os.path.dirname(__file__), "..", "sindy", "karman")
THRESHOLDS = [0.0, 0.001, 0.002, 0.005, 0.01, 0.015, 0.02, 0.03, 0.05, 0.1]
def load_data(scene_name: str, scene_subdir: str = "karman") -> tuple:
"""Load sensors/forces/actions from a scene's controlled.npz."""
data_dir = os.path.join(os.path.dirname(__file__), "..", "data", scene_subdir, scene_name)
npz = np.load(os.path.join(data_dir, "controlled.npz"))
sensors = npz["sensors"].astype(np.float64)
forces = npz["forces"].astype(np.float64)
actions = npz["actions"].astype(np.float64)
return sensors, forces, actions
def run(scene_names: Optional[List[str]] = None):
"""Run SINDy fitting for given scene names."""
if scene_names is None:
scene_names = get_scene_list("karman")
per_scene = {}
for sn in scene_names:
print(f"\n{'='*60}")
print(f"Scene: {sn}")
print(f"{'='*60}")
cfg = get_scene(sn)
sensors, forces, actions_phys = load_data(sn, cfg["scene_id"])
T = sensors.shape[0]
mu = cfg["mu"]
print(f" T={T}, mu={mu:.6f}")
# Build feature matrices
Theta_f, Theta_r, Y, fn_f, fn_r = get_feature_matrix_from_data(
sensors, forces, actions_phys, mu, u0=cfg["u0"],
alpha_mode=False, include_mu=True, n_warmup=2,
)
print(f" Front features: {Theta_f.shape}")
print(f" Rear features: {Theta_r.shape}")
print(f" Y: {Y.shape}")
# Front channel (ci=0, no bias)
print(f"\n --- Front (no bias) ---")
front_results = fit_sindy(Theta_f, Y[:, 0], THRESHOLDS)
best_f = max(front_results, key=lambda r: r["r2"])
print(f" Best: th={best_f['threshold']:.4f} nz={best_f['nz']:2d} R2={best_f['r2']:.6f}")
# Top channel (ci=2, rear shared-head)
print(f"\n --- Top (rear shared-head) ---")
top_results = fit_sindy(Theta_r, Y[:, 2], THRESHOLDS)
best_t = max(top_results, key=lambda r: r["r2"])
print(f" Best: th={best_t['threshold']:.4f} nz={best_t['nz']:2d} R2={best_t['r2']:.6f}")
# Bottom (independent, for comparison)
print(f"\n --- Bottom (independent) ---")
bot_results = fit_sindy(Theta_r, Y[:, 1], THRESHOLDS)
best_b = max(bot_results, key=lambda r: r["r2"])
print(f" Best: th={best_b['threshold']:.4f} nz={best_b['nz']:2d} R2={best_b['r2']:.6f}")
per_scene[sn] = {
"scene": sn,
"re_code": cfg["re_code"],
"mu": mu,
"n_samples": Theta_f.shape[0],
"n_features_front": Theta_f.shape[1],
"n_features_rear": Theta_r.shape[1],
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"front": {
"results": [{k: v for k, v in r.items() if k != "coef"}
for r in front_results],
"best": {k: v for k, v in best_f.items() if k != "coef"},
"best_coef": best_f["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"])
for r in front_results],
},
"top": {
"results": [{k: v for k, v in r.items() if k != "coef"}
for r in top_results],
"best": {k: v for k, v in best_t.items() if k != "coef"},
"best_coef": best_t["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"])
for r in top_results],
},
"bottom": {
"results": [{k: v for k, v in r.items() if k != "coef"}
for r in bot_results],
"best": {k: v for k, v in best_b.items() if k != "coef"},
"best_coef": best_b["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"])
for r in bot_results],
},
}
# Save
os.makedirs(SINDY_DIR, exist_ok=True)
out_path = os.path.join(SINDY_DIR, "sindy_results.json")
result = {
"thresholds": THRESHOLDS,
"all_feature_names_front": fn_f,
"all_feature_names_rear": fn_r,
"per_scene": per_scene,
}
with open(out_path, "w") as f:
json.dump(result, f, indent=2)
print(f"\nSaved: {out_path}")
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--re-codes", type=str, default=None,
help="Comma-separated Re codes (default: all karman)")
ap.add_argument("--scene-names", type=str, default=None,
help="Comma-separated scene names (overrides --re-codes)")
args = ap.parse_args()
if args.scene_names:
names = [s.strip() for s in args.scene_names.split(",")]
elif args.re_codes:
codes = [int(r) for r in args.re_codes.split(",")]
names = [f"karman_re{rc}" for rc in codes]
else:
names = None # all karman
run(names)
if __name__ == "__main__":
main()
src/SR_analysis/sindy/run_pareto.py
-148
View File
@@ -1,148 +0,0 @@
"""Pareto analysis of SINDy threshold grid for any scene.
Loads sindy_results.json and prints Pareto-optimal tradeoffs.
Usage:
python sindy/run_pareto.py --scene karman_re100
python sindy/run_pareto.py --scene illusion_1L
"""
from __future__ import annotations
import argparse
import json
import os
import sys
from typing import List, Tuple
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)
def load_sindy(scene_name: str, sindy_dir: str) -> dict:
"""Load sindy_results.json for a scene."""
path = os.path.join(sindy_dir, scene_name, "sindy_results.json")
if os.path.isfile(path):
return json.load(path)
raise FileNotFoundError(f"Missing {path}")
def pareto(points: List[Tuple[float, float]]) -> List[Tuple[float, float]]:
"""Compute Pareto frontier: lower nz and lower error is better."""
sp = sorted(points, key=lambda x: (x[0], x[1]))
front, best = [], float("inf")
for c, e in sp:
if e < best:
front.append((c, e))
best = e
return front
def fmt(fn: List[str], coef: List[float], threshold: float) -> str:
"""Format a control law string, showing terms above relative threshold."""
ca = np.array(coef, dtype=np.float64)
sc = np.max(np.abs(ca)) if np.max(np.abs(ca)) > 0 else 1.0
mask = np.abs(ca) / sc >= threshold
terms = [f"{ca[i]:+.4f}*{fn[i]}" for i in range(len(fn)) if mask[i]]
return " ".join(terms) if terms else "0"
def analyze(name: str, feat_names: List[str], channel_data: dict) -> dict:
"""Print and return Pareto analysis for one channel."""
grid = channel_data["results"]
pts = [(g["nz"], 1.0 - g["r2"]) for g in grid]
front = pareto(pts)
best = channel_data["best"]
coef = channel_data["best_coef"]
print(f"\n {name}:")
for nz, err in front:
r2 = 1.0 - err
for g in grid:
if g["nz"] == nz and abs(1.0 - g["r2"] - err) < 1e-10:
th = g["threshold"]
print(f" nz={nz:2d} R2={r2:.6f} th={th:.4f}")
if nz <= 8:
s = fmt(feat_names, coef, th)
print(f" {s[:120]}")
print(f" Best: R2={best['r2']:.6f}")
return {
"channel": name,
"best_r2": best["r2"],
"best_nz": sum(1 for c in coef if abs(float(c)) > 1e-8),
"pareto": [{"nz": nz, "r2": 1.0 - e} for nz, e in front],
}
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--scene", type=str, required=True,
help="Scene name (e.g. karman_re100, illusion_1L)")
ap.add_argument("--sindy-dir", type=str, default=None,
help="SINDy results directory")
ap.add_argument("--out", type=str, default=None,
help="Output JSON path")
args = ap.parse_args()
if args.sindy_dir is None:
args.sindy_dir = os.path.join(os.path.dirname(__file__), "..", "sindy")
# Map scene name to subdirectory
# Extract series prefix: karman_* -> karman, illusion_* -> illusion, etc.
first_part = args.scene.split("_")[0]
known_series = {"karman": "karman", "illusion": "illusion", "vortex": "vortex", "steady": "steady"}
series_dir = known_series.get(first_part, first_part)
# Try sindy/{series}/sindy_results.json
json_path = os.path.join(args.sindy_dir, series_dir, "sindy_results.json")
if not os.path.isfile(json_path):
# Fallback: try flat file
json_path = os.path.join(args.sindy_dir, "sindy_results.json")
if not os.path.isfile(json_path):
print(f"ERROR: No sindy results found for {args.scene}")
print(f" Tried: {os.path.join(args.sindy_dir, series_dir, 'sindy_results.json')}")
print(f" Tried: {json_path}")
return 1
with open(json_path) as f:
data = json.load(f)
# Look up the scene in per_scene (multi-scene format)
per = data.get("per_scene", {}).get(args.scene)
if per is not None:
fn_f = per["feature_names_front"]
fn_r = per["feature_names_rear"]
chs = [("front", fn_f, per["front"]),
("top", fn_r, per["top"]),
("bottom", fn_r, per["bottom"])]
else:
# Single-scene format
fn_f = data.get("feature_names_front")
fn_r = data.get("feature_names_rear")
if fn_f is None:
print(f"ERROR: No scene data found for {args.scene} in {json_path}")
return 1
chs = [("front", fn_f, data["front"]),
("top", fn_r, data["top"]),
("bottom", fn_r, data["bottom"])]
print(f"Pareto SR: {args.scene}")
results = {"scene": args.scene,
"channels": [analyze(*c) for c in chs]}
if args.out:
os.makedirs(os.path.dirname(args.out), exist_ok=True)
with open(args.out, "w") as f:
json.dump(results, f, indent=2)
print(f"Saved: {args.out}")
if __name__ == "__main__":
main()
src/SR_analysis/sindy/run_pysr.py
+102 -58
View File
@@ -1,15 +1,11 @@
#!/usr/bin/env python3
"""Restricted PySR on SINDy whitelist features.
Uses the frozen whitelists (whitelist_karman.json, whitelist_illusion.json)
and the controlled CFD data to search for compact closed-form control laws.
Environment: env_sr (conda).
"""Restricted PySR on phase-state features.
Usage:
conda run -n env_sr python src/SR_analysis/sindy/run_pysr.py --scene karman_re100
conda run -n env_sr python src/SR_analysis/sindy/run_pysr.py --scene illusion_1L
conda run -n env_sr python src/SR_analysis/sindy/run_pysr.py --scene all
conda run -n sr_env python src/SR_analysis/sindy/run_pysr.py --scene illusion_1L --phase
conda run -n sr_env python src/SR_analysis/sindy/run_pysr.py --scene karman_re100 --phase
conda run -n sr_env python src/SR_analysis/sindy/run_pysr.py --scene karman_re100 --phase --include-mu
conda run -n sr_env python src/SR_analysis/sindy/run_pysr.py --scene all --phase
"""
from __future__ import annotations
@@ -17,6 +13,7 @@ import argparse
import json
import os
import sys
import traceback
from typing import Dict, List, Optional
import numpy as np
@@ -31,15 +28,19 @@ sys.path.insert(0, _SRC)
from SR_analysis.utils.feature_builder import (
compute_dimensionless, compute_features, build_feature_matrix,
CORE_FEAT_KEYS_V2,
PHASE_STATE_KEYS, ILLUSION_PHASE_KEYS, MU_FEAT_KEYS,
)
from SR_analysis.configs import get_scene
SINDY_DIR = os.path.join(os.path.dirname(__file__))
# All mu keys produced by compute_features(include_mu=True)
# (MU_FEAT_KEYS has 5, but compute_features generates an additional mu_Cl_tot)
ALL_MU_KEYS = list(MU_FEAT_KEYS) + ["mu_Cl_tot"]
def load_controlled_data(scene_name: str):
"""Load controlled data for a scene, returning sensors, forces, actions_phys, target_forces."""
"""Load controlled data for a scene."""
cfg = get_scene(scene_name)
scene_id = cfg["scene_id"]
data_dir = os.path.join(SINDY_DIR, "..", "data", scene_id, scene_name)
@@ -52,34 +53,47 @@ def load_controlled_data(scene_name: str):
return sensors, forces, actions_phys, target_forces, cfg
def run_pysr_scene(scene_name: str, out_dir: str):
def get_feature_keys(scene_id: str, use_phase: bool, include_mu: bool) -> List[str]:
"""Determine feature keys for a scene."""
if not use_phase:
# Fallback: use CORE_FEAT_KEYS_V2 style (old approach)
from SR_analysis.utils.feature_builder import CORE_FEAT_KEYS_V2
keys = list(CORE_FEAT_KEYS_V2)
elif scene_id == "illusion":
keys = list(ILLUSION_PHASE_KEYS)
else:
# karman, steady, vortex
keys = list(PHASE_STATE_KEYS)
if include_mu:
keys += ALL_MU_KEYS
return keys
def run_pysr_scene(scene_name: str, out_dir: str, use_phase: bool, include_mu: bool):
"""Run PySR on a single scene."""
# Select whitelist based on scene family
cfg = get_scene(scene_name)
scene_id = cfg["scene_id"]
if scene_id == "karman":
whitelist_path = os.path.join(SINDY_DIR, "whitelist_karman.json")
elif scene_id == "illusion":
whitelist_path = os.path.join(SINDY_DIR, "whitelist_illusion.json")
else:
print(f" SKIP: no whitelist for scene_id={scene_id}")
return
with open(whitelist_path) as f:
wl = json.load(f)
sensors, forces, actions_phys, target_forces, cfg = load_controlled_data(scene_name)
mu = cfg["mu"]
u0 = cfg["u0"]
sample_interval = cfg.get("sample_interval", 800)
t0_steps = 2000 # T0 = D/U0 = 20/0.01 = 2000 LBM steps
dt_c = sample_interval / t0_steps
front_keys = wl["front_active"]
rear_keys = [k for k in wl["rear_active"] if k != "bias"]
# Determine feature keys
front_keys = get_feature_keys(scene_id, use_phase, include_mu)
rear_keys = list(front_keys) # rear shares same physics features; bias added by build_feature_matrix
# Build features with PROPER lags (matching sindy_fitter.py)
print(f"\n{'='*60}")
print(f"PySR: {scene_name} (scene_id={scene_id}, phase={use_phase}, mu={include_mu})")
print(f" Front keys ({len(front_keys)}): {front_keys}")
print(f" Rear keys ({len(rear_keys)}): {rear_keys}")
print(f"{'='*60}")
# Load data
sensors, forces, actions_phys, target_forces, cfg = load_controlled_data(scene_name)
mu_val = cfg["mu"]
# Build features with PROPER lags and sensors_raw/forces_raw for derivatives
T = sensors.shape[0]
a_prev = np.zeros((T, 3), dtype=np.float64)
a_prev2 = np.zeros((T, 3), dtype=np.float64)
@@ -87,29 +101,41 @@ def run_pysr_scene(scene_name: str, out_dir: str):
a_prev2[2:] = actions_phys[:-2]
dim = compute_dimensionless(sensors, forces, u0=u0, d=20.0)
has_mu = any("mu" in k for k in front_keys)
sym = compute_features(
dim, a_prev, a_prev2, mu,
alpha_mode=False, include_mu=has_mu,
dim, a_prev, a_prev2, mu_val,
alpha_mode=False, include_mu=include_mu,
include_cos_sin=False, u0=u0,
target_forces=target_forces,
dt_c=dt_c,
sensors_raw=sensors, # CRITICAL: enables du_a_dt, dCl_tot_dt, etc.
forces_raw=forces, # CRITICAL: enables derivative features
)
n_warmup = 2
X_front = build_feature_matrix(sym, front_keys, add_bias=False)[n_warmup:]
X_rear = build_feature_matrix(sym, rear_keys, add_bias=True)[n_warmup:]
Y = actions_phys[n_warmup:] # target = omega(t) for t >= n_warmup
Y = actions_phys[n_warmup:] # physical omega (target)
# Also compute non-dimensional alpha for reference
Y_alpha = Y / u0
print(f"\n{'='*60}")
print(f"PySR: {scene_name}")
print(f" Front features ({len(front_keys)}): {front_keys}")
print(f" Rear features ({len(rear_keys)}): {rear_keys}")
print(f" Samples: {X_front.shape[0]}")
print(f"{'='*60}")
print(f" Y (omega) range: [{Y.min():.6f}, {Y.max():.6f}]")
print(f" alpha range: [{Y_alpha.min():.6f}, {Y_alpha.max():.6f}]")
# --- Front channel ---
print(f"\n=== PySR: {scene_name} Front ===")
# Quick sanity check on derivative features before running PySR
for key in front_keys:
val = sym.get(key, None)
if val is not None:
print(f" {key}: mean={np.mean(val):.6f}, std={np.std(val):.6f}")
else:
print(f" {key}: MISSING from sym dict!")
# Suffix for output filename
suffix = "_mu" if include_mu else ""
# --- Front channel (no bias, predict alpha_F) ---
print(f"\n=== PySR: {scene_name} Front (predict alpha_F) ===")
model_f = PySRRegressor(
binary_operators=["+", "-", "*", "/"],
unary_operators=["square"],
@@ -121,24 +147,26 @@ def run_pysr_scene(scene_name: str, out_dir: str):
extra_sympy_mappings={},
batching=False,
)
model_f.fit(X_front, Y[:, 0], variable_names=front_keys)
model_f.fit(X_front, Y_alpha[:, 0], variable_names=front_keys)
results = {
"scene": scene_name,
"scene_id": scene_id,
"channel": "front",
"output": "alpha",
"feature_names": front_keys,
"equations": model_f.equations_.to_dict(orient="records"),
"best_sympy": str(model_f.sympy()),
"best_score": float(model_f.score(X_front, Y[:, 0])),
"best_score": float(model_f.score(X_front, Y_alpha[:, 0])),
}
out_path = os.path.join(out_dir, f"pysr_{scene_name}_front.json")
out_path = os.path.join(out_dir, f"pysr_{scene_name}_front{suffix}.json")
with open(out_path, "w") as f:
json.dump(results, f, indent=2, default=str)
print(f" Best: {model_f.sympy()}")
print(f" Score: {results['best_score']:.4f}")
print(f" Best sympy: {model_f.sympy()}")
print(f" Score (R2): {results['best_score']:.4f}")
print(f" Saved: {out_path}")
# --- Rear/Top channel ---
print(f"\n=== PySR: {scene_name} Top ===")
# --- Rear/Top channel (with bias, predict alpha_T) ---
print(f"\n=== PySR: {scene_name} Top (predict alpha_T, v23 shared-head) ===")
model_t = PySRRegressor(
binary_operators=["+", "-", "*", "/"],
unary_operators=["square"],
@@ -150,45 +178,61 @@ def run_pysr_scene(scene_name: str, out_dir: str):
extra_sympy_mappings={},
batching=False,
)
model_t.fit(X_rear, Y[:, 2], variable_names=["bias"] + rear_keys)
model_t.fit(X_rear, Y_alpha[:, 2], variable_names=["bias"] + rear_keys)
results_t = {
"scene": scene_name,
"scene_id": scene_id,
"channel": "top",
"output": "alpha",
"feature_names": ["bias"] + rear_keys,
"equations": model_t.equations_.to_dict(orient="records"),
"best_sympy": str(model_t.sympy()),
"best_score": float(model_t.score(X_rear, Y[:, 2])),
"best_score": float(model_t.score(X_rear, Y_alpha[:, 2])),
}
out_path = os.path.join(out_dir, f"pysr_{scene_name}_top.json")
out_path = os.path.join(out_dir, f"pysr_{scene_name}_top{suffix}.json")
with open(out_path, "w") as f:
json.dump(results_t, f, indent=2, default=str)
print(f" Best: {model_t.sympy()}")
print(f" Score: {results_t['best_score']:.4f}")
print(f" Best sympy: {model_t.sympy()}")
print(f" Score (R2): {results_t['best_score']:.4f}")
print(f" Saved: {out_path}")
# Verify front formula is no-bias (action ≈ 0 when features ≈ 0)
zero_vec = np.zeros((1, len(front_keys)), dtype=np.float64)
pred_zero = float(model_f.predict(zero_vec)[0])
print(f"\n Sanity check: front prediction at zero-input = {pred_zero:.6f} (should be ~0)")
def main():
ap = argparse.ArgumentParser(description="Restricted PySR on SINDy whitelists")
ap = argparse.ArgumentParser(
description="PySR on phase-state features (uses phase-state by default)"
)
ap.add_argument("--scene", type=str, default="all",
help='Scene name like "karman_re100" or "illusion_1L", or "all"')
ap.add_argument("--out", type=str, default=None)
ap.add_argument("--phase", action="store_true", default=True,
help="Use phase-state features (default: True)")
ap.add_argument("--no-phase", action="store_false", dest="phase",
help="Use CORE_FEAT_KEYS_V2 instead of phase-state")
ap.add_argument("--include-mu", action="store_true", default=False,
help="Include mu (1/Re) modulation features (Karman cross-Re)")
args = ap.parse_args()
if args.out is None:
args.out = os.path.join(SINDY_DIR, "..", "validate", "results")
os.makedirs(args.out, exist_ok=True)
# Default scenes (skip 1.5L — bang-bang regime)
if args.scene == "all":
scenes = ["karman_re100", "illusion_0.75L", "illusion_1L", "illusion_1.5L"]
scenes = ["illusion_0.75L", "illusion_1L", "karman_re100"]
else:
scenes = [args.scene]
for sn in scenes:
try:
run_pysr_scene(sn, args.out)
run_pysr_scene(sn, args.out, args.phase, args.include_mu)
except Exception as e:
print(f" ERROR on {sn}: {e}")
import traceback; traceback.print_exc()
print(f"\n ERROR on {sn}: {e}")
traceback.print_exc()
print("\nDone.")
src/SR_analysis/sindy/run_vortex.py
-133
View File
@@ -1,133 +0,0 @@
"""SINDy fitting for Vortex scenes.
Usage:
conda run -n pycuda_3_10 python sindy/run_vortex.py
conda run -n pycuda_3_10 python sindy/run_vortex.py --vortex-types lamb,taylor
"""
from __future__ import annotations
import argparse
import json
import os
import sys
from typing import List, Optional
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
if _REPO not in sys.path:
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from SR_analysis.utils.sindy_fitter import fit_sindy, get_feature_matrix_from_data
from SR_analysis.configs import get_scene, get_scene_list
SINDY_DIR = os.path.join(os.path.dirname(__file__), "..", "sindy", "vortex")
THRESHOLDS = [0.0, 0.001, 0.002, 0.005, 0.01, 0.015, 0.02, 0.03, 0.05, 0.1]
def load_data(scene_name: str) -> tuple:
data_dir = os.path.join(os.path.dirname(__file__), "..", "data", "vortex", scene_name)
npz = np.load(os.path.join(data_dir, "controlled.npz"))
sensors = npz["sensors"].astype(np.float64)
forces = npz["forces"].astype(np.float64)
actions = npz["actions"].astype(np.float64)
return sensors, forces, actions
def run(scene_names: Optional[List[str]] = None):
if scene_names is None:
scene_names = get_scene_list("vortex")
per_scene = {}
for sn in scene_names:
print(f"\n{'='*60}")
print(f"Scene: {sn}")
print(f"{'='*60}")
cfg = get_scene(sn)
sensors, forces, actions_phys = load_data(sn)
mu = cfg["mu"]
print(f" T={sensors.shape[0]}, mu={mu:.6f}")
Theta_f, Theta_r, Y, fn_f, fn_r = get_feature_matrix_from_data(
sensors, forces, actions_phys, mu, u0=cfg["u0"],
alpha_mode=False, include_mu=True, n_warmup=2,
)
print(f" Front: {Theta_f.shape}, Rear: {Theta_r.shape}")
# Front channel
print(f"\n --- Front (no bias) ---")
front_results = fit_sindy(Theta_f, Y[:, 0], THRESHOLDS)
best_f = max(front_results, key=lambda r: r["r2"])
print(f" Best: th={best_f['threshold']:.4f} nz={best_f['nz']:2d} R2={best_f['r2']:.6f}")
# Top channel
print(f"\n --- Top (rear shared-head) ---")
top_results = fit_sindy(Theta_r, Y[:, 2], THRESHOLDS)
best_t = max(top_results, key=lambda r: r["r2"])
print(f" Best: th={best_t['threshold']:.4f} nz={best_t['nz']:2d} R2={best_t['r2']:.6f}")
# Bottom
print(f"\n --- Bottom (independent) ---")
bot_results = fit_sindy(Theta_r, Y[:, 1], THRESHOLDS)
best_b = max(bot_results, key=lambda r: r["r2"])
print(f" Best: th={best_b['threshold']:.4f} nz={best_b['nz']:2d} R2={best_b['r2']:.6f}")
per_scene[sn] = {
"scene": sn,
"re_code": cfg["re_code"],
"mu": mu,
"n_samples": Theta_f.shape[0],
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"front": {
"results": [{k: v for k, v in r.items() if k != "coef"} for r in front_results],
"best": {k: v for k, v in best_f.items() if k != "coef"},
"best_coef": best_f["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in front_results],
},
"top": {
"results": [{k: v for k, v in r.items() if k != "coef"} for r in top_results],
"best": {k: v for k, v in best_t.items() if k != "coef"},
"best_coef": best_t["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in top_results],
},
"bottom": {
"results": [{k: v for k, v in r.items() if k != "coef"} for r in bot_results],
"best": {k: v for k, v in best_b.items() if k != "coef"},
"best_coef": best_b["coef"],
"sparsity_curve": [(r["threshold"], r["nz"], r["r2"]) for r in bot_results],
},
}
os.makedirs(SINDY_DIR, exist_ok=True)
out_path = os.path.join(SINDY_DIR, "sindy_results.json")
result = {"thresholds": THRESHOLDS, "per_scene": per_scene}
with open(out_path, "w") as f:
json.dump(result, f, indent=2)
print(f"\nSaved: {out_path}")
def main():
ap = argparse.ArgumentParser()
ap.add_argument("--vortex-types", type=str, default=None,
help="Comma-separated vortex types (e.g. lamb,taylor)")
ap.add_argument("--scene-names", type=str, default=None)
args = ap.parse_args()
if args.scene_names:
names = [s.strip() for s in args.scene_names.split(",")]
elif args.vortex_types:
names = [f"vortex_{v.strip()}" for v in args.vortex_types.split(",")]
else:
names = None
run(names)
if __name__ == "__main__":
main()
src/SR_analysis/sindy/vortex/pareto_vortex_lamb.json
-110
View File
@@ -1,110 +0,0 @@
{
"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
}
]
}
]
}
src/SR_analysis/sindy/vortex/pareto_vortex_taylor.json
-58
View File
@@ -1,58 +0,0 @@
{
"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
}
]
}
]
}
src/SR_analysis/sindy/vortex/sindy_results.json
-996
View File
@@ -1,996 +0,0 @@
{
"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
]
]
}
}
}
}
src/SR_analysis/sindy/whitelist_illusion.json
-30
View File
@@ -1,30 +0,0 @@
{
"scene": "illusion_joint",
"threshold": 0.001,
"selection": "Cleaned features: removed aF_lag1/daF etc, replaced with daF_dt (time-normalized). Features from union across 0.75L/1L/1.5L dense fits with target_Cd/target_Cl.",
"front_active": [
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"daF_dt",
"daB_dt",
"daT_dt",
"target_Cd",
"target_Cl"
],
"rear_active": [
"bias",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"daF_dt",
"daB_dt",
"daT_dt",
"target_Cd",
"target_Cl"
],
"feature_count": 16,
"note": "Cleaned Illusion whitelist: no discrete lags, only time-derivatives + physical features. Separate SINDy closed-loop sims: 0.75L=0.908, 1L=0.962, 1.5L=0.926."
}
src/SR_analysis/sindy/whitelist_karman.json
-18
View File
@@ -1,18 +0,0 @@
{
"scene": "karman_re100",
"threshold": 0.003,
"front_active": [
"daF_dt"
],
"rear_active": [
"bias",
"u_a",
"Cd_rear",
"Cl_diff",
"daF_dt",
"daB_dt",
"daT_dt"
],
"feature_count": 14,
"note": "Cleaned Karman whitelist: removed aF_lag1/daF, replaced with daF_dt (time-normalized). Based on th=0.003 support from joint fit. Separate SINDy closed-loop sim=0.901 (94.4% of PPO)."
}
src/SR_analysis/sindy/wrap_joint.py
-66
View File
@@ -1,66 +0,0 @@
#!/usr/bin/env python3
"""Wrap joint model into per_scene format.
Usage: python3 src/SR_analysis/sindy/wrap_joint.py [--scene illusion]
"""
import json, sys, os
sys.path.insert(0, os.path.abspath("."))
sys.path.insert(0, os.path.abspath("src"))
import argparse
from SR_analysis.configs import get_scene_list
def wrap_joint(scene_id):
joint_path = f"src/SR_analysis/sindy/{scene_id}/sindy_joint_v2.json"
out_path = f"src/SR_analysis/sindy/{scene_id}/sindy_joint_wrapped.json"
with open(joint_path) as f:
joint = json.load(f)
fn_f = joint["feature_names_front"]
fn_r = joint["feature_names_rear"]
wrapped = {
"thresholds": joint["thresholds"],
"feature_names_front": fn_f,
"feature_names_rear": fn_r,
"per_scene": {},
}
# Include all scenes in the family (training + generalization)
all_scene_names = get_scene_list(scene_id)
for scene_name in all_scene_names:
wrapped["per_scene"][scene_name] = {
"scene": scene_name, "re_code": "0", "mu": 0.0,
"feature_names_front": fn_f, "feature_names_rear": fn_r,
"front": {
"results": joint["front"]["results"],
"best": joint["front"]["best"],
"best_coef": joint["front"]["best_coef"],
"sparsity_curve": joint["front"]["sparsity_curve"],
},
"top": {
"results": joint["top"]["results"],
"best": joint["top"]["best"],
"best_coef": joint["top"]["best_coef"],
"sparsity_curve": joint["top"]["sparsity_curve"],
},
"bottom": {
"results": joint["bottom"]["results"],
"best": joint["bottom"]["best"],
"best_coef": joint["bottom"]["best_coef"],
"sparsity_curve": joint["bottom"]["sparsity_curve"],
},
}
with open(out_path, "w") as f:
json.dump(wrapped, f, indent=2)
print(f"Saved {out_path}")
print(f" Scenes ({len(all_scene_names)}): {all_scene_names}")
if __name__ == "__main__":
ap = argparse.ArgumentParser()
ap.add_argument("--scene", type=str, default="karman",
choices=["karman", "illusion"])
args = ap.parse_args()
wrap_joint(args.scene)
src/SR_analysis/sindy_sr_knoeledge.md
-339
View File
@@ -1,339 +0,0 @@
# SINDy 与 SR 背景知识
## 文档作用
这份文档只负责一件事：**给正在工作的 coder 提供背景知识、已经确认的经验、已踩过的坑和当前结论强度。**
它不是任务清单，不直接安排“下一步做什么”。凡是执行顺序、阶段划分、最小交付物，统一写在 `sindy_sr_notes`。这份 knowledge 只保留：
- 已确认的技术事实
- 历史错误与纠正
- 结果该如何理解
- 哪些话可以说，哪些话现在还不能说
- 代码和实验上最容易踩的坑
---
## 一、这条线在项目里的位置
SINDy 与 SR 不是独立课题，而是 pinball 后处理主线中的一段工具链。项目真正要解释的是：
\[
\text{obs} \rightarrow \text{act} \rightarrow \text{flow structure} \rightarrow \text{signature}
\]
SINDy 与 SR 当前只直接处理其中的 `obs -> act` 白箱化，但它们的价值在于：
- 检验控制是否真的依赖少数物理量
- 识别不同 cloak 场景中是否复用了同一类反馈结构
- 为后续把控制律与 force、mean wake、observable-related structure 接起来提供接口
因此，任何 SINDy/SR 结果都不应脱离项目总体物理主线单独解读。
---
## 二、当前已经确认的技术事实
### 1. Kármán cloak 的跨 \(Re_D\) 统一骨架存在
这是目前最硬的一批证据之一。跨 Re 的 leave-one-Re-out，尤其 holdout_200，已经显示：
- 用 Re50 + Re100 拟合，可高精度预测 Re200
- 这说明统一骨架不是偶然的特征工程产物，而是真实存在于 PPO 策略中的共享结构
### 2. 对称性问题已经纠偏
最重要的 bug 是镜像变换 \(G\) 对动作的写法错误。
**错误版本**：
\[
[a_F,a_T,a_B] \mapsto [-a_F,a_B,a_T]
\]
**正确版本**：
\[
[a_F,a_T,a_B] \mapsto [-a_F,-a_B,-a_T]
\]
也就是说：
- top / bottom 不仅交换
- 三个动作都要变号
修正后，rear equivariance 误差从约 100% 降到约 10%，原来“PPO 不尊重交换对称性”的结论应正式撤回。
当前应保留的结构关系是：
\[
\alpha_F(Gx) \approx -\alpha_F(x)
\]
\[
\alpha_B(x) \approx -\alpha_T(Gx)
\]
### 3. front no-bias 被数据支持
front 通道不需要常偏置。去掉 bias 后：
- one-step 基本不变
- 关键闭环也基本不变
因此，front odd structure 现在可以作为默认先验，而不是可选修饰。
### 4. rear shared-head 不是纯粹美学约束，而是有效结构
`bottom(x) = -top(Gx)` 的结构不是只让模型更优雅，它在闭环里确实提供了稳健性。v23 的结果说明：
- 结构约束有助于防止 rear 两通道在闭环中各走各路
- 它比无结构的独立 rear 拟合更适合作为解释模型
### 5. 无量纲化不是问题根源
这件事已经确认，不应再反复争论。
- \(u \to u/U_0\)
- \(F \to C_D, C_L\)
- \(\Omega \to \alpha\)
这些都是可逆缩放，不会丢失信息。早期 v3 崩坏来自：
- 错误 \(G\)
- 多项改动一次性叠加
而不是无量纲化本身。
### 6. one-step 与闭环是两回事
这是这条线最重要的工作方法教训之一：
\[
\text{one-step R² 高} \not\Rightarrow \text{闭环好}
\]
早期 v3(old) 就是明确反例。因此：
- one-step 只能说明局部拟合能力
- 闭环验证是核心证据，不是附加项
### 7. PySR 现在可用
之前关于“PySR 不可用”的说法应删除。当前已知：
- `sr_env` 下 PySR 可用
- 后续 SR 主工具应优先考虑 PySR
- threshold Pareto 扫描仍有价值，但不能再混称为完整 SR
---
## 三、哪些结论现在还不能说得太满
### 1. “所有 cloak 已经共享同一骨架”
还不能这么说。当前最强证据只够支持：
- Kármán cloak 跨 \(Re_D\) 统一骨架成立
- steady 初步显示出明显简化版结构
但 all-cloak 统一骨架仍是当前主问题，不是已证结论。
### 2. “steady cloak 已经严格证明是 Kármán 的子模型”
这个说法也太满。更稳的表述是：
- 在当前 steady 数据定义下，steady 的 support 呈现出 Kármán support 的明显简化版
- 这支持 `shared backbone + scene-specific activation` 方向
但 steady 当前的证据强度仍弱于 Kármán，因为 steady 不是同类 DRL 闭环策略数据。
### 3. “高 Re 退化已经证明是采样率问题”
现在还不能这么写。更稳的说法是：
- 这是一个强工作假设
- 需要在时间尺度显式化后重新检验
### 4. “SR 已经做完一轮”
如果实际做的只是 threshold 网格 + Pareto 分析，就不能写成“完整 SR 已完成”。
要区分：
- `threshold Pareto scan`
- `true constrained SR`
---
## 四、统一变量与约束的背景知识
### 1. primitive variables 应统一
后续所有 cloak 场景都应基于同一批 primitive variables：
- \(\hat u, \hat v\)
- \(C_D, C_L\)
- \(\alpha\)
- lagged \(\alpha\)
- \(\Delta\alpha\)
- \(\mu = 1/Re_D\)
- scene metadata 与 \(\Delta t_c\)
### 2. \(G\) 算子必须在 primitive level 定义
不要再直接对压缩特征猜符号。统一规则为：
| 量 | 变换 |
|---|---|
| \((u_B,u_C,u_T)\) | \((u_T,u_C,u_B)\) |
| \((v_B,v_C,v_T)\) | \((-v_T,-v_C,-v_B)\) |
| \((C_{D,F},C_{D,T},C_{D,B})\) | \((C_{D,F},C_{D,B},C_{D,T})\) |
| \((C_{L,F},C_{L,T},C_{L,B})\) | \((-C_{L,F},-C_{L,B},-C_{L,T})\) |
| \((\alpha_F,\alpha_T,\alpha_B)\) | \((-\alpha_F,-\alpha_B,-\alpha_T)\) |
| lag / increment | 同动作规则 |
| \(\mu\) | 不变 |
并且：
\[
G(G(x)) = x
\]
必须作为基本测试。
### 3. 默认结构约束
当前最稳的默认结构是：
- front no-bias
- front odd structure
- rear shared-head
即：
\[
\alpha_T(x)=g_R(x),
\qquad
\alpha_B(x)=-g_R(Gx)
\]
不再把三通道完全独立当默认。
---
## 五、SINDy 与 SR 的正确分工
### SINDy 负责什么
SINDy 的主要价值是：
- 在受限物理库上识别主项
- 给出 support 证据
- 支持跨场景比较
- 给 SR 提供 whitelist
SINDy 是骨架识别器，不是最终公式生成器。
### SR 负责什么
SR 的价值不止是压短公式。它还可以：
- 吸收若干看似分散的 SINDy 项
- 暴露不同场景是否存在同形公式
- 给出比 threshold scan 更强的闭式线索
因此 SR 应该在受限物理库上做，而不是在 raw feature 上自由乱搜。
---
## 六、时间尺度问题的背景知识
### 1. 当前公式混有采样周期信息
lagged action 和 \(\Delta a\) 都隐式绑定了 control interval。也就是说，当前公式里混着：
- 物理骨架
- 离散实现方式
- 采样周期
因此时间尺度问题是结构问题，不是附带工程问题。
### 2. 控制频率测试的严格原则
不能把在 800-step cadence 下拟合出的系数，直接拿到 400-step 或 200-step cadence 下执行，并把结果当正式证据。因为此时：
- 输入分布变了
- memory 项的物理意义变了
- \(\Delta a\) 的尺度也变了
因此，如果要比较不同 control interval，必须：
1. 在目标 \(\Delta t_c\) 下重新采集特征
2. 重新拟合模型
3. 再比较 support、系数结构与闭环
之前的频率扫描结果最多只能当线索，不能当结论。
---
## 七、steady 结果应该怎样理解
当前 steady 的结果有启发性，但需要克制解释。
可以说：
- steady front 全零很合理
- steady rear 比 Kármán 更简单
- steady 当前 support 呈现出 Kármán 的明显简化版
不宜说：
- steady 已经严格证明是 Kármán 的子模型
- steady 与 Kármán 现在证据强度相同
因为 steady 当前的数据来源与 Kármán 不完全对等。
---
## 八、代码与工程层面的已知经验
### 1. 环境分工
- `pycuda_3_10`：CFD、DRL 模型加载、数据采集、SINDy
- `sr_env`：PySR 与 SR 相关工作
### 2. 常见坑
- feature names 与矩阵列顺序不一致
- JSON 保存时未统一处理 numpy 类型
- 不同脚本用不同 channel 命名规则
- 没有统一 validator，导致 \(G\) 与闭环输入错位难以及早发现
### 3. 推荐工程习惯
- 任何场景都先过 validator，再进拟合
- separate fit 的结果按场景 × 方法存储
- support、公式、闭环三类结果必须一起保存
---
## 九、当前最值得牢记的判断
这条线现在最稳的总结是：
\[
\boxed{
\text{Kármán cloak 的跨 }Re_D\text{ 统一骨架已确认，且满足明确的镜像等变结构；v23 是当前最可信的解释模型。}
}
\]
同时必须保留另一句：
\[
\boxed{
\text{all-cloak 的最终 shared backbone 还没有定论；steady、单涡、时间尺度显式化与真正的受限 SR 仍在探索中。}
}
\]
这两句话一起保留，能避免后续工作再次滑向“把局部结果写成全局结论”。
src/SR_analysis/sindy_sr_knowledge.md
+102 -41
View File
@@ -32,15 +32,29 @@ SINDy 与 SR 当前只直接处理其中的 `obs -> act` 白箱化，但它们
---
## 二、当前已经确认的技术事实
## 二、当前已经确认的技术事实（2026-06-23 更新）
### 1. Kármán cloak 的跨 Re_D 统一骨架存在
### 1. Karman 跨 Re 联合公式成立（2026-06-23 新增）
这是目前最硬的一批证据之一。跨 Re 的 leave-one-Re-out，尤其 holdout_200，已经显示：
- 用 Re50 + Re100 拟合，可高精度预测 Re200
- 这说明统一骨架不是偶然的特征工程产物，而是真实存在于 PPO 策略中的共享结构
\[
\alpha_F = \Delta a_F / \Delta t - 14.952 \cdot \mu C_{l,\text{tot}}
\]
### 2. 对称性问题已经纠偏
四个雷诺数（50/100/200/400）CFD 闭环验证全部通过，平均 0.847，远超旧基线 0.735。**这是目前最硬的统一公式证据。** 关键设计：`daF_dt` 用 `dt_c` 归一化消除采样率依赖，`mu * Cl_tot` 项自适应粘度变化。
### 2. 独立 Re PySR 公式形态各异但都成功（2026-06-23 新增）
四个雷诺数的独立公式闭环均 >0.89，但公式形态完全不同：从 Re=50 的阻力阻尼主导，逐步过渡到 Re=400 的升力主导。这反映 pinball 流态在 Re=50~400 区间经历了分岔，PPO 学会了针对性策略。
### 3. 训练分布偏移的假象项（2026-06-23 新增）
PySR 联合搜索初期返回的公式包含 `daB_dt` 项（`alpha_F = daF_dt*0.84 + ...`），但这是 PPO 轨迹中前后端动作共线导致的**虚假相关**。部署时后端恒速旋转，`daB_dt = 0`，该项无贡献。**最终公式必须手动审查，去除此类假象项。**
### 4. Kármán cloak 的跨 Re_D 统一骨架存在
这是旧的事实，与上述联合公式结论一致。
### 5. 对称性问题已经纠偏
最重要的 bug 是镜像变换 G 对动作的写法错误。
@@ -56,9 +70,11 @@ SINDy 与 SR 当前只直接处理其中的 `obs -> act` 白箱化，但它们
修正后，rear equivariance 误差从约 100% 降到约 10%，原来"PPO 不尊重交换对称性"的结论应正式撤回。
### 3. front no-bias 被数据支持，rear shared-head 是有效结构
### 6. front no-bias 被数据支持，rear shared-head 是有效结构
### 4. one-step R² 与闭环是两回事
所有独立和联合公式均验证：前端公式无偏置项（零输入时输出接近零）。
### 7. one-step R² 与闭环是两回事
\[
\text{one-step R}^2 \gg 0.95 \not\Rightarrow \text{闭环好}
@@ -66,81 +82,126 @@ SINDy 与 SR 当前只直接处理其中的 `obs -> act` 白箱化，但它们
**关键证据**：full-lag 16-dim 模型 R²=0.939 但闭环仅 0.619；static 8-dim 模型 R²=0.321 但闭环 0.745。
### 5. 时间特征加剧 rollout mismatch（2026-06-15 确认）
### 8. Illusion 0.75L 和 1L 的 PySR 公式已跑通（2026-06-23 新增）
训练时使用 `x(t-1)` 或 `dx/dt` 特征让 one-step R² 跃升，但闭环验证下降。原因是**分布偏移**：训练时特征来自 PPO 真实轨迹，部署时来自 SINDy 控制轨迹。带时间记忆的模型更容易在自由滚动时暴露分布偏移。
| 场景 | Front 公式 | 闭环 | % 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% |
### 6. phase-state + absolute action 路线已被验证（2026-06-15 核心结论）
两公式结构不同（0.75L 用压升力、1L 用相位），说明**目标尺寸差异导致控制机制不同**。
Illusion 0.75L 和 1L 已证明**不需要动作历史特征**也能达到 PPO 的 96%+ 闭环性能：
- 输入：相位状态 `u_a, du_a/dt, Cl_tot, dCl_tot/dt, Cd_tot, Cd_rear` + error-state
- 输出：直接预测 absolute alpha，不做导数积分
- 结构：v23（front no-bias, rear shared-head）
### 9. Illusion 1.5L 是独立机制（2026-06-25 修正）
### 7. Illusion 1.5L 是独立机制（2026-06-15 确认）
1.5L 之前被描述为"bang-bang（动作饱和）"，但实际分析后发现并非如此：
- 动作使用完整 [-1, 1] 范围但仅在极端值处饱和约 1% 的时间
- 控制是**高频率、高振幅的周期调制**，频率 f=0.24（控制空间），是目标涡脱频率的 5.6 倍
- 所有 10 维 ILLUSION_PHASE 特征与动作之间的线性相关性均 < 0.33
- 动作自相关 lag-2 ≈ -0.9，说明每两步切换一次方向
- 标准 PySR/SINDy 均不适用（已确认 PySR 给出 `aF_lag1 * 0.01` 的平凡解）
1.5L 动作饱和到 [-8,8] 范围，自相关 r=0.07，线性 SINDy 不适用。
### 10. Illusion 联合公式已验证（2026-06-25 新增）
联合公式（0.75L + 1L 数据拼接，无 target_diameter 标记）：
- Front: `Cd_tot - (Cd_err + 5.428) - 0.00978*(du_a_dt + u_a)`，R2=0.907
- Top: `(Cd_err - (Cd_rear - Cl_err))*0.535 + 2.782`，R2=0.828
- CFD 闭环：**0.75L=0.978**，**1L=0.970**
- 联合公式能在两场景上都取得接近最优的效果，说明控制器通用性强
方法 A（target_diameter 标记法）确认：target_diameter **被最优公式使用**，说明两场景公式结构确实不同。
### 11. Karman 联合公式泛化到 Vortex 场景（2026-06-25 新增）
Karman 跨 Re 联合 PySR 公式直接应用于 vortex cloak 场景（未重新训练）：
- vortex_lamb（Lamb 偶极子）：CFD 闭环 **0.949**（超越 PPO 基线 0.942）
- vortex_taylor（Taylor 单极子）：CFD 闭环 **0.905**（接近 PPO 基线 0.916）
- 说明联合公式提取的物理反馈结构具有很强的通用性
### 12. Karman re400 最优控制间隔（2026-06-25 新增）
re400 的 PPO 效果较差（相似度 0.795），怀疑是采样间隔 SI=800 过长。
测试不同 SI 下联合公式的 CFD 闭环表现：
- SI=800: 0.806（基线）
- **SI=400: 0.819（最优）**
- SI=200: 0.794
结论：re400 最优控制间隔为 SI=400，但仍受限于 PPO 本身在高 Re 下的控制质量。
---
## 三、哪些结论现在还不能说得太满
### 1. "所有 cloak 已经共享同一骨架"
还不能这么说。all-cloak 统一骨架仍是当前主问题。
不能。Karman 和 Illusion 的公式结构不同，不能统一。
### 2. "Karman 新路线已经成功"
不能。Karman 无动作历史最佳仅 0.699（旧 v23 为 0.901），说明状态还不充分。
### 2. "Karman 跨 Re 联合公式已完全解决"
不能。后端为常数（`alpha_T = 3.414`），说明后端控制信息未被充分利用。且 re400 独立公式未做闭环验证。
### 3. "高 Re 退化已经证明是采样率问题"
这是一个强工作假设，需要在时间尺度显式化后重新检验。
### 4. "SR 已经做完一轮"
如果实际做的只是 threshold 网格 + Pareto 分析，不能写成"完整 SR 已完成"。
Karman 完整，Illusion 联合公式已验证 + 结构比较已完成，1.5L 已分析为不可 SR 的独立机制。
---
## 四、代码与工程层面的已知经验
### 1. 环境分工
- `pycuda_3_10`：CFD、DRL 模型加载、数据采集、SINDy（pysindy）
- `sr_env`：PySR 与 SR 相关工作
- `pycuda_3_10`：CFD、DRL 模型加载、数据采集、SINDy（pysindy）、CFD 闭环验证
- `sr_env`：PySR 符号回归
### 2. 常见坑
- **actions.npz 是归一化动作 [-1,1]，不是物理 omega**。所有 SINDy 拟合代码需通过 `(norm * scale + bias) * u0` 转换。
- **actions.npz 是归一化动作 [-1,1]，不是物理 omega**。所有拟合代码需通过 `(norm * scale + bias) * u0` 转换。
- **Illusion "2U" 误解**："2U" 在模型名中表示 S_DIM=14（多了2维目标力观测），不是两倍来流速度。u0 始终为 0.01。
- **FIFO bias ≠ DRL action bias**：FIFO bias 用 1U `[0, -U0, U0]` 填充历史；DRL action decoder 用 2U `action*8+[0,-2,2]`。
- **SAMPLE_INTERVAL 因场景而异**：0.75L=400, 1L=600, 1.5L=800, Karman=800。
- **验证器 `n_steps` 不能过短**：需保证控制传播到传感器，S=400→320步, S=600→214步, S=800→160步。
- **保存 `controlled.npz` 时必须包含 `target_forces` 字段**（Illusion 场景的 s_dim=14 推理数据）。
- **单步 validator（predict_v23_deriv）中需要传入 `sensors_raw`/`forces_raw`** 以计算相位特征中的导数项。
- **单步 validator 中需要传入 `sensors_raw`/`forces_raw`** 以计算相位特征中的导数项。
- **PySR 搜索时必须传 `sensors_raw`/`forces_raw`** 给 `compute_features()`，否则所有 `du_a_dt`, `dCl_tot_dt` 等导数特征均为零。这是 run_pysr.py 的 CRITICAL bug（已修复）。
- **输出目标必须是 alpha 而非 omega**：`Y = actions_phys / u0`，否则 coefficients 受 u0 缩放。
- **联合公式中的 `daB_dt` 是训练分布假象**：PPO 轨迹中前后端动作共线导致虚假相关，部署时后端恒速该项为零，应手动移除。
### 3. 当前可信结果（2026-06-15）
### 3. 当前可信结果（2026-06-25 更新）
| 场景 | 路线 | 特征 | 输出 | 是否有动作历史 | 闭环 sim | % of PPO |
|------|------|------|:----:|:-------------:|:--------:|:--------:|
| illusion_0.75L | phase+error | ILLUSION_PHASE (10dim) | alpha | **否** | **0.974** | 100.2% |
| illusion_1L | phase+error | ILLUSION_PHASE (10dim) | alpha | **否** | **0.958** | 98.5% |
| illusion_1.5L | phase+error | ILLUSION_PHASE (10dim) | alpha | **否** | **N/A** | bang-bang |
| karman_re100 | phase-state | PHASE_STATE (6dim) | alpha | **否** | **0.699** | 73.3% |
| karman_re100 | old v23 | CORE_FEAT_KEYS_V2 + a_lag | alpha | **是** | **0.901** | 94.4% |
| 场景 | 路线 | 特征 | 输出 | 动作历史? | 闭环 sim | % of PPO |
|------|------|------|:----:|:---------:|:--------:|:--------:|
| illusion_0.75L | PySR | ILLUSION_PHASE (10dim) | alpha | 否 | **0.979** | **100.7%** |
| illusion_1L | PySR | ILLUSION_PHASE (10dim) | alpha | 否 | **0.957** | 98.4% |
| **illusion joint** | **PySR** | **ILLUSION_PHASE (10dim)** | **alpha** | **否** | **0.978/0.970** | — |
| illusion_1.5L | — | — | — | — | N/A | 高频周期调制 |
| **karman 联合** | **PySR** | **PHYS_DADT+mu (17dim)** | **alpha** | **daF_dt** | **re50=0.847, re100=0.888, re200=0.845, re400=0.806** | **统一公式** |
| karman_re50 | PySR | PHYS_DADT+mu | alpha | daF_dt | **0.895** | **153.8%** |
| karman_re100 | PySR | PHYS_DADT+mu | alpha | daF_dt | **0.888** | 98.6% |
| karman_re200 | PySR | PHYS_DADT+mu | alpha | daF_dt | **0.916** | **115.5%** |
| karman_re400 (SI=400) | PySR | PHYS_DADT+mu | alpha | daF_dt | **0.819** | 最优SI |
| **vortex_lamb** | **Karman joint** | — | alpha | daF_dt | **0.949** | **超越PPO(0.942)** |
| **vortex_taylor** | **Karman joint** | — | alpha | daF_dt | **0.905** | 接近PPO(0.916) |
---
## 五、Bug Audit 补充（2026-06-14~15）
## 五、Bug Audit 全部记录
### 新发现的 Bug
### 2026-06-14~15 发现的 Bug
| # | Bug | 文件 | 严重度 | 修复 |
|---|-----|------|--------|------|
| 13 | `run_pysr.py` 滞后构造错误：`actions_phys` 直接当 `actions_prev` 传入 `compute_features`，导致 `aF_lag1 = alpha(t)` 而非 `alpha(t-1)` — 拟合恒等式 | sindy/run_pysr.py | CRITICAL | 改为正确滞后：`a_prev[1:]=actions_phys[:-1]` |
| 14 | `predict_v23_deriv` 中 `needs_aug` 只检测 `"lag1"` 关键字，忽略 `"_dt"` 结尾的导数特征 → 闭环中 phase-state 的 `du_a_dt` 等始终为零 | validate/run_closed_loop.py | CRITICAL | `needs_aug = any(k.endswith("_dt") or k.endswith("_lag1") for k in ...)` |
| 15 | `compute_features` 中 `target_forces` 未处理 1D 输入（单步验证时 shape 为 (2,) 而非 (1,2)） | utils/feature_builder.py | 中 | `if tf.ndim == 1: tf = tf.reshape(1, -1)` |
| 13 | `run_pysr.py` 滞后构造错误：`actions_phys` 直接当 `actions_prev` 传入 `compute_features`，导致 `aF_lag1 = alpha(t)` 而非 `alpha(t-1)` — 拟合恒等式 | sindy/run_pysr.py | CRITICAL | `a_prev[1:]=actions_phys[:-1]` |
| 14 | `predict_v23_deriv` 中 `needs_aug` 只检测 `"lag1"` 关键字，忽略 `"_dt"` 结尾的导数特征 → 闭环中 phase-state 的 `du_a_dt` 等始终为零 | validate/run_closed_loop.py | CRITICAL | `needs_aug = any(k.endswith("_dt") or k.endswith("_lag1") ...)` |
| 15 | `compute_features` 中 `target_forces` 未处理 1D 输入 | utils/feature_builder.py | 中 | `if tf.ndim == 1: tf = tf.reshape(1, -1)` |
### 关键方法教训（2026-06-15 更新）
### 2026-06-23 修复的关键 Bug
| # | Bug | 文件 | 严重度 | 修复 |
|---|-----|------|--------|------|
| **16** | **`run_pysr.py` 调用 `compute_features()` 时未传 `sensors_raw`/`forces_raw`**，导致所有 `du_a_dt`, `dCl_tot_dt`, `dCd_err_dt`, `dCl_err_dt` 等相位导数特征在 `sym` 字典中缺失 → `build_feature_matrix()` 用全零填充 → PySR 用零向量拟合真实动作 → 公式必然古怪 | **sindy/run_pysr.py** | **CRITICAL** | 在 `compute_features()` 调用中增加 `sensors_raw=sensors, forces_raw=forces` |
### 关键方法教训
1. **PPO 验证优先于 SINDy**。在确认 PPO 推理复现正确之前，SINDy 结果无意义。
2. **控制时长必须 >= NX/U0**。50 步验证结果不可靠。
3. **DDF 保存位置决定验证器起始状态**。正确模式：`stabilize → save_ddf(1) → norm → apply_ddf → bias FIFO → save_ddf(2) → apply_ddf`。
4. **状态消融比输出消融更重要**：`x_n` 不够 → 加 `x_{n-1}` 大幅提升 → 加 `a_{n-1}` 边际仅 3% → 说明缺的是相位状态，不是动作历史。
5. **离线 rollout 好 ≠ CFD 闭环好**。phase-state 模型 offline 50 步漂移仅 11%，但 CFD 闭环 -7% vs 静态。最终判据是 CFD 短闭环。
6. **如遇数值错误，先检查 obs 切片索引、添加顺序、target_forces 维度**。
5. **离线 rollout 好 ≠ CFD 闭环好**。最终判据是 CFD 短闭环。
6. **联合公式最终必须人工审查**：移除训练分布中的假象项（如 `daB_dt` 在后端恒速时恒为零）。
7. **PySR 输出的最佳公式不一定合适**：因为 SR 只优化 loss，不考虑部署分布偏移。检查 Pareto 前沿，理解每个项在闭环中是否真正有贡献。
src/SR_analysis/sindy_sr_notes.md
+107 -157
View File
@@ -1,188 +1,138 @@
# SINDy 与 SR 执行计划（2026-06-15 修订版 v4）
# SINDy-SR Task List & Execution Notes
## 文档作用
这份文档只回答一件事：**接下来要做什么。**
- 只写执行路线、阶段目标、优先级、输出要求
- 不长篇复述历史争论
- 凡是背景知识、已知结论、踩坑记录，统一放到 `sindy_sr_knowledge.md`
> This file is the **active task list** for SINDy / PySR symbolic regression work.
> It is referenced by `sindy_sr_knowledge.md` as the source of truth for execution
> order, phase boundaries, and current task status.
>
> Background knowledge, confirmed facts, bug history, and result interpretation
> belong in `sindy_sr_knowledge.md`. This file only tracks **what to do next**.
---
## 当前总目标
## 1. Current Task Status (2026-06-28)
用 SINDy + SR 在所有 cloak / illusion 场景上产出有物理结论的控制律公式，并验证其跨场景泛化能力。
### Completed
回答三个问题：
1. **哪些物理量是控制的核心变量？**（力、速度、相位、误差之间的取舍）
2. **不同场景是否共享同一公式骨架？**（family 内 -> family 间）
3. **白箱控制律能否在未训练工况工作？**（泛化性 + 失效边界）
- [x] Karman cross-Re joint formula: `α_F = daF_dt − 14.952·μ·Cl_tot`, `α_T = 3.414` (constant)
- [x] Karman deep individual per-Re formulas (niter=120)
- [x] Karman re400 short-SI test (SI=400 optimal at 0.819)
- [x] Illusion 0.75L individual PySR: `α_F = −0.169·(Cl_tot + dCl_tot/dt) − 1.240`
- [x] Illusion 1L individual PySR: `α_F = (du_a/dt + u_a + 26.5)·0.0123`
- [x] Illusion joint formula (Method B): `α_F = Cd_tot − (Cd_err + 5.428) − 0.00978·(du_a/dt + u_a)`, CFD 0.75L=0.978, 1L=0.970
- [x] Illusion generalization CFD (0.5L−2.0L) via joint formula
- [x] Vortex generalization (Karman joint formula): lamb=0.949, taylor=0.905
- [x] Illusion 1.5L characterized as high-frequency periodic modulation (not SR-amenable)
- [x] G-equivariance bug fixed
- [x] Bug #13/#14/#15/#16 all fixed
- [x] `docs/SR_analysis_report.md` written (465 lines)
- [x] `docs/illusion_joint_formula_analysis.md` written
### In Progress (Phase P0 — 2026-06-28)
- [ ] P0.1: Triage 136 JSON result files → classify, archive intermediates
- [ ] P0.2: Create `validate/results/README.md` (reference table for canonical files)
- [ ] P0.3: Update `README.md` directory tree and Key Documentation table
- [ ] P0.4: Update `SR_analysis_report.md` §8 with generalization rows
- [ ] P0.5: Fix doc error in `illusion_joint_formula_analysis.md` line 46 (1.5L value)
### Pending (Phase P2)
- [ ] P2.1: Regenerate PPO rollouts for 0.75L, 1L, 1.5L
- [ ] P2.2: Run SR closed-loop CFD for 1.2L and 2L with full telemetry
- [ ] P2.3: Create `scripts/analyze_illusion_degradation.py`
- [ ] P2.4: Document degradation analysis in `illusion_joint_formula_analysis.md` §2.3
### Pending (Phase P3)
- [ ] P3.1: Run 1L with-target CFD validation
- [ ] P3.2: Document results
---
## 已完成的工作（2026-06-14~15）
## 2. Known Open Items (Cannot Claim Closure Yet)
### 核心成果
From `sindy_sr_knowledge.md` §三:
1. **Phase-state 特征体系**：`PHASE_STATE_KEYS` (6维) + `ILLUSION_PHASE_KEYS` (10维) — 无需动作历史
2. **绝对动作输出模式**：`get_feature_matrix_deriv(output_mode="absolute")` — 无积分累积
3. **闭合验证器支持**：`predict_v23_deriv(output_mode="absolute")` + 闭环 `mode="abs"`
4. **离线 rollout 评估**：`validate/eval_rollout.py` — 多步滚动误差分析
5. **消融实验**：5 种输入配置 × 2 种输出模式 × CFD 闭环验证
### Illusion 新路线成功
| 场景 | 闭环 sim | 动作历史 | 结论 |
|------|:--------:|:--------:|------|
| 0.75L | **0.974** | 无 | 新路线成立 |
| 1L | **0.958** | 无 | 新路线成立 |
| 1.5L | N/A | 无 | bang-bang/不同机制 |
### Karman 需要补状态
| 配置 | 闭环 sim | 对比旧 v23 |
|------|:--------:|:----------:|
| 旧 v23 (a_lag) | 0.901 | 基线 |
| phase-state → abs (最佳新) | **0.699** | -22% |
| phase-state → deriv | 0.656 | -27% |
1. **"All cloak share same skeleton"** — Cannot claim. Karman and Illusion formulas are structurally different.
2. **"Karman cross-Re joint formula fully solved"** — Cannot claim. Rear is constant (α_T = 3.414), rear control information not fully utilized.
3. **"High-Re degradation proven to be sampling rate issue"** — Strong hypothesis, not proven.
4. **"SR is completely done"** — Karman is complete. Illusion has gaps: 1.5L mechanism unexplained, 1L with-target untested, joint + with-target not run.
---
## 当前最值得优先做的实验
## 3. Unfinished Work from Handover (2026-06-25)
### 第一优先级：Illusion 0.75L / 1L → separate PySR → 公式比较
输入特征（10 维）：
```python
ILLUSION_PHASE_KEYS = [
"u_a", "du_a_dt", # 振荡相位
"Cl_tot", "dCl_tot_dt", # 升力动力学
"Cd_tot", "Cd_rear", # 阻力反馈
"Cd_err", "Cl_err", # 力误差
"dCd_err_dt", "dCl_err_dt", # 误差动力学
]
```
输出：absolute alpha（非维动作，不积分）
环境：`conda run -n sr_env`
代码：`sindy/run_pysr.py`
三个产出要求：
- 最短可接受公式
- 闭环可运行公式
- 0.75L 与 1L 的公式形态比较 / 同形骨架判断
### 第二优先级：Karman 状态补强（CCD/OID 进场）
当前 Karman 新路线最佳 0.699，需要在 phase-state 基础上补充缺失状态量：
- 候选：回流区长度、尾迹中心线偏移、POD 模态系数
- 方法：CCD（Correlation-based Decomposition）或 OID（Observer-based Identification）
- 目标：在无动作历史前提下将 Karman 闭环提升到 0.85+
### 第三优先级：Illusion 跨直径联合
条件：0.75L 和 1L 的 PySR 公式确认同形骨架后
方法：error-state 特征 + 目标 signature 描述符（target St、amplitude）
| # | Item | Priority |
|---|------|----------|
| 1 | Rebuild `sindy_sr_notes.md` | P0 |
| 2 | Write generalization results table into `SR_analysis_report.md` §8 | P0 |
| 3 | Fix `illusion_joint_formula_analysis.md` 1.5L line 46 value | P0 |
| 4 | README directory tree: add `gen_illusion_target.py`, `batch_illusion_generalization.sh` | P0 |
| 5 | README Key Documentation table: add `illusion_joint_formula_analysis.md` | P0 |
| 6 | Create `validate/results/README.md` | P0 |
| 7 | Illusion large-diameter deep analysis (1.5L/2.0L) | P2 |
| 8 | 1L with-target CFD validation | P3 |
| 9 | Illusion joint + with-target search | Future |
| 10 | re400 independent PySR formula CFD | Future |
| 11 | Erase scene (no PPO data yet) | Blocked |
---
## 当前暂不建议做的事
## 4. Recommended Formula Files
- 不回退到动作历史主导（Illusion 已经证明不需要）
- 不先做跨 Re Karman 联合拟合（单 Re 还没站稳）
- 不先做 vortex 扩展（不是当前瓶颈）
- 不先上 PySR 到 Karman（状态还有问题，公式不会更物理）
### Canonical (use these for all papers/reports)
| Scene | Front JSON | Top JSON | CFD Similarity |
|-------|-----------|----------|:---:|
| Karman cross-Re (joint) | `validate/results/karman_joint_deep_front.json` | `validate/results/karman_joint_deep_top.json` | 0.847 avg |
| Illusion joint (0.75L+1L) | `validate/results/pysr_illusion_joint_front.json` | `validate/results/pysr_illusion_joint_top.json` | 0.978/0.970 |
| Karman re50 (individual) | `validate/results/pysr_karman_re50_front_deep.json` | `validate/results/pysr_karman_re50_top_deep.json` | 0.895 |
| Karman re100 (individual) | `validate/results/pysr_karman_re100_front_deep.json` | `validate/results/pysr_karman_re100_top_deep.json` | 0.888 |
| Karman re200 (individual) | `validate/results/pysr_karman_re200_front_deep.json` | `validate/results/pysr_karman_re200_top_deep.json` | 0.916 |
| Illusion 0.75L (individual) | `validate/results/pysr_illusion_0.75L_front.json` | `validate/results/pysr_illusion_0.75L_top.json` | 0.979 |
| Illusion 1L (individual) | `validate/results/pysr_illusion_1L_front.json` | `validate/results/pysr_illusion_1L_top.json` | 0.957 |
---
## 常用命令
## 5. Priority Markers
| Marker | Meaning |
|--------|---------|
| **P0** | Must complete for paper-ready state |
| **P1** | Analysis deepening (defer until P0 done) |
| **P2** | Research extension (1.5L/2.0L mechanism) |
| **P3** | Optional micro-experiments |
| **Future** | Nice-to-have, not urgent |
| **Blocked** | Prerequisites not met |
---
## 6. Quick Commands Reference
### SINDy 拟合
```bash
# Illusion phase-state + absolute action
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes illusion_0.75L,illusion_1L --deriv --phase --output-mode absolute
# PPO inference (pycuda_3_10, GPU 2)
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
# Karman phase-state + absolute action
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re100 --deriv --phase --output-mode absolute
# Target generation for generalization scenes
conda run -n pycuda_3_10 python src/SR_analysis/scripts/gen_illusion_target.py --scene all --device 2
# Karman expanded features
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re100 --deriv --karman-expand --output-mode absolute
# PySR search (sr_env)
conda run -n sr_env python src/SR_analysis/sindy/run_pysr_deep.py --both
conda run -n sr_env python src/SR_analysis/sindy/run_pysr_deep_illusion.py --both
# Karman phase-state + mu modulation
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes karman_re100 --deriv --karman-mu --output-mode absolute
# CFD closed-loop validation (pycuda_3_10, GPU 2)
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
# Illusion separate with error-state (old style for comparison)
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes illusion_0.75L,illusion_1L,illusion_1.5L
# Karman cross-Re joint
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_all_v2.py \
--scenes "karman_re50,karman_re100,karman_re200,karman_re400" --joint
```
### CFD 闭环验证
```bash
# Illusion phase-state + absolute action
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop_illusion.py \
--scene illusion_0.75L --device 0 --steps 320 \
--sindy-results src/SR_analysis/sindy/illusion/sindy_results_deriv.json
--scene illusion_1L --device 2 --steps 214 --mode pysr \
--pysr-front validate/results/pysr_illusion_joint_front.json \
--pysr-top validate/results/pysr_illusion_joint_top.json
# Karman phase-state + absolute action
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \
--scene karman_re100 --device 1 --steps 200 --mode abs \
--sindy-results src/SR_analysis/sindy/karman/sindy_results_deriv.json
# Karman old v23 for comparison
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \
--scene karman_re100 --device 1 --steps 200 --mode v23 \
--sindy-results src/SR_analysis/sindy/karman/sindy_joint_wrapped.json
# Batch generalization validation
bash src/SR_analysis/validate/batch_illusion_generalization.sh
```
### 离线 rollout 评估
```bash
python3 src/SR_analysis/validate/eval_rollout.py \
--sindy-results src/SR_analysis/sindy/karman/sindy_results_deriv.json \
--scene karman_re100
```
### PySR 符号回归
```bash
conda run -n sr_env python src/SR_analysis/sindy/run_pysr.py \
--scene illusion_1L
```
---
## 关键文件索引
| 文件 | 用途 |
|------|------|
| `configs.py` | 统一场景元数据（10+场景） |
| `utils/feature_builder.py` | 特征工程：PHYSICS_FEAT_KEYS, PHASE_STATE_KEYS, ILLUSION_PHASE_KEYS, 导数/滞后特征 |
| `utils/sindy_fitter.py` | STLSQ拟合：get_feature_matrix_v2, get_feature_matrix_deriv, compute_action_deriv |
| `utils/cfd_interface.py` | LegacyCelerisLab封装（GPU, 推理, norm） |
| `utils/g_operator.py` | 等变性诊断 |
| `scripts/infer_*.py` | Karman/Illusion/Vortex 推理管线 |
| `sindy/run_all_v2.py` | 统一 SINDy 拟合入口（支持 --deriv, --phase, --output-mode, --karman-expand 等） |
| `sindy/run_pysr.py` | 受限 PySR 符号回归 |
| `sindy/wrap_joint.py` | 联合模型 → wrapped 格式 |
| `validate/run_closed_loop.py` | Karman 闭环验证器（v23, deriv, abs 模式） |
| `validate/run_closed_loop_illusion.py` | Illusion 闭环验证器 |
| `validate/eval_rollout.py` | **新**：离线多步 rollout 评估 |
| `compare/support_overlap.py` | 跨场景 support 比较 |
| `compare/shared_core.py` | 多场景 shared core 检测 |
| `docs/SR_analysis_results.md` | **新**：完整分析报告 |
| `docs/figures/SR_analysis/fig*.png` | **新**：结果图表（6张） |
## 环境
| 环境 | conda 名 | 用途 |
|------|---------|------|
| CFD + DRL + SINDy | `pycuda_3_10` | infer, run_all_v2, closed-loop |
| PySR 符号回归 | `sr_env` | run_pysr.py |
| GPU | device 0, device 1 | |
src/SR_analysis/validate/run_closed_loop.py
+42 -7
View File
@@ -436,6 +436,13 @@ def run_validation(
feat_keys_front=coefs["feat_keys_front"],
feat_keys_rear=coefs["feat_keys_rear"],
output_mode="absolute")
elif mode == "pysr":
from SR_analysis.validate.predict_pysr import pysr_predict_v23
omega_phys = pysr_predict_v23(
obs, a_prev_phys, a_prev2_phys, mu, u0, dt_c,
coefs["front_func"], coefs["top_func"],
feat_keys_front=coefs["feat_keys_front"],
feat_keys_rear=coefs["feat_keys_rear"])
elif mode == "unstructured":
omega_phys = predict_unstructured(
obs, a_prev_phys, a_prev2_phys, mu, u0,
@@ -487,21 +494,49 @@ def main():
ap.add_argument("--device", type=int, default=2, help="GPU device")
ap.add_argument("--steps", type=int, default=100)
ap.add_argument("--mode", type=str, default="v23",
choices=["v23", "unstructured", "deriv", "abs"],
help="Control law mode: v23 (default), unstructured, deriv, or abs")
choices=["v23", "unstructured", "deriv", "abs", "pysr"],
help="Control law mode: v23 (default), unstructured, deriv, abs, or pysr")
ap.add_argument("--sindy-results", type=str, default=None,
help="Path to sindy_results.json")
help="Path to sindy_results.json (not used in pysr mode)")
ap.add_argument("--threshold", type=float, default=None,
help="SINDy threshold for sparsity (default: best_R2)")
ap.add_argument("--pysr-front", type=str, default=None,
help="PySR front JSON (required for --mode pysr)")
ap.add_argument("--pysr-top", type=str, default=None,
help="PySR top JSON (required for --mode pysr)")
ap.add_argument("--out", type=str, default=None,
help="Output directory for result JSON")
args = ap.parse_args()
if args.sindy_results is None:
args.sindy_results = os.path.join(
os.path.dirname(__file__), "..", "sindy", "karman", "sindy_results_v2.json")
if args.mode == "pysr":
# Load PySR formulas directly
if args.pysr_front is None or args.pysr_top is None:
raise ValueError("--pysr-front and --pysr-top required for --mode pysr")
from SR_analysis.validate.predict_pysr import (
load_pysr_formulas, build_pysr_functions,
)
front_expr, top_expr, fn_f, fn_r = load_pysr_formulas(
args.pysr_front, args.pysr_top)
front_func, top_func = build_pysr_functions(
front_expr, top_expr, fn_f, fn_r)
# Package as pseudo-coefs dict for run_validation
feat_keys_front = [k for k in fn_f if k != "bias"]
feat_keys_rear = [k for k in fn_r if k != "bias"]
coefs = {
"mode": "pysr",
"front_func": front_func,
"top_func": top_func,
"feat_keys_front": feat_keys_front,
"feat_keys_rear": feat_keys_rear,
"feat_names_front": fn_f,
"feat_names_rear": fn_r,
}
else:
if args.sindy_results is None:
args.sindy_results = os.path.join(
os.path.dirname(__file__), "..", "sindy", "karman", "sindy_results_v2.json")
coefs = load_sindy_coefs(args.sindy_results, args.scene, threshold=args.threshold)
coefs = load_sindy_coefs(args.sindy_results, args.scene, threshold=args.threshold)
result = run_validation(args.scene, coefs, args.device,
n_steps=args.steps, mode=args.mode)
src/SR_analysis/validate/run_closed_loop_illusion.py
+202 -13
View File
@@ -54,12 +54,15 @@ def run_validation_illusion(
n_steps: int = 0,
threshold: Optional[float] = None,
out_dir: Optional[str] = None,
mode_override: Optional[str] = None,
) -> dict:
"""Run closed-loop validation for an Illusion scene.
Parameters
----------
n_steps : int, default=0 (auto: >= 2*NX/U0 / sample_interval)
mode_override : str, optional
Override auto-detected mode: "v23", "absolute", or None for auto-detect.
"""
cfg = get_scene(scene_name)
u0 = cfg["u0"]
@@ -175,8 +178,10 @@ def run_validation_illusion(
a_prev = bias_norm.copy()
a_prev2 = a_prev.copy()
# Determine mode from SINDy results metadata (default: v23)
sindy_mode = coefs.get("mode", "v23")
# Determine mode from SINDy results metadata (default: v23), unless overridden
sindy_mode = mode_override if mode_override else coefs.get("mode", "v23")
if mode_override:
print(f" mode overridden: {mode_override}")
for step in range(n_steps):
obs = fifo[-1] if fifo else np.zeros(12, dtype=np.float32)
@@ -190,7 +195,10 @@ def run_validation_illusion(
if target_harmonics is not None:
target_forces_step = gen_target_states_at(step, target_harmonics)
if sindy_mode == "absolute":
if sindy_mode == "pysr":
# PySR mode -- should not reach here (uses run_validation_pysr_illusion instead)
raise ValueError("PySR mode should use run_validation_pysr_illusion()")
elif sindy_mode == "absolute":
omega_phys = predict_v23_deriv(
obs, obs_prev, a_prev_phys, a_prev2_phys, mu, u0, sample_interval/2000.0,
coefs["front_coef"], coefs["top_coef"],
@@ -296,6 +304,162 @@ def run_validation_illusion(
}
def run_validation_pysr_illusion(
scene_name: str,
device_id: int,
front_func,
top_func,
feat_names_front: list,
feat_names_rear: list,
n_steps: int = 0,
out_dir: Optional[str] = None,
) -> dict:
"""Run closed-loop validation using PySR formulas for an Illusion scene."""
from SR_analysis.validate.predict_pysr import pysr_predict_v23
cfg = get_scene(scene_name)
u0 = cfg["u0"]
mu = cfg["mu"]
l0 = 20.0
sample_interval = cfg["sample_interval"]
conv_len = cfg.get("conv_len", 36)
action_scale = cfg["action_scale"]
action_bias = cfg["action_bias"]
n_obj_total = cfg["n_objects_env"]
sensor_x = cfg["sensor_x"]
front_x = cfg["pinball_front_x"]
rear_x = cfg["pinball_rear_x"]
min_steps = int(1 * NX / u0 / sample_interval)
if n_steps == 0 or n_steps < min_steps:
n_steps = max(min_steps, 200)
print(f" auto-set steps={n_steps} (min_steps={min_steps})")
print(f"\n=== Validating {scene_name} (mode=pysr, device={device_id}, steps={n_steps}) ===")
if out_dir is None:
out_dir = os.path.join(os.path.dirname(__file__), "results")
os.makedirs(out_dir, exist_ok=True)
# Load target data
data_dir = os.path.join(
os.path.dirname(__file__), "..", "data", "illusion", scene_name,
)
target_npz = np.load(os.path.join(data_dir, "target.npz"))
target_states = target_npz.get("target_sensors", target_npz["target_states"][:, 2:8])
print(f" target loaded: {target_states.shape}")
target_harmonics = None
harm_path = os.path.join(data_dir, "target_harmonics.json")
if os.path.isfile(harm_path):
with open(harm_path) as f:
target_harmonics = json.load(f)
print(f" target harmonics loaded: {len(target_harmonics)} channels")
# Build environment
cuda_cfg, field_cfg = load_legacy_configs(LEGACY_CFG_DIR)
field_cfg = field_cfg._replace(viscosity=float(cfg["nu"]))
ff = FlowField(field_cfg, cuda_cfg, device_id=device_id)
ny = ff.FIELD_SHAPE[1]
for y_off in [2.0, 0.0, -2.0]:
sc = (sensor_x * l0, (ny - 1) / 2 + y_off * l0, 0.0)
ff.add_sensor(sc, l0 / 4.0)
ff.add_cylinder((front_x * l0, (ny - 1) / 2, 0.0), l0 / 2.0)
ff.add_cylinder((rear_x * l0, (ny - 1) / 2 + 0.75 * l0, 0.0), l0 / 2.0)
ff.add_cylinder((rear_x * l0, (ny - 1) / 2 - 0.75 * l0, 0.0), l0 / 2.0)
n_obj = ff.obs.size // 2
assert n_obj == 6, f"Expected 6 objects, got {n_obj}"
stabilize_steps = int(4 * NX / u0)
print(f" stabilising pinball ({stabilize_steps} steps)...")
ff.run(stabilize_steps, np.zeros(n_obj, dtype=DATA_TYPE))
ff.get_ddf()
ff.save_ddf()
# Norm + bias FIFO (same as run_validation_illusion)
print(f" norm collection ({FIFO_LEN} steps)...")
fifo = deque(maxlen=FIFO_LEN)
for _ in range(FIFO_LEN):
ff.run(sample_interval, np.zeros(n_obj, dtype=DATA_TYPE))
fifo.append(ff.obs.copy()[0:12])
ff.apply_ddf()
bias_arr = np.zeros(n_obj, dtype=DATA_TYPE)
bias_arr[3] = 0.0
bias_arr[4] = -1.0 * u0
bias_arr[5] = 1.0 * u0
fifo.clear()
for _ in range(FIFO_LEN):
ff.run(sample_interval, bias_arr)
fifo.append(ff.obs.copy()[0:12])
ff.get_ddf()
ff.save_ddf()
ff.apply_ddf()
# Closed-loop with PySR
feat_keys_front = [k for k in feat_names_front if k != "bias"]
feat_keys_rear = [k for k in feat_names_rear if k != "bias"]
sens_list, actions_list = [], []
bias_norm = np.array([0.0, 0.125, -0.125], dtype=np.float64)
a_prev = bias_norm.copy()
a_prev2 = a_prev.copy()
for step in range(n_steps):
obs = fifo[-1] if fifo else np.zeros(12, dtype=np.float32)
a_prev_phys = (a_prev * action_scale + np.array(action_bias, dtype=np.float64)) * u0
a_prev2_phys = (a_prev2 * action_scale + np.array(action_bias, dtype=np.float64)) * u0
target_forces_step = None
if target_harmonics is not None:
target_forces_step = gen_target_states_at(step, target_harmonics)
omega_phys = pysr_predict_v23(
obs, a_prev_phys, a_prev2_phys, mu, u0, sample_interval / 2000.0,
front_func, top_func,
feat_keys_front, feat_keys_rear,
target_forces=target_forces_step,
)
norm_a = (omega_phys / u0 - np.array(action_bias, dtype=np.float64)) / action_scale
norm_a = np.clip(norm_a, -1.0, 1.0).astype(np.float32)
action_arr = scale_action(norm_a, scale=action_scale, bias=action_bias,
u0=u0, n_total_bodies=n_obj_total)
ff.context.push()
ff.run(sample_interval, action_arr)
ff.context.pop()
obs_new = ff.obs.copy()[0:12]
fifo.append(obs_new)
sens_list.append(obs_new[0:6])
actions_list.append(omega_phys.copy())
a_prev2 = a_prev.copy()
a_prev = norm_a.copy().astype(np.float64)
sens_arr = np.array(sens_list, dtype=np.float32)
from SR_analysis.utils.cfd_interface import compute_similarity
sim_full = compute_similarity(target_states, sens_arr, conv_len)
sim_tail = compute_similarity(target_states[-conv_len:], sens_arr[-conv_len:], conv_len // 2)
action_range = float(np.max(np.abs(actions_list)))
print(f" similarity (full)={sim_full:.4f} (tail)={sim_tail:.4f} action_range={action_range:.4f}")
del ff
return {
"scene": scene_name,
"mode": "pysr",
"similarity_full": sim_full,
"similarity_tail": sim_tail,
"action_range": action_range,
"n_steps": n_steps,
}
def main():
ap = argparse.ArgumentParser(description="Illusion closed-loop SINDy validation")
ap.add_argument("--scene", type=str, required=True)
@@ -304,22 +468,47 @@ def main():
help="Steps (default: auto-set to cover 1*NX/U0)")
ap.add_argument("--threshold", type=float, default=None)
ap.add_argument("--sindy-results", type=str, default=None)
ap.add_argument("--mode", type=str, default=None,
help='Override mode: "v23", "absolute", "pysr" (default: auto-detect)')
ap.add_argument("--pysr-front", type=str, default=None,
help="PySR front JSON (required for --mode pysr)")
ap.add_argument("--pysr-top", type=str, default=None,
help="PySR top JSON (required for --mode pysr)")
ap.add_argument("--out", type=str, default=None)
args = ap.parse_args()
if args.sindy_results is None:
args.sindy_results = os.path.join(
os.path.dirname(__file__), "..", "sindy", "illusion", "sindy_results_v2.json")
if args.mode == "pysr":
if args.pysr_front is None or args.pysr_top is None:
raise ValueError("--pysr-front and --pysr-top required for --mode pysr")
from SR_analysis.validate.predict_pysr import (
load_pysr_formulas, build_pysr_functions, pysr_predict_v23,
)
front_expr, top_expr, fn_f, fn_r = load_pysr_formulas(
args.pysr_front, args.pysr_top)
front_func, top_func = build_pysr_functions(
front_expr, top_expr, fn_f, fn_r)
result = run_validation_pysr_illusion(
args.scene, args.device, front_func, top_func, fn_f, fn_r,
n_steps=args.steps, out_dir=args.out,
)
else:
if args.sindy_results is None:
args.sindy_results = os.path.join(
os.path.dirname(__file__), "..", "sindy", "illusion", "sindy_results_v2.json")
result = run_validation_illusion(
args.scene, args.sindy_results, args.device,
n_steps=args.steps, threshold=args.threshold, out_dir=args.out,
mode_override=args.mode,
)
result = run_validation_illusion(
args.scene, args.sindy_results, args.device,
n_steps=args.steps, threshold=args.threshold, out_dir=args.out,
)
th_str = f"_th{args.threshold}" if args.threshold is not None else ""
if args.mode == "pysr":
out_name = f"{args.scene}_pysr"
else:
th_str = f"_th{args.threshold}" if args.threshold is not None else ""
out_name = f"{args.scene}_v23{th_str}"
out_dir = args.out or os.path.join(os.path.dirname(__file__), "results")
os.makedirs(out_dir, exist_ok=True)
out_path = os.path.join(out_dir, f"{args.scene}_v23{th_str}.json")
out_path = os.path.join(out_dir, f"{out_name}.json")
with open(out_path, "w") as f:
json.dump(result, f, indent=2)
print(f"Saved: {out_path}")