import json
from pathlib import Path
import numpy as np

# Toy flow: dz/dt = v0 + L z, z(0)=x.  The pseudo-trajectory is an
# exact trajectory with elapsed time (1+a)t, hence it has a pure timing defect.
def flow(t, x, v0, L):
    t = np.asarray(t, dtype=float)
    if abs(L) < 1e-12:
        return x + v0 * t
    return (x + v0 / L) * np.exp(L * t) - v0 / L

def vector_field(z, v0, L):
    return v0 + L * z

def sup_speed(T, x, v0, L, n=20001):
    t = np.linspace(0.0, T, n)
    return float(np.max(np.abs(vector_field(flow(t, x, v0, L), v0, L))))

def standard_error(T, eps, a, x, v0, L, n=2001):
    # For a>=0 and an increasing affine flow, the closest admissible h is
    # h(t)=(1+min(eps,a))t. Its secant slopes are exactly in Rep(eps).
    t = np.linspace(0.0, T, n)
    h_slope = 1.0 + min(float(eps), float(a))
    return float(np.max(np.abs(flow((1.0 + a) * t, x, v0, L)
                         - flow(h_slope * t, x, v0, L))))

def fixed_time_error(T, a, x, v0, L):
    return standard_error(T, 0.0, a, x, v0, L)

def oriented_error(T, a, x, v0, L):
    # h(t)=(1+a)t is an unrestricted increasing homeomorphism.
    return 0.0

def empirical_threshold(T, rho, a, x, v0, L):
    lo, hi = 0.0, max(a, 1e-12)
    for _ in range(70):
        mid = (lo + hi) / 2
        if standard_error(T, mid, a, x, v0, L) <= rho:
            hi = mid
        else:
            lo = mid
    return hi

def run():
    x, v0, a, rho = 0.0, 1.0, 0.08, 0.02
    rows = []
    for L in (0.2, 0.5, 1.0):
        for T in (0.5, 1.0, 2.0):
            speed = sup_speed(T, x, v0, L)
            pred = rho / (T * L * speed)
            obs = empirical_threshold(T, rho, a, x, v0, L)
            rows.append({"L": L, "T": T, "sup_speed": speed,
                         "predicted_epsilon_crit": pred,
                         "observed_epsilon_crit": obs,
                         "ratio_observed_to_predicted": obs / pred})

    # Direct verification of the mechanism over epsilon: bounded standard
    # tracking transitions to zero at eps=a; oriented tracking is always zero.
    sweep = []
    T = 1.0; L = 0.5
    for eps in np.linspace(0.0, 0.12, 13):
        sweep.append({"epsilon": float(eps),
                      "fixed_h_error": fixed_time_error(T, a, x, v0, L),
                      "standard_error": standard_error(T, eps, a, x, v0, L),
                      "oriented_error": oriented_error(T, a, x, v0, L)})

    # Numerical checks of the definitions: h is monotone and every sampled
    # secant slope lies in [1-eps,1+eps]. Also verify pseudo-path defect.
    eps = 0.05; T = 1.0
    t = np.linspace(0, T, 1001)
    h = (1 + eps) * t
    slopes = np.diff(h) / np.diff(t)
    ztilde = flow((1+a)*t, x, v0, L)
    # Since z_tilde(t)=phi_{(1+a)t}(x), its exact defect is a*f(z_tilde).
    defect = a * sup_speed(T, x, v0, L)
    checks = {"h_min_slope": float(slopes.min()), "h_max_slope": float(slopes.max()),
              "secant_constraint_holds": bool(np.all((slopes >= 1-eps-1e-10) & (slopes <= 1+eps+1e-10))),
              "pseudo_defect_numeric": float(defect),
              "pseudo_defect_bound_delta": float(a * sup_speed(T, x, v0, L))}
    out = {"setup": {"x": x, "v0": v0, "timing_distortion_a": a, "rho": rho},
           "threshold_sweep": rows, "tracking_sweep": sweep, "math_checks": checks}
    Path("toy_results.json").write_text(json.dumps(out, indent=2))
    print(json.dumps(out, indent=2))

if __name__ == "__main__":
    run()