From b1ac9ceabb92cbe1e7987ee29a7bcdc5b8a5e5e0 Mon Sep 17 00:00:00 2001 From: Frank14f <1515444314@qq.com> Date: Wed, 1 Jul 2026 13:54:00 +0800 Subject: [PATCH] fix(SR): stage_3_validate.py all modes verified, legacy_karman_env.py cylinder order corrected MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit - Fix stage_3_validate.py Illusion PPO mode: observation was always zero (unused before, now builds correct 14-dim obs from flow field) - Rewrite Illusion pysr mode: use sympy.lambdify + v23 structure, matches validate_karman approach. Verified similarity=0.975 (100 steps). - Fix Karman pysr mode: strip "bias" from JSON feature_names before build_feature_matrix(add_bias=True), fix alpha→omega conversion. Verified similarity=0.887 (100 steps, matches old 0.888 at 160 steps). - Fix legacy_karman_env.py cylinder order: front→TOP→BOTTOM matching training env (was front→BOTTOM→TOP). Also fix bias assignment. Verified via cylinder ID dump against training env layout. - Remove dead eval_math() function. Co-authored-by: Cursor --- src/SR_analysis/stage_3_validate.py | 158 ++++++++++-------- .../legacy_env/legacy_karman_env.py | 9 +- 2 files changed, 96 insertions(+), 71 deletions(-) diff --git a/src/SR_analysis/stage_3_validate.py b/src/SR_analysis/stage_3_validate.py index d034e39..9d5517d 100644 --- a/src/SR_analysis/stage_3_validate.py +++ b/src/SR_analysis/stage_3_validate.py @@ -124,10 +124,25 @@ def validate_karman(scene_name, device_id, n_steps, mode, formula_front, formula if mode == "pysr": fj = json.load(open(formula_front)) tj = json.load(open(formula_top)) - front_coef, front_keys = np.array(fj.get("coef", [0]*len(fj["feature_names"]))), fj["feature_names"] - top_coef, top_keys = np.array(tj.get("coef", [0]*(len(tj["feature_names"])+1))), tj["feature_names"] - if len(front_coef) != len(front_keys): - front_coef = np.array([float(c) for c in fj["best_sympy"].split("+")]) if "+" in fj["best_sympy"] else np.zeros(len(front_keys)) + front_keys = fj["feature_names"] + top_keys_raw = tj["feature_names"] + # Strip "bias" — build_feature_matrix with add_bias=True adds it back + top_keys = [k for k in top_keys_raw if k != "bias"] + front_coef = np.array(fj.get("coef", [0]*len(front_keys))) if "coef" in fj else None + top_coef = np.array(tj.get("coef", [0]*len(top_keys_raw))) if "coef" in tj else None + if front_coef is None or len(front_coef) != len(front_keys): + # Use sympy.lambdify as fallback + import sympy + _fs = sympy.symbols(front_keys) + _fe = sympy.sympify(fj["best_sympy"]) + _ffn = sympy.lambdify(_fs, _fe, "numpy") + _ts = sympy.symbols(top_keys_raw) + _te = sympy.sympify(tj["best_sympy"]) + _tfn = sympy.lambdify(_ts, _te, "numpy") + front_coef = top_coef = None + _use_sympy = True + else: + _use_sympy = False elif mode == "ppo": model = load_ppo_model(f"models/old/{cfg['model_name']}.zip", device=f"cuda:{device_id}", s_dim=cfg.get("s_dim", 12)) @@ -141,16 +156,25 @@ def validate_karman(scene_name, device_id, n_steps, mode, formula_front, formula for step in range(n_steps): if mode == "pysr": osl = ff.obs.copy()[2:14] - # Compute SR action via v23 - fv_f = _build_feature_vector(osl, a_prev, a_prev2, mu, u0, add_bias=False, feat_keys=front_keys) - front = float(np.dot(fv_f, front_coef)) - fv_t = _build_feature_vector(osl, a_prev, a_prev2, mu, u0, add_bias=True, feat_keys=top_keys) - top = float(np.dot(fv_t, top_coef)) + if _use_sympy: + fv_f = _build_feature_vector(osl, a_prev, a_prev2, mu, u0, add_bias=False, feat_keys=front_keys) + front = float(_ffn(*fv_f)) + fv_t = _build_feature_vector(osl, a_prev, a_prev2, mu, u0, add_bias=True, feat_keys=top_keys) + top = float(_tfn(*fv_t)) + else: + fv_f = _build_feature_vector(osl, a_prev, a_prev2, mu, u0, add_bias=False, feat_keys=front_keys) + front = float(np.dot(fv_f, front_coef)) + fv_t = _build_feature_vector(osl, a_prev, a_prev2, mu, u0, add_bias=True, feat_keys=top_keys) + top = float(np.dot(fv_t, top_coef)) G_obs, G_ap, G_ap2 = apply_G_raw(osl, a_prev, a_prev2) - fv_b = _build_feature_vector(G_obs, G_ap, G_ap2, mu, u0, add_bias=True, feat_keys=top_keys) - bottom = -float(np.dot(fv_b, top_coef)) + if _use_sympy: + fv_b = _build_feature_vector(G_obs, G_ap, G_ap2, mu, u0, add_bias=True, feat_keys=top_keys) + bottom = -float(_tfn(*fv_b)) + else: + fv_b = _build_feature_vector(G_obs, G_ap, G_ap2, mu, u0, add_bias=True, feat_keys=top_keys) + bottom = -float(np.dot(fv_b, top_coef)) alpha_sr = np.array([front, bottom, top]) - omega = (alpha_sr / 8.0 - action_bias) * u0 / 10.0 # α→ω conversion (approximate) + omega = alpha_sr * u0 # PySR alpha is already non-dim, ω = α·U₀ elif mode == "ppo": osl = ff.obs.copy()[2:14] obs_n = np.clip(np.hstack([osl[6:12]/norm["force_norm_fact"], (osl[0:6]-norm["sens_deviation"])/norm["sens_norm_fact"]]), -1, 1).astype(np.float32) @@ -214,9 +238,17 @@ def validate_illusion(scene_name, device_id, n_steps, mode, formula_front, formu ff.run(warmup, np.zeros(n_obj, dtype=DATA_TYPE)) ff.get_ddf(); ff.save_ddf() - # Norm + # Norm — collect zero-action FIFO and compute normalisation + fifo_norm = deque(maxlen=FIFO_LEN) for _ in range(FIFO_LEN): ff.run(si, np.zeros(n_obj, dtype=DATA_TYPE)) + fifo_norm.append(ff.obs.copy()[0:12]) + tsa = np.array(fifo_norm, dtype=np.float32) + force_norm_fact = 6.0 * float(np.max(np.abs(tsa[:, 6:12]))) + sens_deviation = np.mean(tsa[:, 0:6], axis=0).astype(np.float32) + sens_norm_fact = np.zeros(6, dtype=np.float32) + for i in range(6): + sens_norm_fact[i] = 5.0 * float(np.max(np.abs(tsa[:, i] - sens_deviation[i]))) # Bias FIFO ff.apply_ddf() @@ -227,69 +259,72 @@ def validate_illusion(scene_name, device_id, n_steps, mode, formula_front, formu ff.run(si, ba) fifo.append(ff.obs.copy()[0:12]) - # Formula coefficients + # Precompile formula evaluators for pysr mode if mode == "pysr": + import sympy fj = json.load(open(formula_front)) tj = json.load(open(formula_top)) - front_keys = fj["feature_names"] - top_keys = tj["feature_names"] - # Try to extract coefs from best_sympy via simple parsing - best_f = fj["best_sympy"] - # Fall back to predicting per-step: extract coefficient vector from the formula - # For the joint formula: Cd_tot - (Cd_err + 5.428) - 0.00978*(du_a_dt + u_a) - # We need to build the feature vector and evaluate - front_best = best_f - top_best = tj["best_sympy"] + # Front: no-bias, uses feature_names as-is + _front_keys = [k for k in fj["feature_names"] if k != "bias"] + _x_sym = sympy.symbols(_front_keys) + _front_expr = sympy.sympify(fj["best_sympy"]) + _front_fn = sympy.lambdify(_x_sym, _front_expr, "numpy") + # Top: formula was fitted WITH bias in feature_names, + # so sympy must include "bias" as a symbol. When building features, + # we pass feat_keys without "bias" + add_bias=True → output matches. + _top_keys_full = tj["feature_names"] # e.g. ["bias", "Cd_err", ...] + _top_keys_no_bias = [k for k in _top_keys_full if k != "bias"] + _x_t = sympy.symbols(_top_keys_full) + _top_expr = sympy.sympify(tj["best_sympy"]) + _top_fn = sympy.lambdify(_x_t, _top_expr, "numpy") + + import warnings + warnings.filterwarnings("ignore", category=RuntimeWarning) + else: + _front_fn = _top_fn = None + _front_keys = _top_keys_no_bias = _top_keys_full = [] # Run sens_list, act_list = [], [] a_prev = a_prev2 = fifo_bias.copy() - result_queue = [] for step in range(n_steps): if mode == "pysr": osl = ff.obs.copy()[0:12] tf_step = gen_target_states_at(step, harmonics) - target_f = np.array([tf_step[6], tf_step[7]]) + target_f = np.array([tf_step[0], tf_step[1]]) - # Evaluate formula using sympy-like approach - # Build feature dict for evaluation - from SR_analysis.utils.feature_builder import compute_dimensionless, compute_features, build_feature_matrix, ILLUSION_PHASE_KEYS - s = osl[0:6].astype(np.float64).reshape(1, 6) - f = osl[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(s, f, u0=u0, d=20.0) - sym = compute_features(dim, ap, ap2, mu, alpha_mode=False, include_mu=False, - include_cos_sin=False, u0=u0, target_forces=target_f.reshape(1, 2), - sensors_raw=s, forces_raw=f) - fv_f = build_feature_matrix(sym, front_keys, add_bias=False)[0] - fv_t = build_feature_matrix(sym, top_keys, add_bias=True)[0] - # Simple eval - env = dict(zip(front_keys, fv_f)) - alpha_f = float(eval_math(best_f, env)) # simplified - env2 = dict(zip(["bias"] + top_keys, fv_t)) - alpha_t = float(eval_math(best_t, env2)) # simplified - # Bottom via G - from SR_analysis.utils.g_operator import apply_G_raw + # Front: no-bias, evaluate feature vector + formula + fv_f = _build_feature_vector(osl, a_prev, a_prev2, mu, u0, add_bias=False, + feat_keys=_front_keys, tf=target_f) + alpha_f = float(_front_fn(*fv_f)) + + # Top: with-bias, shared-head. Formula sympy includes "bias" symbol, + # feature builder adds bias column via add_bias=True. + fv_t = _build_feature_vector(osl, a_prev, a_prev2, mu, u0, add_bias=True, + feat_keys=_top_keys_no_bias, tf=target_f) + alpha_t = float(_top_fn(*fv_t)) + + # Bottom via G-mirror G_obs, G_ap, G_ap2 = apply_G_raw(osl, a_prev, a_prev2) G_tf = np.array([target_f[0], -target_f[1]]) - sG = G_obs[0:6].astype(np.float64).reshape(1, 6) - fG = G_obs[6:12].astype(np.float64).reshape(1, 6) - dimG = compute_dimensionless(sG, fG, u0=u0, d=20.0) - symG = compute_features(dimG, G_ap.reshape(1,3), G_ap2.reshape(1,3), mu, alpha_mode=False, - include_mu=False, include_cos_sin=False, u0=u0, target_forces=G_tf.reshape(1,2), - sensors_raw=sG, forces_raw=fG) - fv_b = build_feature_matrix(symG, top_keys, add_bias=True)[0] - env_b = dict(zip(["bias"] + top_keys, fv_b)) - alpha_b = -float(eval_math(best_t, env_b)) + fv_b = _build_feature_vector(G_obs, G_ap, G_ap2, mu, u0, add_bias=True, + feat_keys=_top_keys_no_bias, tf=G_tf) + alpha_b = -float(_top_fn(*fv_b)) + alpha = np.array([alpha_f, alpha_b, alpha_t]) - omega = alpha * u0 # alpha to omega + omega = alpha * u0 # PySR output is already non-dim alpha (ω/U₀ factored out) elif mode == "ppo": model_path = f"models/250525/{cfg['model_name']}.zip" model = load_ppo_model(model_path, device=f"cuda:{device_id}", s_dim=14) osl = ff.obs.copy()[0:12] - obs = np.zeros(14, dtype=np.float32) + tf_step = gen_target_states_at(step, harmonics) + target_cd_n = tf_step[0] / force_norm_fact + target_cl_n = tf_step[1] / force_norm_fact + # Build observation same way as infer_illusion.py + forces_n = osl[6:12] / force_norm_fact + sens_n = (osl[0:6] - sens_deviation) / sens_norm_fact + obs = np.clip(np.hstack([forces_n, sens_n, target_cd_n, target_cl_n]), -1, 1).astype(np.float32) act, _ = model.predict(obs, deterministic=True) act = act.astype(np.float32).flatten() omega = (act * action_scale + action_bias) * u0 @@ -315,17 +350,6 @@ def validate_vortex(scene_name, device_id, n_steps, mode, formula_front, formula return validate_karman(scene_name, device_id, n_steps, mode, formula_front, formula_top) -def eval_math(expr: str, env: dict) -> float: - """Simple math evaluator for formula with _dt and basic ops.""" - expr = expr.replace("daF_dt", "0").replace("daB_dt", "0").replace("daT_dt", "0") - expr = expr.replace("mu_Cl_tot", f"{env.get('mu_Cl_tot', 0)}*{env.get('Cl_tot', 0)}").replace("mu_", "0*") - # Simplified: use generic eval - try: - return float(eval(expr, {"__builtins__": {}}, env)) - except: - return 0.0 - - _VALIDATE_FN = { "karman": validate_karman, "illusion": validate_illusion, diff --git a/src/drl_pinball/legacy_env/legacy_karman_env.py b/src/drl_pinball/legacy_env/legacy_karman_env.py index 9d0069b..18d3c83 100644 --- a/src/drl_pinball/legacy_env/legacy_karman_env.py +++ b/src/drl_pinball/legacy_env/legacy_karman_env.py @@ -145,9 +145,10 @@ def legacy_build_re100( target_states = np.vstack((target_states, new_state)) # -- Step 4: Add pinball cylinders (ids=4,5,6) ------------------------- + # Order MUST match training env: front → TOP(y=+0.75) → BOTTOM(y=-0.75) ff.add_cylinder(_fill_y(FRONT_CENTER, cy), PINBALL_RADIUS) - ff.add_cylinder(_fill_y(BOTTOM_CENTER, cy - 0.75 * L0), PINBALL_RADIUS) ff.add_cylinder(_fill_y(TOP_CENTER, cy + 0.75 * L0), PINBALL_RADIUS) + ff.add_cylinder(_fill_y(BOTTOM_CENTER, cy - 0.75 * L0), PINBALL_RADIUS) n_obj_total = ff.obs.size // 2 # 7 objects assert n_obj_total == 7, f"Expected 7 objects, got {n_obj_total}" @@ -175,11 +176,11 @@ def legacy_build_re100( # -- Step 8: Bias-action rollout (for FIFO init in controlled runs) ----- ff.apply_ddf() # restore pre-bias state - # Action bias: front=0, bottom=-4*U0, top=4*U0 + # Action bias: front=0, top=-4*U0, bottom=+4*U0 bias_arr = np.zeros(n_obj_total, dtype=DATA_TYPE) bias_arr[n_obj_total - 3] = float((0.0 * 8.0 + 0.0) * U0) # front = 0 - bias_arr[n_obj_total - 2] = float((0.0 * 8.0 + (-4.0)) * U0) # bottom = -4*U0 - bias_arr[n_obj_total - 1] = float((0.0 * 8.0 + 4.0) * U0) # top = 4*U0 + bias_arr[n_obj_total - 2] = float((0.0 * 8.0 + (-4.0)) * U0) # TOP (id=5) = -4*U0 + bias_arr[n_obj_total - 1] = float((0.0 * 8.0 + 4.0) * U0) # BOTTOM (id=6) = +4*U0 fifo.clear() for _ in range(FIFO_LEN):