"""Unified UTFANSWF v23 Row-19 -> Row-23 gravity/scalar/GW research block.

Integrated v23 reproducibility-side closure. Protected main.tex/main.pdf remain unchanged.

Purpose:
  1) consume the audited Row-19 EH/Litim fixed point;
  2) solve the explicit printed scalar beta functions on that gravity background;
  3) correct the R^2 metric-variation normalization directly from f'(R);
  4) derive the canonical tensor propagation map for a stationary constant-scalar
     black-hole background;
  5) hand that map to the existing Row-23 Kerr-QNM executable.

The result is deliberately narrower than a full general curvature-squared/nonminimal
FRG truncation. It closes the canonical black-hole tensor branch needed by Row 23.
"""
from __future__ import annotations

import json
import math
from pathlib import Path
from typing import Any

import numpy as np
from scipy.optimize import root

from tools.row19.frg_eh_litim import benchmark_result as eh_benchmark
from harness.gate_row23_ringdown import gate_ringdown_echoes
from harness.gate_row22_nz import gate_nz

PI = math.pi
ROOT = Path(__file__).resolve().parents[2]


def scalar_beta(y: np.ndarray, g: float, lam: float) -> np.ndarray:
    """Printed v23 scalar beta functions evaluated on audited EH fixed point.

    Ordering y = [lambda_phi_tilde, m2_tilde].
    """
    lphi, m2 = y[0], y[1]
    d = 1.0 - 2.0 * lam
    if abs(1.0 + m2) < 1e-12 or abs(d) < 1e-12:
        return np.asarray([1e30, 1e30], dtype=np.result_type(y, float))

    eta_phi = (
        lphi / (6.0 * PI**2 * (1.0 + m2) ** 2)
        - g / (6.0 * PI * d**2)
    )
    beta_m2 = (
        (-2.0 + eta_phi) * m2
        + lphi / (6.0 * PI**2 * (1.0 + m2))
        - g / (PI * d)
    )
    beta_lphi = (
        2.0 * eta_phi * lphi
        + 3.0 * lphi**2 / (2.0 * PI**2 * (1.0 + m2) ** 2)
        - g * lphi / (2.0 * PI * d**2)
    )
    return np.asarray([beta_lphi, beta_m2])


def scalar_jacobian_complex_step(x: np.ndarray, g: float, lam: float, eps: float = 1e-30) -> np.ndarray:
    x = np.asarray(x, dtype=float)
    J = np.empty((2, 2), dtype=float)
    for j in range(2):
        xc = x.astype(complex)
        xc[j] += 1j * eps
        J[:, j] = np.imag(scalar_beta(xc, g, lam)) / eps
    return J


def find_scalar_fixed_points(g: float, lam: float) -> list[np.ndarray]:
    sols: list[np.ndarray] = []
    seeds_l = np.linspace(-20.0, 20.0, 25)
    seeds_m = np.linspace(-0.95, 4.0, 25)
    for lp0 in seeds_l:
        for m0 in seeds_m:
            rr = root(lambda z: np.real(scalar_beta(z, g, lam)).astype(float), [lp0, m0], method="hybr")
            if not rr.success:
                continue
            x = np.asarray(rr.x, dtype=float)
            if not np.all(np.isfinite(x)):
                continue
            if abs(1.0 + x[1]) < 1e-8:
                continue
            if np.linalg.norm(np.real(scalar_beta(x, g, lam))) > 1e-9:
                continue
            if not any(np.linalg.norm(x - z) < 1e-7 for z in sols):
                sols.append(x)
    return sorted(sols, key=lambda z: (float(z[0]), float(z[1])))


def classify_scalar_fp(x: np.ndarray, g: float, lam: float) -> dict[str, Any]:
    lphi, m2 = map(float, x)
    b = np.real(scalar_beta(x, g, lam)).astype(float)
    J = scalar_jacobian_complex_step(x, g, lam)
    eig = np.linalg.eigvals(J)
    bfb = bool(lphi > 0.0 or (abs(lphi) < 1e-10 and m2 >= 0.0))
    out: dict[str, Any] = {
        "lambda_phi_tilde": lphi,
        "m2_tilde": m2,
        "beta_norm": float(np.linalg.norm(b)),
        "bounded_below_quartic_truncation": bfb,
        "stability_eigenvalues": [
            {"re": float(e.real), "im": float(e.imag)} for e in eig
        ],
    }
    if lphi > 0.0 and m2 < 0.0:
        v2 = -6.0 * m2 / lphi
        out.update({
            "stationary_vacuum_type": "Z2_broken_pair",
            "phi0_over_k_abs": math.sqrt(v2),
            "fluctuation_m2_over_k2": -2.0 * m2,
            "vacuum_energy_over_k4_without_constant_shift": -1.5 * m2 * m2 / lphi,
            "note": "Vacuum-energy constant renormalizes the cosmological constant; UV fixed-point values are not treated as the astrophysical IR vacuum values.",
        })
    elif abs(lphi) < 1e-10 and m2 < 0.0:
        out["stationary_vacuum_type"] = "UNBOUNDED_AT_QUARTIC_ORDER"
    else:
        out["stationary_vacuum_type"] = "OTHER"
    return out


def r2_variation_audit() -> dict[str, Any]:
    """Audit coefficients for L_R2 = alpha_R/(12 M^2) R^2.

    If B = alpha_R/(12 M^2), variation of B R^2 is
      B [2 R Rmn - 1/2 R^2 gmn - 2 nabla_m nabla_n R + 2 gmn Box R].

    The legacy appendix prints the same bracket with alpha_R/(6 M^2), i.e. 2B.
    """
    correct_prefactor_units = 1.0 / 12.0
    legacy_prefactor_units = 1.0 / 6.0
    ratio = legacy_prefactor_units / correct_prefactor_units
    return {
        "lagrangian_coefficient": "B = alpha_R/(12 M^2)",
        "correct_metric_variation": "H_R2 = alpha_R/(12 M^2) [2 R R_mn - (1/2)R^2 g_mn - 2 nabla_m nabla_n R + 2 g_mn Box R]",
        "equivalent_factored_form": "H_R2 = alpha_R/(6 M^2) [R R_mn - (1/4)R^2 g_mn - nabla_m nabla_n R + g_mn Box R]",
        "legacy_appendix_prefactor_relative_to_correct": ratio,
        "legacy_factor_of_two_correction_required": bool(abs(ratio - 2.0) < 1e-15),
        "scalaron_mass_squared_convention": "m_s^2 = M_Pl^2 M^2 / alpha_R (vacuum trace equation, using canonical normalization)",
    }



def global_response_operator_map() -> dict[str, Any]:
    """Replace legacy scalar channel equality by operator universality.

    The universal object is a positive/passive response operator on the direct
    sum of physical channels. Observable recycling efficiencies are channel
    projections. Equality of two projected scalar outputs requires an
    independently derived intertwiner/symmetry; it is not automatic.
    """
    # Algebra-only sanity example: a PSD response can vanish in the TT channel
    # while remaining nonzero in another channel. 0.37 is not a physics value.
    X = np.diag([0.0, 0.37])
    P_tt = np.diag([1.0, 0.0])
    P_ax = np.diag([0.0, 1.0])
    evals = np.linalg.eigvalsh(X)
    xi_tt = float(np.trace(P_tt @ X @ P_tt))
    xi_ax = float(np.trace(P_ax @ X @ P_ax))
    return {
        "universal_object": "dimensionless passive HB/NZH response operator X_Sigma >= 0 on direct-sum channel space",
        "channel_law": "Xi_c = <u_c| P_c X_Sigma P_c |u_c>/<u_c|u_c>",
        "entropy_law": "Sdot_c = (A_Sigma,c kappa_eff,c / 2pi) Xi_c",
        "legacy_scalar_equality": "R_a = R_echo = Xi is NOT canonical unless an inter-channel intertwiner is independently derived",
        "cross_channel_equality_condition": "Need U_ab with P_b X P_b = U_ab (P_a X P_a) U_ab^dagger and mapped normalized source states; absent this, projected scalars need not be equal.",
        "canonical_TT_projection": "P_TT X_Sigma P_TT = 0 on stationary Kerr constant-scalar branch",
        "axion_projection": "May be nonzero through the axion-photon response; Row 21 remains an independent channel gate.",
        "checks": {
            "algebra_example_operator_is_positive_semidefinite": bool(float(evals.min()) >= -1e-15),
            "zero_TT_projection_can_coexist_with_nonzero_other_projection": bool(abs(xi_tt) < 1e-15 and xi_ax > 0.0),
            "no_new_channel_projection_fit_parameter_introduced": True,
        },
        "algebra_example_only": {"Xi_TT": xi_tt, "Xi_other": xi_ax, "not_a_physics_calibration": True},
        "disposition": "REPLACE_DIRECT_SCALAR_CHANNEL_EQUALITY_WITH_OPERATOR_PROJECTION_UNIVERSALITY",
    }


def canonical_flrw_tensor_map() -> dict[str, Any]:
    """Canonical nonstationary FLRW tensor map for f(R)=EH+R^2.

    For f(R) gravity there is no extra massive spin-2 pole and c_T=1. The
    effective tensor kinetic normalization is M_*^2=2 f_R. This fixes the GW
    friction law symbolically. Numerical damping remains unfrozen until the
    canonical action combination alpha_R/M^2 is fixed in the IR.
    """
    return {
        "branch": "canonical_EH_plus_R2__FLRW_linear_TT",
        "f_R": "M_Pl^2/2 + alpha_R R/(6 M^2)",
        "Mstar_squared": "M_Pl^2 + alpha_R R(z)/(3 M^2) = M_Pl^2[1+R(z)/(3 m_s^2)]",
        "cT_squared": 1.0,
        "mT_squared": 0.0,
        "nu_general": "nu(z)=d ln Mstar^2/d ln a=-(1+z)R'(z)/[3 m_s^2+R(z)]",
        "gw_luminosity_distance": "dL_GW/dL_EM=sqrt([1+R(0)/(3m_s^2)]/[1+R(z)/(3m_s^2)])",
        "scalaron_ownership": "m_s^2=M_Pl^2 M^2/alpha_R; Mstar and nu are derived from this same R^2-sector mass and are not independent GW knobs",
        "flat_LCDM_R": "R(z)=3 H0^2 [Omega_m(1+z)^3 + 4 Omega_Lambda]",
        "flat_LCDM_dimensionless_x": "x_R = alpha_R H0^2/(M^2 M_Pl^2)=H0^2/m_s^2",
        "flat_LCDM_Mstar_ratio": "Mstar^2/M_Pl^2 = 1 + x_R[Omega_m(1+z)^3+4 Omega_Lambda]",
        "flat_LCDM_nu": "nu(z)=-3 x_R Omega_m(1+z)^3/[1+x_R(Omega_m(1+z)^3+4 Omega_Lambda)]",
        "stability_condition": "Mstar^2(z)>0 over the propagated redshift interval",
        "stationary_or_constant_R_limit": {"nu": 0.0, "dL_GW_over_dL_EM": 1.0},
        "closure": "FUNCTIONAL_FORM_CLOSED__NUMERICAL_x_R_NOT_FROZEN",
        "claim_boundary": "This is the complete canonical linear-TT propagation map for the EH+R^2 branch. Its only unfrozen propagation scale is the already-owned scalaron mass m_s (equivalently x_R=H0^2/m_s^2); R_mn^2/nonminimal phi R remain separate extensions.",
    }

def canonical_tensor_map() -> dict[str, Any]:
    """Canonical R + R^2 tensor map on stationary constant-scalar branch.

    For a Ricci-flat Kerr background R0=0 and TT perturbation delta R=0:
      delta H_R2^TT = 0.
    More generally f(R) carries an extra scalaron but no extra massive spin-2 pole.
    """
    return {
        "branch": "stationary_constant_scalar__canonical_R_plus_R2__no_Ricci_tensor_squared",
        "assumptions": {
            "scalar_background_stationary": True,
            "grad_phi0": 0.0,
            "Vprime_phi0": 0.0,
            "vacuum_energy_absorbed_into_Lambda_eff": True,
            "astrophysical_row23_limit": "Lambda_eff M_BH^2 << 1 so Kerr is the active benchmark",
            "preferred_frame_background": False,
        },
        "tensor_quantities": {
            "alpha1_preferred_frame": 0.0,
            "cT_squared": 1.0,
            "mT_squared": 0.0,
            "Mstar_squared_general_constant_R0": "M_Pl^2 + alpha_R R0/(3 M^2)",
            "Mstar_squared_Ricci_flat": "M_Pl^2",
            "nu_stationary": 0.0,
            "spin2_relative_propagation_delay": 0.0,
            "TT_equation_Ricci_flat": "delta G_mn^TT = 0",
            "extra_mode": "scalaron only: (Box - m_s^2) delta R = 0",
        },
        "row23_boundary": {
            "future_horizon_condition": "purely ingoing",
            "DtN": "D_HB = -i omega_hat",
            "R_NZH_amplitude": 0.0,
            "epsilon_1": 0.0,
            "Upsilon": 0.0,
        },
        "nonstationary_FLRW": canonical_flrw_tensor_map(),
        "global_response_universality": global_response_operator_map(),
        "claim_boundary": "This closes the canonical EH+R^2 linear tensor map for stationary Kerr and gives the exact nonstationary FLRW friction form. It does not claim R_mn^2/nonminimal phi R extensions or a frozen numerical IR x_R.",
    }



def row19_core_projection() -> dict[str, Any]:
    """Return only the upstream Row-19-owned portion of the unified block."""
    eh = eh_benchmark()
    g = float(eh["metrics"]["gstar"])
    lam = float(eh["metrics"]["lambdastar"])
    fps = find_scalar_fixed_points(g, lam)
    scalar = [classify_scalar_fp(x, g, lam) for x in fps]
    scalar_bfb = [s for s in scalar if s["bounded_below_quartic_truncation"]]
    r2 = r2_variation_audit()
    tensor = canonical_tensor_map()
    checks = {
        "eh_audited_fixed_point_passes": bool(eh["pass"]),
        "scalar_fixed_points_solved": len(scalar) >= 1,
        "exactly_one_bfb_scalar_branch": len(scalar_bfb) == 1,
        "scalar_bfb_branch_has_stationary_minimum": bool(
            len(scalar_bfb) == 1
            and scalar_bfb[0].get("stationary_vacuum_type") == "Z2_broken_pair"
            and scalar_bfb[0].get("fluctuation_m2_over_k2", -1.0) > 0.0
        ),
        "r2_metric_variation_factor_two_issue_identified": bool(r2["legacy_factor_of_two_correction_required"]),
        "canonical_tensor_speed_exactly_luminal": tensor["tensor_quantities"]["cT_squared"] == 1.0,
        "canonical_no_massive_spin2_mode": tensor["tensor_quantities"]["mT_squared"] == 0.0,
        "canonical_no_preferred_frame_alpha1": tensor["tensor_quantities"]["alpha1_preferred_frame"] == 0.0,
        "canonical_stationary_tensor_friction_zero": tensor["tensor_quantities"]["nu_stationary"] == 0.0,
        "canonical_horizon_reflectivity_zero": tensor["row23_boundary"]["R_NZH_amplitude"] == 0.0,
        "canonical_FLRW_tensor_friction_functional_form_closed": tensor["nonstationary_FLRW"]["closure"].startswith("FUNCTIONAL_FORM_CLOSED"),
        "global_response_projection_repair_passes": all(tensor["global_response_universality"]["checks"].values()),
    }
    return {
        "eh_frg": eh,
        "scalar_portal": {
            "fixed_points": scalar,
            "admissible_bfb": scalar_bfb,
            "interpretation": "UV fixed-point branch establishes an admissible stable scalar sector; it is not identified with a numerical astrophysical IR vacuum.",
        },
        "r2_metric_variation_audit": r2,
        "canonical_tensor_map": tensor,
        "checks": checks,
        "pass": all(checks.values()),
    }

def run() -> dict[str, Any]:
    eh = eh_benchmark()
    g = float(eh["metrics"]["gstar"])
    lam = float(eh["metrics"]["lambdastar"])

    fps = find_scalar_fixed_points(g, lam)
    scalar = [classify_scalar_fp(x, g, lam) for x in fps]
    scalar_bfb = [s for s in scalar if s["bounded_below_quartic_truncation"]]

    r2 = r2_variation_audit()
    tensor = canonical_tensor_map()

    nz = gate_nz(str(ROOT / "datasets"))
    row23 = gate_ringdown_echoes(str(ROOT / "datasets"))

    checks = {
        "eh_audited_fixed_point_passes": bool(eh["pass"]),
        "scalar_fixed_points_solved": len(scalar) >= 1,
        "exactly_one_bfb_scalar_branch": len(scalar_bfb) == 1,
        "scalar_bfb_branch_has_stationary_minimum": bool(
            len(scalar_bfb) == 1
            and scalar_bfb[0].get("stationary_vacuum_type") == "Z2_broken_pair"
            and scalar_bfb[0].get("fluctuation_m2_over_k2", -1.0) > 0.0
        ),
        "r2_metric_variation_factor_two_issue_identified": bool(r2["legacy_factor_of_two_correction_required"]),
        "canonical_tensor_speed_exactly_luminal": tensor["tensor_quantities"]["cT_squared"] == 1.0,
        "canonical_no_massive_spin2_mode": tensor["tensor_quantities"]["mT_squared"] == 0.0,
        "canonical_no_preferred_frame_alpha1": tensor["tensor_quantities"]["alpha1_preferred_frame"] == 0.0,
        "canonical_stationary_tensor_friction_zero": tensor["tensor_quantities"]["nu_stationary"] == 0.0,
        "canonical_horizon_reflectivity_zero": tensor["row23_boundary"]["R_NZH_amplitude"] == 0.0,
        "canonical_FLRW_tensor_friction_functional_form_closed": tensor["nonstationary_FLRW"]["closure"].startswith("FUNCTIONAL_FORM_CLOSED"),
        "global_response_projection_repair_passes": all(tensor["global_response_universality"]["checks"].values()),
        "row22_nz_gate_still_passes": nz.get("verdict") == "PASS",
        "row23_kerr_qnm_benchmark_still_passes": row23.get("verdict") == "PASS",
    }

    verdict = "PASS_ROW19_TO_ROW23_CANONICAL_GRAVITY_SCALAR_GW_BLOCK" if all(checks.values()) else "FAIL"
    return {
        "artifact": "UTFANSWF_V23_ROW19_TO_ROW23_CANONICAL_GRAVITY_SCALAR_GW_BLOCK_v1",
        "role": "integrated pre-PDF reproducibility hardening; protected main.tex/main.pdf unchanged",
        "eh_frg": eh,
        "scalar_portal": {
            "fixed_points": scalar,
            "admissible_bfb": scalar_bfb,
            "interpretation": "UV fixed-point branch establishes an admissible stable scalar sector; it is not identified with a numerical astrophysical IR vacuum.",
        },
        "r2_metric_variation_audit": r2,
        "canonical_tensor_map": tensor,
        "row22_handoff": nz,
        "row23_handoff": row23,
        "checks": checks,
        "verdict": verdict,
        "remaining_general_boundary": [
            "Full arbitrary-background curvature-squared FRG including R_mn R^mn is not needed for the canonical Row-23 branch and remains future work.",
            "A nonminimal phi R response is not numerically fixed and is not used in the canonical black-hole tensor branch.",
            "A nonzero reflective NZH branch would require additional derived dynamics; it is not part of the canonical v23 prediction.",
        ],
    }


def main() -> int:
    result = run()
    out = ROOT / "results" / "row19_to_row23_gravity_scalar_gw_v1.json"
    out.write_text(json.dumps(result, indent=2, sort_keys=True) + "\n", encoding="utf-8")
    print(json.dumps({
        "verdict": result["verdict"],
        "g_star": result["eh_frg"]["metrics"]["gstar"],
        "lambda_star": result["eh_frg"]["metrics"]["lambdastar"],
        "n_scalar_fixed_points": len(result["scalar_portal"]["fixed_points"]),
        "n_bfb_scalar_fixed_points": len(result["scalar_portal"]["admissible_bfb"]),
        "r2_legacy_prefactor_ratio": result["r2_metric_variation_audit"]["legacy_appendix_prefactor_relative_to_correct"],
        "cT_squared": result["canonical_tensor_map"]["tensor_quantities"]["cT_squared"],
        "R_NZH": result["canonical_tensor_map"]["row23_boundary"]["R_NZH_amplitude"],
        "row22": result["row22_handoff"]["verdict"],
        "row23": result["row23_handoff"]["verdict"],
    }, indent=2))
    print("Wrote", out)
    return 0 if result["verdict"].startswith("PASS") else 1


if __name__ == "__main__":
    raise SystemExit(main())
