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

第二轮：整理两个工作目录

Browse Source
This commit is contained in:
Frank14f
2026-06-10 15:59:52 +08:00
parent d1b9922c6b
commit 096d9dcd0f
130 changed files with 13171 additions and 7263 deletions
Download Patch File Download Diff File Expand all files Collapse all files
src/SR_analysis/README.md
+289
View File
@@ -0,0 +1,289 @@
# SR_analysis: Unified SINDy-SR Analysis Pipeline
## Overview
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.
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.
For background, see:
- `src/sindy_sr_notes.md` -- execution plan
- `src/sindy_sr_knoeledge.md` -- confirmed facts and known pitfalls
## Directory Structure
```
SR_analysis/
configs.py # Unified scene metadata (all 10 scenes)
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 (from analysis_cloak)
sindy_fitter.py # STLSQ threshold grid, feature matrix builder
cfd_interface.py # LegacyCelerisLab wrapper (requires pycuda_3_10)
g_operator.py # Equivariance diagnostics
data/
karman/ # Karman cloak: karman_re50, re100, re200, re400
steady/ # Steady cloak: steady_data.npz
illusion/ # Illusion: illusion_0.75L, illusion_1L, illusion_1.5L
vortex/ # Vortex cloak: vortex_lamb, vortex_taylor
scripts/
infer_karman.py # Inference: LegacyCFD + PPO -> controlled.npz
infer_illusion.py # Inference: for 0.75L, 1L, 1.5L diameters
infer_vortex.py # Inference: for Lamb dipole + Taylor monopole
sindy/
run_karman.py # SINDy fitting for Karman scenes
run_illusion.py # SINDy fitting for Illusion scenes
run_vortex.py # SINDy fitting for Vortex scenes
run_pareto.py # Pareto-front analysis from SINDy results
karman/ # Output: sindy_results.json, pareto_*.json
illusion/ # Output: sindy_results.json, pareto_*.json
vortex/ # Output: sindy_results.json, pareto_*.json
validate/
run_closed_loop.py # Unified closed-loop validator (v23 + unstructured modes)
compare/
support_overlap.py # Pairwise support set comparison
shared_core.py # Multi-scene shared-core detection
```
## Key Design Decisions
### 1. Scene Metadata Driven
All scene parameters (Re, action scaling, geometry, model paths) are defined
once in `configs.py`, not hard-coded in scripts. Adding a new scene means
adding one dict to `configs.py`.
### 2. Data / Features / Models Separation
- `data/` -- raw sensor/force/action arrays (.npz), one-time generation
- `sindy/` -- SINDy fitting results (JSON), reusable for comparison
- `scripts/` -- inference pipelines that produce `data/`
### 3. Unified Feature Builder
Every scene uses the same `utils/feature_builder.py`, which produces
21 dimensionless features from raw lattice-unit sensor/force data:
**Sensor features (nondim):**
- `u_m`, `u_a`, `u_c` -- streamwise: mean, antisymmetric, centre
- `v_a` -- antisymmetric cross-stream
- `sin_ua`, `cos_ua` -- phase encoding via u_a
**Force features (Cd/Cl):**
- `Cd_tot`, `Cd_rear` -- total and rear-cylinder drag
- `Cl_tot`, `Cl_diff` -- total and differential lift
**Memory features (nondim alpha):**
- `aF_lag1`, `aB_lag1`, `aT_lag1` -- lagged actions (t-1)
- `daF`, `daB`, `daT` -- action increments (t-1)-(t-2)
**Reynolds modulation:**
- `mu` (= 1/Re_D), `mu_u_a`, `mu_v_a`, `mu_Cd_tot`, `mu_Cl_diff`
### 4. G-Equivariant Structure (v23)
Default control law structure (confirmed as the best v23 model):
```
Front(t) = f_front(x(t)) # no bias, odd under G
Top(t) = f_rear(x(t)) # with bias
Bottom(t) = -f_rear(G[x(t)]) # shared-head: bottom = -top(Gx)
```
Where G is the mirror operator (y -> -y) with corrected sign rules:
- `[aF, aT, aB] -> [-aF, -aB, -aT]`
- Sensor swap: top <-> bottom, negate v
- Force swap: front unchanged, bottom <-> top, negate Cl
### 5. STLSQ Threshold Grid
Default thresholds: `[0, 0.001, 0.002, 0.005, 0.01, 0.015, 0.02, 0.03, 0.05, 0.1]`
Per-channel: front (no bias), top (shared-head), bottom (independent, for comparison)
## Scene Inventory
| Scene Name | Description | Re_code | Sample Interval | Action | U0 |
|---|---|---|---|---|---|
| karman_re50 | Karman cloak at low Re | 50 | 800 | 8x + [0,-4,4] | 0.01 |
| karman_re100 | Karman cloak (default) | 100 | 800 | 8x + [0,-4,4] | 0.01 |
| karman_re200 | Karman cloak at high Re | 200 | 800 | 8x + [0,-4,4] | 0.01 |
| karman_re400 | Karman cloak at highest Re | 400 | 800 | 8x + [0,-4,4] | 0.01 |
| steady | Open-loop constant rotation | 100 | 800 | 8x + [0,-5.1,5.1] | 0.01 |
| illusion_0.75L | Imitate 0.75D cylinder | 100 | 600 | 8x + [0,-2,2] | 0.01 |
| illusion_1L | Imitate 1.0D cylinder | 100 | 600 | 8x + [0,-2,2] | 0.01 |
| illusion_1.5L | Imitate 1.5D cylinder | 100 | 600 | 8x + [0,-2,2] | 0.02 |
| vortex_lamb | Cloak Lamb dipole | 100 | 800 | 4x + [0,-4,4] | 0.01 |
| vortex_taylor | Cloak Taylor monopole | 100 | 800 | 4x + [0,-4,4] | 0.01 |
Note: "Re_code" uses reference length 2*D (code convention).
Physical Re_D = Re_code / 2. E.g. Re_code=100 -> Re_D=50.
## Re-generation Commands
All commands run from repo root (`/home/frank14f/DynamisLab`).
### Data Generation (requires GPU, pycuda_3_10 env)
```bash
# Karman cloak -- all 4 training Re
conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_karman.py --re all --device 0
# Karman cloak -- single Re
conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_karman.py --re 100 --device 0 --steps 200
# Illusion -- all 3 diameters
conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_illusion.py --diameter all --device 0
# Vortex -- both types
conda run -n pycuda_3_10 python src/SR_analysis/scripts/infer_vortex.py --type all --device 0
```
### SINDy Fitting (no GPU needed, pycuda_3_10 env for pysindy)
```bash
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_karman.py
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_illusion.py
conda run -n pycuda_3_10 python src/SR_analysis/sindy/run_vortex.py
```
### Pareto Analysis (no GPU, no conda needed)
```bash
python3 src/SR_analysis/sindy/run_pareto.py --scene karman_re100
python3 src/SR_analysis/sindy/run_pareto.py --scene illusion_1L
```
### Closed-loop Validation (requires GPU)
```bash
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \
--scene karman_re70 --device 2 \
--sindy-results src/SR_analysis/sindy/karman/sindy_results.json
# With custom mode
conda run -n pycuda_3_10 python src/SR_analysis/validate/run_closed_loop.py \
--scene karman_re70 --device 2 --mode unstructured
```
### Cross-scene Comparison (no GPU)
```bash
# Pairwise support overlap
python3 src/SR_analysis/compare/support_overlap.py \
--sindy-results src/SR_analysis/sindy/karman/sindy_results.json \
--scenes karman_re100 illusion_1L
# Multi-scene shared core
python3 src/SR_analysis/compare/shared_core.py \
--sindy-results src/SR_analysis/sindy/karman/sindy_results.json \
--scenes karman_re50 karman_re100 karman_re200 karman_re400
```
## Key Results Summary
### Data Quality (similarity scores)
| Scene | PPO Similarity |
|---|---|
| karman_re50 | 0.962 |
| karman_re100 | 0.954 |
| karman_re200 | 0.884 |
| karman_re400 | 0.795 (inferred, not verified) |
| vortex_lamb | 0.942 |
| vortex_taylor | 0.916 |
| illusion_1L | ~0.55 (metric not directly comparable) |
### SINDy Fit Quality (R2 scores for one-step prediction)
| Scene | Front | Top (shared) | Bottom |
|---|---|---|---|
| karman_re50 | 0.998 | 0.989 | 0.996 |
| karman_re100 | 0.995 | 0.993 | 0.997 |
| karman_re200 | 0.957 | 0.914 | 0.918 |
| karman_re400 | 0.991 | 0.979 | 0.969 |
| illusion_0.75L | 0.991 | 0.989 | 0.990 |
| illusion_1L | 0.979 | 0.984 | 0.984 |
| illusion_1.5L | 0.959 | 0.928 | 0.932 |
| vortex_lamb | 0.904 | 0.980 | 0.933 |
| vortex_taylor | 0.960 | 0.810 | 0.643 |
### Shared Core Features
**Karman cross-Re (active in all re50/100/200):**
- Front core: `mu`, `mu_Cd_tot`, `mu_Cl_diff`, `mu_v_a` (mu-modulated terms dominate)
- Top core: `Cl_tot`, `bias`, `mu_Cd_tot`, `mu_Cl_diff`, `mu_u_a`, `mu_v_a`
- Scene-specific: lower-Re scenes have additional `Cd_tot`, `Cl_diff`, `aT_lag1` etc.
**Illusion cross-diameter (active in all 0.75L/1L/1.5L):**
- Front core: `mu`, `mu_Cd_tot`, `mu_Cl_diff` (same structure as Karman front!)
- Top core: `Cd_rear`, `Cl_tot`, `bias`, `mu_Cd_tot`, `mu_Cl_diff`
- This suggests a **shared mu-modulated feedback structure** exists across both scenes
## Known Issues and Caveats
1. **Vortex Taylor rear channels** have low R2 (0.64-0.81). The weak monopole
produces near-zero rear action, making SINDy fitting noisy. Use Lamb as the
primary vortex reference.
2. **Closed-loop validator** (`validate/run_closed_loop.py`) has been ported but
NOT yet tested end-to-end. The original `validate_v23.py` verified Karman
but the new unified version has not been run.
3. **Illusion similarity scores** use the Karman CONV_LEN=30 metric, giving
lower raw numbers. The controlled.npz data itself is valid for SINDy.
4. **Steady cloak** is open-loop constant rotation, not PPO-derived. It serves
as a physical consistency check, not a primary comparison scene.
5. **SINDy one-step R2 is not sufficient** -- a high R2 does not guarantee good
closed-loop performance. Always validate via `validate/run_closed_loop.py`.
6. **Scene key naming**: keys like `illusion_1L`, `illusion_1.5L` use the short
float format from Python (1.0 -> "1L", 1.5 -> "1.5L", 0.75 -> "0.75L").
## Next Steps (Future Work)
1. **PySR symbolic regression** -- Run PySR on the SINDy-identified active
features (in `sr_env` conda env) to find closed-form formulas. Essential
reading: `src/pysr.md`.
2. **Closed-loop validation of all new scenes** -- Run
`validate/run_closed_loop.py` for illusion and vortex scenes using their
SINDy coefficients.
3. **Cross-scene shared backbone test** -- Fit a single SINDy model on merged
Karman + Illusion data, test if it performs on both.
4. **Time-scale explicit formulation** -- Make the sample interval an explicit
feature to compare control laws across different frequencies.
5. **Steady as consistency check** -- Validate that Karman-derived control laws
can reproduce the steady cloak result as a sanity check.
## File Reference
| File | Lines | Purpose |
|---|---|---|
| configs.py | ~205 | Unified scene metadata |
| utils/feature_builder.py | ~212 | Dimensionless features + G-op |
| utils/sindy_fitter.py | ~175 | STLSQ fitting, feature matrix builder |
| utils/cfd_interface.py | ~370 | LegacyCelerisLab wrapper |
| utils/g_operator.py | ~170 | Equivariance diagnostics |
| utils/__init__.py | ~10 | Selective exports |
| scripts/infer_karman.py | ~250 | Karman inference pipeline |
| scripts/infer_illusion.py | ~270 | Illusion inference pipeline |
| scripts/infer_vortex.py | ~280 | Vortex inference pipeline |
| sindy/run_karman.py | ~160 | Karman SINDy fitting |
| sindy/run_illusion.py | ~110 | Illusion SINDy fitting |
| sindy/run_vortex.py | ~110 | Vortex SINDy fitting |
| sindy/run_pareto.py | ~140 | Pareto analysis |
| validate/run_closed_loop.py | ~270 | Closed-loop validator |
| compare/support_overlap.py | ~150 | Pairwise support comparison |
| compare/shared_core.py | ~140 | Multi-scene shared core detection |
src/SR_analysis/compare/illusion_shared_core.json
+185
View File
@@ -0,0 +1,185 @@
{
"scenes": [
"illusion_0.75L",
"illusion_1L",
"illusion_1.5L"
],
"threshold": 0.02,
"channels": {
"front": {
"n_scenes": 3,
"n_core": 3,
"core_features": {
"mu": {
"group": "mu_mod",
"coef": {
"mean": 6.124328107212669,
"std": 12.593845646937403,
"per_scene": {
"illusion_0.75L": -5.325014552992803,
"illusion_1L": 23.663898368547343,
"illusion_1.5L": 0.03410050608346959
}
}
},
"mu_Cd_tot": {
"group": "mu_mod",
"coef": {
"mean": 1.2528987158986185,
"std": 0.6780019496130162,
"per_scene": {
"illusion_0.75L": 0.852203136137665,
"illusion_1L": 2.2076417965839896,
"illusion_1.5L": 0.6988512149742004
}
}
},
"mu_Cl_diff": {
"group": "mu_mod",
"coef": {
"mean": -0.7229594207792337,
"std": 0.4999554272357686,
"per_scene": {
"illusion_0.75L": -0.23779821359027198,
"illusion_1L": -0.5201221350686689,
"illusion_1.5L": -1.4109579136787602
}
}
}
},
"scene_specific": {
"illusion_0.75L": [
"mu_v_a"
],
"illusion_1.5L": [
"Cd_rear",
"Cl_diff",
"Cl_tot",
"mu_u_a"
]
}
},
"top": {
"n_scenes": 3,
"n_core": 5,
"core_features": {
"Cd_rear": {
"group": "force",
"coef": {
"mean": -0.012196232003454469,
"std": 0.04820789952138249,
"per_scene": {
"illusion_0.75L": -0.06463454014637952,
"illusion_1L": 0.05175447708136916,
"illusion_1.5L": -0.02370863294535305
}
}
},
"Cl_tot": {
"group": "force",
"coef": {
"mean": 0.06996842377218454,
"std": 0.01752531021573054,
"per_scene": {
"illusion_0.75L": 0.05319377957208298,
"illusion_1L": 0.09415648101989169,
"illusion_1.5L": 0.06255501072457896
}
}
},
"bias": {
"group": "bias",
"coef": {
"mean": 0.50722496890139,
"std": 0.7941369244496769,
"per_scene": {
"illusion_0.75L": 1.5684596958510442,
"illusion_1L": 0.2949090708624881,
"illusion_1.5L": -0.3416938600093622
}
}
},
"mu_Cd_tot": {
"group": "mu_mod",
"coef": {
"mean": -0.5748401414253728,
"std": 1.3356390970379715,
"per_scene": {
"illusion_0.75L": -0.034745982702969365,
"illusion_1L": -2.4124080088752256,
"illusion_1.5L": 0.7226335673020766
}
}
},
"mu_Cl_diff": {
"group": "mu_mod",
"coef": {
"mean": -0.14655148885046754,
"std": 0.42192581541357527,
"per_scene": {
"illusion_0.75L": -0.7427999157338147,
"illusion_1L": 0.13162402418325947,
"illusion_1.5L": 0.1715214249991526
}
}
}
},
"scene_specific": {
"illusion_1.5L": [
"Cd_tot"
]
}
},
"bottom": {
"n_scenes": 3,
"n_core": 3,
"core_features": {
"bias": {
"group": "bias",
"coef": {
"mean": -0.07180166971995111,
"std": 0.08175050418248041,
"per_scene": {
"illusion_0.75L": 0.03345586870310438,
"illusion_1L": -0.16584728761620943,
"illusion_1.5L": -0.08301359024674829
}
}
},
"mu_Cd_tot": {
"group": "mu_mod",
"coef": {
"mean": 0.5765450851901482,
"std": 0.4820935221835714,
"per_scene": {
"illusion_0.75L": -0.03489998074530567,
"illusion_1L": 1.1434618843744861,
"illusion_1.5L": 0.6210733519412643
}
}
},
"mu_Cl_diff": {
"group": "mu_mod",
"coef": {
"mean": 0.1661185617369393,
"std": 0.34042296056983895,
"per_scene": {
"illusion_0.75L": 0.6070993919053495,
"illusion_1L": -0.22165490297479395,
"illusion_1.5L": 0.11291119628026249
}
}
}
},
"scene_specific": {
"illusion_0.75L": [
"aF_lag1",
"mu_u_a"
],
"illusion_1.5L": [
"daB"
]
}
}
}
}
src/SR_analysis/compare/karman_shared_core.json
+218
View File
@@ -0,0 +1,218 @@
{
"scenes": [
"karman_re50",
"karman_re100",
"karman_re200"
],
"threshold": 0.02,
"channels": {
"front": {
"n_scenes": 3,
"n_core": 4,
"core_features": {
"mu": {
"group": "mu_mod",
"coef": {
"mean": -0.2118470353646019,
"std": 1.3418924147536682,
"per_scene": {
"karman_re50": 0.5373930579297568,
"karman_re100": 0.9234972689284461,
"karman_re200": -2.0964314329520084
}
}
},
"mu_Cd_tot": {
"group": "mu_mod",
"coef": {
"mean": 0.40446115366058594,
"std": 0.2882660447342232,
"per_scene": {
"karman_re50": 0.21424605058608626,
"karman_re100": 0.18730338573081926,
"karman_re200": 0.8118340246648522
}
}
},
"mu_Cl_diff": {
"group": "mu_mod",
"coef": {
"mean": 0.49780940740664015,
"std": 0.4399793858420079,
"per_scene": {
"karman_re50": -0.12436843773164996,
"karman_re100": 0.8155192599082862,
"karman_re200": 0.8022774000432843
}
}
},
"mu_v_a": {
"group": "mu_mod",
"coef": {
"mean": -0.015356732494797251,
"std": 0.0636302488843807,
"per_scene": {
"karman_re50": -0.01474954032019323,
"karman_re100": 0.0622687182995332,
"karman_re200": -0.09358937546373172
}
}
}
},
"scene_specific": {}
},
"top": {
"n_scenes": 3,
"n_core": 6,
"core_features": {
"Cl_tot": {
"group": "force",
"coef": {
"mean": 0.052870132404981444,
"std": 0.03551999603386398,
"per_scene": {
"karman_re50": 0.005437814717321779,
"karman_re100": 0.09090888202182665,
"karman_re200": 0.0622637004757959
}
}
},
"bias": {
"group": "bias",
"coef": {
"mean": 0.5261511607510746,
"std": 0.3945210247055603,
"per_scene": {
"karman_re50": 0.009006304945500938,
"karman_re100": 0.9660827819346605,
"karman_re200": 0.6033643953730624
}
}
},
"mu_Cd_tot": {
"group": "mu_mod",
"coef": {
"mean": -1.0229551996731832,
"std": 0.6119923552005874,
"per_scene": {
"karman_re50": -0.25178195059980246,
"karman_re100": -1.0682906975825734,
"karman_re200": -1.7487929508371736
}
}
},
"mu_Cl_diff": {
"group": "mu_mod",
"coef": {
"mean": -0.3025183164705127,
"std": 1.1363636819618512,
"per_scene": {
"karman_re50": -0.25926759753886264,
"karman_re100": 1.0671077959431223,
"karman_re200": -1.7153951478157978
}
}
},
"mu_u_a": {
"group": "mu_mod",
"coef": {
"mean": -0.08031494383668124,
"std": 0.08982478076029217,
"per_scene": {
"karman_re50": 0.023534741132533694,
"karman_re100": -0.19559725673163528,
"karman_re200": -0.06888231591094213
}
}
},
"mu_v_a": {
"group": "mu_mod",
"coef": {
"mean": 0.030282322344902513,
"std": 0.13310118373865282,
"per_scene": {
"karman_re50": -0.08783423089354239,
"karman_re100": -0.03758555104936961,
"karman_re200": 0.21626674897761955
}
}
}
},
"scene_specific": {
"karman_re50": [
"Cd_tot",
"Cl_diff",
"aT_lag1"
],
"karman_re200": [
"cos_ua"
]
}
},
"bottom": {
"n_scenes": 3,
"n_core": 4,
"core_features": {
"bias": {
"group": "bias",
"coef": {
"mean": 0.07937251506166355,
"std": 0.3195026042260621,
"per_scene": {
"karman_re50": -0.05202015913990901,
"karman_re100": -0.22933046011719774,
"karman_re200": 0.5194681644420974
}
}
},
"mu_Cd_tot": {
"group": "mu_mod",
"coef": {
"mean": -0.8264381469826733,
"std": 0.9571803403036669,
"per_scene": {
"karman_re50": -0.03551514238315014,
"karman_re100": -0.2705206993407136,
"karman_re200": -2.1732785992241563
}
}
},
"mu_Cl_diff": {
"group": "mu_mod",
"coef": {
"mean": -0.18340550520310042,
"std": 0.23298894865543504,
"per_scene": {
"karman_re50": 0.08969580568091838,
"karman_re100": -0.16030801380832282,
"karman_re200": -0.4796043074818968
}
}
},
"mu_v_a": {
"group": "mu_mod",
"coef": {
"mean": 0.09749842111867539,
"std": 0.07771330943469408,
"per_scene": {
"karman_re50": 0.08216046763939658,
"karman_re100": 0.010919861645215943,
"karman_re200": 0.19941493407141364
}
}
}
},
"scene_specific": {
"karman_re50": [
"Cl_diff",
"daB",
"daF",
"mu",
"u_c",
"u_m",
"v_a"
]
}
}
}
}
src/SR_analysis/compare/shared_core.py
+152
View File
@@ -0,0 +1,152 @@
"""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
@@ -0,0 +1,158 @@
"""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
+219
View File
@@ -0,0 +1,219 @@
"""Unified scene configuration for SR_analysis.
All scene metadata in one place. Each scene dict contains all parameters
needed for data generation, SINDy fitting, and validation.
Re convention:
- "re_code" uses reference length 2*D (matching model file naming).
- mu = 1/Re_D = 2/re_code.
- Re_D = re_code / 2 is the true physical Reynolds number.
"""
from __future__ import annotations
import os
from typing import Any, Dict, List, Optional, Tuple
# -- Root paths (resolved when configs.py is imported) -----------------------
_PROJ = os.path.abspath(os.path.join(os.path.dirname(__file__), ".."))
MODEL_DIR = os.path.join(_PROJ, "..", "models")
LEGACY_CFG_DIR = os.path.join(os.path.dirname(__file__), "configs", "legacy")
# -- Physics constants -------------------------------------------------------
U0 = 0.01 # default inlet velocity (lattice)
D_CYL = 20.0
D_REF = 40.0
L0 = 20.0
NX = 1280
NY = 512
CENTER_Y = (NY - 1) / 2.0
FIFO_LEN = 150
CONV_LEN = 30
def nu_from_re(re_code: float, u0: float = U0) -> float:
"""Viscosity from code Reynolds number (reference length = 2*D)."""
return u0 * D_REF / re_code
# -- Scene definitions -------------------------------------------------------
# Each scene dict has fields:
# scene_id, re_code, mu, nu, has_disturbance, sample_interval,
# action_scale, action_bias (tuple), source ("PPO_inference"|"open_loop"),
# model_name (str or None), n_objects_env, obs_slice [start,end],
# sensor_x, pinball_front_x, pinball_rear_x,
# target_type ("periodic"|"steady"|"transient"),
# s_dim (DRL observation dim, 12 or 14)
SCENES: Dict[str, Any] = {}
# -- Karman Cloak (cross-Re) ------------------------------------------------
for rc, mn in [(50, "d1a3o12_re50"), (100, "d1a3o12_re100"),
(200, "d1a3o12_re200"), (400, "d1a3o12_re400")]:
key = f"karman_re{rc}"
SCENES[key] = {
"scene_id": "karman",
"re_code": rc,
"mu": 2.0 / rc,
"nu": nu_from_re(rc),
"has_disturbance": True,
"sample_interval": 800,
"action_scale": 8.0,
"action_bias": (0.0, -4.0, 4.0),
"source": "PPO_inference",
"model_name": mn,
"n_objects_env": 7,
"obs_slice": (2, 14),
"sensor_x": 40.0,
"pinball_front_x": 30.0,
"pinball_rear_x": 31.3,
"target_type": "periodic",
"s_dim": 12,
"u0": U0,
}
# -- Steady Cloak (open-loop constant rotation) -----------------------------
SCENES["steady"] = {
"scene_id": "steady",
"re_code": 100,
"mu": 2.0 / 100,
"nu": nu_from_re(100),
"has_disturbance": False,
"sample_interval": 800,
"action_scale": 8.0,
"action_bias": (0.0, -5.1, 5.1), # from gen_steady_data.py defaults
"source": "open_loop",
"model_name": None,
"n_objects_env": 6,
"obs_slice": (0, 12),
"sensor_x": 40.0,
"pinball_front_x": 30.0,
"pinball_rear_x": 31.3,
"target_type": "steady",
"s_dim": 12,
"u0": U0,
}
# -- Illusion (cylinder imitation, 3 diameters, 1U=0.01) --------------------
def _illusion_key(diam: float) -> str:
"""Generate clean illusion scene key."""
s = f"{diam:.3f}".rstrip("0").rstrip(".")
return f"illusion_{s}L"
_ILLUSION_1U = [
(0.75, "d1a3o12_250525_imit_075L_1U"),
(1.0, "d1a3o12_250525_imit_1L_1U"),
]
for diam, mn in _ILLUSION_1U:
key = _illusion_key(diam)
SCENES[key] = {
"scene_id": "illusion",
"target_diameter": diam,
"re_code": 100,
"mu": 2.0 / 100,
"nu": nu_from_re(100),
"has_disturbance": False,
"sample_interval": 600,
"action_scale": 8.0,
"action_bias": (0.0, -2.0, 2.0),
"source": "PPO_inference",
"model_name": mn,
"n_objects_env": 6,
"obs_slice": (0, 12),
"sensor_x": 30.0,
"pinball_front_x": 19.0,
"pinball_rear_x": 20.3,
"target_type": "periodic",
"s_dim": 12,
"u0": U0,
}
# 1.5L Illusion (2U=0.02 model)
SCENES[_illusion_key(1.5)] = {
"scene_id": "illusion",
"target_diameter": 1.5,
"re_code": 100,
"mu": 2.0 / 100,
"nu": nu_from_re(100, u0=0.02),
"has_disturbance": False,
"sample_interval": 600,
"action_scale": 8.0,
"action_bias": (0.0, -2.0, 2.0),
"source": "PPO_inference",
"model_name": "d1a3o14_250525_imit_15L_2U",
"n_objects_env": 6,
"obs_slice": (0, 12),
"sensor_x": 30.0,
"pinball_front_x": 19.0,
"pinball_rear_x": 20.3,
"target_type": "periodic",
"s_dim": 14,
"u0": 0.02,
}
# -- Vortex Cloak (Lamb dipole + Taylor monopole) --------------------------
_SCENES_VORTEX = [
("lamb", "vortex_lamb", 0.5),
("taylor", "vortex_taylor", 0.03),
]
for vtype, mn, strength in _SCENES_VORTEX:
key = f"vortex_{vtype}"
SCENES[key] = {
"scene_id": "vortex",
"vortex_type": vtype,
"vortex_strength": strength,
"re_code": 100,
"mu": 2.0 / 100,
"nu": nu_from_re(100),
"has_disturbance": False,
"sample_interval": 800,
"max_steps": 150,
"action_scale": 4.0,
"action_bias": (0.0, -4.0, 4.0),
"source": "PPO_inference",
"model_name": mn,
"n_objects_env": 6,
"obs_slice": (0, 12),
"sensor_x": 40.0,
"pinball_front_x": 30.0,
"pinball_rear_x": 31.3,
"target_type": "transient",
"s_dim": 12,
"u0": U0,
}
# -- Utility helpers ---------------------------------------------------------
def get_scene(name: str) -> dict:
"""Return scene config dict by name. Raises KeyError if not found."""
if name not in SCENES:
raise KeyError(f"Unknown scene: {name}. Available: {list(SCENES.keys())}")
return dict(SCENES[name])
def get_scene_list(scene_id: Optional[str] = None) -> List[str]:
"""Return list of scene names, optionally filtered by scene_id."""
if scene_id is None:
return list(SCENES.keys())
return [k for k, v in SCENES.items() if v["scene_id"] == scene_id]
def model_path_for_scene(scene_name: str) -> Optional[str]:
"""Return absolute path to PPO model .zip file, or None."""
cfg = get_scene(scene_name)
mn = cfg.get("model_name")
if mn is None:
return None
# Check model directories in priority order
candidate_dirs = [
os.path.join(MODEL_DIR, "old"),
os.path.join(MODEL_DIR, "250525"),
os.path.join(MODEL_DIR, "250729"),
os.path.join(MODEL_DIR, "250326"),
]
for d in candidate_dirs:
p = os.path.join(d, f"{mn}.zip")
if os.path.isfile(p):
return p
return None
src/SR_analysis/configs/legacy/config_cuda.json
+9
View File
@@ -0,0 +1,9 @@
{
"multi_gpu": false,
"gpu_connection": "NVLink",
"required_cuda_capability": "7.0",
"threads_per_block": 128,
"X_1U": 128,
"Y_1U": 32,
"Z_1U": 1
}
src/SR_analysis/configs/legacy/config_flowfield.json
+13
View File
@@ -0,0 +1,13 @@
{
"data_type": "FP32",
"dimensionality": 2,
"lattice": 9,
"field_dim_in_U": [10, 16, 1],
"viscosity": 0.004,
"velocity": 0.01,
"boundary_conditions": {
"x": ["parabolic", "outflow"],
"y": ["noslip", "noslip"],
"z": ["none", "none"]
}
}
src/SR_analysis/data/illusion/illusion_0.75L/config.json
+21
View File
@@ -0,0 +1,21 @@
{
"scene_id": "illusion",
"target_diameter": 0.75,
"re_code": 100,
"mu": 0.02,
"nu": 0.004,
"has_disturbance": false,
"sample_interval": 600,
"action_scale": 8.0,
"action_bias": "(0.0, -2.0, 2.0)",
"source": "PPO_inference",
"model_name": "d1a3o12_250525_imit_075L_1U",
"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,
"u0": 0.01
}
src/SR_analysis/data/illusion/illusion_0.75L/norm.json
+24
View File
@@ -0,0 +1,24 @@
{
"force_norm_fact": 0.013476977124810219,
"sens_deviation": [
0.962617814540863,
-0.12039308249950409,
0.6415857672691345,
0.011103342287242413,
0.9339056611061096,
0.11935960501432419
],
"sens_norm_fact": [
2.0483264923095703,
2.5809006690979004,
0.7443606853485107,
3.4969263076782227,
2.1811583042144775,
2.5745153427124023
],
"action_bias": [
0.0,
-2.0,
2.0
]
}
src/SR_analysis/data/illusion/illusion_0.75L/result.json
+5
View File
@@ -0,0 +1,5 @@
{
"scene": "illusion_0.75L",
"controlled": true,
"similarity": 0.18393706196948187
}
src/SR_analysis/data/illusion/illusion_1.5L/config.json
+21
View File
@@ -0,0 +1,21 @@
{
"scene_id": "illusion",
"target_diameter": 1.5,
"re_code": 100,
"mu": 0.02,
"nu": 0.008,
"has_disturbance": false,
"sample_interval": 600,
"action_scale": 8.0,
"action_bias": "(0.0, -2.0, 2.0)",
"source": "PPO_inference",
"model_name": "d1a3o14_250525_imit_15L_2U",
"n_objects_env": 6,
"obs_slice": "(0, 12)",
"sensor_x": 30.0,
"pinball_front_x": 19.0,
"pinball_rear_x": 20.3,
"target_type": "periodic",
"s_dim": 14,
"u0": 0.02
}
src/SR_analysis/data/illusion/illusion_1.5L/norm.json
+24
View File
@@ -0,0 +1,24 @@
{
"force_norm_fact": 0.054184265434741974,
"sens_deviation": [
1.9146838188171387,
-0.23843951523303986,
1.305143117904663,
0.0009254463366232812,
1.8709181547164917,
0.2255263477563858
],
"sens_norm_fact": [
4.165492057800293,
5.171088695526123,
1.5217405557632446,
6.904999732971191,
4.384937763214111,
5.106513023376465
],
"action_bias": [
0.0,
-2.0,
2.0
]
}
src/SR_analysis/data/illusion/illusion_1.5L/result.json
+5
View File
@@ -0,0 +1,5 @@
{
"scene": "illusion_15L",
"controlled": true,
"similarity": 0.30179914651024675
}
src/SR_analysis/data/illusion/illusion_1L/config.json
+21
View File
@@ -0,0 +1,21 @@
{
"scene_id": "illusion",
"target_diameter": 1.0,
"re_code": 100,
"mu": 0.02,
"nu": 0.004,
"has_disturbance": false,
"sample_interval": 600,
"action_scale": 8.0,
"action_bias": "(0.0, -2.0, 2.0)",
"source": "PPO_inference",
"model_name": "d1a3o12_250525_imit_1L_1U",
"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,
"u0": 0.01
}
src/SR_analysis/data/illusion/illusion_1L/norm.json
+24
View File
@@ -0,0 +1,24 @@
{
"force_norm_fact": 0.013487594202160835,
"sens_deviation": [
0.9391862154006958,
-0.09573575109243393,
0.6525679230690002,
0.03420928493142128,
0.949753999710083,
0.13210755586624146
],
"sens_norm_fact": [
2.1654911041259766,
2.459106683731079,
0.7431825995445251,
3.613541603088379,
2.1102607250213623,
2.6434059143066406
],
"action_bias": [
0.0,
-2.0,
2.0
]
}
src/SR_analysis/data/illusion/illusion_1L/result.json
+5
View File
@@ -0,0 +1,5 @@
{
"scene": "illusion_1.0L",
"controlled": true,
"similarity": 0.5543174409436081
}
src/SR_analysis/data/karman/karman_re100/config.json
+20
View File
@@ -0,0 +1,20 @@
{
"scene_id": "karman",
"re_code": 100,
"mu": 0.02,
"nu": 0.004,
"has_disturbance": true,
"sample_interval": 800,
"action_scale": 8.0,
"action_bias": "(0.0, -4.0, 4.0)",
"source": "PPO_inference",
"model_name": "d1a3o12_re100",
"n_objects_env": 7,
"obs_slice": "(2, 14)",
"sensor_x": 40.0,
"pinball_front_x": 30.0,
"pinball_rear_x": 31.3,
"target_type": "periodic",
"s_dim": 12,
"u0": 0.01
}
src/SR_analysis/data/karman/karman_re100/norm.json
+24
View File
@@ -0,0 +1,24 @@
{
"force_norm_fact": 0.019199129194021225,
"sens_deviation": [
0.8231719732284546,
-0.12661591172218323,
0.24832786619663239,
-0.01064519677311182,
0.7844515442848206,
0.1161285787820816
],
"sens_norm_fact": [
3.3014419078826904,
3.2062995433807373,
1.8544995784759521,
3.4928226470947266,
3.1099960803985596,
2.815072774887085
],
"action_bias": [
0.0,
-4.0,
4.0
]
}
src/SR_analysis/data/karman/karman_re100/result.json
+6
View File
@@ -0,0 +1,6 @@
{
"scene": "karman_re100",
"controlled": true,
"avg_reward_last100": 0.6665451352155793,
"similarity": 0.9538050162761162
}
src/SR_analysis/data/karman/karman_re100/vorticity_controlled.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 206 KiB

src/SR_analysis/data/karman/karman_re100/vorticity_uncontrolled.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 266 KiB

src/SR_analysis/data/karman/karman_re200/config.json
+20
View File
@@ -0,0 +1,20 @@
{
"scene_id": "karman",
"re_code": 200,
"mu": 0.01,
"nu": 0.002,
"has_disturbance": true,
"sample_interval": 800,
"action_scale": 8.0,
"action_bias": "(0.0, -4.0, 4.0)",
"source": "PPO_inference",
"model_name": "d1a3o12_re200",
"n_objects_env": 7,
"obs_slice": "(2, 14)",
"sensor_x": 40.0,
"pinball_front_x": 30.0,
"pinball_rear_x": 31.3,
"target_type": "periodic",
"s_dim": 12,
"u0": 0.01
}
src/SR_analysis/data/karman/karman_re200/norm.json
+24
View File
@@ -0,0 +1,24 @@
{
"force_norm_fact": 0.02405486349016428,
"sens_deviation": [
0.761349618434906,
-0.10393908619880676,
-0.0060332342982292175,
-0.01062991376966238,
0.7603892087936401,
0.08925710618495941
],
"sens_norm_fact": [
2.458379030227661,
2.4950430393218994,
0.986889123916626,
2.259662389755249,
2.737121820449829,
2.521576404571533
],
"action_bias": [
0.0,
-4.0,
4.0
]
}
src/SR_analysis/data/karman/karman_re200/result.json
+6
View File
@@ -0,0 +1,6 @@
{
"scene": "karman_re200",
"controlled": true,
"avg_reward_last100": 0.3298547239579881,
"similarity": 0.8842106151498026
}
src/SR_analysis/data/karman/karman_re200/vorticity_controlled.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 232 KiB

src/SR_analysis/data/karman/karman_re200/vorticity_uncontrolled.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 270 KiB

src/SR_analysis/data/karman/karman_re400/config.json
+10
View File
@@ -0,0 +1,10 @@
{
"re_code": 400,
"nu": 0.001,
"u0": 0.01,
"sample_interval": 800,
"fifo_len": 150,
"conv_len": 30,
"device_id": 2,
"model_path": "/home/frank14f/DynamisLab/models/old/d1a3o12_re400.zip"
}
src/SR_analysis/data/karman/karman_re400/norm.json
+24
View File
@@ -0,0 +1,24 @@
{
"force_norm_fact": 0.030264523811638355,
"sens_deviation": [
0.8747888207435608,
-0.024021463468670845,
0.5912007689476013,
0.017280835658311844,
0.9475194215774536,
0.07682034373283386
],
"sens_norm_fact": [
5.08115291595459,
5.131664276123047,
3.3446834087371826,
5.320921897888184,
4.4046711921691895,
5.202882766723633
],
"action_bias": [
0.0,
-4.0,
4.0
]
}
src/SR_analysis/data/karman/karman_re400/result.json
+7
View File
@@ -0,0 +1,7 @@
{
"re_code": 400,
"uncontrolled": true,
"controlled": true,
"avg_reward_last100": 0.4389868174760177,
"similarity": 0.7950085552241137
}
src/SR_analysis/data/karman/karman_re400/vorticity_controlled.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 250 KiB

src/SR_analysis/data/karman/karman_re400/vorticity_uncontrolled.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 274 KiB

src/SR_analysis/data/karman/karman_re50/config.json
+20
View File
@@ -0,0 +1,20 @@
{
"scene_id": "karman",
"re_code": 50,
"mu": 0.04,
"nu": 0.008,
"has_disturbance": true,
"sample_interval": 800,
"action_scale": 8.0,
"action_bias": "(0.0, -4.0, 4.0)",
"source": "PPO_inference",
"model_name": "d1a3o12_re50",
"n_objects_env": 7,
"obs_slice": "(2, 14)",
"sensor_x": 40.0,
"pinball_front_x": 30.0,
"pinball_rear_x": 31.3,
"target_type": "periodic",
"s_dim": 12,
"u0": 0.01
}
src/SR_analysis/data/karman/karman_re50/norm.json
+24
View File
@@ -0,0 +1,24 @@
{
"force_norm_fact": 0.015911692287772894,
"sens_deviation": [
0.6078279614448547,
-0.04162348061800003,
0.07135022431612015,
-0.0008823801181279123,
0.6027681827545166,
0.0418982058763504
],
"sens_norm_fact": [
0.789494514465332,
1.1795930862426758,
0.18662318587303162,
1.1806247234344482,
0.8472481369972229,
1.2091511487960815
],
"action_bias": [
0.0,
-4.0,
4.0
]
}
src/SR_analysis/data/karman/karman_re50/result.json
+6
View File
@@ -0,0 +1,6 @@
{
"scene": "karman_re50",
"controlled": true,
"avg_reward_last100": 0.5020737935656798,
"similarity": 0.9614716421942122
}
src/SR_analysis/data/karman/karman_re50/vorticity_controlled.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 188 KiB

src/SR_analysis/data/karman/karman_re50/vorticity_uncontrolled.png
BIN
View File
Binary file not shown.
Side by Side

After

Width:  |  Height:  |  Size: 264 KiB

src/SR_analysis/data/vortex/vortex_lamb/config.json
+23
View File
@@ -0,0 +1,23 @@
{
"scene_id": "vortex",
"vortex_type": "lamb",
"vortex_strength": 0.5,
"re_code": 100,
"mu": 0.02,
"nu": 0.004,
"has_disturbance": false,
"sample_interval": 800,
"max_steps": 150,
"action_scale": 4.0,
"action_bias": "(0.0, -4.0, 4.0)",
"source": "PPO_inference",
"model_name": "vortex_lamb",
"n_objects_env": 6,
"obs_slice": "(0, 12)",
"sensor_x": 40.0,
"pinball_front_x": 30.0,
"pinball_rear_x": 31.3,
"target_type": "transient",
"s_dim": 12,
"u0": 0.01
}
src/SR_analysis/data/vortex/vortex_lamb/norm.json
+24
View File
@@ -0,0 +1,24 @@
{
"force_norm_fact": 0.032297031953930855,
"sens_deviation": [
0.8963788747787476,
-0.10795877873897552,
-0.014931724406778812,
1.4363668924488593e-05,
0.8963659405708313,
0.10797087103128433
],
"sens_norm_fact": [
1.4594829082489014,
1.0411686897277832,
5.520665168762207,
0.0022848353255540133,
1.4571585655212402,
1.0411614179611206
],
"action_bias": [
0.0,
-4.0,
4.0
]
}
src/SR_analysis/data/vortex/vortex_lamb/result.json
+5
View File
@@ -0,0 +1,5 @@
{
"scene": "vortex_lamb",
"controlled": true,
"similarity": 0.9421189523977421
}
src/SR_analysis/data/vortex/vortex_taylor/config.json
+23
View File
@@ -0,0 +1,23 @@
{
"scene_id": "vortex",
"vortex_type": "taylor",
"vortex_strength": 0.03,
"re_code": 100,
"mu": 0.02,
"nu": 0.004,
"has_disturbance": false,
"sample_interval": 800,
"max_steps": 150,
"action_scale": 4.0,
"action_bias": "(0.0, -4.0, 4.0)",
"source": "PPO_inference",
"model_name": "vortex_taylor",
"n_objects_env": 6,
"obs_slice": "(0, 12)",
"sensor_x": 40.0,
"pinball_front_x": 30.0,
"pinball_rear_x": 31.3,
"target_type": "transient",
"s_dim": 12,
"u0": 0.01
}
src/SR_analysis/data/vortex/vortex_taylor/norm.json
+24
View File
@@ -0,0 +1,24 @@
{
"force_norm_fact": 0.03310442715883255,
"sens_deviation": [
0.9501951336860657,
-0.15176598727703094,
0.49544230103492737,
-0.009104072116315365,
1.0117771625518799,
0.1462564468383789
],
"sens_norm_fact": [
4.190938949584961,
2.9906840324401855,
3.5891337394714355,
4.104613780975342,
3.012289047241211,
3.3520655632019043
],
"action_bias": [
0.0,
-4.0,
4.0
]
}
src/SR_analysis/data/vortex/vortex_taylor/result.json
+5
View File
@@ -0,0 +1,5 @@
{
"scene": "vortex_taylor",
"controlled": true,
"similarity": 0.9158487825490536
}
src/SR_analysis/pysr.md
+842
View File
@@ -0,0 +1,842 @@
# 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/scripts/infer_illusion.py
+309
View File
@@ -0,0 +1,309 @@
"""Inference pipeline for Illusion (cylinder imitation) scenes.
Generates controlled data for a given target cylinder diameter using
LegacyCelerisLab + trained PPO model.
Usage:
conda run -n pycuda_3_10 python scripts/infer_illusion.py \\
--diameter 1.0 --device 0
conda run -n pycuda_3_10 python scripts/infer_illusion.py \\
--diameter all --device 2
"""
from __future__ import annotations
import argparse
import json
import os
import sys
import time
from collections import deque
from typing import Optional
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
if _REPO not in sys.path:
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from LegacyCelerisLab import FlowField # noqa: E402
from LegacyCelerisLab import utils as legacy_utils # noqa: E402
from SR_analysis.utils.cfd_interface import (
load_legacy_configs, build_observation,
scale_action, load_ppo_model, compute_similarity,
)
from SR_analysis.configs import (
get_scene, get_scene_list, model_path_for_scene,
LEGACY_CFG_DIR, FIFO_LEN, CONV_LEN,
)
DATA_TYPE = np.float32
def run_single_illusion(
scene_name: str,
device_id: int,
output_root: str,
n_infer_steps: int = 200,
) -> dict:
"""Run full inference pipeline for one Illusion scene."""
cfg = get_scene(scene_name)
nu = cfg["nu"]
u0 = cfg["u0"]
l0 = 20.0
sample_interval = cfg["sample_interval"]
action_scale = cfg["action_scale"]
action_bias = cfg["action_bias"]
n_obj_total = cfg["n_objects_env"]
sensor_x = cfg["sensor_x"] # 30.0 for illusion
front_x = cfg["pinball_front_x"] # 19.0
rear_x = cfg["pinball_rear_x"] # 20.3
target_diam = cfg["target_diameter"]
os.makedirs(output_root, exist_ok=True)
print(f"\n{'='*60}")
print(f"Scene: {scene_name} Diam={target_diam}L u0={u0} device={device_id}")
print(f"{'='*60}")
# Save config
with open(os.path.join(output_root, "config.json"), "w") as f:
json.dump({k: str(v) if not isinstance(v, (int, float, list, bool))
else v for k, v in cfg.items()}, f, indent=2)
# Load legacy CFD configs with overridden viscosity and velocity
cuda_cfg, field_cfg = load_legacy_configs(LEGACY_CFG_DIR)
field_cfg = field_cfg._replace(viscosity=float(nu))
if u0 != 0.01:
field_cfg = field_cfg._replace(velocity=float(u0))
# -- Phase 1: Target recording (target cylinder + 3 sensors) ------------
ff = FlowField(field_cfg, cuda_cfg, device_id=device_id)
ny = ff.FIELD_SHAPE[1]
# Add target cylinder at x=20*L0, with radius = target_diam * L0
print(f" Adding target cylinder: diam={target_diam}L, pos=({20*l0:.0f}, {ny/2:.0f})")
ff.add_cylinder((20.0 * l0, (ny - 1) / 2, 0.0), target_diam * l0)
# Add 3 sensors at x = sensor_x * L0
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)
n_obj_phase1 = ff.obs.size // 2
print(f" Phase 1 objects: {n_obj_phase1}")
# Stabilize
stabilize_steps = int(4 * ff.FIELD_SHAPE[0] / u0)
print(f" Stabilising ({stabilize_steps} steps)...")
ff.run(stabilize_steps, np.zeros(n_obj_phase1, dtype=DATA_TYPE))
# Record target
target_states = np.empty((0, 8), dtype=DATA_TYPE)
for _ in range(FIFO_LEN):
ff.run(sample_interval, np.zeros(n_obj_phase1, dtype=DATA_TYPE))
new_state = ff.obs.copy()[0:8] # sensor[6] + cylinder force[2]
target_states = np.vstack((target_states, new_state))
print(f" Target recorded: {target_states.shape}")
# Save target
np.savez(os.path.join(output_root, "target.npz"), target_states=target_states)
# Clean up and create pinball env
del ff
ff = FlowField(field_cfg, cuda_cfg, device_id=device_id)
# -- Phase 2: Pinball env -----------------------------------------------
# Add 3 sensors (same positions as target phase)
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)
# Add 3 pinball cylinders (illusion positions)
# Front at x=front_x*L0, rear at x=rear_x*L0
ff.add_cylinder((front_x * l0, (ny - 1) / 2, 0.0), l0 / 2.0)
ff.add_cylinder((rear_x * l0, (ny - 1) / 2 + 0.75 * l0, 0.0), l0 / 2.0)
ff.add_cylinder((rear_x * l0, (ny - 1) / 2 - 0.75 * l0, 0.0), l0 / 2.0)
n_obj = ff.obs.size // 2
print(f" Pinball env objects: {n_obj}")
assert n_obj == 6, f"Expected 6 objects, got {n_obj}"
# Stabilize with zero action
print(f" Stabilising pinball ({stabilize_steps} steps)...")
ff.run(stabilize_steps, np.zeros(n_obj, dtype=DATA_TYPE))
# Checkpoint
ff.get_ddf()
ff.save_ddf()
# Norm collection (zero action)
fifo = deque(maxlen=FIFO_LEN)
for _ in range(FIFO_LEN):
ff.run(sample_interval, np.zeros(n_obj, dtype=DATA_TYPE))
fifo.append(ff.obs.copy()[0:12])
temp_states = np.array(fifo, dtype=DATA_TYPE)
force_norm_fact = 6.0 * float(np.max(np.abs(temp_states[:, 6:12])))
sens_deviation = np.mean(temp_states[:, 0:6], axis=0).astype(DATA_TYPE)
sens_norm_fact = np.zeros(6, dtype=DATA_TYPE)
for i in range(6):
sens_norm_fact[i] = 5.0 * float(np.max(np.abs(temp_states[:, i] - sens_deviation[i])))
norm = {
"force_norm_fact": force_norm_fact,
"sens_deviation": sens_deviation.tolist(),
"sens_norm_fact": sens_norm_fact.tolist(),
"action_bias": list(action_bias),
}
print(f" norm: force_norm_fact={force_norm_fact:.6f}")
# Bias-action rollout
ff.apply_ddf()
bias_arr = np.zeros(n_obj, dtype=DATA_TYPE)
bias_arr[3] = float(action_bias[0] * u0)
bias_arr[4] = float(action_bias[1] * u0)
bias_arr[5] = float(action_bias[2] * u0)
print(f" bias action: {bias_arr}")
fifo.clear()
for _ in range(FIFO_LEN):
ff.run(sample_interval, bias_arr)
fifo.append(ff.obs.copy()[0:12])
save_states = np.array(list(fifo), dtype=DATA_TYPE)
norm["save_states"] = save_states
ff.apply_ddf()
# Save norm
norm_json = {k: v for k, v in norm.items() if not isinstance(v, np.ndarray)}
with open(os.path.join(output_root, "norm.json"), "w") as f:
json.dump(norm_json, f, indent=2)
# -- Phase 3: Controlled inference ---------------------------------------
result = {"scene": scene_name, "controlled": False}
model_path = model_path_for_scene(scene_name)
if model_path is not None:
s_dim = cfg.get("s_dim", 12)
print(f" loading model: {model_path} (s_dim={s_dim})")
model = load_ppo_model(model_path, device=f"cuda:{device_id}", s_dim=s_dim)
model.set_random_seed(0)
print(f" controlled rollout ({n_infer_steps} steps) ...")
ff.restore_ddf()
ff.apply_ddf()
# Re-bias FIFO
fifo = deque(maxlen=FIFO_LEN)
for _ in range(FIFO_LEN):
ff.context.push()
ff.run(sample_interval, bias_arr)
ff.context.pop()
fifo.append(ff.obs.copy()[0:12])
sens_list, forc_list, action_list = [], [], []
obs = np.zeros(s_dim, dtype=np.float32)
for step in range(n_infer_steps):
action, _states = model.predict(obs, deterministic=True)
action = action.astype(np.float32).flatten()
action_list.append(action.copy())
# Convert to legacy action array (6 objects: sensors[3] + pinball[3])
temp = np.zeros(n_obj, dtype=DATA_TYPE)
temp[3:6] = np.array(
(action * action_scale + list(action_bias)) * u0,
dtype=DATA_TYPE)
ff.context.push()
ff.run(sample_interval, temp)
ff.context.pop()
obs_slice = ff.obs.copy()[0:12]
fifo.append(obs_slice)
sens_list.append(obs_slice[0:6])
forc_list.append(obs_slice[6:12])
# Build normalized obs (just forces + sens for S_DIM=12)
forces_norm = obs_slice[6:12] / force_norm_fact
sens_norm = (obs_slice[0:6] - sens_deviation) / sens_norm_fact
obs12 = np.clip(np.hstack([forces_norm, sens_norm]), -1.0, 1.0).astype(np.float32)
if s_dim == 14:
# Need target values -- for inference we zero-pad
obs = np.zeros(14, dtype=np.float32)
obs[:12] = obs12
else:
obs = obs12
np.savez(os.path.join(output_root, "controlled.npz"),
sensors=np.array(sens_list, dtype=np.float32),
forces=np.array(forc_list, dtype=np.float32),
actions=np.array(action_list, dtype=np.float32))
# Compute similarity (use the target cylinder's sensor-only signals)
# For comparison, compute similarity between controlled sensors and target
target_sensors = target_states[:, 0:6]
sim = compute_similarity(target_sensors,
np.array(sens_list, dtype=np.float32), CONV_LEN)
print(f" similarity (vs target cylinder) = {sim:.4f}")
result["controlled"] = True
result["similarity"] = sim
else:
print(f" WARNING: no model for {scene_name}")
del ff
with open(os.path.join(output_root, "result.json"), "w") as f:
json.dump(result, f, indent=2)
return result
def main():
ap = argparse.ArgumentParser(description="Illusion inference")
ap.add_argument("--diameter", type=str, default="1.0",
help='Diameter: 0.75, 1.0, 1.5, or "all"')
ap.add_argument("--device", type=int, default=0, help="GPU device ID")
ap.add_argument("--steps", type=int, default=200)
ap.add_argument("--out-root", type=str, default=None)
args = ap.parse_args()
# Determine scene names
if args.diameter.lower() == "all":
scene_names = get_scene_list("illusion")
else:
d = float(args.diameter)
# Match by target_diameter field
scene_names = []
for sn in get_scene_list("illusion"):
cfg = get_scene(sn)
if abs(cfg["target_diameter"] - d) < 0.01:
scene_names.append(sn)
if not scene_names:
print(f"ERROR: no illusion scene found for diameter={d}")
return 1
if args.out_root is None:
out_root = os.path.join(os.path.dirname(__file__), "..", "data", "illusion")
else:
out_root = args.out_root
t_start = time.time()
for sn in scene_names:
case_dir = os.path.join(out_root, sn)
result = run_single_illusion(sn, args.device, case_dir,
n_infer_steps=args.steps)
print(f" Done: {sn} -> {case_dir}")
elapsed = time.time() - t_start
print(f"\nTotal time: {elapsed:.1f}s")
if __name__ == "__main__":
main()
src/SR_analysis/scripts/infer_karman.py
+253
View File
@@ -0,0 +1,253 @@
"""Inference pipeline for Karman cloak across Re.
Generates controlled/uncontrolled data for a given Re case using
LegacyCelerisLab + trained PPO model.
Usage:
conda run -n pycuda_3_10 python scripts/infer_karman.py --re 100 --device 0
conda run -n pycuda_3_10 python scripts/infer_karman.py --re all --device 2
"""
from __future__ import annotations
import argparse
import json
import os
import sys
import time
from collections import deque
import numpy as np
# Add repo root for LegacyCelerisLab and src/ for SR_analysis
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
if _REPO not in sys.path:
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from LegacyCelerisLab import FlowField # noqa: E402
from SR_analysis.utils.cfd_interface import (
nu_from_re, load_legacy_configs,
build_karman_cloak_env, add_pinball, build_observation,
scale_action, load_ppo_model, save_vorticity_png,
vorticity_from_ddf, compute_similarity, ACTION_SMOOTH_WEIGHT,
)
from SR_analysis.configs import (
SCENES, get_scene, get_scene_list, model_path_for_scene,
LEGACY_CFG_DIR, FIFO_LEN, CONV_LEN,
)
DATA_TYPE = np.float32
def run_single_re(
scene_name: str,
device_id: int,
output_root: str,
n_infer_steps: int = 200,
) -> dict:
"""Run full inference pipeline for one Karman Re case."""
cfg = get_scene(scene_name)
re_code = cfg["re_code"]
nu = cfg["nu"]
u0 = cfg["u0"]
l0 = 20.0
sample_interval = cfg["sample_interval"]
action_scale = cfg["action_scale"]
action_bias = cfg["action_bias"]
n_obj_total = cfg["n_objects_env"]
os.makedirs(output_root, exist_ok=True)
print(f"\n{'='*60}")
print(f"Scene: {scene_name} Re_code={re_code} nu={nu:.6f} device={device_id}")
print(f"{'='*60}")
# Save config
with open(os.path.join(output_root, "config.json"), "w") as f:
json.dump({k: str(v) if not isinstance(v, (int, float, list, bool))
else v for k, v in cfg.items()}, f, indent=2)
# Load legacy CFD configs with overridden viscosity
cuda_cfg, field_cfg = load_legacy_configs(LEGACY_CFG_DIR)
field_cfg = field_cfg._replace(viscosity=float(nu))
# Build env: dist cylinder + sensors, record target
ff = FlowField(field_cfg, cuda_cfg, device_id=device_id)
target_states, env_info = build_karman_cloak_env(
ff, u0=u0, l0=l0, sample_interval=sample_interval,
fifo_len=FIFO_LEN, data_type=DATA_TYPE,
)
np.savez(os.path.join(output_root, "target.npz"), target_states=target_states)
# Add pinball, compute norm
norm = add_pinball(
ff, l0=l0, u0=u0, sample_interval=sample_interval,
fifo_len=FIFO_LEN, data_type=DATA_TYPE,
action_bias=action_bias, pinball_front_x=cfg["pinball_front_x"],
pinball_rear_x=cfg["pinball_rear_x"],
obs_slice_start=cfg["obs_slice"][0], obs_slice_end=cfg["obs_slice"][1],
)
norm_for_json = {k: v for k, v in norm.items()
if not isinstance(v, np.ndarray)}
with open(os.path.join(output_root, "norm.json"), "w") as f:
json.dump(norm_for_json, f, indent=2)
# Uncontrolled rollout
print(" uncontrolled rollout ...")
ff.restore_ddf()
ff.apply_ddf()
sens_list, forc_list = [], []
for _ in range(n_infer_steps):
ff.run(sample_interval, np.zeros(n_obj_total, dtype=DATA_TYPE))
obs_slice = ff.obs.copy()[2:14]
sens_list.append(obs_slice[0:6])
forc_list.append(obs_slice[6:12])
np.savez(os.path.join(output_root, "uncontrolled.npz"),
sensors=np.array(sens_list, dtype=np.float32),
forces=np.array(forc_list, dtype=np.float32))
omega_unc = vorticity_from_ddf(ff, u0=u0)
save_vorticity_png(os.path.join(output_root, "vorticity_uncontrolled.png"),
omega_unc, title=f"{scene_name} uncontrolled")
# Controlled rollout
result = {"scene": scene_name, "controlled": False}
model_path = model_path_for_scene(scene_name)
if model_path is not None:
s_dim = cfg.get("s_dim", 12)
print(f" loading model: {model_path} (s_dim={s_dim})")
model = load_ppo_model(model_path, device=f"cuda:{device_id}", s_dim=s_dim)
model.set_random_seed(0)
print(f" controlled rollout ({n_infer_steps} steps) ...")
ff.restore_ddf()
ff.apply_ddf()
# Bias FIFO init
fifo = deque(maxlen=FIFO_LEN)
bias_action = scale_action(
np.zeros(3, dtype=np.float32),
scale=action_scale, bias=action_bias, u0=u0,
n_total_bodies=n_obj_total,
)
for _ in range(FIFO_LEN):
ff.context.push()
ff.run(sample_interval, bias_action)
ff.context.pop()
fifo.append(ff.obs.copy()[2:14])
sens_list_c, forc_list_c, action_list_c = [], [], []
reward_list_c = []
obs = np.zeros(12, dtype=np.float32)
for step in range(n_infer_steps):
action, _states = model.predict(obs, deterministic=True)
action = action.astype(np.float32).flatten()
action_list_c.append(action.copy())
action_arr = scale_action(
action, scale=action_scale, bias=action_bias,
u0=u0, n_total_bodies=n_obj_total,
)
ff.context.push()
ff.run(sample_interval, action_arr)
ff.context.pop()
obs_slice = ff.obs.copy()[2:14]
fifo.append(obs_slice)
sens_list_c.append(obs_slice[0:6])
forc_list_c.append(obs_slice[6:12])
obs = build_observation(obs_slice, norm)
# Compute reward
states_arr = np.array(list(fifo), dtype=np.float32)
if len(states_arr) >= CONV_LEN:
forces = states_arr[-1, 6:12] / norm["force_norm_fact"]
cd = float((forces[0] + forces[2] + forces[4]) / 3.0)
cl = float((forces[1] + forces[3] + forces[5]) / 3.0)
sim = compute_similarity(target_states, states_arr[:, 0:6], CONV_LEN)
r_cd = np.exp(-abs(cd * 20.0))
r_cl = np.exp(-abs(cl * 80.0))
r_sim = np.exp(-10.0 * abs(sim - 1.0))
reward = min(0.3 * r_cd + 0.4 * r_cl + 0.3 * r_sim, 1.0)
reward_list_c.append(float(reward))
np.savez(os.path.join(output_root, "controlled.npz"),
sensors=np.array(sens_list_c, dtype=np.float32),
forces=np.array(forc_list_c, dtype=np.float32),
actions=np.array(action_list_c, dtype=np.float32),
rewards=np.array(reward_list_c, dtype=np.float32))
omega_con = vorticity_from_ddf(ff, u0=u0)
save_vorticity_png(os.path.join(output_root, "vorticity_controlled.png"),
omega_con, title=f"{scene_name} controlled")
avg_reward = (float(np.mean(reward_list_c[-100:]))
if len(reward_list_c) >= 100
else float(np.mean(reward_list_c)))
sim_score = compute_similarity(
target_states, np.array(sens_list_c, dtype=np.float32), CONV_LEN)
result["controlled"] = True
result["avg_reward_last100"] = avg_reward
result["similarity"] = sim_score
print(f" avg_reward(last100)={avg_reward:.4f} similarity={sim_score:.4f}")
else:
print(f" no model for {scene_name}, skipping controlled rollout")
del ff
with open(os.path.join(output_root, "result.json"), "w") as f:
json.dump(result, f, indent=2)
return result
def main():
ap = argparse.ArgumentParser(description="Karman cloak inference")
ap.add_argument("--re", type=str, default="100",
help='Re case: 50,100,200,400, or "all", or "validation"')
ap.add_argument("--device", type=int, default=0, help="GPU device ID")
ap.add_argument("--steps", type=int, default=200,
help="Number of inference steps per rollout")
ap.add_argument("--out-root", type=str, default=None,
help="Output root (default: SR_analysis/data/karman)")
args = ap.parse_args()
# Determine scene names
selection = args.re.lower()
if selection == "all":
scene_names = get_scene_list("karman")
elif selection == "validation":
scene_names = [f"karman_re{rc}" for rc in [35, 70, 150]]
# Check which are defined
scene_names = [s for s in scene_names if s in SCENES]
else:
rc = int(selection)
scene_names = [f"karman_re{rc}"]
if args.out_root is None:
out_root = os.path.join(os.path.dirname(__file__), "..", "data", "karman")
else:
out_root = args.out_root
t_start = time.time()
for sn in scene_names:
case_dir = os.path.join(out_root, sn)
result = run_single_re(sn, args.device, case_dir,
n_infer_steps=args.steps)
print(f" Done: {sn} -> {case_dir}")
elapsed = time.time() - t_start
print(f"\nTotal time: {elapsed:.1f}s")
if __name__ == "__main__":
main()
src/SR_analysis/scripts/infer_vortex.py
+300
View File
@@ -0,0 +1,300 @@
"""Inference pipeline for Vortex cloak (Lamb dipole + Taylor monopole).
Generates controlled data for vortex cloak scenes using
LegacyCelerisLab + trained PPO model.
Vortex env characteristics:
- No disturbance cylinder
- Vortex is added AFTER DDF checkpoint
- Transient: MAX_STEPS=150
- Action scaling: action*4 + [0,-4,4]
Usage:
conda run -n pycuda_3_10 python scripts/infer_vortex.py \\
--type lamb --device 0
conda run -n pycuda_3_10 python scripts/infer_vortex.py \\
--type all --device 2
"""
from __future__ import annotations
import argparse
import json
import os
import sys
import time
from collections import deque
from typing import Optional
import numpy as np
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
if _REPO not in sys.path:
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from LegacyCelerisLab import FlowField # noqa: E402
from SR_analysis.utils.cfd_interface import (
load_legacy_configs, build_observation,
scale_action, load_ppo_model, compute_similarity,
)
from SR_analysis.configs import (
get_scene, get_scene_list, model_path_for_scene,
LEGACY_CFG_DIR, FIFO_LEN, CONV_LEN,
)
DATA_TYPE = np.float32
def run_single_vortex(
scene_name: str,
device_id: int,
output_root: str,
n_infer_steps: Optional[int] = None,
) -> dict:
"""Run full inference pipeline for one Vortex scene."""
cfg = get_scene(scene_name)
nu = cfg["nu"]
u0 = cfg["u0"]
l0 = 20.0
sample_interval = cfg["sample_interval"]
action_scale = cfg["action_scale"]
action_bias = cfg["action_bias"]
n_obj_pinball = cfg["n_objects_env"]
max_steps = cfg.get("max_steps", 150)
vtype = cfg["vortex_type"]
vstrength = cfg["vortex_strength"]
if n_infer_steps is None:
n_infer_steps = max_steps # transient -- use full episode
os.makedirs(output_root, exist_ok=True)
print(f"\n{'='*60}")
print(f"Scene: {scene_name} Vortex={vtype} strength={vstrength} device={device_id}")
print(f"{'='*60}")
# Save config
with open(os.path.join(output_root, "config.json"), "w") as f:
json.dump({k: str(v) if not isinstance(v, (int, float, list, bool))
else v for k, v in cfg.items()}, f, indent=2)
# Load legacy CFD configs
cuda_cfg, field_cfg = load_legacy_configs(LEGACY_CFG_DIR)
field_cfg = field_cfg._replace(viscosity=float(nu))
ff = FlowField(field_cfg, cuda_cfg, device_id=device_id)
ny = ff.FIELD_SHAPE[1]
# -- Phase 1: Sensors only env, record target with vortex -----------------
# Add 3 sensors at x=40*L0
for y_off in [2.0, 0.0, -2.0]:
sc = (40.0 * l0, (ny - 1) / 2 + y_off * l0, 0.0)
ff.add_sensor(sc, l0 / 4.0)
n_obj_sensors = ff.obs.size // 2
print(f" Sensor-only objects: {n_obj_sensors}")
# Short stabilize (1*NX/U0 instead of 4* for vortex)
stabilize_steps_short = int(1 * ff.FIELD_SHAPE[0] / u0)
ff.run(stabilize_steps_short, np.zeros(n_obj_sensors, dtype=DATA_TYPE))
# Save clean flow DDF (for later restore)
ff.get_ddf()
ff.save_ddf()
# Add vortex
print(f" Adding vortex: type={vtype}, center=({10*l0:.0f}, {ny/2:.0f})")
ff.add_vortex((10.0 * l0, (ny - 1) / 2, 0.0),
2.0 * l0, vstrength * u0, 0, vtype)
# Record target (vortex evolving through sensor-only env)
target_states = np.empty((0, 6), dtype=DATA_TYPE)
for _ in range(max_steps):
ff.run(sample_interval, np.zeros(n_obj_sensors, dtype=DATA_TYPE))
target_states = np.vstack((target_states, ff.obs.copy()))
print(f" Target recorded: {target_states.shape}")
np.savez(os.path.join(output_root, "target.npz"), target_states=target_states)
# -- Phase 2: Restore clean flow, add pinball, add vortex, record norm ----
ff.restore_ddf()
ff.apply_ddf()
# Add 3 pinball cylinders
ff.add_cylinder((30.0 * l0, (ny - 1) / 2, 0.0), l0 / 2.0)
ff.add_cylinder((31.3 * l0, (ny - 1) / 2 + 0.75 * l0, 0.0), l0 / 2.0)
ff.add_cylinder((31.3 * l0, (ny - 1) / 2 - 0.75 * l0, 0.0), l0 / 2.0)
n_obj = ff.obs.size // 2
print(f" Pinball env objects: {n_obj}")
assert n_obj == 6, f"Expected 6, got {n_obj}"
# Stabilize with zero action
ff.run(stabilize_steps_short, np.zeros(n_obj, dtype=DATA_TYPE))
# Stabilize with bias action (following vortex env code)
bias_arr = np.zeros(n_obj, dtype=DATA_TYPE)
bias_arr[3] = float(action_bias[0] * u0)
bias_arr[4] = float(action_bias[1] * u0)
bias_arr[5] = float(action_bias[2] * u0)
ff.run(stabilize_steps_short, bias_arr)
# Add vortex at x=15*L0 (pinball env) and save DDF
print(f" Adding vortex for pinball env: type={vtype}")
ff.add_vortex((15.0 * l0, (ny - 1) / 2, 0.0),
2.0 * l0, vstrength * u0, 0, vtype)
# SAVE DDF AFTER vortex (vortex env reset restores this mid-transient state)
ff.get_ddf()
ff.save_ddf()
print(" DDF saved with vortex active")
# Norm collection (zero action)
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())
temp_states = np.array(fifo, dtype=DATA_TYPE)
force_norm_fact = 6.0 * float(np.max(np.abs(temp_states[:, 6:12])))
sens_deviation = np.mean(temp_states[:, 0:6], axis=0).astype(DATA_TYPE)
sens_norm_fact = np.zeros(6, dtype=DATA_TYPE)
for i in range(6):
sens_norm_fact[i] = 5.0 * float(np.max(np.abs(temp_states[:, i] - sens_deviation[i])))
norm = {
"force_norm_fact": force_norm_fact,
"sens_deviation": sens_deviation.tolist(),
"sens_norm_fact": sens_norm_fact.tolist(),
"action_bias": list(action_bias),
}
print(f" norm: force_norm_fact={force_norm_fact:.6f}")
# Bias rollout (restore before vortex was added, then add vortex again)
ff.apply_ddf() # restore to DDF with vortex active
fifo.clear()
for _ in range(FIFO_LEN):
ff.run(sample_interval, bias_arr)
fifo.append(ff.obs.copy())
save_states = np.array(list(fifo), dtype=DATA_TYPE)
norm["save_states"] = save_states
ff.apply_ddf()
norm_json = {k: v for k, v in norm.items() if not isinstance(v, np.ndarray)}
with open(os.path.join(output_root, "norm.json"), "w") as f:
json.dump(norm_json, f, indent=2)
# -- Phase 3: Controlled inference ----------------------------------------
result = {"scene": scene_name, "controlled": False}
model_path = model_path_for_scene(scene_name)
if model_path is not None:
print(f" loading model: {model_path}")
model = load_ppo_model(model_path, device=f"cuda:{device_id}", s_dim=12)
model.set_random_seed(0)
print(f" controlled rollout ({n_infer_steps} steps) ...")
ff.restore_ddf()
ff.apply_ddf()
# Bias FIFO init (restore with vortex active)
fifo = deque(maxlen=FIFO_LEN)
for _ in range(FIFO_LEN):
ff.context.push()
ff.run(sample_interval, bias_arr)
ff.context.pop()
fifo.append(ff.obs.copy())
sens_list, forc_list, action_list = [], [], []
obs = np.zeros(12, dtype=np.float32)
for step in range(n_infer_steps):
action, _states = model.predict(obs, deterministic=True)
action = action.astype(np.float32).flatten()
action_list.append(action.copy())
# Action: action*4 + [0,-4,4]
temp = np.zeros(n_obj, dtype=DATA_TYPE)
temp[3:6] = np.array(
(action * action_scale + list(action_bias)) * u0,
dtype=DATA_TYPE)
ff.context.push()
ff.run(sample_interval, temp)
ff.context.pop()
obs_slice = ff.obs.copy()
fifo.append(obs_slice)
sens_list.append(obs_slice[0:6])
forc_list.append(obs_slice[6:12])
forces_norm = obs_slice[6:12] / force_norm_fact
sens_norm = (obs_slice[0:6] - sens_deviation) / sens_norm_fact
obs = np.clip(np.hstack([forces_norm, sens_norm]), -1.0, 1.0).astype(np.float32)
np.savez(os.path.join(output_root, "controlled.npz"),
sensors=np.array(sens_list, dtype=np.float32),
forces=np.array(forc_list, dtype=np.float32),
actions=np.array(action_list, dtype=np.float32))
# Similarity: align by step index (no lag for transient)
states_arr = np.array(sens_list, dtype=np.float32)
target_arr = target_states[:n_infer_steps, :] if n_infer_steps <= max_steps else target_states
n_align = min(states_arr.shape[0], target_arr.shape[0])
if n_align >= CONV_LEN:
sim = compute_similarity(target_arr, states_arr[:n_align], CONV_LEN)
else:
sim = 0.0
print(f" similarity (vs target) = {sim:.4f}")
result["controlled"] = True
result["similarity"] = sim
else:
print(f" WARNING: no model for {scene_name}")
del ff
with open(os.path.join(output_root, "result.json"), "w") as f:
json.dump(result, f, indent=2)
return result
def main():
ap = argparse.ArgumentParser(description="Vortex cloak inference")
ap.add_argument("--type", type=str, default="lamb",
help='Vortex type: lamb, taylor, or "all"')
ap.add_argument("--device", type=int, default=0, help="GPU device ID")
ap.add_argument("--steps", type=int, default=None,
help="Inference steps (default: max_steps for scene)")
ap.add_argument("--out-root", type=str, default=None)
args = ap.parse_args()
if args.type.lower() == "all":
scene_names = get_scene_list("vortex")
else:
scene_names = [f"vortex_{args.type.lower()}"]
if args.out_root is None:
out_root = os.path.join(os.path.dirname(__file__), "..", "data", "vortex")
else:
out_root = args.out_root
t_start = time.time()
for sn in scene_names:
case_dir = os.path.join(out_root, sn)
result = run_single_vortex(sn, args.device, case_dir,
n_infer_steps=args.steps)
print(f" Done: {sn} -> {case_dir}")
elapsed = time.time() - t_start
print(f"\nTotal time: {elapsed:.1f}s")
if __name__ == "__main__":
main()
src/SR_analysis/sindy/illusion/pareto_illusion_1.5L.json
+102
View File
@@ -0,0 +1,102 @@
{
"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
@@ -0,0 +1,110 @@
{
"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
@@ -0,0 +1,118 @@
{
"scene": "karman_re100",
"channels": [
{
"channel": "front",
"best_r2": 0.9950013415086292,
"best_nz": 21,
"pareto": [
{
"nz": 2,
"r2": 0.8681678005649112
},
{
"nz": 3,
"r2": 0.9692414821446799
},
{
"nz": 4,
"r2": 0.989378747581318
},
{
"nz": 6,
"r2": 0.9908017426802015
},
{
"nz": 12,
"r2": 0.9939890287801123
},
{
"nz": 13,
"r2": 0.9947660566975605
},
{
"nz": 18,
"r2": 0.9949963098276752
},
{
"nz": 21,
"r2": 0.9950013415086292
}
]
},
{
"channel": "top",
"best_r2": 0.9928439394510484,
"best_nz": 22,
"pareto": [
{
"nz": 1,
"r2": 0.44188145385324007
},
{
"nz": 2,
"r2": 0.9285296099294669
},
{
"nz": 6,
"r2": 0.9812946866555053
},
{
"nz": 12,
"r2": 0.9909377708950154
},
{
"nz": 17,
"r2": 0.9927843560231949
},
{
"nz": 20,
"r2": 0.992824536423302
},
{
"nz": 22,
"r2": 0.9928439394510484
}
]
},
{
"channel": "bottom",
"best_r2": 0.9965335484991993,
"best_nz": 22,
"pareto": [
{
"nz": 1,
"r2": 0.8456588970543057
},
{
"nz": 2,
"r2": 0.9030386794798415
},
{
"nz": 8,
"r2": 0.9947803148283133
},
{
"nz": 10,
"r2": 0.9962551509064745
},
{
"nz": 13,
"r2": 0.9963877299663628
},
{
"nz": 16,
"r2": 0.9964335809851655
},
{
"nz": 18,
"r2": 0.9965311727162228
},
{
"nz": 22,
"r2": 0.9965335484991993
}
]
}
]
}
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_illusion.py
+133
View File
@@ -0,0 +1,133 @@
"""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
@@ -0,0 +1,159 @@
"""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
@@ -0,0 +1,148 @@
"""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_vortex.py
+133
View File
@@ -0,0 +1,133 @@
"""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
@@ -0,0 +1,110 @@
{
"scene": "vortex_lamb",
"channels": [
{
"channel": "front",
"best_r2": 0.9035051679446613,
"best_nz": 21,
"pareto": [
{
"nz": 3,
"r2": 0.8536388609680176
},
{
"nz": 8,
"r2": 0.8957625139671737
},
{
"nz": 9,
"r2": 0.8985847902462778
},
{
"nz": 16,
"r2": 0.9017891237504362
},
{
"nz": 20,
"r2": 0.9035028044287374
},
{
"nz": 21,
"r2": 0.9035051679446613
}
]
},
{
"channel": "top",
"best_r2": 0.979810012198325,
"best_nz": 22,
"pareto": [
{
"nz": 0,
"r2": -0.2926478447353347
},
{
"nz": 1,
"r2": 0.9262664550660489
},
{
"nz": 2,
"r2": 0.9602153300731789
},
{
"nz": 5,
"r2": 0.9685659511773673
},
{
"nz": 6,
"r2": 0.9752589669369754
},
{
"nz": 19,
"r2": 0.9796029724451923
},
{
"nz": 20,
"r2": 0.9797408009698567
},
{
"nz": 22,
"r2": 0.979810012198325
}
]
},
{
"channel": "bottom",
"best_r2": 0.9334422885170565,
"best_nz": 22,
"pareto": [
{
"nz": 0,
"r2": -5.4349952892154265
},
{
"nz": 1,
"r2": 0.6935730348615337
},
{
"nz": 5,
"r2": 0.8721441125489685
},
{
"nz": 7,
"r2": 0.8963211646824344
},
{
"nz": 11,
"r2": 0.9294742132351927
},
{
"nz": 15,
"r2": 0.9318911728011019
},
{
"nz": 22,
"r2": 0.9334422885170565
}
]
}
]
}
src/SR_analysis/sindy/vortex/pareto_vortex_taylor.json
+58
View File
@@ -0,0 +1,58 @@
{
"scene": "vortex_taylor",
"channels": [
{
"channel": "front",
"best_r2": 0.9603622630700551,
"best_nz": 21,
"pareto": [
{
"nz": 0,
"r2": -1762.7026716770156
},
{
"nz": 21,
"r2": 0.9603622630700551
}
]
},
{
"channel": "top",
"best_r2": 0.809824114603052,
"best_nz": 22,
"pareto": [
{
"nz": 0,
"r2": -3813.3981416722418
},
{
"nz": 1,
"r2": 4.909409545561516e-09
},
{
"nz": 22,
"r2": 0.809824114603052
}
]
},
{
"channel": "bottom",
"best_r2": 0.6431303693566448,
"best_nz": 22,
"pareto": [
{
"nz": 0,
"r2": -11389.175969107222
},
{
"nz": 1,
"r2": 2.5346330034814457e-08
},
{
"nz": 22,
"r2": 0.6431303693566448
}
]
}
]
}
src/SR_analysis/sindy/vortex/sindy_results.json
+996
View File
@@ -0,0 +1,996 @@
{
"thresholds": [
0.0,
0.001,
0.002,
0.005,
0.01,
0.015,
0.02,
0.03,
0.05,
0.1
],
"per_scene": {
"vortex_lamb": {
"scene": "vortex_lamb",
"re_code": 100,
"mu": 0.02,
"n_samples": 148,
"feature_names_front": [
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff"
],
"feature_names_rear": [
"bias",
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff"
],
"front": {
"results": [
{
"threshold": 0.0,
"nz": 21,
"r2": 0.9035051679446613,
"mae": 0.05502827225008196
},
{
"threshold": 0.001,
"nz": 20,
"r2": 0.9035028044287374,
"mae": 0.054966606007031876
},
{
"threshold": 0.002,
"nz": 20,
"r2": 0.9035028044287374,
"mae": 0.054966606007031876
},
{
"threshold": 0.005,
"nz": 16,
"r2": 0.9017891237504362,
"mae": 0.053710513734883655
},
{
"threshold": 0.01,
"nz": 9,
"r2": 0.8985847902462778,
"mae": 0.054491226319598914
},
{
"threshold": 0.015,
"nz": 8,
"r2": 0.8957625139671737,
"mae": 0.05306021520916354
},
{
"threshold": 0.02,
"nz": 8,
"r2": 0.8957625139671737,
"mae": 0.05306021520916354
},
{
"threshold": 0.03,
"nz": 8,
"r2": 0.8957625139671737,
"mae": 0.05306021520916354
},
{
"threshold": 0.05,
"nz": 3,
"r2": 0.8536388609680176,
"mae": 0.048637995761707624
},
{
"threshold": 0.1,
"nz": 3,
"r2": 0.8536388609680176,
"mae": 0.048637995761707624
}
],
"best": {
"threshold": 0.0,
"nz": 21,
"r2": 0.9035051679446613,
"mae": 0.05502827225008196
},
"best_coef": [
0.0072499249196507284,
-0.013406647625693383,
-0.00038101504054493135,
-0.01577043426582744,
-0.021357793291808133,
0.09538149016633536,
-0.17597430070013861,
0.01105460767691196,
0.013467531810937875,
-0.008447877114831191,
-0.009842767133851753,
0.019209968175028087,
-0.036801729600496,
0.0036358387699529284,
0.006079481319894673,
0.003342294121453103,
-0.2971343137386559,
-0.6703323813926846,
-0.7885217073102587,
-1.0678896759007668,
0.5527301613931733
],
"sparsity_curve": [
[
0.0,
21,
0.9035051679446613
],
[
0.001,
20,
0.9035028044287374
],
[
0.002,
20,
0.9035028044287374
],
[
0.005,
16,
0.9017891237504362
],
[
0.01,
9,
0.8985847902462778
],
[
0.015,
8,
0.8957625139671737
],
[
0.02,
8,
0.8957625139671737
],
[
0.03,
8,
0.8957625139671737
],
[
0.05,
3,
0.8536388609680176
],
[
0.1,
3,
0.8536388609680176
]
]
},
"top": {
"results": [
{
"threshold": 0.0,
"nz": 22,
"r2": 0.979810012198325,
"mae": 0.009405479490894075
},
{
"threshold": 0.001,
"nz": 20,
"r2": 0.9797408009698567,
"mae": 0.009464924709826631
},
{
"threshold": 0.002,
"nz": 19,
"r2": 0.9796029724451923,
"mae": 0.009503048001896178
},
{
"threshold": 0.005,
"nz": 6,
"r2": 0.9752589669369754,
"mae": 0.009083428201966606
},
{
"threshold": 0.01,
"nz": 5,
"r2": 0.9685659511773673,
"mae": 0.01027038394304665
},
{
"threshold": 0.015,
"nz": 2,
"r2": 0.9602153300731789,
"mae": 0.012152564325363832
},
{
"threshold": 0.02,
"nz": 2,
"r2": 0.9602153300731789,
"mae": 0.012152564325363832
},
{
"threshold": 0.03,
"nz": 2,
"r2": 0.9602153300731789,
"mae": 0.012152564325363832
},
{
"threshold": 0.05,
"nz": 1,
"r2": 0.9262664550660489,
"mae": 0.013270022046144844
},
{
"threshold": 0.1,
"nz": 0,
"r2": -0.2926478447353347,
"mae": 0.10377883511859723
}
],
"best": {
"threshold": 0.0,
"nz": 22,
"r2": 0.979810012198325,
"mae": 0.009405479490894075
},
"best_coef": [
1.7326565320927485,
-0.021510750075876994,
0.0007162636287450304,
0.003690011764588065,
-0.001255780396435901,
0.008069876064011276,
-0.027165643005235728,
0.019732103606457083,
-0.007023240411992188,
-0.00413143280625655,
0.0016423995129504244,
0.0008079534031552331,
-0.005050487340309005,
0.011582526033145867,
0.0010455718346845003,
0.0012590379370475096,
0.0019097817325384912,
0.0346531306245988,
0.03581318435812541,
-0.06278853147020641,
0.40349380438892807,
-0.3511614171629877
],
"sparsity_curve": [
[
0.0,
22,
0.979810012198325
],
[
0.001,
20,
0.9797408009698567
],
[
0.002,
19,
0.9796029724451923
],
[
0.005,
6,
0.9752589669369754
],
[
0.01,
5,
0.9685659511773673
],
[
0.015,
2,
0.9602153300731789
],
[
0.02,
2,
0.9602153300731789
],
[
0.03,
2,
0.9602153300731789
],
[
0.05,
1,
0.9262664550660489
],
[
0.1,
0,
-0.2926478447353347
]
]
},
"bottom": {
"results": [
{
"threshold": 0.0,
"nz": 22,
"r2": 0.9334422885170565,
"mae": 0.0062053698726041735
},
{
"threshold": 0.001,
"nz": 15,
"r2": 0.9318911728011019,
"mae": 0.00612106433470266
},
{
"threshold": 0.002,
"nz": 11,
"r2": 0.9294742132351927,
"mae": 0.006263040564511156
},
{
"threshold": 0.005,
"nz": 7,
"r2": 0.8963211646824344,
"mae": 0.006779203944819282
},
{
"threshold": 0.01,
"nz": 5,
"r2": 0.8721441125489685,
"mae": 0.007291293800356046
},
{
"threshold": 0.015,
"nz": 1,
"r2": 0.6935730348615337,
"mae": 0.009771975563999018
},
{
"threshold": 0.02,
"nz": 1,
"r2": 0.6935730348615337,
"mae": 0.009771975563999018
},
{
"threshold": 0.03,
"nz": 1,
"r2": 3.300804074513053e-12,
"mae": 0.029445380091027484
},
{
"threshold": 0.05,
"nz": 1,
"r2": 3.300804074513053e-12,
"mae": 0.029445380091027484
},
{
"threshold": 0.1,
"nz": 0,
"r2": -5.4349952892154265,
"mae": 0.0796928089615461
}
],
"best": {
"threshold": 0.0,
"nz": 22,
"r2": 0.9334422885170565,
"mae": 0.0062053698726041735
},
"best_coef": [
0.6365299396095444,
-0.00680014903595817,
-0.0014763845224131282,
0.0008331648138976364,
-0.0013963333722837737,
0.0026865900945637236,
-0.014425589645304323,
-0.008472370593784707,
-6.004332296000475e-05,
-0.0014926889561804521,
-0.0004258139447154254,
-0.0018644416734531306,
0.006364757842372448,
-0.003199874412701356,
0.001370932025212518,
0.002115447057480394,
0.0019992859663057654,
0.012730598783987558,
-0.07381922521108229,
-0.069816526655473,
0.13432950141437106,
-0.0030019765287261236
],
"sparsity_curve": [
[
0.0,
22,
0.9334422885170565
],
[
0.001,
15,
0.9318911728011019
],
[
0.002,
11,
0.9294742132351927
],
[
0.005,
7,
0.8963211646824344
],
[
0.01,
5,
0.8721441125489685
],
[
0.015,
1,
0.6935730348615337
],
[
0.02,
1,
0.6935730348615337
],
[
0.03,
1,
3.300804074513053e-12
],
[
0.05,
1,
3.300804074513053e-12
],
[
0.1,
0,
-5.4349952892154265
]
]
}
},
"vortex_taylor": {
"scene": "vortex_taylor",
"re_code": 100,
"mu": 0.02,
"n_samples": 148,
"feature_names_front": [
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff"
],
"feature_names_rear": [
"bias",
"u_m",
"u_a",
"u_c",
"v_a",
"Cd_tot",
"Cd_rear",
"Cl_tot",
"Cl_diff",
"sin_ua",
"cos_ua",
"aF_lag1",
"aB_lag1",
"aT_lag1",
"daF",
"daB",
"daT",
"mu",
"mu_u_a",
"mu_v_a",
"mu_Cd_tot",
"mu_Cl_diff"
],
"front": {
"results": [
{
"threshold": 0.0,
"nz": 21,
"r2": 0.9603622630700551,
"mae": 5.5722956333111334e-05
},
{
"threshold": 0.001,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.002,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.005,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.01,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.015,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.02,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.03,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.05,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
},
{
"threshold": 0.1,
"nz": 0,
"r2": -1762.7026716770156,
"mae": 0.019129430850011273
}
],
"best": {
"threshold": 0.0,
"nz": 21,
"r2": 0.9603622630700551,
"mae": 5.5722956333111334e-05
},
"best_coef": [
0.00013185241334928363,
0.0013869486093058166,
-2.0689454061431898e-05,
0.0003650495654675631,
-7.646260510377518e-05,
0.0016834710794748663,
0.002138870811614627,
0.0009672545588015969,
-0.0011958663539923654,
-0.00015017143903531285,
0.007403929744140764,
-0.0015484862269750206,
-0.0003240821985992223,
-0.0010978795962499714,
-0.0002774795287867582,
-0.00041159658606479107,
0.00011190538695352297,
0.06934742736821846,
0.018252511507361284,
-0.0038231311262574906,
0.04836265570630274
],
"sparsity_curve": [
[
0.0,
21,
0.9603622630700551
],
[
0.001,
0,
-1762.7026716770156
],
[
0.002,
0,
-1762.7026716770156
],
[
0.005,
0,
-1762.7026716770156
],
[
0.01,
0,
-1762.7026716770156
],
[
0.015,
0,
-1762.7026716770156
],
[
0.02,
0,
-1762.7026716770156
],
[
0.03,
0,
-1762.7026716770156
],
[
0.05,
0,
-1762.7026716770156
],
[
0.1,
0,
-1762.7026716770156
]
]
},
"top": {
"results": [
{
"threshold": 0.0,
"nz": 22,
"r2": 0.809824114603052,
"mae": 0.00019280734797529785
},
{
"threshold": 0.001,
"nz": 1,
"r2": 4.909409545561516e-09,
"mae": 0.0004627442799109909
},
{
"threshold": 0.002,
"nz": 1,
"r2": 4.909409545561516e-09,
"mae": 0.0004627442799109909
},
{
"threshold": 0.005,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
},
{
"threshold": 0.01,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
},
{
"threshold": 0.015,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
},
{
"threshold": 0.02,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
},
{
"threshold": 0.03,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
},
{
"threshold": 0.05,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
},
{
"threshold": 0.1,
"nz": 0,
"r2": -3813.3981416722418,
"mae": 0.06224470555379584
}
],
"best": {
"threshold": 0.0,
"nz": 22,
"r2": 0.809824114603052,
"mae": 0.00019280734797529785
},
"best_coef": [
-0.003855648809544881,
-0.00016402917649747364,
-0.00930469925984968,
0.0003725286141146794,
-0.0018143741612594137,
-0.0006228672933056743,
-0.004123657830294623,
0.009617379353839065,
-0.0020693580870050163,
0.0064472830467429834,
-0.0009363904689749538,
0.002848747579355727,
0.008718629599353982,
0.006611618826853361,
0.0026442625498493562,
-0.000822158875902874,
-0.00019526846815097482,
-7.71129751901901e-05,
-0.4652349584088788,
-0.09071890000133181,
-0.031143453311218216,
-0.10346632929440941
],
"sparsity_curve": [
[
0.0,
22,
0.809824114603052
],
[
0.001,
1,
4.909409545561516e-09
],
[
0.002,
1,
4.909409545561516e-09
],
[
0.005,
0,
-3813.3981416722418
],
[
0.01,
0,
-3813.3981416722418
],
[
0.015,
0,
-3813.3981416722418
],
[
0.02,
0,
-3813.3981416722418
],
[
0.03,
0,
-3813.3981416722418
],
[
0.05,
0,
-3813.3981416722418
],
[
0.1,
0,
-3813.3981416722418
]
]
},
"bottom": {
"results": [
{
"threshold": 0.0,
"nz": 22,
"r2": 0.6431303693566448,
"mae": 0.00012115305583366485
},
{
"threshold": 0.001,
"nz": 1,
"r2": 2.5346330034814457e-08,
"mae": 0.00029202565622359403
},
{
"threshold": 0.002,
"nz": 1,
"r2": 2.5346330034814457e-08,
"mae": 0.00029202565622359403
},
{
"threshold": 0.005,
"nz": 1,
"r2": 2.5346330034814457e-08,
"mae": 0.00029202565622359403
},
{
"threshold": 0.01,
"nz": 1,
"r2": 2.5346330034814457e-08,
"mae": 0.00029202565622359403
},
{
"threshold": 0.015,
"nz": 1,
"r2": 2.5346330034814457e-08,
"mae": 0.00029202565622359403
},
{
"threshold": 0.02,
"nz": 1,
"r2": 2.5346330034814457e-08,
"mae": 0.00029202565622359403
},
{
"threshold": 0.03,
"nz": 0,
"r2": -11389.175969107222,
"mae": 0.05019249195686063
},
{
"threshold": 0.05,
"nz": 0,
"r2": -11389.175969107222,
"mae": 0.05019249195686063
},
{
"threshold": 0.1,
"nz": 0,
"r2": -11389.175969107222,
"mae": 0.05019249195686063
}
],
"best": {
"threshold": 0.0,
"nz": 22,
"r2": 0.6431303693566448,
"mae": 0.00012115305583366485
},
"best_coef": [
0.024327554221764875,
-0.0005249278989068122,
-0.0012993577286068104,
2.2680966036791235e-07,
-0.00010725143675978186,
-6.83864821520398e-05,
0.0010096871646746155,
0.010716542986185212,
0.0009908500661735154,
0.0011374615793514392,
0.0001758268322788275,
-0.0008997304730805852,
-0.0002053415257431248,
0.0003868365622255378,
0.005272105474899691,
0.0033010584355117824,
0.0027208004849191268,
0.0004865510823514911,
-0.06496786884052097,
-0.005364174292120229,
-0.003419346453559459,
0.04954205911186436
],
"sparsity_curve": [
[
0.0,
22,
0.6431303693566448
],
[
0.001,
1,
2.5346330034814457e-08
],
[
0.002,
1,
2.5346330034814457e-08
],
[
0.005,
1,
2.5346330034814457e-08
],
[
0.01,
1,
2.5346330034814457e-08
],
[
0.015,
1,
2.5346330034814457e-08
],
[
0.02,
1,
2.5346330034814457e-08
],
[
0.03,
0,
-11389.175969107222
],
[
0.05,
0,
-11389.175969107222
],
[
0.1,
0,
-11389.175969107222
]
]
}
}
}
}
src/SR_analysis/sindy_sr_knoeledge.md
+339
View File
@@ -0,0 +1,339 @@
# 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_notes.md
+347
View File
@@ -0,0 +1,347 @@
# SINDy 与 SR 执行计划
## 文档作用
这份文档只回答一件事：**接下来要做什么。**
- 只写执行路线、阶段目标、优先级、输出要求
- 不长篇复述历史争论
- 不把背景知识和任务顺序混在一起
- 凡是历史 bug、经验教训、哪些结论已经成立、哪些还不能说，统一放到 `sindy_sr_knowladge`
---
## 当前总目标
把 **所有 cloak 场景** 纳入同一受约束分析框架，使用 **SINDy + SR 并行** 探索控制骨架，先分别拟合，再横向比较，从而判断：
- 哪些项是 shared core
- 哪些项是 scene-specific activation
- 哪些差异来自时间尺度写法，而不是物理骨架本身
当前默认路线不是“先强行做统一总公式”，而是：
\[
\text{separate fit} \rightarrow \text{compare} \rightarrow \text{shared-backbone test}
\]
---
## 当前工作原则
后续默认遵守：
1. **所有场景共用同一套 primitive variables 与同一套 \(G\) 规则**
2. **front 默认 no-bias + odd structure**
3. **rear 默认 shared-head**
4. **SINDy 与 SR 并行推进**，SR 不是替换 SINDy，而是并行工具
5. **先分场景拟合，再做横向比较**
6. **闭环验证不可省略**，但不要求每个场景第一轮都做全套闭环
7. **时间尺度问题必须显式进入特征定义**
8. **不再只围绕跨 Re 的 Kármán 单线深挖**；跨 Re 现在是已站住的第一证据，不是全部主线
9. **当前不讨论功率与能量分析**，这不是这条线的优先任务
10. **当前不把时延作为主要矛盾**；对稳定周期控制，当前第一矛盾是变量骨架与时间尺度写法
---
## 当前场景优先级
### 第一层：本轮重点场景
- Kármán cloak
- steady cloak
这两个场景本轮必须做完整输出，因为它们最适合先建立场景间比较模板。
### 第二层：下一轮扩展场景
- 单涡 cloak（monopole / taylor / lamb 等已有单涡场景）
- erase
- 其他已有 cloak 场景
这批先做轻量版 separate fit，再决定哪些值得补完整闭环与深挖。
---
## 阶段 0
## 统一接口
### 0.1 统一 primitive variables
所有场景统一输出：
- 无量纲 sensor：\(\hat u, \hat v\)
- 力系数：\(C_D, C_L\)
- 无量纲动作：\(\alpha\)
- lagged \(\alpha\)
- \(\Delta \alpha\) 或显式含 \(\Delta t_c\) 的版本
- \(\mu = 1/Re_D\)
- scene metadata：scene id、\(Re_D\)、control interval \(\Delta t_c\)、target type、采样设置
### 0.2 固定 \(G\) 算子
统一使用同一套 \(G\) 规则，不允许不同场景临时改写。动作必须满足：
\[
(\alpha_F,\alpha_T,\alpha_B) \mapsto (-\alpha_F,-\alpha_B,-\alpha_T)
\]
### 0.3 统一 feature builder
统一生成三层特征：
| 层级 | 内容 | 用途 |
|---|---|---|
| core | \(\hat u, \hat v, C_D, C_L, \alpha^-, \Delta\alpha^-, \mu\) | 所有场景共用 |
| derived | 对称/反对称组合、总量/差量 | SINDy 主库 |
| time-scale | 显式含 \(\Delta t_c\) 的版本 | 时间尺度分析 |
### 0.4 必做测试
任何场景进入拟合前，先过：
- \(G(G(x)) = x\)
- feature names 与矩阵列顺序一致
- 闭环预测器输入维度一致
- front / rear 的结构约束在数据接口层正确落地
### 阶段 0 输出
- 统一 `feature_builder`
- 统一 `symmetry`
- 统一 `time_scale`
- 统一 `validators`
- 场景级最小数据摘要（样本数、变量范围、control interval）
---
## 阶段 1
## Kármán 与 steady 的第一轮 separate SINDy
### 目标
在统一变量与统一约束下，先得到两个重点场景各自最可信的稀疏 support。
### 默认约束
- front no-bias
- front odd structure
- rear shared-head
- correct \(G\) consistency
### 每个场景必须输出
| 输出 | 说明 |
|---|---|
| best support | 最优 support 列表 |
| sparsity curve | 稀疏度-误差曲线 |
| front / rear 主项表 | 主导项与系数 |
| one-step metrics | R²、RMSE |
| key closed-loop result | 至少一个关键闭环指标 |
| support stability | threshold / window / bootstrap 稳定性 |
| contribution table | 主要项贡献度，不只看是否出现 |
### 本阶段比较项
Kármán 与 steady 做第一轮 support overlap：
- 哪些项共同出现
- 哪些项只在 Kármán 激活
- steady 是否表现为明显简化版
- overlap 不能只看布尔出现，还要看贡献量级与稳定性
---
## 阶段 2
## Kármán 与 steady 的第一轮受限 SR
### 目标
不是追求最终公式，而是看：
- 在受限物理库上，是否能压出更短闭式
- SINDy 中看似不同的项，是否会被更统一的表达吸收
- Kármán 与 steady 是否出现同形公式
### SR 输入规则
SR 只能使用：
- 阶段 0 的统一变量
- 阶段 1 的 SINDy 已筛出主项及其邻近项
- 受限运算集合
### 当前允许的运算
- 加减乘
- protected divide
- 少量 square
- 必要时有限放开 tanh
### 当前不允许的运算
- raw trig 乱搜
- 高次幂
- 指数
- 深层嵌套
### 每个场景必须输出
| 输出 | 说明 |
|---|---|
| shortest acceptable formula | 最短可接受公式 |
| complexity-error pareto | 复杂度 vs 误差 |
| 和 SINDy 的关系 | 压缩了什么、保留了什么、吸收了什么 |
| key closed-loop result | 最佳 SR 公式至少一版关键闭环 |
| formula family notes | 同一场景内是否存在多种同等可接受闭式 |
---
## 阶段 3
## 横向比较
当 Kármán 与 steady 的 SINDy + SR 都出来后，立即做横向比较。
### 3.1 support 比较
不要只看“是否出现”，必须同时看：
- 是否出现
- 系数或贡献量级
- 稳定性
- SR 是否把它吸收到更高层表达里
### 3.2 公式形态比较
至少检查：
- 是否都包含同类 force feedback 核心
- 是否都包含同类 memory 核心
- steady 是否只是删掉了周期相关项
- 是否存在同形结构 + 少量场景激活项
### 3.3 第一轮 shared-backbone 判断
这一轮只回答：
1. 是否存在 Kármán 与 steady 的 shared core
2. steady 是否可以视为 Kármán 的明显简化版
3. 哪些项更像 scene-specific activation，而不是 backbone 本身
本阶段不急着拟合 all-cloak 联合总公式。
---
## 阶段 4
## 扩展到单涡与其他 cloak
在 Kármán 与 steady 的第一轮比较完成后，再把单涡 cloak 与其他场景纳入。
### 轻量版输出要求
- separate SINDy
- separate SR
- one-step metrics
- support 形态
- 必要时补关键闭环
### 比较目标
重点看：
- 单涡是否保留 shared core
- 单涡是否主要新增 history / transient 项
- 是否开始出现子家族结构
---
## 阶段 5
## 时间尺度显式化
这条线并行推进，但先不抢在全部场景前面。
### 当前目标
- 让 \(\Delta t_c\) 显式进入特征
- 不再把“1 个采样步”默认当物理可比量
- 比较显式化前后，support 是否更收敛
- 为后面严肃讨论采样率影响扫清接口问题
### 第一批测试场景
- Kármán cloak
- steady cloak
### 第一批对比内容
- 旧的 discrete lag / \(\Delta a\)
- 显式含 \(\Delta t_c\) 的版本
### 关注结果
- support 是否变化
- SR 公式是否更统一
- 跨场景比较是否更干净
注意：当前阶段的目标是**时间尺度显式化**，不是立刻给出高采样率优于 PPO 的最终结论。
---
## 本轮必须完成的最小任务
1. 统一 Kármán 与 steady 的 feature builder
2. 统一 Kármán 与 steady 的 \(G\) / time-scale / validators
3. 跑 Kármán 与 steady 的第一轮 separate SINDy
4. 跑 Kármán 与 steady 的第一轮受限 SR（PySR）
5. 输出 Kármán vs steady 的：
- support overlap
- 公式形态比较
- shared core / scene-specific 的初步分类
---
## 本轮结束时应交付的结果包
每个重点场景至少有一张 summary 表：
| method | sparsity | one-step | closed-loop | key terms | notes |
|---|---:|---:|---:|---|---|
| SINDy | | | | | |
| SR | | | | | |
以及一个跨场景比较表：
| comparison | shared core | scene-enhanced | scene-specific | notes |
|---|---|---|---|---|
| Kármán vs steady | | | | |
---
## 当前不该做的事
- 不继续围绕单一 Kármán across Re 版本号升级
- 不把 threshold 网格或简单 Pareto 扫描直接当成完整 SR
- 不在 raw feature 上做自由 SR
- 不在不同采样间隔下直接复用旧系数并据此下正式结论
- 不在场景比较证据还弱时，提前宣布 all-cloak 统一总公式
- 不把功率、能量、时延这些非当前主矛盾问题拉入本轮 SINDy/SR 主线
---
## 当前这份 notes 的直接收束
接下来核心任务不是“继续优化某个跨 Re 模型”，而是：
\[
\boxed{
\text{把 Kármán 与 steady 先做成可比较的 separate SINDy + separate SR 结果包，然后以它们为模板扩到单涡与其他 cloak。}
}
\]
这一步完成后，才进入更严肃的 all-cloak shared-backbone 判断。
src/SR_analysis/utils/__init__.py
+9
View File
@@ -0,0 +1,9 @@
from .feature_builder import (
compute_dimensionless, compute_features, build_feature_matrix,
get_feature_names, apply_G_alpha, apply_G_x,
CORE_FEAT_KEYS, MU_FEAT_KEYS, ALL_FEAT_KEYS,
)
from .sindy_fitter import (
fit_channel, fit_sindy, print_control_law,
get_active_support, get_feature_matrix_from_data,
)
src/SR_analysis/utils/cfd_interface.py
+397
View File
@@ -0,0 +1,397 @@
"""CFD interface for LegacyCelerisLab (pycuda_3_10 env).
All functions use the LegacyCelerisLab (old) CFD API via:
from LegacyCelerisLab import FlowField
Must be run inside: conda run -n pycuda_3_10
NOTE: This module should be imported directly, not through SR_analysis.utils
because it requires pycuda. Other utils (sindy_fitter, feature_builder, g_operator)
do NOT require pycuda and can be imported from the __init__.
"""
from __future__ import annotations
import json
import os
import sys
from collections import deque
from typing import Any, Dict, List, Optional, Tuple
import numpy as np
# -- Import legacy CFD -------------------------------------------------------
# LegacyCelerisLab lives at the repo root; SR_analysis is at repo_root/src/SR_analysis.
_REPO = os.path.abspath(os.path.join(os.path.dirname(__file__), "..", "..", ".."))
if _REPO not in sys.path:
sys.path.insert(0, _REPO)
_SRC = os.path.join(_REPO, "src")
if _SRC not in sys.path:
sys.path.insert(0, _SRC)
from LegacyCelerisLab import FlowField # noqa: E402
from LegacyCelerisLab import utils as legacy_utils # noqa: E402
# ---------------------------------------------------------------------------
# Action-smoothing constant (legacy run() internal)
# ---------------------------------------------------------------------------
ACTION_SMOOTH_WEIGHT = 0.1 # used by FlowField.run() internally
def nu_from_re(re_code: float, u0: float = 0.01, d_ref: float = 40.0) -> float:
"""Return kinematic viscosity for a given code Reynolds number.
``re_code`` uses reference length *2*D* = 40.0 (matching model file naming).
"""
return u0 * d_ref / re_code
def load_legacy_configs(config_dir: str) -> Tuple[Any, Any]:
"""Load and return legacy (cuda_config, field_config) from *config_dir*."""
cuda_cfg = legacy_utils.load_cuda_config(
os.path.join(config_dir, "config_cuda.json")
)
field_cfg = legacy_utils.load_flow_field_config(
os.path.join(config_dir, "config_flowfield.json")
)
return cuda_cfg, field_cfg
# ---------------------------------------------------------------------------
# Environment helpers -- Karman cloak (disturbance cylinder + pinball)
# ---------------------------------------------------------------------------
def build_karman_cloak_env(
flow_field: FlowField,
*,
u0: float,
l0: float,
sample_interval: int,
fifo_len: int,
data_type: type,
) -> Tuple[np.ndarray, dict]:
"""Phase 0-1: add dist-cylinder & 3 sensors, stabilize, record target.
Steps (mirrors env_karman_cloak_standard.__init__):
1. add dist_cylinder (id=0)
2. add 3 sensors (id=1,2,3)
3. stabilize run(4*NX/U0, zero-action[4])
4. record FIFO_LEN x run(SAMPLE_INTERVAL, zero[4]), collect obs[2:8]
Returns
-------
target_states : ndarray (FIFO_LEN, 6) -- sensor0/1/2 ux,uy
info : dict with n_objects, NX, NY
"""
# dist cylinder
center = (10.0 * l0, (flow_field.FIELD_SHAPE[1] - 1) / 2, 0.0)
flow_field.add_cylinder(center, l0)
# sensors
for y_off in [2.0, 0.0, -2.0]:
sc = (40.0 * l0, (flow_field.FIELD_SHAPE[1] - 1) / 2 + y_off * l0, 0.0)
flow_field.add_sensor(sc, l0 / 4.0)
n_obj = flow_field.obs.size // 2
# stabilize
stabilize_steps = int(4 * flow_field.FIELD_SHAPE[0] / u0)
print(f" stabilising ({stabilize_steps} steps)...")
flow_field.run(stabilize_steps, np.zeros(n_obj, dtype=data_type))
# record target (only sensor signals = obs[2:8])
target_states = np.empty((0, 6), dtype=data_type)
for _ in range(fifo_len):
flow_field.run(sample_interval, np.zeros(n_obj, dtype=data_type))
new_state = flow_field.obs.copy()[2:8]
target_states = np.vstack((target_states, new_state))
print(f" target recorded: {target_states.shape}")
return target_states, {"n_objects": n_obj, "NX": flow_field.FIELD_SHAPE[0],
"NY": flow_field.FIELD_SHAPE[1]}
def add_pinball(
flow_field: FlowField,
*,
l0: float,
u0: float,
sample_interval: int,
fifo_len: int,
data_type: type,
action_bias: Optional[Tuple[float, float, float]] = None,
pinball_front_x: float = 30.0,
pinball_rear_x: float = 31.3,
obs_slice_start: int = 2,
obs_slice_end: int = 14,
n_objects_total: Optional[int] = None,
) -> dict:
"""Add pinball cylinders, stabilize, compute norm, bias rollout.
Steps:
1. add front, bottom, top cylinders
2. stabilize run(4*NX/U0, zero-action)
3. get_ddf() + save_ddf() (checkpoint)
4. FIFO_LEN x run(SAMPLE_INTERVAL, zero) -> compute norm
5. apply_ddf() (restore pre-bias state)
6. FIFO_LEN x run(SAMPLE_INTERVAL, bias-action) -> save_states
7. apply_ddf()
Parameters
----------
pinball_front_x, pinball_rear_x : pinball geometry (L0 units).
Default 30.0/31.3 for Karman; 19.0/20.3 for Illusion.
obs_slice_start, obs_slice_end : slice of obs for norm.
Default [2:14] for Karman (7 objects); [0:12] for Illusion (6 objects).
n_objects_total : if provided, used for bias array length.
Default: inferred from flow_field after adding cylinders.
Returns dict with norm values.
"""
if action_bias is None:
action_bias = (0.0, -4.0, 4.0)
u0_float = float(u0)
# add 3 pinball cylinders
ny = flow_field.FIELD_SHAPE[1]
centers = [
(pinball_front_x * l0, (ny - 1) / 2, 0.0),
(pinball_rear_x * l0, (ny - 1) / 2 + 0.75 * l0, 0.0),
(pinball_rear_x * l0, (ny - 1) / 2 - 0.75 * l0, 0.0),
]
for c in centers:
flow_field.add_cylinder(c, l0 / 2.0)
n_obj = flow_field.obs.size // 2 if n_objects_total is None else n_objects_total
print(f" bodies after pinball: {n_obj}")
# stabilize
stabilize_steps = int(4 * flow_field.FIELD_SHAPE[0] / u0_float)
print(f" stabilising pinball ({stabilize_steps} steps)...")
flow_field.run(stabilize_steps, np.zeros(n_obj, dtype=data_type))
# checkpoint DDF
flow_field.get_ddf()
flow_field.save_ddf()
# --- norm phase (zero-action) ---
fifo = deque(maxlen=fifo_len)
for _ in range(fifo_len):
flow_field.run(sample_interval, np.zeros(n_obj, dtype=data_type))
fifo.append(flow_field.obs.copy()[obs_slice_start:obs_slice_end])
temp_states = np.array(fifo, dtype=data_type)
# forces are at indices [6:12] relative to the slice end
force_start = obs_slice_end - obs_slice_start - 6
force_end = force_start + 6
force_norm_fact = 6.0 * float(np.max(np.abs(temp_states[:, force_start:force_end])))
sens_deviation = np.mean(temp_states[:, 0:6], axis=0).astype(data_type)
sens_norm_fact = np.zeros(6, dtype=data_type)
for i in range(6):
sens_norm_fact[i] = 5.0 * float(np.max(np.abs(temp_states[:, i] - sens_deviation[i])))
print(f" norm: force_norm_fact={force_norm_fact:.6f}")
print(f" norm: sens_deviation={sens_deviation}")
print(f" norm: sens_norm_fact={sens_norm_fact}")
# --- bias-action rollout ---
flow_field.apply_ddf()
bias = np.zeros(n_obj, dtype=data_type)
bias[n_obj - 3] = float(action_bias[0] * u0_float)
bias[n_obj - 2] = float(action_bias[1] * u0_float)
bias[n_obj - 1] = float(action_bias[2] * u0_float)
print(f" bias action: {bias}")
fifo.clear()
for _ in range(fifo_len):
flow_field.run(sample_interval, bias)
fifo.append(flow_field.obs.copy()[obs_slice_start:obs_slice_end])
save_states = np.array(list(fifo), dtype=data_type)
flow_field.apply_ddf()
return {
"force_norm_fact": force_norm_fact,
"sens_deviation": sens_deviation.tolist(),
"sens_norm_fact": sens_norm_fact.tolist(),
"action_bias": list(action_bias),
"save_states": save_states,
}
def build_observation(
obs_slice: np.ndarray,
norm: dict,
) -> np.ndarray:
"""Assemble normalised DRL observation (12-dim) from a single obs slice.
``obs_slice`` is 12-element: sensor[0:6] + force[6:12].
Returns clipped 12-dim array in [-1, 1].
"""
forces = obs_slice[6:12] / norm["force_norm_fact"]
sens = (obs_slice[0:6] - norm["sens_deviation"]) / norm["sens_norm_fact"]
obs = np.clip(np.hstack([forces, sens]), -1.0, 1.0).astype(np.float32)
return obs
def action_to_physical(
action_norm: np.ndarray,
*,
scale: float = 8.0,
bias: Tuple[float, float, float] = (0.0, -4.0, 4.0),
u0: float = 0.01,
) -> np.ndarray:
"""Convert normalized action [-1,1] to physical omega (lattice units).
physical_omega[i] = (action_norm[i] * scale + bias[i]) * u0
"""
action_norm = np.asarray(action_norm, dtype=np.float64).reshape(-1, 3)
bias_arr = np.array(bias, dtype=np.float64)
return (action_norm * scale + bias_arr) * u0
def scale_action(
action_norm: np.ndarray,
*,
scale: float = 8.0,
bias: Tuple[float, float, float] = (0.0, -4.0, 4.0),
u0: float = 0.01,
n_total_bodies: int = 7,
) -> np.ndarray:
"""Convert normalised action ([-1,1]^3) to legacy CFD action array.
Returns array of length *n_total_bodies* with cylinders' omegas at the
last 3 slots.
"""
a = np.zeros(n_total_bodies, dtype=np.float32)
omega = (np.array(action_norm, dtype=np.float32) * scale + np.array(bias, dtype=np.float32)) * u0
a[n_total_bodies - 3:] = omega
return a
# ---------------------------------------------------------------------------
# Vorticity & field export
# ---------------------------------------------------------------------------
def vorticity_from_ddf(flow_field: FlowField, u0: float) -> np.ndarray:
"""Compute z-vorticity from current DDF on host."""
flow_field.get_ddf()
ddf = flow_field.ddf.copy().reshape((9, flow_field.FIELD_SHAPE[1],
flow_field.FIELD_SHAPE[0])).transpose(2, 1, 0)
ux = (ddf[:, :, 1] + ddf[:, :, 5] + ddf[:, :, 8]
- ddf[:, :, 3] - ddf[:, :, 6] - ddf[:, :, 7]) / u0
uy = (ddf[:, :, 2] + ddf[:, :, 5] + ddf[:, :, 6]
- ddf[:, :, 4] - ddf[:, :, 7] - ddf[:, :, 8]) / u0
omega = np.gradient(uy, axis=1) - np.gradient(ux, axis=0)
return omega.astype(np.float64)
def save_vorticity_png(path: str, omega: np.ndarray, title: str = ""):
"""Save vorticity field as a PNG with symmetric colour bar."""
import matplotlib
matplotlib.use("Agg")
import matplotlib.pyplot as plt
abs_o = np.abs(omega[np.isfinite(omega)])
vmax = float(np.percentile(abs_o, 99.5)) if abs_o.size > 0 else 1.0
if vmax <= 0:
vmax = 1.0
ny, nx = omega.shape
fig, ax = plt.subplots(figsize=(min(18, max(8, nx / 60)), min(10, max(3, ny / 40))))
im = ax.imshow(omega, origin="lower", aspect="equal", cmap="RdBu_r",
vmin=-vmax, vmax=vmax, extent=(0, nx - 1, 0, ny - 1))
ax.set_xlabel("x (lattice)")
ax.set_ylabel("y (lattice)")
if title:
ax.set_title(title)
fig.colorbar(im, ax=ax, fraction=0.046, pad=0.04, label=r"$\omega_z$")
fig.tight_layout()
fig.savefig(path, dpi=150, bbox_inches="tight")
plt.close(fig)
# ---------------------------------------------------------------------------
# DTW similarity
# ---------------------------------------------------------------------------
def calc_lag(target: np.ndarray, state: np.ndarray) -> int:
"""Find lag that maximises cross-correlation between two 1-D signals."""
t = target - np.mean(target)
s = state - np.mean(state)
corr = np.correlate(t, s, mode="full")
lags = np.arange(-len(target) + 1, len(target))
return int(lags[np.argmax(corr)])
def calc_dtw_sim(target: np.ndarray, state: np.ndarray) -> float:
"""DTW-based similarity: 1 - (DTW distance / len(target))."""
n, m = len(target), len(state)
dtw = np.full((n + 1, m + 1), np.inf)
dtw[0, 0] = 0.0
for i in range(1, n + 1):
for j in range(1, m + 1):
cost = abs(float(target[i - 1]) - float(state[j - 1]))
dtw[i, j] = cost + min(dtw[i - 1, j], dtw[i, j - 1], dtw[i - 1, j - 1])
return float(1.0 - dtw[n, m] / n)
def compute_similarity(
target_states: np.ndarray,
state_series: np.ndarray,
conv_len: int,
) -> float:
"""Compute lag-compensated DTW similarity over *conv_len* window."""
ref = target_states[conv_len:2 * conv_len, 1]
cur = state_series[-conv_len:, 1]
lag = calc_lag(ref, cur)
sim_sum = 0.0
for i in range(6):
target_seq = np.roll(target_states[:, i], -lag)[conv_len:2 * conv_len]
state_seq = state_series[-conv_len:, i]
sim_sum += calc_dtw_sim(target_seq, state_seq) / 6.0
return float(sim_sum)
# ---------------------------------------------------------------------------
# Dummy env for loading SB3 models
# ---------------------------------------------------------------------------
def create_dummy_env(s_dim: int = 12, a_dim: int = 3):
"""Return a gym.Env with correct observation/action spaces for model loading."""
import gymnasium as gym
from gymnasium import spaces
class DummyEnv(gym.Env):
def __init__(self):
super().__init__()
self.observation_space = spaces.Box(low=-1, high=1, shape=(s_dim,), dtype=np.float32)
self.action_space = spaces.Box(low=-1, high=1, shape=(a_dim,), dtype=np.float32)
def reset(self, seed=None):
return np.zeros(s_dim, dtype=np.float32), {}
def step(self, action):
return np.zeros(s_dim, dtype=np.float32), 0.0, False, False, {}
def render(self):
pass
return DummyEnv()
def load_ppo_model(model_path: str, device: str = "cuda:0", s_dim: int = 12, a_dim: int = 3):
"""Load a PPO model with Sin activation."""
import torch
from torch.nn import Module
from stable_baselines3 import PPO
class Sin(Module):
def forward(self, x):
return torch.sin(x)
dummy_env = create_dummy_env(s_dim, a_dim)
model = PPO.load(model_path, env=dummy_env, device=device)
return model
src/SR_analysis/utils/feature_builder.py
+204
View File
@@ -0,0 +1,204 @@
"""Unified feature builder for all cloak scenes.
Produces dimensionless features with consistent G-equivariant structure.
All scenes (Karman, steady, vortex, illusion) use this same builder.
Copy of analysis_cloak/common/feature_builder.py -- kept as canonical source.
"""
from __future__ import annotations
from typing import Dict, List, Tuple
import numpy as np
# -- Physical constants ------------------------------------------------------
U0 = 0.01 # inlet velocity (lattice units)
D_CYL = 20.0 # cylinder diameter (lattice)
# -- Dimensionless conversion ------------------------------------------------
def compute_dimensionless(
sensors: np.ndarray, # (T, 6) raw lattice [s0_ux,s0_uy, s1_ux,s1_uy, s2_ux,s2_uy]
forces: np.ndarray, # (T, 6) raw lattice [f0_fx,f0_fy, f1_fx,f1_fy, f2_fx,f2_fy]
u0: float = U0,
d: float = D_CYL,
rho: float = 1.0,
) -> Dict[str, np.ndarray]:
"""Convert raw lattice CFD data to dimensionless quantities.
Sensor order: [s0_ux,s0_uy, s1_ux,s1_uy, s2_ux,s2_uy]
where s0=top(y=+2L0), s1=mid(y=0), s2=bottom(y=-2L0)
Force order: [front_fx,front_fy, bottom_fx,bottom_fy, top_fx,top_fy]
Returns:
u_hat_B, u_hat_C, u_hat_T: nondim streamwise velocity (bottom/centre/top)
v_hat_B, v_hat_C, v_hat_T: nondim crosswise velocity
Cd_F, Cd_T, Cd_B: drag coefficient per cylinder
Cl_F, Cl_T, Cl_B: lift coefficient per cylinder
"""
s = np.asarray(sensors, dtype=np.float64)
f = np.asarray(forces, dtype=np.float64)
# Sensor positions: s0=top, s1=centre, s2=bottom
# Convention: B=bottom=s2, C=centre=s1, T=top=s0
return {
"u_hat_T": s[:, 0] / u0,
"v_hat_T": s[:, 1] / u0,
"u_hat_C": s[:, 2] / u0,
"v_hat_C": s[:, 3] / u0,
"u_hat_B": s[:, 4] / u0,
"v_hat_B": s[:, 5] / u0,
"Cd_F": 2.0 * f[:, 0] / (rho * u0**2 * d),
"Cl_F": 2.0 * f[:, 1] / (rho * u0**2 * d),
"Cd_B": 2.0 * f[:, 2] / (rho * u0**2 * d),
"Cl_B": 2.0 * f[:, 3] / (rho * u0**2 * d),
"Cd_T": 2.0 * f[:, 4] / (rho * u0**2 * d),
"Cl_T": 2.0 * f[:, 5] / (rho * u0**2 * d),
}
# -- G operator (corrected) --------------------------------------------------
def apply_G_alpha(alpha: np.ndarray) -> np.ndarray:
"""Apply mirror G to action: [aF, aT, aB] -> [-aF, -aB, -aT]."""
return np.array([-alpha[0], -alpha[2], -alpha[1]], dtype=alpha.dtype)
def apply_G_x(dim: Dict[str, np.ndarray],
a_prev: np.ndarray,
a_prev2: np.ndarray) -> Tuple[Dict, np.ndarray, np.ndarray]:
"""Apply G to dimensionless state.
Returns (G_dim, G_a_prev, G_a_prev2) with corrected sign rules.
"""
G_dim = {
"u_hat_B": dim["u_hat_T"], "u_hat_C": dim["u_hat_C"], "u_hat_T": dim["u_hat_B"],
"v_hat_B": -dim["v_hat_T"], "v_hat_C": -dim["v_hat_C"], "v_hat_T": -dim["v_hat_B"],
"Cd_F": dim["Cd_F"], "Cd_T": dim["Cd_B"], "Cd_B": dim["Cd_T"],
"Cl_F": -dim["Cl_F"], "Cl_T": -dim["Cl_B"], "Cl_B": -dim["Cl_T"],
}
G_a_prev = np.column_stack([-a_prev[:, 0], -a_prev[:, 2], -a_prev[:, 1]])
G_a_prev2 = np.column_stack([-a_prev2[:, 0], -a_prev2[:, 2], -a_prev2[:, 1]])
return G_dim, G_a_prev, G_a_prev2
# -- Feature key definitions -------------------------------------------------
CORE_FEAT_KEYS = [
"u_m", "u_a", "u_c",
"v_a",
"Cd_tot", "Cd_rear",
"Cl_tot", "Cl_diff",
"sin_ua", "cos_ua",
"aF_lag1", "aB_lag1", "aT_lag1",
"daF", "daB", "daT",
]
MU_FEAT_KEYS = ["mu", "mu_u_a", "mu_v_a", "mu_Cd_tot", "mu_Cl_diff"]
ALL_FEAT_KEYS = CORE_FEAT_KEYS + MU_FEAT_KEYS
# -- Feature computation -----------------------------------------------------
def compute_features(
dim: Dict[str, np.ndarray],
actions_prev: np.ndarray, # (T, 3) physical omega(t-1) or nondim alpha(t-1)
actions_prev2: np.ndarray, # (T, 3) physical omega(t-2)
mu: float,
alpha_mode: bool = False, # if True, actions_prev are already nondim alpha
include_mu: bool = True,
u0: float = U0, # inlet velocity for omega->alpha conversion
) -> Dict[str, np.ndarray]:
"""Compute unified feature dictionary from dimensionless primitives.
Args:
dim: from compute_dimensionless()
actions_prev: lagged actions (physical omega or nondim alpha)
actions_prev2: twice-lagged actions
mu: 1/Re_D
alpha_mode: if True, actions are already nondim; else convert
include_mu: include mu modulation terms
u0: inlet velocity (lattice), used only when alpha_mode=False
Returns dict with all features as (T,) or (T,3) arrays.
"""
T = actions_prev.shape[0]
u_B, u_C, u_T = dim["u_hat_B"], dim["u_hat_C"], dim["u_hat_T"]
v_B, v_C, v_T = dim["v_hat_B"], dim["v_hat_C"], dim["v_hat_T"]
Cd_F, Cd_T, Cd_B = dim["Cd_F"], dim["Cd_T"], dim["Cd_B"]
Cl_F, Cl_T, Cl_B = dim["Cl_F"], dim["Cl_T"], dim["Cl_B"]
# If actions are in physical omega, convert to nondim alpha
if alpha_mode:
a = actions_prev.astype(np.float64)
a2 = actions_prev2.astype(np.float64)
else:
a = actions_prev.astype(np.float64) / u0
a2 = actions_prev2.astype(np.float64) / u0
sym = {}
# Sensor combinations (nondim)
sym["u_m"] = (u_B + u_C + u_T) / 3.0
sym["u_a"] = (u_T - u_B) / 2.0
sym["u_c"] = u_C.copy()
sym["v_a"] = (v_T - v_B) / 2.0
# Force combinations (dimensionless Cd/Cl)
sym["Cd_tot"] = Cd_F + Cd_T + Cd_B
sym["Cd_rear"] = Cd_T + Cd_B
sym["Cl_tot"] = Cl_F + Cl_T + Cl_B
sym["Cl_diff"] = Cl_T - Cl_B
# Phase
sym["sin_ua"] = np.sin(np.pi * sym["u_a"])
sym["cos_ua"] = np.cos(np.pi * sym["u_a"])
# Memory (nondim alpha)
sym["aF_lag1"] = a[:, 0]
sym["aB_lag1"] = a[:, 1]
sym["aT_lag1"] = a[:, 2]
sym["daF"] = a[:, 0] - a2[:, 0]
sym["daB"] = a[:, 1] - a2[:, 1]
sym["daT"] = a[:, 2] - a2[:, 2]
# Mu modulation
if include_mu:
sym["mu"] = np.full(T, mu, dtype=np.float64)
sym["mu_u_a"] = sym["u_a"] * mu
sym["mu_v_a"] = sym["v_a"] * mu
sym["mu_Cd_tot"] = sym["Cd_tot"] * mu
sym["mu_Cl_diff"] = sym["Cl_diff"] * mu
return sym
def build_feature_matrix(
sym: Dict[str, np.ndarray],
feat_keys: List[str],
add_bias: bool = True,
) -> np.ndarray:
"""Build feature matrix (T, N) from symbol dict."""
cols = []
if add_bias:
cols.append(np.ones(sym[feat_keys[0]].shape[0], dtype=np.float64))
for k in feat_keys:
if k in sym:
cols.append(np.asarray(sym[k], dtype=np.float64))
else:
# Missing key -> zero
T = sym.get("u_m", np.ones(1)).shape[0]
cols.append(np.zeros(T, dtype=np.float64))
return np.column_stack(cols)
def get_feature_names(feat_keys: List[str], add_bias: bool = True) -> List[str]:
"""Get feature names matching build_feature_matrix output."""
names = []
if add_bias:
names.append("bias")
names.extend(feat_keys)
return names
src/SR_analysis/utils/g_operator.py
+191
View File
@@ -0,0 +1,191 @@
"""G-operator and equivariance tools.
Provides G-operator transformations, dimensionless conversion,
and equivariance diagnostics for PPO control laws.
"""
from __future__ import annotations
from typing import Any, Dict, Optional, Tuple
import numpy as np
from .feature_builder import compute_dimensionless as _compute_dimless
def apply_G_alpha(alpha: np.ndarray) -> np.ndarray:
"""Apply mirror G to action: [aF, aT, aB] -> [-aF, -aB, -aT]."""
return np.array([-alpha[0], -alpha[2], -alpha[1]], dtype=alpha.dtype)
def apply_G_raw(obs_slice: np.ndarray,
a_prev: np.ndarray,
a_prev2: np.ndarray) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
"""Apply G to raw obs slice [sensor(6)+force(6)] and action arrays.
Parameters
----------
obs_slice : (12,) raw obs [s0_ux,uy, s1_ux,uy, s2_ux,uy, f0_fx,fy, f1_fx,fy, f2_fx,fy]
a_prev : (3,) physical omega at t-1
a_prev2 : (3,) physical omega at t-2
Returns
-------
G_obs : (12,) transformed obs slice
G_a_prev : (3,) transformed a_prev
G_a_prev2 : (3,) transformed a_prev2
"""
G_obs = np.zeros(12, dtype=np.float64)
# sensors: swap top(0,1) <-> bottom(4,5), negate v
G_obs[0] = obs_slice[4]
G_obs[1] = -obs_slice[5]
G_obs[2] = obs_slice[2]
G_obs[3] = -obs_slice[3]
G_obs[4] = obs_slice[0]
G_obs[5] = -obs_slice[1]
# forces: swap bottom(2,3) <-> top(4,5), negate fy
G_obs[6] = obs_slice[6]
G_obs[7] = -obs_slice[7]
G_obs[8] = obs_slice[10]
G_obs[9] = -obs_slice[11]
G_obs[10] = obs_slice[8]
G_obs[11] = -obs_slice[9]
G_a_prev = np.array([-a_prev[0], -a_prev[2], -a_prev[1]], dtype=np.float64)
G_a_prev2 = np.array([-a_prev2[0], -a_prev2[2], -a_prev2[1]], dtype=np.float64)
return G_obs, G_a_prev, G_a_prev2
def check_equivariance(
model: Any,
obs_slice_series: np.ndarray, # (T, 12) raw obs
actions_phys: np.ndarray, # (T, 3) physical omega
norm: dict,
action_scale: float = 8.0,
action_bias: Tuple[float, float, float] = (0.0, -4.0, 4.0),
u0: float = 0.01,
) -> Dict[str, float]:
"""Check G-equivariance of a PPO model over a time series.
Returns dict with front/rear equivariance errors.
"""
from .cfd_interface import build_observation, action_to_physical
T = min(obs_slice_series.shape[0], actions_phys.shape[0])
ef, eb, et = [], [], []
for t in range(2, T):
# Get current obs
osl = obs_slice_series[t]
a_prev = actions_phys[t - 1] if t > 0 else actions_phys[0]
a_prev2 = actions_phys[t - 2] if t > 1 else actions_phys[0]
# Predict action for current state
obs = build_observation(osl, norm)
act, _ = model.predict(obs, deterministic=True)
act = act.astype(np.float32).flatten()
alpha = action_to_physical(act.reshape(1, 3),
scale=action_scale, bias=action_bias, u0=u0).flatten()
# Apply G to state
G_obs, _, _ = apply_G_raw(osl, a_prev, a_prev2)
obs_G = build_observation(G_obs, norm)
act_G, _ = model.predict(obs_G, deterministic=True)
act_G = act_G.astype(np.float32).flatten()
alpha_G = action_to_physical(act_G.reshape(1, 3),
scale=action_scale, bias=action_bias, u0=u0).flatten()
# Expected: G(alpha) = [-aF, -aB, -aT]
expected = apply_G_alpha(alpha)
ef.append(abs(float(alpha_G[0]) - float(expected[0])))
eb.append(abs(float(alpha_G[1]) - float(expected[1])))
et.append(abs(float(alpha_G[2]) - float(expected[2])))
ef_arr = np.array(ef)
eb_arr = np.array(eb)
et_arr = np.array(et)
alpha_range = float(np.max(np.abs(actions_phys[2:])))
return {
"front_mean_abs_error": float(np.mean(ef_arr)),
"front_rel_error": float(np.mean(ef_arr) / (alpha_range + 1e-12)),
"rear_bottom_rel_error": float(np.mean(eb_arr) / (alpha_range + 1e-12)),
"rear_top_rel_error": float(np.mean(et_arr) / (alpha_range + 1e-12)),
"alpha_range": alpha_range,
}
def diagnose_one_re(model, ff, target_states, norm, config, n_steps=150) -> dict:
"""Run PPO inference and check equivariance.
Parameters
----------
model : loaded PPO model
ff : FlowField instance (must be at saved checkpoint state)
target_states : (FIFO_LEN, 6) target sensor signals
norm : norm dict
config : scene config dict with action_scale, action_bias, u0, etc.
Returns
-------
dict with equivariance metrics.
"""
from collections import deque
from .cfd_interface import (build_observation, scale_action,
action_to_physical, compute_similarity)
action_scale = config.get("action_scale", 8.0)
action_bias = config.get("action_bias", (0.0, -4.0, 4.0))
u0 = config.get("u0", 0.01)
sample_interval = config.get("sample_interval", 800)
fifo_len = config.get("fifo_len", 150)
n_obj_total = config.get("n_objects_total", 7)
ff.restore_ddf()
ff.apply_ddf()
# Bias FIFO init
fifo = deque(maxlen=fifo_len)
bias_arr = scale_action(np.zeros(3, dtype=np.float32),
scale=action_scale, bias=action_bias,
u0=u0, n_total_bodies=n_obj_total)
for _ in range(fifo_len):
ff.run(sample_interval, bias_arr)
fifo.append(ff.obs.copy()[2:14])
# Inference
obs_array = []
action_array = []
obs = np.zeros(12, dtype=np.float32)
for _ in range(n_steps):
act, _ = model.predict(obs, deterministic=True)
act = act.astype(np.float32).flatten()
action_array.append(act.copy())
action_arr = scale_action(act, scale=action_scale, bias=action_bias,
u0=u0, n_total_bodies=n_obj_total)
ff.context.push()
ff.run(sample_interval, action_arr)
ff.context.pop()
obs_slice = ff.obs.copy()[2:14]
fifo.append(obs_slice)
obs_array.append(obs_slice)
obs = build_observation(obs_slice, norm)
obs_series = np.array(obs_array, dtype=np.float64)
actions_phys = action_to_physical(np.array(action_array),
scale=action_scale, bias=action_bias, u0=u0)
states_arr = np.array(list(fifo), dtype=np.float32)
sim = compute_similarity(target_states, states_arr[:, 0:6],
config.get("conv_len", 30))
# Equivariance check
eq = check_equivariance(model, obs_series, actions_phys, norm,
action_scale, action_bias, u0)
return {
"similarity": sim,
"equivariance": eq,
}
src/SR_analysis/utils/sindy_fitter.py
+181
View File
@@ -0,0 +1,181 @@
"""SINDy fitting utilities: STLSQ threshold grid, feature matrix building.
All features are built using the unified feature_builder module.
"""
from __future__ import annotations
from typing import Dict, List, Optional, Tuple
import numpy as np
from .feature_builder import (
compute_dimensionless, compute_features, build_feature_matrix,
get_feature_names, ALL_FEAT_KEYS, U0,
)
# Default thresholds used across all scenes
DEFAULT_THRESHOLDS = [0.0, 0.001, 0.002, 0.005, 0.01, 0.015, 0.02, 0.03, 0.05, 0.1]
def fit_channel(
Theta: np.ndarray,
y: np.ndarray,
thresholds: Optional[List[float]] = None,
alpha: float = 1e-4,
max_iter: int = 25,
) -> Tuple[List[dict], dict]:
"""Fit a single channel (one cylinder) with STLSQ threshold grid.
Returns
-------
rows : list of dict per threshold
best : dict with best threshold entry (highest R2)
"""
import pysindy as ps
if thresholds is None:
thresholds = DEFAULT_THRESHOLDS
# Normalise features for thresholding stability
std = np.std(Theta, axis=0)
std = np.where(std < 1e-8, 1.0, std)
Theta_s = Theta / std
best = None
rows = []
for th in thresholds:
opt = ps.STLSQ(threshold=th, alpha=alpha, max_iter=max_iter)
opt.fit(Theta_s, y)
coef = np.asarray(opt.coef_, dtype=np.float64).flatten() / std
y_pred = Theta @ coef
ssr = float(np.sum((y - y_pred) ** 2))
sst = float(np.sum((y - np.mean(y)) ** 2) + 1e-12)
r2 = 1.0 - ssr / sst
mae = float(np.mean(np.abs(y - y_pred)))
nz = int(np.sum(np.abs(coef) > 1e-8))
entry = {"threshold": float(th), "nz": nz, "r2": r2, "mae": mae, "coef": coef}
rows.append(entry)
if best is None or r2 > best["r2"]:
best = entry
return rows, best
def fit_sindy(
Theta: np.ndarray,
y: np.ndarray,
thresholds: Optional[List[float]] = None,
) -> List[dict]:
"""Run SINDy with threshold grid, return results list.
Each result dict has keys: threshold, nz, r2, mae, coef.
"""
if thresholds is None:
thresholds = DEFAULT_THRESHOLDS
std = np.std(Theta, axis=0)
std = np.where(std < 1e-8, 1.0, std)
Theta_s = Theta / std
results = []
for th in thresholds:
import pysindy as ps
opt = ps.STLSQ(threshold=th, alpha=1e-4, max_iter=25)
opt.fit(Theta_s, y)
coef = np.asarray(opt.coef_, dtype=np.float64).flatten() / std
y_pred = Theta @ coef
ssr = float(np.sum((y - y_pred) ** 2))
sst = float(np.sum((y - np.mean(y)) ** 2) + 1e-12)
r2 = 1.0 - ssr / sst
mae = float(np.mean(np.abs(y - y_pred)))
nz = int(np.sum(np.abs(coef) > 1e-8))
results.append({
"threshold": float(th), "nz": nz, "r2": r2,
"mae": mae, "coef": [float(c) for c in coef],
})
return results
def print_control_law(feature_names: List[str], coef: np.ndarray, channel_label: str = "ch"):
"""Pretty-print a sparse control law."""
terms = []
for i, c in enumerate(coef):
if abs(c) > 1e-8:
terms.append(f"{c:.6f} * {feature_names[i]}")
print(f" {channel_label}: {' + '.join(terms)}")
nz = sum(1 for c in coef if abs(c) > 1e-8)
print(f" non-zero terms: {nz}")
def get_active_support(
coef: np.ndarray,
feat_names: List[str],
relative_threshold: float = 0.02,
) -> Dict[str, float]:
"""Extract active features from coefficient vector.
Features with |coef| / max(|coef|) >= relative_threshold are considered active.
"""
max_c = np.max(np.abs(coef))
if max_c < 1e-12:
return {}
active = {}
for name, c in zip(feat_names, coef):
if abs(c) / max_c >= relative_threshold:
active[name] = float(c)
return active
def get_feature_matrix_from_data(
sensors: np.ndarray, # (T, 6)
forces: np.ndarray, # (T, 6)
actions_phys: np.ndarray, # (T, 3) physical omega
mu: float,
u0: float = U0,
alpha_mode: bool = False,
include_mu: bool = True,
n_warmup: int = 2,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray, List[str], List[str]]:
"""Build feature matrices from raw CFD data.
Constructs dimensionless features via feature_builder, creates front (no bias)
and rear (with bias) feature matrices, and returns them aligned with Y.
Parameters
----------
sensors, forces, actions_phys : raw data arrays.
mu : 1/Re_D.
u0 : inlet velocity (lattice units).
alpha_mode : if True, actions_phys are already nondim alpha.
include_mu : include mu modulation features.
n_warmup : number of warmup steps to discard (default 2 for lag/da).
Returns
-------
Theta_front : (T-warmup, N_front) feature matrix, NO bias column
Theta_rear : (T-warmup, N_rear) feature matrix, WITH bias column
Y : (T-warmup, 3) target action matrix
feat_names_front : list of feature names for front
feat_names_rear : list of feature names for rear
"""
T = sensors.shape[0]
a_prev = np.zeros((T, 3), dtype=np.float64)
a_prev2 = np.zeros((T, 3), dtype=np.float64)
a_prev[1:] = actions_phys[:-1]
a_prev2[2:] = actions_phys[:-2]
dim = compute_dimensionless(sensors, forces, u0=u0, d=20.0)
sym = compute_features(dim, a_prev, a_prev2, mu,
alpha_mode=alpha_mode, include_mu=include_mu, u0=u0)
Theta_f = build_feature_matrix(sym, ALL_FEAT_KEYS, add_bias=False)
Theta_r = build_feature_matrix(sym, ALL_FEAT_KEYS, add_bias=True)
feat_names_front = get_feature_names(ALL_FEAT_KEYS, add_bias=False)
feat_names_rear = get_feature_names(ALL_FEAT_KEYS, add_bias=True)
return (Theta_f[n_warmup:], Theta_r[n_warmup:],
actions_phys[n_warmup:],
feat_names_front, feat_names_rear)
src/SR_analysis/validate/run_closed_loop.py
+320
View File
@@ -0,0 +1,320 @@
"""Unified closed-loop validator for SINDy control laws.
Validates a SINDy-derived control law by running it in a closed-loop CFD
environment and measuring similarity to the target flow.
Two modes:
- v23 (default): front no-bias + rear shared-head [bottom=-top(Gx)]
- unstructured: front with bias + rear independent
Usage:
conda run -n pycuda_3_10 python validate/run_closed_loop.py \\
--scene karman_re70 --device 2 \\
--sindy-results sindy/karman/sindy_results.json
"""
from __future__ import annotations
import argparse
import json
import os
import sys
from collections import deque
from typing import Any, Dict, List, Optional, 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 LegacyCelerisLab import FlowField # noqa: E402
from SR_analysis.utils.cfd_interface import (
nu_from_re, load_legacy_configs, build_karman_cloak_env, add_pinball,
scale_action, action_to_physical, compute_similarity,
load_ppo_model,
)
from SR_analysis.utils.sindy_fitter import get_feature_matrix_from_data
from SR_analysis.utils.feature_builder import (
compute_dimensionless, compute_features, build_feature_matrix,
apply_G_x, ALL_FEAT_KEYS, U0,
)
from SR_analysis.utils.g_operator import apply_G_raw
from SR_analysis.configs import (
get_scene, SCENES, LEGACY_CFG_DIR, FIFO_LEN, CONV_LEN,
)
DATA_TYPE = np.float32
def build_feature_vector(
obs_slice: np.ndarray,
a_prev: np.ndarray,
a_prev2: np.ndarray,
mu: float,
u0: float,
add_bias: bool,
) -> np.ndarray:
"""Build a single-row feature vector from raw obs and action state.
Matches the feature_builder logic but for a single time step.
"""
sensors = obs_slice[0:6].astype(np.float64).reshape(1, 6)
forces = obs_slice[6:12].astype(np.float64).reshape(1, 6)
ap = a_prev.astype(np.float64).reshape(1, 3)
ap2 = a_prev2.astype(np.float64).reshape(1, 3)
dim = compute_dimensionless(sensors, forces, u0=u0, d=20.0)
sym = compute_features(dim, ap, ap2, mu, alpha_mode=False, include_mu=True, u0=u0)
feat = build_feature_matrix(sym, ALL_FEAT_KEYS, add_bias=add_bias)
return feat[0] # single row
def predict_v23(
obs_slice: np.ndarray,
a_prev: np.ndarray,
a_prev2: np.ndarray,
mu: float,
u0: float,
front_coef: np.ndarray,
top_coef: np.ndarray,
feat_names_front: List[str],
feat_names_rear: List[str],
) -> np.ndarray:
"""Predict actions using v23: front no-bias + rear shared-head.
Returns (3,) physical omega array: [front, bottom, top].
"""
# Front channel: no bias
front = float(np.dot(
build_feature_vector(obs_slice, a_prev, a_prev2, mu, u0, add_bias=False),
front_coef))
# Top channel: with bias
top = float(np.dot(
build_feature_vector(obs_slice, a_prev, a_prev2, mu, u0, add_bias=True),
top_coef))
# Bottom = -top(Gx) using shared-head
G_obs, G_a_prev, G_a_prev2 = apply_G_raw(obs_slice, a_prev, a_prev2)
bottom = -float(np.dot(
build_feature_vector(G_obs, G_a_prev, G_a_prev2, mu, u0, add_bias=True),
top_coef))
return np.array([front, bottom, top], dtype=np.float64)
def predict_unstructured(
obs_slice: np.ndarray,
a_prev: np.ndarray,
a_prev2: np.ndarray,
mu: float,
u0: float,
front_coef: np.ndarray,
bottom_coef: np.ndarray,
top_coef: np.ndarray,
feat_names: List[str],
) -> np.ndarray:
"""Predict actions using unstructured: each channel independent with bias."""
feat = build_feature_vector(obs_slice, a_prev, a_prev2, mu, u0, add_bias=True)
front = float(np.dot(feat, front_coef))
bottom = float(np.dot(feat, bottom_coef))
top = float(np.dot(feat, top_coef))
return np.array([front, bottom, top], dtype=np.float64)
def load_sindy_coefs(sindy_path: str, scene_name: str) -> Dict[str, Any]:
"""Load SINDy coefficients for a scene from results JSON.
Returns dict with keys: front_coef, top_coef, bottom_coef,
feat_names_front, feat_names_rear, front_bias_mode.
"""
with open(sindy_path) as f:
data = json.load(f)
per = data["per_scene"].get(scene_name)
if per is None:
raise KeyError(f"Scene {scene_name} not found in {sindy_path}")
fn_f = per["feature_names_front"]
fn_r = per["feature_names_rear"]
front_coef = np.array(per["front"]["best_coef"], dtype=np.float64)
top_coef = np.array(per["top"]["best_coef"], dtype=np.float64)
bottom_coef = np.array(per["bottom"]["best_coef"], dtype=np.float64)
# Detect if front was fitted with bias (fn_f has "bias") or without
front_has_bias = fn_f[0] == "bias" if len(fn_f) > 0 else False
return {
"front_coef": front_coef,
"top_coef": top_coef,
"bottom_coef": bottom_coef,
"feat_names_front": fn_f,
"feat_names_rear": fn_r,
"front_has_bias": front_has_bias,
}
def run_validation(
scene_name: str,
coefs: Dict[str, Any],
device_id: int,
n_steps: int = 100,
mode: str = "v23",
) -> dict:
"""Run closed-loop validation using a SINDy control law.
Parameters
----------
scene_name : e.g. "karman_re70"
coefs : dict from load_sindy_coefs()
device_id : GPU device
n_steps : number of closed-loop steps
mode : "v23" (rear shared-head) or "unstructured"
Returns dict with similarity, actions range, etc.
"""
cfg = get_scene(scene_name)
re_code = cfg["re_code"]
nu = cfg["nu"]
u0 = cfg["u0"]
mu = cfg["mu"]
l0 = 20.0
sample_interval = cfg["sample_interval"]
action_scale = cfg["action_scale"]
action_bias = cfg["action_bias"]
n_obj_total = cfg["n_objects_env"]
print(f"\n=== Validating {scene_name} (mode={mode}, device={device_id}) ===")
# Build environment
cuda_cfg, field_cfg = load_legacy_configs(LEGACY_CFG_DIR)
field_cfg = field_cfg._replace(viscosity=float(nu))
ff = FlowField(field_cfg, cuda_cfg, device_id=device_id)
# Record target, add pinball, compute norm
target_states, _ = build_karman_cloak_env(
ff, u0=u0, l0=l0, sample_interval=sample_interval,
fifo_len=FIFO_LEN, data_type=DATA_TYPE,
)
norm = add_pinball(
ff, l0=l0, u0=u0, sample_interval=sample_interval,
fifo_len=FIFO_LEN, data_type=DATA_TYPE,
action_bias=action_bias, pinball_front_x=cfg["pinball_front_x"],
pinball_rear_x=cfg["pinball_rear_x"],
obs_slice_start=cfg["obs_slice"][0], obs_slice_end=cfg["obs_slice"][1],
)
# Reset to checkpoint
ff.restore_ddf()
ff.apply_ddf()
# Bias FIFO init
fifo = deque(maxlen=FIFO_LEN)
bias_arr = scale_action(np.zeros(3, dtype=np.float32), scale=action_scale,
bias=action_bias, u0=u0, n_total_bodies=n_obj_total)
for _ in range(FIFO_LEN):
ff.run(sample_interval, bias_arr)
fifo.append(ff.obs.copy()[2:14])
# Closed-loop with SINDy law
sens_list = []
actions_list = []
a_prev = action_to_physical(np.zeros((1, 3), dtype=np.float32),
scale=action_scale, bias=action_bias, u0=u0).flatten()
a_prev2 = a_prev.copy()
for _ in range(n_steps):
obs = fifo[-1] if fifo else np.zeros(12, dtype=np.float32)
# Predict using SINDy law
if mode == "v23":
omega = predict_v23(
obs, a_prev, a_prev2, mu, u0,
coefs["front_coef"], coefs["top_coef"],
coefs["feat_names_front"], coefs["feat_names_rear"])
elif mode == "unstructured":
omega = predict_unstructured(
obs, a_prev, a_prev2, mu, u0,
coefs["front_coef"], coefs["bottom_coef"], coefs["top_coef"],
coefs["feat_names_rear"])
else:
raise ValueError(f"Unknown mode: {mode}")
# Clip to valid action range
norm_a = (omega / u0 - np.array(action_bias, dtype=np.float64)) / action_scale
norm_a = np.clip(norm_a, -1.0, 1.0).astype(np.float32)
# Apply to CFD
action_arr = scale_action(norm_a, scale=action_scale, bias=action_bias,
u0=u0, n_total_bodies=n_obj_total)
ff.run(sample_interval, action_arr)
obs_new = ff.obs.copy()[2:14]
fifo.append(obs_new)
sens_list.append(obs_new[0:6])
actions_list.append(omega.copy())
a_prev2 = a_prev.copy()
a_prev = omega.copy()
# Evaluate
sens_arr = np.array(sens_list, dtype=np.float32)
actions_arr = np.array(actions_list, dtype=np.float64)
sim = compute_similarity(target_states, sens_arr, CONV_LEN)
action_range = float(np.max(np.abs(actions_arr)))
print(f" similarity={sim:.4f} action_range={action_range:.4f}")
del ff
return {
"scene": scene_name,
"mode": mode,
"similarity": sim,
"action_range": action_range,
"n_steps": n_steps,
}
def main():
ap = argparse.ArgumentParser(description="Closed-loop SINDy validation")
ap.add_argument("--scene", type=str, required=True, help="Scene name")
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"])
ap.add_argument("--sindy-results", type=str, default=None,
help="Path to sindy_results.json")
ap.add_argument("--out", type=str, default=None,
help="Output directory for result JSON")
args = ap.parse_args()
if args.sindy_results is None:
args.sindy_results = os.path.join(
os.path.dirname(__file__), "..", "sindy", "karman", "sindy_results.json")
coefs = load_sindy_coefs(args.sindy_results, args.scene)
result = run_validation(args.scene, coefs, args.device,
n_steps=args.steps, mode=args.mode)
if args.out is None:
out_dir = os.path.join(os.path.dirname(__file__), "..",
"validate", "results")
else:
out_dir = args.out
os.makedirs(out_dir, exist_ok=True)
out_path = os.path.join(out_dir, f"{args.scene}_{args.mode}.json")
with open(out_path, "w") as f:
json.dump(result, f, indent=2)
print(f"Saved: {out_path}")
if __name__ == "__main__":
main()