1"""Toy verification of Rankine--Hugoniot front tokens on Burgers' equation.
  2
  3Run with the prescribed Python interpreter. Outputs a concise JSON report and
  4an optional CSV-like result file in this directory.
  5"""
  6import json, math
  7from pathlib import Path
  8import numpy as np
  9
 10SEED = 7
 11rng = np.random.default_rng(SEED)
 12
 13
 14def flux(u):
 15    return 0.5 * u * u
 16
 17
 18def rh_speed(ul, ur):
 19    den = ul - ur
 20    return (flux(ul) - flux(ur)) / den if abs(den) > 1e-12 else ul
 21
 22
 23def godunov_flux(ul, ur):
 24    # Exact scalar Burgers Godunov flux for the Riemann problem.
 25    if ul <= ur:  # rarefaction
 26        if ul >= 0: return flux(ul)
 27        if ur <= 0: return flux(ur)
 28        return 0.0
 29    # shock
 30    s = rh_speed(ul, ur)
 31    return flux(ul) if s >= 0 else flux(ur)
 32
 33
 34def godunov_step(u, dx, dt):
 35    # Fixed exterior states equal to the nearest interior value (shock stays away).
 36    ext = np.concatenate(([u[0]], u, [u[-1]]))
 37    f = np.array([godunov_flux(ext[i], ext[i+1]) for i in range(len(ext)-1)])
 38    return u - dt / dx * (f[1:] - f[:-1])
 39
 40
 41def exact_shock(x, xfront, ul=1.0, ur=0.0):
 42    return np.where(x < xfront, ul, ur)
 43
 44
 45def detect_front(u, x, threshold=0.35):
 46    jumps = np.abs(np.diff(u))
 47    ids = np.flatnonzero(jumps >= threshold)
 48    if len(ids) == 0:
 49        return []
 50    # For this toy, adjacent flagged edges are one cluster; code also handles many.
 51    clusters = []
 52    start = prev = int(ids[0])
 53    for q in ids[1:]:
 54        q = int(q)
 55        if q > prev + 1:
 56            clusters.append((start, prev)); start = q
 57        prev = q
 58    clusters.append((start, prev))
 59    out = []
 60    for a, b in clusters:
 61        # traces are pooled from the cells immediately outside the flagged cluster
 62        li = max(0, a)
 63        ri = min(len(u)-1, b+1)
 64        out.append({"x": float(0.5*(x[a]+x[b+1])),
 65                    "ul": float(u[li]), "ur": float(u[ri])})
 66    return out
 67
 68
 69def token_rollout(u0, x, dt, nsteps):
 70    """Propagate detected fronts analytically and render the piecewise field."""
 71    tokens = detect_front(u0, x)
 72    # In the intended module this is where attention/residual prediction enters;
 73    # for a clean shock the exact traces are already the sufficient local state.
 74    xs = []
 75    fields = []
 76    u = u0.copy()
 77    for _ in range(nsteps + 1):
 78        if tokens:
 79            t = tokens[0]
 80            fields.append(exact_shock(x, t["x"], t["ul"], t["ur"]))
 81            xs.append(t["x"])
 82        else:
 83            fields.append(u.copy()); xs.append(float("nan"))
 84        if tokens:
 85            for t in tokens:
 86                t["x"] += dt * rh_speed(t["ul"], t["ur"])
 87    return np.asarray(fields), np.asarray(xs)
 88
 89
 90def crossing_position(u, x, level=0.5):
 91    # Linear interpolation of the numerical shock crossing.
 92    k = np.flatnonzero((u[:-1] >= level) & (u[1:] < level))
 93    if len(k) == 0: return float("nan")
 94    i = int(k[0]); den = u[i] - u[i+1]
 95    return float(x[i] + (level-u[i]) * (x[i+1]-x[i]) / den) if den else float(x[i])
 96
 97
 98def main():
 99    # Core mathematical check over random traces and x-independent flux.
100    pairs = rng.uniform(-1.0, 1.0, (100, 2))
101    pairs = pairs[np.abs(pairs[:,0]-pairs[:,1]) > .1]
102    speed_err = []
103    for ul, ur in pairs:
104        s = rh_speed(ul, ur)
105        speed_err.append(abs((flux(ul)-flux(ur))/(ul-ur)-s))
106    math_check = {"max_rh_identity_error": float(max(speed_err)),
107                  "mean_rh_identity_error": float(np.mean(speed_err)),
108                  "identity_pass": bool(max(speed_err) < 1e-12)}
109
110    # Clean shock benchmark at three resolutions, same physical horizon.
111    records = []
112    for N in (64, 128, 256):
113        dx = 1.0/N; x = (np.arange(N)+.5)*dx
114        x0 = .25; ul, ur = 1., 0.; s = rh_speed(ul, ur)
115        u0 = exact_shock(x, x0, ul, ur)
116        # Conservative baseline has a standard CFL timestep.
117        dt = .35*dx/max(abs(ul), abs(ur)); steps = int(.20/dt)
118        dt = .20/steps
119        ub = u0.copy(); base_pos=[]
120        for _ in range(steps):
121            ub = godunov_step(ub, dx, dt)
122            base_pos.append(crossing_position(ub, x))
123        tf, token_pos = token_rollout(u0, x, dt, steps)
124        times = np.arange(steps+1)*dt
125        true_pos = x0 + s*times
126        # compare only positions available at the grid output times
127        bpos = np.asarray(base_pos)
128        bpos_err = float(np.nanmean(np.abs(bpos-true_pos[1:])))
129        tpos_err = float(np.mean(np.abs(token_pos-true_pos)))
130        bfield = float(np.mean(np.abs(ub-exact_shock(x,true_pos[-1],ul,ur))))
131        tfield = float(np.mean(np.abs(tf[-1]-exact_shock(x,true_pos[-1],ul,ur))))
132        records.append({"N":N,"dt":dt,"steps":steps,
133                        "baseline_front_L1":bpos_err,"idea_front_L1":tpos_err,
134                        "baseline_field_L1":bfield,"idea_field_L1":tfield})
135
136    # One-step propagation error vs dt: discretization of the initial location is
137    # held fixed, so the analytic kinematic update should have no time integration error.
138    N=128; dx=1/N; x=(np.arange(N)+.5)*dx; x0=.25; u0=exact_shock(x,x0)
139    tok=detect_front(u0,x)[0]; xdet=tok['x']; s=rh_speed(tok['ul'],tok['ur'])
140    scales=[]
141    for dt in (0.002, 0.004, 0.008):
142        pred=xdet+dt*s; truth=x0+dt*s
143        scales.append({"dt":dt,"abs_position_error":abs(pred-truth)})
144    report={"seed":SEED,"math_check":math_check,"records":records,
145            "one_step_scaling":scales,
146            "note":"Front-token field uses exact local traces on a clean Riemann shock; no learned residual is needed in this deliberately minimal case."}
147    Path("results.json").write_text(json.dumps(report, indent=2))
148    print(json.dumps(report, indent=2))
149
150if __name__ == '__main__': main()