"""Analytic consistency checks for the partial-reduction BTZ phase diagram.

This module uses only Python's standard library. It does not generate research
data or fit parameters; it independently verifies the algebraic identities,
channel inequalities, transition domain, and cutoff expansion in the manuscript.
"""

from __future__ import annotations

import math


def _positive(name: str, value: float) -> None:
    if not math.isfinite(value) or value <= 0.0:
        raise ValueError(f"{name} must be finite and positive")


def fg_radius(z: float, ell: float, r_plus: float) -> float:
    """Return the BTZ radius r(z) in Fefferman-Graham gauge."""
    _positive("z", z)
    _positive("ell", ell)
    _positive("r_plus", r_plus)
    return ell**2 / z + r_plus**2 * z / (4.0 * ell**2)


def fg_cutoff(r_c: float, ell: float, r_plus: float) -> float:
    """Return the boundary cutoff z_c corresponding to r=r_c."""
    _positive("r_c", r_c)
    _positive("ell", ell)
    _positive("r_plus", r_plus)
    if r_c < r_plus:
        raise ValueError("r_c must lie on or outside the outer horizon")
    return 2.0 * ell**2 / (r_c + math.sqrt(r_c**2 - r_plus**2))


def surface_gravity(ell: float, r_plus: float) -> float:
    _positive("ell", ell)
    _positive("r_plus", r_plus)
    return r_plus / ell**2


def reduced_newton_constant(g3: float, ell: float) -> float:
    """Return G2=G3/ell in the partial angular reduction."""
    _positive("g3", g3)
    _positive("ell", ell)
    return g3 / ell


def reduced_dilaton(radius: float, alpha: float, ell: float) -> float:
    """Return Phi=2*pi*alpha*r/ell for a sector of size 2*pi*alpha."""
    _positive("radius", radius)
    _positive("alpha", alpha)
    _positive("ell", ell)
    if alpha > 1.0:
        raise ValueError("alpha must satisfy 0 < alpha <= 1")
    return 2.0 * math.pi * alpha * radius / ell


def dilaton_entropy_one_endpoint(
    radius: float, alpha: float, ell: float, g3: float
) -> float:
    """Return Phi/(4*G2) for one reduced QES endpoint."""
    phi = reduced_dilaton(radius, alpha, ell)
    g2 = reduced_newton_constant(g3, ell)
    return phi / (4.0 * g2)


def parent_arc_entropy(radius: float, alpha: float, g3: float) -> float:
    """Return the parent arc length 2*pi*alpha*r divided by 4*G3."""
    _positive("radius", radius)
    _positive("alpha", alpha)
    _positive("g3", g3)
    if alpha > 1.0:
        raise ValueError("alpha must satisfy 0 < alpha <= 1")
    return 2.0 * math.pi * alpha * radius / (4.0 * g3)


def angular_image(delta_phi: float, winding: int) -> float:
    """Return the global-BTZ angular image Delta phi + 2*pi*k."""
    if not isinstance(winding, int):
        raise TypeError("winding must be an integer")
    return delta_phi + 2.0 * math.pi * winding


def symmetric_interface_images(alpha: float) -> tuple[float, float]:
    """Return the two short images joining endpoints at phi=+/-pi*alpha.

    For the ordered endpoints (p_minus, p_plus), the k=0 image crosses the
    reduced sector and the k=+1 image crosses the
    complementary bath sector.
    """
    _positive("alpha", alpha)
    if alpha > 1.0:
        raise ValueError("alpha must satisfy 0 < alpha <= 1")
    return -2.0 * math.pi * alpha, 2.0 * math.pi * (1.0 - alpha)


def transition_tau_finite_cutoff(
    delta_phi: float, r_c: float, ell: float, r_plus: float
) -> float | None:
    """Return tau=kappa*t for the connected/disconnected crossing.

    Returns None when the disconnected branch is already shorter at t=0,
    so no nonnegative crossing exists.  This implements Eq. (exacttau) of
    the manuscript without fitted parameters.
    """
    for name, value in (("r_c", r_c), ("ell", ell), ("r_plus", r_plus)):
        _positive(name, value)
    if r_c <= r_plus:
        raise ValueError("r_c must lie strictly outside the outer horizon")
    a_ratio = (r_c / r_plus) ** 2
    x = r_plus * abs(delta_phi) / (2.0 * ell)
    xi = -1.0 + (2.0 * a_ratio / (a_ratio - 1.0)) * math.sinh(x) ** 2
    tolerance = 32.0 * math.ulp(1.0)
    if xi < 1.0 - tolerance:
        return None
    return 0.5 * math.acosh(max(1.0, xi))


def transition_tau_asymptotic(
    delta_phi: float, ell: float, r_plus: float
) -> float | None:
    """Return the r_c -> infinity limit of the dimensionless crossing time."""
    _positive("ell", ell)
    _positive("r_plus", r_plus)
    sinh_x = math.sinh(r_plus * abs(delta_phi) / (2.0 * ell))
    tolerance = 32.0 * math.ulp(1.0)
    if sinh_x < 1.0 - tolerance:
        return None
    return math.acosh(max(1.0, sinh_x))


def saturation_fraction(alpha: float) -> float:
    """Return the shorter complementary angular fraction min(alpha, 1-alpha)."""
    _positive("alpha", alpha)
    if alpha >= 1.0:
        raise ValueError("alpha must satisfy 0 < alpha < 1")
    return min(alpha, 1.0 - alpha)


def saturation_angle(alpha: float) -> float:
    """Return the angular size of the physical disconnected network."""
    return 2.0 * math.pi * saturation_fraction(alpha)


def physical_transition_tau(
    alpha: float, r_c: float, ell: float, r_plus: float
) -> float | None:
    """Return the exact finite-cutoff crossing of the physical lower envelope."""
    return transition_tau_finite_cutoff(
        saturation_angle(alpha), r_c, ell, r_plus
    )


def transition_alpha_critical(r_c: float, ell: float, r_plus: float) -> float:
    """Return the finite-cutoff lower boundary alpha_c of the transition region."""
    for name, value in (("r_c", r_c), ("ell", ell), ("r_plus", r_plus)):
        _positive(name, value)
    if r_c <= r_plus:
        raise ValueError("r_c must lie strictly outside the outer horizon")
    return (ell / (math.pi * r_plus)) * math.asinh(
        math.sqrt(1.0 - (r_plus / r_c) ** 2)
    )


def transition_invariants(
    alpha: float, tau: float, r_c: float, ell: float, r_plus: float
) -> tuple[float, float]:
    """Return (C_connected, C_saturation) for the two surviving networks."""
    if not math.isfinite(tau) or tau < 0.0:
        raise ValueError("tau must be finite and nonnegative")
    for name, value in (("r_c", r_c), ("ell", ell), ("r_plus", r_plus)):
        _positive(name, value)
    if r_c <= r_plus:
        raise ValueError("r_c must lie strictly outside the outer horizon")
    a_ratio = (r_c / r_plus) ** 2
    x = math.pi * r_plus * saturation_fraction(alpha) / ell
    c_connected = a_ratio + (a_ratio - 1.0) * math.cosh(2.0 * tau)
    c_saturation = 1.0 + 2.0 * a_ratio * math.sinh(x) ** 2
    return c_connected, c_saturation


def cutoff_shift_coefficient(alpha: float, ell: float, r_plus: float) -> float:
    """Return the O((r_plus/r_c)^2) coefficient of the crossing time."""
    _positive("ell", ell)
    _positive("r_plus", r_plus)
    s_value = math.sinh(
        math.pi * r_plus * saturation_fraction(alpha) / ell
    )
    if s_value <= 1.0:
        raise ValueError("the continuum transition must lie inside its open domain")
    return s_value / (2.0 * math.sqrt(s_value**2 - 1.0))


def exterior_horizon_coordinate(radius: float, r_plus: float) -> float:
    """Return q=(r-r_plus)/(r+r_plus) for the exterior domain."""
    _positive("radius", radius)
    _positive("r_plus", r_plus)
    if radius < r_plus:
        raise ValueError("the exterior coordinate requires radius >= r_plus")
    return (radius - r_plus) / (radius + r_plus)


def radius_from_horizon_coordinate(q: float, r_plus: float) -> float:
    """Invert the exterior compact horizon coordinate, 0 <= q < 1."""
    _positive("r_plus", r_plus)
    if not math.isfinite(q) or q < 0.0 or q >= 1.0:
        raise ValueError("q must satisfy 0 <= q < 1")
    return r_plus * (1.0 + q) / (1.0 - q)


def proper_distance_from_horizon_coordinate(q: float) -> float:
    """Return chi=arcosh(r/r_plus)=2 artanh(sqrt(q))."""
    if not math.isfinite(q) or q < 0.0 or q >= 1.0:
        raise ValueError("q must satisfy 0 <= q < 1")
    return 2.0 * math.atanh(math.sqrt(q))


def same_exterior_cosh_length(
    r1: float,
    r2: float,
    delta_phi: float,
    delta_t: float,
    ell: float,
    r_plus: float,
) -> float:
    """Return cosh(L/ell) for a spacelike same-exterior BTZ geodesic."""
    for name, value in (("r1", r1), ("r2", r2), ("ell", ell), ("r_plus", r_plus)):
        _positive(name, value)
    if r1 < r_plus or r2 < r_plus:
        raise ValueError("endpoints must be on or outside the outer horizon")
    angular = (r1 * r2 / r_plus**2) * math.cosh(r_plus * delta_phi / ell)
    radial = math.sqrt((r1**2 - r_plus**2) * (r2**2 - r_plus**2)) / r_plus**2
    temporal = radial * math.cosh(surface_gravity(ell, r_plus) * delta_t)
    return angular - temporal


def opposite_exterior_cosh_length(
    r1: float,
    r2: float,
    delta_phi: float,
    t_left: float,
    t_right: float,
    ell: float,
    r_plus: float,
) -> float:
    """Return cosh(L/ell) for an opposite-exterior BTZ geodesic."""
    for name, value in (("r1", r1), ("r2", r2), ("ell", ell), ("r_plus", r_plus)):
        _positive(name, value)
    if r1 < r_plus or r2 < r_plus:
        raise ValueError("endpoints must be on or outside the outer horizon")
    angular = (r1 * r2 / r_plus**2) * math.cosh(r_plus * delta_phi / ell)
    radial = math.sqrt((r1**2 - r_plus**2) * (r2**2 - r_plus**2)) / r_plus**2
    temporal = radial * math.cosh(surface_gravity(ell, r_plus) * (t_left + t_right))
    return angular + temporal


def geodesic_length_from_cosh(cosh_length: float, ell: float) -> float:
    _positive("ell", ell)
    if cosh_length < 1.0:
        raise ValueError("the selected endpoints are not spacelike in this branch")
    return ell * math.acosh(cosh_length)


def equal_time_cosh_length(r_c: float, delta_phi: float, ell: float, r_plus: float) -> float:
    """Closed form of the equal-time, equal-radius same-exterior result."""
    return 1.0 + 2.0 * (r_c / r_plus) ** 2 * math.sinh(
        r_plus * delta_phi / (2.0 * ell)
    ) ** 2


def run_self_checks() -> None:
    ell = 2.3
    r_plus = 1.1
    r_c = 50.0
    delta_phi = 0.8
    alpha = 0.37
    g3 = 0.09

    # Exact local equality between the reduced dilaton entropy and
    # the Bekenstein-Hawking entropy of the parent angular arc.
    s_dilaton = dilaton_entropy_one_endpoint(r_c, alpha, ell, g3)
    s_arc = parent_arc_entropy(r_c, alpha, g3)
    assert math.isclose(s_dilaton, s_arc, rel_tol=2e-14, abs_tol=2e-14)

    z_c = fg_cutoff(r_c, ell, r_plus)
    reconstructed = fg_radius(z_c, ell, r_plus)
    assert math.isclose(reconstructed, r_c, rel_tol=2e-14, abs_tol=2e-14)

    general = same_exterior_cosh_length(
        r_c, r_c, delta_phi, 0.0, ell, r_plus
    )
    reduced = equal_time_cosh_length(r_c, delta_phi, ell, r_plus)
    assert math.isclose(general, reduced, rel_tol=2e-14, abs_tol=2e-14)

    # k and -k describe images whose lengths agree after reversing the
    # interval orientation.
    image_plus = same_exterior_cosh_length(
        r_c, r_c, angular_image(delta_phi, 1), 0.0, ell, r_plus
    )
    image_reversed = same_exterior_cosh_length(
        r_c, r_c, angular_image(-delta_phi, -1), 0.0, ell, r_plus
    )
    assert math.isclose(image_plus, image_reversed, rel_tol=2e-14, abs_tol=2e-14)

    delta_reduced, delta_bath = symmetric_interface_images(alpha)
    assert math.isclose(delta_reduced, -2.0 * math.pi * alpha)
    assert math.isclose(delta_bath, 2.0 * math.pi * (1.0 - alpha))

    t = 0.7
    opposite = opposite_exterior_cosh_length(
        r_c, r_c, 0.0, t, t, ell, r_plus
    )
    expected = (r_c / r_plus) ** 2 + (
        (r_c**2 - r_plus**2) / r_plus**2
    ) * math.cosh(2.0 * surface_gravity(ell, r_plus) * t)
    assert math.isclose(opposite, expected, rel_tol=2e-14, abs_tol=2e-14)

    q = exterior_horizon_coordinate(r_c, r_plus)
    assert math.isclose(
        radius_from_horizon_coordinate(q, r_plus),
        r_c,
        rel_tol=2e-14,
        abs_tol=2e-14,
    )
    chi_direct = math.acosh(r_c / r_plus)
    chi_from_q = proper_distance_from_horizon_coordinate(q)
    assert math.isclose(chi_direct, chi_from_q, rel_tol=2e-14, abs_tol=2e-14)

    exact_length = geodesic_length_from_cosh(reduced, ell)
    asymptotic_length = 2.0 * ell * math.log(
        (2.0 * r_c / r_plus)
        * math.sinh(r_plus * delta_phi / (2.0 * ell))
    )
    relative_difference = abs(exact_length - asymptotic_length) / exact_length
    assert relative_difference < 5e-3

    # Verify the exact finite-cutoff crossing formula against independent
    # evaluations of the two BTZ geodesic expressions.
    ell_transition = 1.0
    r_plus_transition = 2.0
    r_c_transition = 50.0
    delta_transition = math.pi
    tau_transition = transition_tau_finite_cutoff(
        delta_transition, r_c_transition, ell_transition, r_plus_transition
    )
    assert tau_transition is not None
    kappa_transition = surface_gravity(ell_transition, r_plus_transition)
    time_transition = tau_transition / kappa_transition
    cosh_connected = opposite_exterior_cosh_length(
        r_c_transition,
        r_c_transition,
        0.0,
        time_transition,
        time_transition,
        ell_transition,
        r_plus_transition,
    )
    cosh_disconnected = same_exterior_cosh_length(
        r_c_transition,
        r_c_transition,
        delta_transition,
        0.0,
        ell_transition,
        r_plus_transition,
    )
    assert math.isclose(
        cosh_connected, cosh_disconnected, rel_tol=3e-14, abs_tol=3e-14
    )

    # A sufficiently small angular image has no connected-to-disconnected
    # crossing because the disconnected branch wins already at t=0.
    assert (
        transition_tau_finite_cutoff(
            0.1, r_c_transition, ell_transition, r_plus_transition
        )
        is None
    )

    tau_large_cutoff = transition_tau_finite_cutoff(
        delta_transition, 1.0e7, ell_transition, r_plus_transition
    )
    tau_continuum = transition_tau_asymptotic(
        delta_transition, ell_transition, r_plus_transition
    )
    assert tau_large_cutoff is not None and tau_continuum is not None
    assert math.isclose(tau_large_cutoff, tau_continuum, rel_tol=2e-14)

    # Test the complete lower-envelope formulas on a deterministic grid.
    grid_count = 0
    for ratio in (0.6, 0.8, 1.0, 1.5, 2.5):
        ell_grid = 1.0
        r_plus_grid = ratio * ell_grid
        for cutoff_ratio in (2.0, 5.0, 20.0, 100.0):
            r_c_grid = cutoff_ratio * r_plus_grid
            alpha_critical = transition_alpha_critical(
                r_c_grid, ell_grid, r_plus_grid
            )
            for alpha_grid in (0.1, 0.25, 0.4, 0.5, 0.6, 0.75, 0.9):
                tau_grid = physical_transition_tau(
                    alpha_grid, r_c_grid, ell_grid, r_plus_grid
                )
                expected_exists = (
                    alpha_critical <= 0.5
                    and alpha_critical <= alpha_grid <= 1.0 - alpha_critical
                )
                assert (tau_grid is not None) == expected_exists
                if tau_grid is not None:
                    c_connected, c_saturation = transition_invariants(
                        alpha_grid,
                        tau_grid,
                        r_c_grid,
                        ell_grid,
                        r_plus_grid,
                    )
                    assert math.isclose(
                        c_connected,
                        c_saturation,
                        rel_tol=8e-13,
                        abs_tol=8e-13,
                    )
                    grid_count += 1

                # Purity symmetry under exchange of reduced sector and bath.
                assert math.isclose(
                    saturation_angle(alpha_grid),
                    saturation_angle(1.0 - alpha_grid),
                    rel_tol=2e-14,
                    abs_tol=2e-14,
                )

    # The crossed left-right pairing has a strictly larger invariant whenever
    # its shortest angular image is nonzero.
    for alpha_grid in (0.1, 0.25, 0.4, 0.6, 0.75, 0.9):
        shortest_crossed = saturation_angle(alpha_grid)
        parallel = opposite_exterior_cosh_length(
            r_c, r_c, 0.0, t, t, ell, r_plus
        )
        crossed = opposite_exterior_cosh_length(
            r_c, r_c, shortest_crossed, t, t, ell, r_plus
        )
        assert crossed > parallel

        # Every nonprimitive same-boundary image has greater angular magnitude
        # than the winning primitive image and therefore greater length.
        base_delta = -2.0 * math.pi * alpha_grid
        primitive_magnitude = saturation_angle(alpha_grid)
        primitive_invariant = equal_time_cosh_length(
            r_c, primitive_magnitude, ell, r_plus
        )
        for winding in range(-5, 6):
            image_magnitude = abs(angular_image(base_delta, winding))
            if image_magnitude > primitive_magnitude + 1.0e-12:
                assert (
                    equal_time_cosh_length(
                        r_c, image_magnitude, ell, r_plus
                    )
                    > primitive_invariant
                )

        # A mixed primitive orientation on the two boundaries carries the
        # nontrivial cylinder cycle and must be completed by a horizon section.
        # Its total length is strictly above the shorter equal-orientation pair.
        length_reduced = geodesic_length_from_cosh(
            equal_time_cosh_length(
                r_c, 2.0 * math.pi * alpha_grid, ell, r_plus
            ),
            ell,
        )
        length_bath = geodesic_length_from_cosh(
            equal_time_cosh_length(
                r_c, 2.0 * math.pi * (1.0 - alpha_grid), ell, r_plus
            ),
            ell,
        )
        horizon_length = 2.0 * math.pi * r_plus
        mixed_with_horizon = length_reduced + length_bath + horizon_length
        shortest_equal_orientation = 2.0 * min(length_reduced, length_bath)
        assert mixed_with_horizon > shortest_equal_orientation

    # Verify the leading finite-cutoff shift away from the critical line.
    alpha_shift = 0.4
    ell_shift = 1.0
    r_plus_shift = 1.0
    tau_infinity_shift = transition_tau_asymptotic(
        saturation_angle(alpha_shift), ell_shift, r_plus_shift
    )
    assert tau_infinity_shift is not None
    coefficient = cutoff_shift_coefficient(
        alpha_shift, ell_shift, r_plus_shift
    )
    r_c_shift = 1.0e4 * r_plus_shift
    tau_finite_shift = physical_transition_tau(
        alpha_shift, r_c_shift, ell_shift, r_plus_shift
    )
    assert tau_finite_shift is not None
    u_shift = (r_plus_shift / r_c_shift) ** 2
    numerical_coefficient = (tau_finite_shift - tau_infinity_shift) / u_shift
    assert math.isclose(numerical_coefficient, coefficient, rel_tol=2e-7)

    print("All BTZ analytic consistency checks passed.")
    print(f"Partial-reduction entropy relative error: {abs(s_dilaton-s_arc)/s_arc:.3e}")
    print(f"FG inverse-map relative error: {abs(reconstructed-r_c)/r_c:.3e}")
    print(f"Finite-cutoff/asymptotic length difference: {relative_difference:.3e}")
    print(f"Exact finite-cutoff crossing tau: {tau_transition:.12f}")
    print(f"Branch equality relative error: {abs(cosh_connected-cosh_disconnected)/cosh_disconnected:.3e}")
    print(f"Admissible lower-envelope grid crossings checked: {grid_count}")
    print("Mixed-orientation homology bound: passed")
    print(f"Cutoff-shift coefficient relative error: {abs(numerical_coefficient-coefficient)/coefficient:.3e}")


if __name__ == "__main__":
    run_self_checks()