|
| 1 | +"""Numerical-correctness check of VaryingCoefSetup on the REAL 14-dataset |
| 2 | +inputs (the seeded-truth features/grids/g_func weights each system feeds the |
| 3 | +estimator). Verifies, per (system, variable): |
| 4 | +
|
| 5 | + A. G symmetric & finite (direct + super-Gram). |
| 6 | + B. super/from_full == direct construction (G, Phiy, yWy, score) -- the |
| 7 | + EqRPS fast path used in the live search. |
| 8 | + C. gamma_0 solves the weighted normal equations (|X^T W (y - X g0)| ~ 0) -- |
| 9 | + robust to collinearity, unlike comparing to a separate WLS solve. |
| 10 | + D. Var(gamma_0) matches sigma^2 * diag((X^T W X)^-1) (reported; the |
| 11 | + equilibration floor makes it approximate on near-singular blocks). |
| 12 | + E. score finite, >= 0, no NaN/Inf anywhere. |
| 13 | + F. Parseval (exact, basis property): var(beta(x)) == NC_raw per feature. |
| 14 | +
|
| 15 | +Usage: python _numcheck.py [system ...] (default: all 14) |
| 16 | +Run, then delete.""" |
| 17 | +from __future__ import annotations |
| 18 | +import os, sys, traceback |
| 19 | + |
| 20 | +_THIS = os.path.dirname(os.path.abspath(__file__)) |
| 21 | +_ROOT = os.path.abspath(os.path.join(_THIS, '..', '..')) |
| 22 | +for _p in (_ROOT, _THIS): |
| 23 | + if _p not in sys.path: |
| 24 | + sys.path.insert(0, _p) |
| 25 | + |
| 26 | +import numpy as np |
| 27 | +import yaml |
| 28 | +import epde.globals as gv |
| 29 | +from epde.operators.common.stability import VaryingCoefSetup as VC |
| 30 | +from kdv_sindy_test import build_pool_only, _normalize_grid_labels |
| 31 | +from thesis_runner import load_config, pipeline_settings, _set_seeds |
| 32 | +from vcoef_stat_compare import _ALL |
| 33 | +from epde.interface.equation_translator import translate_equation |
| 34 | + |
| 35 | + |
| 36 | +def _inputs(system): |
| 37 | + """[(var, Z (N,n_terms), target_idx, w (N,), grid_shape), ...] for truth.""" |
| 38 | + cfg = load_config(system) |
| 39 | + _set_seeds(0) |
| 40 | + gv.set_gram_config('vcoef') |
| 41 | + search = build_pool_only(cfg, pipeline_settings('new')) |
| 42 | + coords, data, variable_names, dim = cfg.load_data() |
| 43 | + all_vars = list(variable_names) |
| 44 | + truth = yaml.safe_load(open(os.path.join(_THIS, 'configs', f'{system}.yaml'))) |
| 45 | + teqs = truth.get('truth_equations') or [] |
| 46 | + seeded = (teqs[0] if len(all_vars) == 1 |
| 47 | + else {v: teqs[i] for i, v in enumerate(all_vars)}) |
| 48 | + seeded = _normalize_grid_labels(seeded) |
| 49 | + soeq = translate_equation(seeded, search.pool, all_vars=all_vars) |
| 50 | + out = [] |
| 51 | + for v in all_vars: |
| 52 | + eq = soeq.vals[v] |
| 53 | + eq.main_var_to_explain = v |
| 54 | + eq.weights_internal = np.ones(len(eq.structure) - 1) |
| 55 | + eq.weights_internal_evald = True |
| 56 | + eq.weights_final_evald = True |
| 57 | + eq.evaluate(normalize=False, return_val=False) # populate grid cache |
| 58 | + Z = np.vstack([t.evaluate(False, grids=None) |
| 59 | + for t in eq.structure]).T.astype(float) |
| 60 | + w = np.asarray(gv.grid_cache.g_func[gv.grid_cache.g_func_mask], float).reshape(-1) |
| 61 | + gshape = tuple(int(n) for n in gv.grid_cache.inner_shape) |
| 62 | + out.append((v, Z, int(eq.target_idx), w, gshape)) |
| 63 | + return out |
| 64 | + |
| 65 | + |
| 66 | +def _check(v, Z, tgt, w, gshape): |
| 67 | + N, n_terms = Z.shape |
| 68 | + feat_idx = [i for i in range(n_terms) if i != tgt] |
| 69 | + Xf = Z[:, feat_idx] |
| 70 | + yt = Z[:, tgt] |
| 71 | + direct = VC(Xf, yt, w, gshape, main_var=v, fit_intercept=True) |
| 72 | + sup = VC.precompute_super(Z, w, gshape, main_var=v) |
| 73 | + ff = VC.from_full(sup, tgt) |
| 74 | + |
| 75 | + m = {} |
| 76 | + # A. symmetric & finite |
| 77 | + G = direct.G |
| 78 | + m['Gsym'] = float(np.abs(G - G.T).max() / (np.abs(G).max() + 1e-30)) |
| 79 | + m['finite'] = bool(np.all(np.isfinite(G)) and np.all(np.isfinite(direct.Phiy)) |
| 80 | + and np.all(np.isfinite(sup['G_super']))) |
| 81 | + # B. super/from_full == direct |
| 82 | + m['dG'] = float(np.abs(ff.G - direct.G).max() / (np.abs(direct.G).max() + 1e-30)) |
| 83 | + m['dPhiy'] = float(np.abs(ff.Phiy - direct.Phiy).max() / (np.abs(direct.Phiy).max() + 1e-30)) |
| 84 | + m['dyWy'] = float(abs(ff.yWy - direct.yWy) / (abs(direct.yWy) + 1e-30)) |
| 85 | + sc_d = direct.score(None) |
| 86 | + sc_f = ff.score(None) |
| 87 | + m['dscore'] = float(np.abs(sc_f - sc_d).max()) |
| 88 | + # C. gamma_0 weighted normal equations |
| 89 | + sol = direct._solve_gammas(None) |
| 90 | + B = sol['B']; nf = sol['nf'] |
| 91 | + g0 = sol['gamma'][np.arange(nf) * B] # const per feature (incl intercept) |
| 92 | + Xa = np.column_stack([Xf, np.ones(N)]) |
| 93 | + XtWy = Xa.T @ (w * yt) |
| 94 | + resid = Xa.T @ (w * (yt - Xa @ g0)) |
| 95 | + m['normeq'] = float(np.abs(resid).max() / (np.abs(XtWy).max() + 1e-30)) |
| 96 | + # D. Var(gamma_0) vs sigma^2 diag((X^T W X)^-1) |
| 97 | + A = Xa.T @ (w[:, None] * Xa) |
| 98 | + Neff = float(w.sum()) |
| 99 | + rss = max(float(yt @ (w * yt) - g0 @ XtWy), 0.0) |
| 100 | + sigma2 = rss / max(Neff - nf, 1.0) |
| 101 | + try: |
| 102 | + var_ref = sigma2 * np.diag(np.linalg.inv(A)) |
| 103 | + var_vc = sol['var'][np.arange(nf) * B] |
| 104 | + rel = np.abs(var_vc - var_ref) / (np.abs(var_ref) + 1e-30) |
| 105 | + m['dVar'] = float(np.nanmax(rel)) |
| 106 | + except np.linalg.LinAlgError: |
| 107 | + m['dVar'] = float('nan') |
| 108 | + m['condA'] = float(np.linalg.cond(A)) |
| 109 | + # E. score sane |
| 110 | + m['score_ok'] = bool(np.all(np.isfinite(sc_d)) and np.all(sc_d >= -1e-9)) |
| 111 | + # F. Parseval exact |
| 112 | + st = direct.beta_field_stats(None) |
| 113 | + g = sol['gamma'] |
| 114 | + par = 0.0 |
| 115 | + for i in range(nf): |
| 116 | + nc_raw = float(np.sum(g[i * B + 1:(i + 1) * B] ** 2)) |
| 117 | + par = max(par, abs(st['std'][i] ** 2 - nc_raw) / (nc_raw + 1e-12)) |
| 118 | + m['parseval'] = float(par) |
| 119 | + return m |
| 120 | + |
| 121 | + |
| 122 | +# tolerances |
| 123 | +TOL = dict(Gsym=1e-10, dG=1e-7, dPhiy=1e-7, dyWy=1e-9, dscore=1e-6, |
| 124 | + normeq=1e-6, parseval=1e-6) |
| 125 | + |
| 126 | + |
| 127 | +def verdict(m): |
| 128 | + bad = [] |
| 129 | + if not m['finite']: |
| 130 | + bad.append('NONFINITE') |
| 131 | + if not m['score_ok']: |
| 132 | + bad.append('score') |
| 133 | + for k, t in TOL.items(): |
| 134 | + if not (m[k] <= t): |
| 135 | + bad.append(f'{k}={m[k]:.1e}') |
| 136 | + return bad |
| 137 | + |
| 138 | + |
| 139 | +def main(): |
| 140 | + systems = sys.argv[1:] or list(_ALL) |
| 141 | + n_ok = n_tot = 0 |
| 142 | + for s in systems: |
| 143 | + try: |
| 144 | + rows = _inputs(s) |
| 145 | + except Exception as e: |
| 146 | + print(f"{s:18s} INPUT-ERROR {type(e).__name__}: {str(e)[:60]}") |
| 147 | + continue |
| 148 | + for (v, Z, tgt, w, gshape) in rows: |
| 149 | + n_tot += 1 |
| 150 | + try: |
| 151 | + m = _check(v, Z, tgt, w, gshape) |
| 152 | + except Exception as e: |
| 153 | + print(f"{s:14s}/{v:3s} CHECK-ERROR {type(e).__name__}: {str(e)[:50]}") |
| 154 | + traceback.print_exc() |
| 155 | + continue |
| 156 | + bad = verdict(m) |
| 157 | + tag = f"{s}/{v}" if len(rows) > 1 else s |
| 158 | + if not bad: |
| 159 | + n_ok += 1 |
| 160 | + print(f"{tag:18s} OK superGmax={m['dG']:.0e} normeq={m['normeq']:.0e} " |
| 161 | + f"parseval={m['parseval']:.0e} dVar={m['dVar']:.0e} condA={m['condA']:.0e}") |
| 162 | + else: |
| 163 | + print(f"{tag:18s} FAIL {', '.join(bad)} (condA={m['condA']:.0e})") |
| 164 | + print(f"\n==== numerically correct: {n_ok}/{n_tot} (system,var) blocks ====") |
| 165 | + |
| 166 | + |
| 167 | +if __name__ == '__main__': |
| 168 | + sys.exit(main()) |
0 commit comments