1import json
 2from pathlib import Path
 3import numpy as np
 4
 5# Toy flow: dz/dt = v0 + L z, z(0)=x.  The pseudo-trajectory is an
 6# exact trajectory with elapsed time (1+a)t, hence it has a pure timing defect.
 7def flow(t, x, v0, L):
 8    t = np.asarray(t, dtype=float)
 9    if abs(L) < 1e-12:
10        return x + v0 * t
11    return (x + v0 / L) * np.exp(L * t) - v0 / L
12
13def vector_field(z, v0, L):
14    return v0 + L * z
15
16def sup_speed(T, x, v0, L, n=20001):
17    t = np.linspace(0.0, T, n)
18    return float(np.max(np.abs(vector_field(flow(t, x, v0, L), v0, L))))
19
20def standard_error(T, eps, a, x, v0, L, n=2001):
21    # For a>=0 and an increasing affine flow, the closest admissible h is
22    # h(t)=(1+min(eps,a))t. Its secant slopes are exactly in Rep(eps).
23    t = np.linspace(0.0, T, n)
24    h_slope = 1.0 + min(float(eps), float(a))
25    return float(np.max(np.abs(flow((1.0 + a) * t, x, v0, L)
26                         - flow(h_slope * t, x, v0, L))))
27
28def fixed_time_error(T, a, x, v0, L):
29    return standard_error(T, 0.0, a, x, v0, L)
30
31def oriented_error(T, a, x, v0, L):
32    # h(t)=(1+a)t is an unrestricted increasing homeomorphism.
33    return 0.0
34
35def empirical_threshold(T, rho, a, x, v0, L):
36    lo, hi = 0.0, max(a, 1e-12)
37    for _ in range(70):
38        mid = (lo + hi) / 2
39        if standard_error(T, mid, a, x, v0, L) <= rho:
40            hi = mid
41        else:
42            lo = mid
43    return hi
44
45def run():
46    x, v0, a, rho = 0.0, 1.0, 0.08, 0.02
47    rows = []
48    for L in (0.2, 0.5, 1.0):
49        for T in (0.5, 1.0, 2.0):
50            speed = sup_speed(T, x, v0, L)
51            pred = rho / (T * L * speed)
52            obs = empirical_threshold(T, rho, a, x, v0, L)
53            rows.append({"L": L, "T": T, "sup_speed": speed,
54                         "predicted_epsilon_crit": pred,
55                         "observed_epsilon_crit": obs,
56                         "ratio_observed_to_predicted": obs / pred})
57
58    # Direct verification of the mechanism over epsilon: bounded standard
59    # tracking transitions to zero at eps=a; oriented tracking is always zero.
60    sweep = []
61    T = 1.0; L = 0.5
62    for eps in np.linspace(0.0, 0.12, 13):
63        sweep.append({"epsilon": float(eps),
64                      "fixed_h_error": fixed_time_error(T, a, x, v0, L),
65                      "standard_error": standard_error(T, eps, a, x, v0, L),
66                      "oriented_error": oriented_error(T, a, x, v0, L)})
67
68    # Numerical checks of the definitions: h is monotone and every sampled
69    # secant slope lies in [1-eps,1+eps]. Also verify pseudo-path defect.
70    eps = 0.05; T = 1.0
71    t = np.linspace(0, T, 1001)
72    h = (1 + eps) * t
73    slopes = np.diff(h) / np.diff(t)
74    ztilde = flow((1+a)*t, x, v0, L)
75    # Since z_tilde(t)=phi_{(1+a)t}(x), its exact defect is a*f(z_tilde).
76    defect = a * sup_speed(T, x, v0, L)
77    checks = {"h_min_slope": float(slopes.min()), "h_max_slope": float(slopes.max()),
78              "secant_constraint_holds": bool(np.all((slopes >= 1-eps-1e-10) & (slopes <= 1+eps+1e-10))),
79              "pseudo_defect_numeric": float(defect),
80              "pseudo_defect_bound_delta": float(a * sup_speed(T, x, v0, L))}
81    out = {"setup": {"x": x, "v0": v0, "timing_distortion_a": a, "rho": rho},
82           "threshold_sweep": rows, "tracking_sweep": sweep, "math_checks": checks}
83    Path("toy_results.json").write_text(json.dumps(out, indent=2))
84    print(json.dumps(out, indent=2))
85
86if __name__ == "__main__":
87    run()