1import json
 2import numpy as np
 3from invariant_reconstruction import pressure, admissible, limited_state, density_theta_bound, GAMMA
 4
 5
 6def state(rho, vx, vy, p):
 7    return np.array([rho, rho*vx, rho*vy,
 8                     p/(GAMMA-1.0) + 0.5*rho*(vx*vx+vy*vy)])
 9
10
11def pressure_root(w0, wr, pf):
12    # The pressure-floor condition after multiplying by rho is quadratic in theta.
13    d = wr - w0
14    def f(t):
15        w = w0 + t*d
16        return w[3]*w[0] - .5*np.dot(w[1:3], w[1:3]) - pf/(GAMMA-1.0)*w[0]
17    y0, y1, y2 = f(0.), f(1.), f(2.)
18    a = .5*(y2 - 2*y1 + y0)
19    b = y1 - y0 - a
20    if abs(a) < 1e-13:
21        roots = [-y0/b] if abs(b) > 1e-13 else []
22    else:
23        roots = np.roots([a, b, y0])
24    good = [float(np.real(x)) for x in roots
25            if abs(np.imag(x)) < 1e-8 and np.real(x) >= 0]
26    return min(good) if good else np.inf
27
28
29def main():
30    rng = np.random.default_rng(1178)
31    rf, pf = 1e-6, 1e-6
32    w0 = state(1.0, .7, -.2, 1.0)
33    out = {"predictions": {}, "comparison": {}}
34
35    # Prediction 1: if pressure is safely high, the exact affine density bound is used.
36    density_rows = []
37    for rr in np.geomspace(1e-9, 0.9*rf, 8):
38        wr = state(rr, 0., 0., 1.0)
39        _, got = limited_state(w0, wr, rf, pf)
40        expected = density_theta_bound(w0, wr, rf)
41        density_rows.append([float(rr), float(got), float(expected), abs(got-expected)])
42    out["predictions"]["density_affine_bound"] = {
43        "prediction": "safe high-pressure endpoints satisfy theta=(rho0-rhofloor)/(rho0-rhoraw)",
44        "max_abs_error": float(max(r[-1] for r in density_rows)), "rows": density_rows}
45
46    # Prediction 2: pressure positivity is enforced at the first analytic quadratic root.
47    pressure_rows = []
48    for mach_momentum in [2., 4., 8., 16., 32.]:
49        # Negative thermodynamic pressure makes the raw endpoint inadmissible.
50        wr = state(1.0, mach_momentum, 0., -0.1)
51        _, got = limited_state(w0, wr, rf, pf)
52        expected = min(1.0, pressure_root(w0, wr, pf))
53        pressure_rows.append([mach_momentum, float(got), float(expected), abs(got-expected),
54                              float(pressure(wr))])
55    out["predictions"]["pressure_root"] = {
56        "prediction": "accepted theta equals the first pressure-floor crossing on a bad segment",
57        "max_abs_error": float(max(r[3] for r in pressure_rows)), "rows": pressure_rows}
58
59    # Prediction 3: scaling a fixed bad ray by Lambda gives theta proportional to 1/Lambda.
60    direction = state(0.2, 18., 0., -0.1) - w0
61    scales = np.array([1., 2., 4., 8., 16., 32.])
62    theta = np.array([limited_state(w0, w0 + lam*direction, rf, pf)[1]
63                      for lam in scales])
64    scaled = scales * theta
65    out["predictions"]["inverse_strength_scaling"] = {
66        "prediction": "for Lambda beyond the first crossing, Lambda*theta is constant",
67        "scales": scales.tolist(), "theta": theta.tolist(),
68        "lambda_theta": scaled.tolist(),
69        "tail_ratio": float(max(scaled[-3:])/min(scaled[-3:]))}
70
71    # Secondary fixed-seed comparison against first-order (cell-average) reconstruction.
72    n = 3000
73    cells, raws = [], []
74    for _ in range(n):
75        c = state(rng.uniform(.5, 1.5), rng.uniform(-1, 1), rng.uniform(-1, 1),
76                  rng.uniform(.3, 2.))
77        raws.append(c + rng.normal(0, 1.8, 4)); cells.append(c)
78    raw_invalid = sum(not admissible(r, rf, pf) for r in raws)
79    safe, ts = [], []
80    for c, r in zip(cells, raws):
81        s, t = limited_state(c, r, rf, pf); safe.append(s); ts.append(t)
82    out["comparison"] = {
83        "samples": n,
84        "raw_invalid_rate": raw_invalid/n,
85        "invariant_invalid_rate": sum(not admissible(s,rf,pf) for s in safe)/n,
86        "baseline_first_order_l2_to_raw": float(np.mean([np.linalg.norm(c-r) for c,r in zip(cells,raws)])),
87        "idea_l2_to_raw": float(np.mean([np.linalg.norm(s-r) for s,r in zip(safe,raws)])),
88        "idea_l2_to_cell": float(np.mean([np.linalg.norm(s-c) for s,c in zip(safe,cells)])),
89        "activation_rate": float(np.mean(np.array(ts) < .999999))}
90    print(json.dumps(out, indent=2))
91    with open("results.json", "w") as f: json.dump(out, f, indent=2)
92
93if __name__ == "__main__":
94    main()