#!/usr/bin/env python3
"""Recalculate a derivative-mixing Hamiltonian; a=1, spatial gradients omitted.

Run with Python, SymPy and NumPy. An optional --output writes the evidence JSON.
The numerical implementation is handwritten separately from the symbolic one.
"""
import argparse
import json
from pathlib import Path

import numpy as np
import sympy as sp


def verify():
    u, v, sigma, rho, mass, pp, ps = sp.symbols('u v sigma rho m p_pi p_sigma', real=True)
    lagrangian = (u**2 + v**2 - mass**2 * sigma**2) / 2 + rho * u * sigma
    momenta = [sp.diff(lagrangian, w) for w in (u, v)]
    velocities = sp.solve([pp - momenta[0], ps - momenta[1]], [u, v])
    legendre = (pp * u + ps * v - lagrangian).subs(velocities)
    hamiltonian = (pp**2 + ps**2) / 2 - rho * sigma * pp + (mass**2 + rho**2) * sigma**2 / 2
    residual = sp.simplify(legendre - hamiltonian)
    assert residual == 0
    decoupled = sp.simplify(hamiltonian.subs(rho, 0) - (pp**2 + ps**2 + mass**2 * sigma**2) / 2)
    assert decoupled == 0
    incorrect_h = hamiltonian - rho**2 * sigma**2 / 2
    incorrect_residual = sp.simplify(legendre - incorrect_h)
    assert incorrect_residual == rho**2 * sigma**2 / 2

    # This path uses the original L and explicit momenta, without lambdify or
    # generated SymPy expressions. Numeric checks include rho=0 and both signs.
    rng = np.random.default_rng(20261007)
    energy_error, velocity_error, detected_negative_cases = 0.0, 0.0, 0
    samples = []
    for n in range(64):
        mixing = 0.0 if n == 0 else float(rng.uniform(-3, 3))
        m = float(rng.uniform(0.2, 3))
        s, du, dv = (float(x) for x in rng.uniform(-2, 2, size=3))
        p, q = du + mixing * s, dv
        original_l = 0.5 * (du * du + dv * dv - m * m * s * s) + mixing * du * s
        direct_energy = p * du + q * dv - original_l

        def numeric_h(p_arg, q_arg):
            return 0.5 * (p_arg * p_arg + q_arg * q_arg) - mixing * s * p_arg + 0.5 * (m * m + mixing * mixing) * s * s

        candidate_energy = numeric_h(p, q)
        relative_energy = abs(candidate_energy - direct_energy) / max(1, abs(direct_energy))
        eps = 1e-5
        numeric_velocities = [(numeric_h(p + eps, q) - numeric_h(p - eps, q)) / (2 * eps),
                              (numeric_h(p, q + eps) - numeric_h(p, q - eps)) / (2 * eps)]
        relative_velocity = max(abs(x - y) / max(1, abs(y)) for x, y in zip(numeric_velocities, [du, dv]))
        assert relative_energy < 1e-12, (n, relative_energy)
        assert relative_velocity < 1e-8, (n, relative_velocity)
        energy_error = max(energy_error, relative_energy)
        velocity_error = max(velocity_error, relative_velocity)
        wrong_energy = candidate_energy - 0.5 * mixing * mixing * s * s
        if abs(wrong_energy - direct_energy) > 1e-10:
            detected_negative_cases += 1
        samples.append({'rho': mixing, 'mass': m, 'sigma': s, 'pi_velocity': du, 'sigma_velocity': dv,
                        'energy_relative_error': relative_energy, 'velocity_relative_error': relative_velocity})
    assert detected_negative_cases == 63
    return {'scope': 'Homogeneous quadratic model; a=1; no spatial gradients or cubic interactions',
            'symbolic_legendre_residual': str(residual), 'zero_mixing_residual': str(decoupled),
            'negative_control_residual': str(incorrect_residual), 'numerical_samples': 64,
            'maximum_energy_relative_error': energy_error, 'maximum_velocity_relative_error': velocity_error,
            'negative_control_detected_cases': detected_negative_cases, 'status': 'passed', 'inputs': samples}


if __name__ == '__main__':
    parser = argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--output', type=Path)
    args = parser.parse_args()
    result = verify()
    if args.output:
        args.output.write_text(json.dumps(result, indent=2) + '\n')
    print(json.dumps({k: v for k, v in result.items() if k != 'inputs'}, indent=2))