import json
import numpy as np


def transverse_matrix(a=1.0, lambda_p=2.0, lambda_i=2.0, k_p=0.0, k_i=0.0):
    """Scalar n=1 instance of the paper's transverse matrix."""
    return np.array([[-a - k_p * lambda_p, -1.0],
                     [k_i * lambda_i, 0.0]], dtype=float)


def euler_radius(M, h):
    return float(np.max(np.abs(np.linalg.eigvals(np.eye(M.shape[0]) + h * M))))


def simulate(a=1.0, bias_difference=1.0, k_p=0.0, k_i=0.0,
             h=0.01, steps=5000):
    """Explicit Euler difference dynamics for two copies and L transverse mode 2."""
    d, z = 0.0, 0.0
    ds = np.empty(steps + 1)
    ds[0] = d
    for t in range(steps):
        dd = -a * d + bias_difference - 2.0 * k_p * d - z
        dz = 2.0 * k_i * d
        d, z = d + h * dd, z + h * dz
        ds[t + 1] = d
    return ds


def eig_pairs(M):
    return [[float(x.real), float(x.imag)] for x in np.linalg.eigvals(M)]


def main():
    a, h = 1.0, 0.01
    cases = {
        "uncoupled": (0.0, 0.0),
        "proportional": (1.0, 0.0),
        "PI_stable": (1.0, 0.2),
        "PI_unstable": (1.0, 200.0),
    }
    results = {}
    for name, (kp, ki) in cases.items():
        d = simulate(k_p=kp, k_i=ki, h=h)
        M = transverse_matrix(k_p=kp, k_i=ki)
        results[name] = {
            "kP": kp, "kI": ki,
            "spectral_radius": euler_radius(M, h),
            "max_continuous_real_part": float(np.max(np.linalg.eigvals(M).real)),
            "continuous_eigenvalues": eig_pairs(M),
            "final_disagreement": float(d[-1]),
            "mean_abs_last_500": float(np.mean(np.abs(d[-500:]))),
            "max_abs": float(np.max(np.abs(d))),
        }

    scan = []
    for ki in np.linspace(0.0, 200.0, 201):
        M = transverse_matrix(k_p=1.0, k_i=float(ki))
        rho = euler_radius(M, h)
        d = simulate(k_p=1.0, k_i=float(ki), h=h, steps=3000)
        scan.append({"kI": float(ki), "rho": rho,
                     "strictly_stable_by_rho": bool(rho < 1.0),
                     "tail_abs": float(np.mean(np.abs(d[-300:])))})
    positive = [x for x in scan if x["kI"] > 0]
    stable = [x for x in positive if x["strictly_stable_by_rho"]]
    unstable = [x for x in positive if not x["strictly_stable_by_rho"]]
    results["scan"] = {
        "largest_positive_stable_kI": max(x["kI"] for x in stable),
        "first_positive_unstable_kI": min(x["kI"] for x in unstable),
        "points": scan,
    }

    # Reproducible checks of both parts of the claimed signature.
    assert results["proportional"]["mean_abs_last_500"] > 0.3
    assert results["PI_stable"]["mean_abs_last_500"] < 0.001
    assert results["PI_stable"]["spectral_radius"] < 1.0
    assert results["PI_unstable"]["spectral_radius"] > 1.0
    assert results["PI_unstable"]["max_abs"] > 1e6

    with open("pi_results.json", "w") as f:
        json.dump(results, f, indent=2)
    summary = {k: v for k, v in results.items() if k != "scan"}
    print(json.dumps(summary, indent=2))
    print("positive-kI Euler boundary:",
          results["scan"]["largest_positive_stable_kI"],
          "to", results["scan"]["first_positive_unstable_kI"])


if __name__ == "__main__":
    main()