"""New AP1-M1 short candidate with inherited high-precision chi sources.

This code never imports the old v7.13 solver or Physical Map.  The frozen
HARD-PASS background is read only as the mode-history prefix.  A new short
N-segment is evolved with minimally coupled standard fluids and iterated
renormalized chi sources.  It is a local candidate gate, not H(z).
"""
from __future__ import annotations

from dataclasses import asdict, dataclass
import json
from pathlib import Path

import numpy as np
from scipy.integrate import solve_ivp
from scipy.interpolate import CubicSpline
from scipy.optimize import brentq

from ap1_m1_background_preflight import decode_hard_pass_state
from ap1_r2c_high_precision_modes import (
    completed_history_decimal,
    propagate_endpoint_decimal,
)
from ap1_r2c_mode_inheritance import inherited_trajectory
from chi_background_closure import (
    ChiParameters,
    ChiTrajectory,
    curvature_uv_tail_n,
    frozen_renormalization_constants,
    order0_uv_tail_n,
    physical_shells,
    portal_terms,
)
from planck2018_neutrino_closure import Planck2018Pilot, shared_standard_background


VS_COEFF = np.array([1.1971212004620583e-11, 9.87022611239063e-14,
                     -1.0336424891583054e-13, 1.189490499106138e-10,
                     1.2236365189926054e-10])


@dataclass(frozen=True)
class CandidateConfig:
    delta_N: float = 1.0e-7
    candidate_nodes: int = 17
    inherited_nodes: int = 1025
    fixed_point_steps: int = 16
    relaxation: float = 0.5
    mode_nodes: int = 8
    digits: int = 60
    rtol: float = 2.0e-9
    atol: float = 2.0e-12

    def validate(self) -> None:
        if not (0.0 < self.delta_N <= 0.01) or self.candidate_nodes < 17:
            raise ValueError("bounded short candidate grid required")
        if self.inherited_nodes < 129 or self.fixed_point_steps < 1:
            raise ValueError("insufficient inheritance or fixed-point depth")
        if not (0.0 < self.relaxation <= 1.0) or self.digits < 32:
            raise ValueError("invalid relaxation or precision")


def _F(z: float, b: float) -> float:
    return 0.25 * z**4 - 0.5 * z**2 + b * z


def _Fz(z: float, b: float) -> float:
    return z**3 - z + b


def _vs(theta: float) -> float:
    return float(sum(VS_COEFF[i] * theta**i for i in range(5)))


def _vs_prime(theta: float) -> float:
    return float(sum(i * VS_COEFF[i] * theta ** (i - 1) for i in range(1, 5)))


def classical_terms(state: np.ndarray) -> dict[str, float]:
    sigma, sd, theta, td, _H = map(float, state)
    Lambda = 2.3e-3; f_sigma = 2.5e-2; Ec = 0.04 * 3.0 * (2.0e-6)**2
    X = 0.5 * td**2; z = sigma / f_sigma; l4 = Lambda**4
    b = 0.35 - 4.78e-2 * theta - 5.0 * X / l4
    kinetic = X + 5.0 * (z - 2.0 / np.sqrt(3.0))**2 * X**2 / (3.0 * l4)
    P = kinetic - Ec * _F(z, b) - _vs(theta)
    PX = 1.0 + 10.0 * (z - 2.0 / np.sqrt(3.0))**2 * X / (3.0 * l4) + 5.0 * Ec * z / l4
    PXX = 10.0 * (z - 2.0 / np.sqrt(3.0))**2 / (3.0 * l4)
    P_sigma = (10.0 * (z - 2.0 / np.sqrt(3.0)) * X**2 / (3.0 * l4) - Ec * _Fz(z, b)) / f_sigma
    P_theta = -Ec * z * (-4.78e-2) - _vs_prime(theta)
    PX_sigma = (20.0 * (z - 2.0 / np.sqrt(3.0)) * X / (3.0 * l4) + 5.0 * Ec / l4) / f_sigma
    rho = 0.5 * sd**2 + 2.0 * X * PX - P
    pressure = 0.5 * sd**2 + P
    return {"rho": rho, "pressure": pressure, "P_sigma": P_sigma,
            "P_theta": P_theta, "PX": PX, "PXX": PXX,
            "PX_sigma": PX_sigma, "K": PX + 2.0 * X * PXX}


def standard_stress(N: np.ndarray, pilot: Planck2018Pilot) -> tuple[np.ndarray, np.ndarray]:
    f = shared_standard_background(np.asarray(N), pilot)
    return (f["rho_b"] + f["rho_c"] + f["rho_gamma"] + f["rho_nu"],
            f["rho_gamma"] / 3.0 + f["p_nu"])


def match_seed_N(seed: np.ndarray, rho_q: float, pilot: Planck2018Pilot) -> float:
    allowance = 3.0 * seed[4]**2 - classical_terms(seed)["rho"] - rho_q
    if allowance <= 0.0:
        raise RuntimeError("non-positive M1 standard-density allowance at seed")
    residual = lambda N: float(standard_stress(np.array([N]), pilot)[0][0] - allowance)
    return float(brentq(residual, -100.0, 0.0, xtol=1.0e-14, rtol=1.0e-14))


def integrate_segment(seed: np.ndarray, Nrel: np.ndarray, N_seed: float,
                      q: dict[str, np.ndarray], pilot: Planck2018Pilot,
                      cfg: CandidateConfig) -> tuple[np.ndarray, np.ndarray]:
    chi_p = ChiParameters()
    def qi(key: str, n: float) -> float:
        return float(np.interp(n, Nrel, q[key]))
    def rhs(n: float, y: np.ndarray) -> np.ndarray:
        sigma, sd, theta, td, H = y
        if H <= 0.0:
            raise RuntimeError("candidate left positive-H branch")
        d = classical_terms(y)
        rho_s, p_s = standard_stress(np.array([N_seed + n]), pilot)
        _m2, dm_s, dm_t = portal_terms(sigma, theta, chi_p)
        chi2 = qi("chi2", n)
        sdd = -3.0 * H * sd + d["P_sigma"] - 0.5 * float(dm_s) * chi2
        tdd = (d["P_theta"] - 0.5 * float(dm_t) * chi2 - 3.0 * H * d["PX"] * td - d["PX_sigma"] * sd * td) / d["K"]
        Hdot = -0.5 * (d["rho"] + d["pressure"] + qi("rho", n) + qi("pressure", n) + float(rho_s[0] + p_s[0]))
        return np.array([sd / H, sdd / H, td / H, tdd / H, Hdot / H])
    sol = solve_ivp(rhs, (0.0, float(Nrel[-1])), seed, t_eval=Nrel,
                    method="DOP853", rtol=cfg.rtol, atol=cfg.atol,
                    max_step=cfg.delta_N / (cfg.candidate_nodes - 1))
    if not sol.success or sol.y.shape[1] != len(Nrel):
        raise RuntimeError(f"candidate integration failed: {sol.message}")
    y = sol.y.T
    Hdot = np.empty(len(Nrel))
    for i, (n, row) in enumerate(zip(Nrel, y)):
        d = classical_terms(row); rho_s, p_s = standard_stress(np.array([N_seed + n]), pilot)
        Hdot[i] = -0.5 * (d["rho"] + d["pressure"] + q["rho"][i] + q["pressure"][i] + rho_s[0] + p_s[0])
    return y, Hdot


def completed_candidate_sources_decimal(
    full: ChiTrajectory,
    candidate: ChiTrajectory,
    k: np.ndarray,
    weights: np.ndarray,
    digits: int,
) -> tuple[dict[str, np.ndarray], dict]:
    """Combine inherited modes with local UV terms on the candidate mesh.

    The resolved modes require the entire inherited trajectory.  In contrast,
    the adiabatic UV orders are local differential functionals and must not be
    differentiated across the strongly nonuniform inheritance junction.
    """
    p = ChiParameters()
    mode_result = propagate_endpoint_decimal(
        full.validated(), k, weights, digits, record_history=True
    )
    count = len(candidate.N)
    resolved = {
        key: np.asarray(mode_result["resolved_pv_history"][key][-count:])
        for key in ("rho", "pressure", "chi2")
    }
    tail0 = order0_uv_tail_n(candidate, p)
    tail24 = curvature_uv_tail_n(candidate, p, nodes_per_octave=6)
    constants = frozen_renormalization_constants(p)
    mass2 = np.asarray(portal_terms(candidate.sigma, candidate.theta, p)[0])
    C = constants["C_chi2"]
    A = constants["A_g"]
    B = constants["B_G_rhs"]
    completed = {
        "rho": resolved["rho"] + tail0["rho"] + tail24["rho"]
        - 0.5 * C * mass2 + A + 3.0 * B * candidate.H**2,
        "pressure": resolved["pressure"] + tail0["pressure"]
        + tail24["pressure"] + 0.5 * C * mass2 - A
        - B * (2.0 * candidate.Hdot + 3.0 * candidate.H**2),
        "chi2": resolved["chi2"] + tail0["chi2"] + tail24["chi2"] - C,
    }
    return completed, mode_result


def run_candidate(state_path: Path, cfg: CandidateConfig | None = None) -> dict:
    cfg = CandidateConfig() if cfg is None else cfg; cfg.validate()
    pilot = Planck2018Pilot(); chi_p = ChiParameters()
    frozen = decode_hard_pass_state(state_path); seed = frozen[-1, :5].copy()
    prefix = inherited_trajectory(state_path, cfg.inherited_nodes)
    k, weights = physical_shells(chi_p.Lambda, nodes=cfg.mode_nodes)
    inherited_sources = completed_history_decimal(prefix, k, weights, cfg.digits)["completed_history"]
    Nrel = np.linspace(0.0, cfg.delta_N, cfg.candidate_nodes)
    q = {key: np.full(cfg.candidate_nodes, float(inherited_sources[key][-1])) for key in ("rho", "pressure", "chi2")}
    iteration_log = []
    for iteration in range(cfg.fixed_point_steps):
        q_used = {key: value.copy() for key, value in q.items()}
        N_seed = match_seed_N(seed, float(q["rho"][0]), pilot)
        y, Hdot = integrate_segment(seed, Nrel, N_seed, q_used, pilot, cfg)
        full = ChiTrajectory(
            N=np.concatenate([prefix.N, prefix.N[-1] + Nrel[1:]]),
            H=np.concatenate([prefix.H, y[1:, 4]]),
            Hdot=np.concatenate([prefix.Hdot, Hdot[1:]]),
            sigma=np.concatenate([prefix.sigma, y[1:, 0]]),
            theta=np.concatenate([prefix.theta, y[1:, 2]]),
        )
        candidate = ChiTrajectory(
            N=prefix.N[-1] + Nrel, H=y[:, 4], Hdot=Hdot,
            sigma=y[:, 0], theta=y[:, 2],
        )
        fresh, mode_result = completed_candidate_sources_decimal(
            full, candidate, k, weights, cfg.digits
        )
        change = max(float(np.max(np.abs(fresh[key] - q_used[key])) / max(np.max(np.abs(fresh[key])), 1.0e-300)) for key in q)
        q = {key: cfg.relaxation * fresh[key] + (1.0 - cfg.relaxation) * q_used[key] for key in q}
        iteration_log.append({"iteration": iteration + 1, "N_seed": N_seed,
                              "source_relative_change": change,
                              "wronskian_relative_error": mode_result["wronskian_relative_error"]})

    # One unrelaxed validation evaluation: background and Friedmann use the
    # same q input; Ward uses the freshly returned source on that background.
    q_background = {key: value.copy() for key, value in q.items()}
    N_seed = match_seed_N(seed, float(q_background["rho"][0]), pilot)
    y, Hdot = integrate_segment(seed, Nrel, N_seed, q_background, pilot, cfg)
    full = ChiTrajectory(
        N=np.concatenate([prefix.N, prefix.N[-1] + Nrel[1:]]),
        H=np.concatenate([prefix.H, y[1:, 4]]),
        Hdot=np.concatenate([prefix.Hdot, Hdot[1:]]),
        sigma=np.concatenate([prefix.sigma, y[1:, 0]]),
        theta=np.concatenate([prefix.theta, y[1:, 2]]),
    )
    candidate = ChiTrajectory(
        N=prefix.N[-1] + Nrel, H=y[:, 4], Hdot=Hdot,
        sigma=y[:, 0], theta=y[:, 2],
    )
    q_fresh, validation_modes = completed_candidate_sources_decimal(
        full, candidate, k, weights, cfg.digits
    )
    validation_change = max(
        float(np.max(np.abs(q_fresh[key] - q_background[key])) /
              max(np.max(np.abs(q_fresh[key])), 1.0e-300))
        for key in q_background
    )
    dlist = [classical_terms(row) for row in y]
    rho_s, p_s = standard_stress(N_seed + Nrel, pilot)
    total = np.array([d["rho"] for d in dlist]) + rho_s + q_background["rho"]
    friedmann = (3.0 * y[:, 4]**2 - total) / max(float(np.max(np.abs(total))), 1.0e-300)
    mass2 = np.asarray(portal_terms(y[:, 0], y[:, 2], chi_p)[0])
    rho_dot = y[:, 4] * CubicSpline(Nrel, q_fresh["rho"])(Nrel, 1)
    mass2_dot = y[:, 4] * CubicSpline(Nrel, mass2)(Nrel, 1)
    ward = rho_dot + 3.0 * y[:, 4] * (q_fresh["rho"] + q_fresh["pressure"]) - 0.5 * mass2_dot * q_fresh["chi2"]
    ward_scale = max(
        float(np.max(np.abs(rho_dot))),
        float(np.max(np.abs(3.0 * y[:, 4] * (q_fresh["rho"] + q_fresh["pressure"])))),
        float(np.max(np.abs(0.5 * mass2_dot * q_fresh["chi2"]))),
        1.0e-300,
    )
    return {
        "schema": "apeiron-ap1-r2c-self-consistent-short-candidate-v1.0",
        "classification": "SHORT_CANDIDATE_RUN_G22_NOT_YET_PASS",
        "config": asdict(cfg), "old_solver_or_physical_map_called": False,
        "seed": {"N_physical": N_seed, "implied_z": float(np.exp(-N_seed) - 1.0)},
        "iterations": iteration_log,
        "diagnostics": {
            "finite": bool(np.all(np.isfinite(y))), "positive_H": bool(np.all(y[:, 4] > 0.0)),
            "min_PX": float(min(d["PX"] for d in dlist)), "min_K": float(min(d["K"] for d in dlist)),
            "max_abs_friedmann_residual": float(np.max(np.abs(friedmann))),
            "ward_normalized": float(np.max(np.abs(ward)) / ward_scale),
            "validation_source_relative_change": validation_change,
            "validation_wronskian_relative_error": validation_modes["wronskian_relative_error"],
            "H_endpoint_Mpl": float(y[-1, 4]), "delta_H_over_H": float(y[-1, 4] / y[0, 4] - 1.0),
        },
        "claim_boundary": "short local candidate only; no z=0 reach and no H(z), growth, lensing, fit, or empirical claim",
    }


def main() -> None:
    import argparse
    parser = argparse.ArgumentParser(); parser.add_argument("state", type=Path); parser.add_argument("--output", type=Path)
    args = parser.parse_args(); rendered = json.dumps(run_candidate(args.state), indent=2) + "\n"
    if args.output: args.output.write_text(rendered, encoding="utf-8")
    else: print(rendered, end="")


if __name__ == "__main__":
    main()