"""Fail-closed AP1-M1 Friedmann--Raychaudhuri junction-match pilot.

The frozen v7.13 trajectory is read only.  This module does not import or run
the old solver or Physical Map.  It evaluates the inherited, renormalized
homogeneous chi source and asks whether a standard M-1 density can satisfy the
Friedmann and Raychaudhuri junction conditions simultaneously.

All roots produced here are pilot brackets.  They are not released as seeds:
matching tolerances and a finer independent validation level must be frozen
before a junction state can be accepted for the new long-background solve.
"""
from __future__ import annotations

from dataclasses import asdict, dataclass
from datetime import datetime, timezone
from hashlib import sha256
import json
from pathlib import Path
from typing import Iterable

import numpy as np
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
from ap1_r2c_mode_inheritance import inherited_trajectory
from ap1_r2c_self_consistent_candidate import classical_terms, standard_stress
from chi_background_closure import ChiParameters, physical_shells, portal_terms
from planck2018_neutrino_closure import Planck2018Pilot


EXPECTED_AUTHORITY_SHA256 = {
    "HANDOFF/CURRENT_STATE_LATEST.npz": (
        "26473f1ce826f7272af22bea2ce92d313f914b7b8ca9601d5f7047e79de9c81f"
    ),
    "AP1/APEIRON_AP1_M1_BACKGROUND_SPEC_LATEST.md": (
        "9a565aef648f4a7229941bf9efa5c51a754fd4763dd11e59e8032701d6ed0259"
    ),
    "AP1/APEIRON_AP1_R2C_G22_PASS_CHECKPOINT_20260902T151240Z.json": (
        "13cb421bc457a59c149ae68840829b867f3637f92acfdff533ce665ff6f26c6e"
    ),
    "AP1/CODE/ap1_m1_background_preflight.py": (
        "8b5754bd0c89f6b6ec77fa7a8e042ac153168ac445ba932d862b142254718475"
    ),
    "AP1/CODE/planck2018_neutrino_closure.py": (
        "ccd3158930dcb9ab74d7634ed06f2eaa00745ac95e2883f9fe02d9039a056acc"
    ),
    "AP1/CODE/chi_background_closure.py": (
        "1097888c72fbe9152e1905a2a3ab37c15df9ca681639571b87d76de6ddf072a7"
    ),
    "AP1/CODE/ap1_r2c_mode_inheritance.py": (
        "221d3f111eb0eb1435d64c62cc4de117a5464052a9345281d98c42f7fccff1cb"
    ),
    "AP1/CODE/ap1_r2c_high_precision_modes.py": (
        "3b76dc1bba339fc7576381661b7a0099501231154f45cf9bb2eba3db16f0d484"
    ),
    "AP1/CODE/ap1_r2c_self_consistent_candidate.py": (
        "2c22454f69817ecdc40d38dacad0efaa72fdddcfa50c9c00a48478d610ed8cb0"
    ),
}


class BoundaryMatchError(RuntimeError):
    """A junction-match prerequisite or numerical domain check failed."""


@dataclass(frozen=True)
class MatchPilotConfig:
    levels: tuple[int, ...] = (1025, 2049, 4097)
    mode_nodes: int = 8
    decimal_digits: int = 60
    physical_N_lower: float = -100.0
    physical_N_upper: float = 0.0

    def validate(self, frozen_nodes: int) -> None:
        if len(self.levels) < 1 or tuple(sorted(set(self.levels))) != self.levels:
            raise ValueError("strictly increasing unique pilot levels required")
        if min(self.levels) < 17 or max(self.levels) > frozen_nodes:
            raise ValueError("pilot levels must lie inside the frozen trajectory")
        if self.mode_nodes < 4 or self.decimal_digits < 32:
            raise ValueError("insufficient momentum or Decimal resolution")
        if not self.physical_N_lower < self.physical_N_upper <= 0.0:
            raise ValueError("physical N bracket must end no later than today")


def file_sha256(path: Path) -> str:
    digest = sha256()
    with path.open("rb") as stream:
        for block in iter(lambda: stream.read(1024 * 1024), b""):
            digest.update(block)
    return digest.hexdigest()


def audit_authorities(apeiron_root: Path) -> dict[str, str]:
    observed: dict[str, str] = {}
    for relative, expected in EXPECTED_AUTHORITY_SHA256.items():
        path = apeiron_root / relative
        if not path.is_file():
            raise BoundaryMatchError(f"missing authority: {relative}")
        value = file_sha256(path)
        if value != expected:
            raise BoundaryMatchError(f"authority hash mismatch: {relative}")
        observed[relative] = value
    return observed


def sampled_rows(frozen: np.ndarray, nodes: int) -> tuple[np.ndarray, np.ndarray]:
    if not 17 <= nodes <= len(frozen):
        raise ValueError("sample count outside frozen trajectory")
    indices = np.unique(
        np.round(np.linspace(0, len(frozen) - 1, nodes)).astype(int)
    )
    if len(indices) != nodes:
        raise BoundaryMatchError("requested sampling did not remain one-to-one")
    return indices, frozen[indices]


def friedmann_matched_physical_N(
    required_standard_density: float,
    pilot: Planck2018Pilot,
    lower: float = -100.0,
    upper: float = 0.0,
) -> float:
    if not np.isfinite(required_standard_density) or required_standard_density <= 0.0:
        raise BoundaryMatchError("positive finite standard density is required")

    def residual(N_value: float) -> float:
        rho, _pressure = standard_stress(np.array([N_value]), pilot)
        return float(rho[0] - required_standard_density)

    left = residual(lower)
    right = residual(upper)
    if left * right > 0.0:
        raise BoundaryMatchError("required density is outside the shared M-1 bracket")
    return float(brentq(residual, lower, upper, xtol=1.0e-13, rtol=1.0e-14))


def interpolate_sign_change_roots(
    N_local: np.ndarray,
    N_physical: np.ndarray,
    raychaudhuri_residual: np.ndarray,
    normalized_residual: np.ndarray,
    frequency_margin: np.ndarray,
    H: np.ndarray,
    required_density: np.ndarray,
) -> list[dict[str, float | int | bool]]:
    arrays = [
        np.asarray(value, dtype=float)
        for value in (
            N_local,
            N_physical,
            raychaudhuri_residual,
            normalized_residual,
            frequency_margin,
            H,
            required_density,
        )
    ]
    if any(value.ndim != 1 or value.shape != arrays[0].shape for value in arrays):
        raise ValueError("matching arrays must be one-dimensional and aligned")
    roots: list[dict[str, float | int | bool]] = []
    for index in range(len(arrays[0]) - 1):
        left = arrays[2][index]
        right = arrays[2][index + 1]
        if not (np.isfinite(left) and np.isfinite(right)):
            continue
        if left == 0.0:
            fraction = 0.0
        elif left * right < 0.0:
            fraction = float(-left / (right - left))
        else:
            continue

        def interp(values: np.ndarray) -> float:
            return float(values[index] + fraction * (values[index + 1] - values[index]))

        margin = interp(arrays[4])
        roots.append(
            {
                "left_index": index,
                "right_index": index + 1,
                "fraction": fraction,
                "N_local": interp(arrays[0]),
                "N_physical": interp(arrays[1]),
                "H_Mpl": interp(arrays[5]),
                "required_standard_density_Mpl4": interp(arrays[6]),
                "vacuum_frequency_margin_Mpl2": margin,
                "positive_vacuum_frequency_margin": bool(margin > 0.0),
                "left_raychaudhuri_residual_Mpl2": float(left),
                "right_raychaudhuri_residual_Mpl2": float(right),
                "left_normalized_residual": float(arrays[3][index]),
                "right_normalized_residual": float(arrays[3][index + 1]),
                "bracket_width_N_local": float(arrays[0][index + 1] - arrays[0][index]),
            }
        )
    return roots


def select_latest_positive_frequency_root(
    roots: Iterable[dict[str, float | int | bool]],
) -> dict[str, float | int | bool] | None:
    admissible = [
        root for root in roots if bool(root["positive_vacuum_frequency_margin"])
    ]
    if not admissible:
        return None
    return max(admissible, key=lambda item: float(item["N_local"]))


def evaluate_level(
    state_path: Path,
    frozen: np.ndarray,
    nodes: int,
    config: MatchPilotConfig,
) -> dict:
    indices, rows = sampled_rows(frozen, nodes)
    trajectory = inherited_trajectory(state_path, nodes)
    if not np.array_equal(trajectory.N, rows[:, 5]):
        raise BoundaryMatchError("trajectory sampling and frozen rows diverged")
    parameters = ChiParameters()
    pilot = Planck2018Pilot()
    momenta, weights = physical_shells(parameters.Lambda, nodes=config.mode_nodes)
    mode_result = completed_history_decimal(
        trajectory,
        momenta,
        weights,
        digits=config.decimal_digits,
    )
    quantum = mode_result["completed_history"]

    N_physical = np.full(nodes, np.nan)
    raychaudhuri = np.full(nodes, np.nan)
    normalized = np.full(nodes, np.nan)
    required_density = np.full(nodes, np.nan)
    frequency_margin = np.asarray(
        portal_terms(rows[:, 0], rows[:, 2], parameters)[0], dtype=float
    ) - 2.25 * rows[:, 4] ** 2

    for index, row in enumerate(rows):
        classical = classical_terms(row[:5])
        required = (
            3.0 * row[4] ** 2 - classical["rho"] - quantum["rho"][index]
        )
        required_density[index] = required
        try:
            matched_N = friedmann_matched_physical_N(
                float(required),
                pilot,
                config.physical_N_lower,
                config.physical_N_upper,
            )
        except BoundaryMatchError:
            continue
        rho_standard, pressure_standard = standard_stress(
            np.array([matched_N]), pilot
        )
        enthalpy = (
            classical["rho"]
            + classical["pressure"]
            + quantum["rho"][index]
            + quantum["pressure"][index]
            + rho_standard[0]
            + pressure_standard[0]
        )
        residual = float(trajectory.Hdot[index] + 0.5 * enthalpy)
        scale = max(
            abs(float(trajectory.Hdot[index])),
            abs(float(0.5 * enthalpy)),
            float(row[4] ** 2),
            1.0e-300,
        )
        N_physical[index] = matched_N
        raychaudhuri[index] = residual
        normalized[index] = residual / scale

    roots = interpolate_sign_change_roots(
        rows[:, 5],
        N_physical,
        raychaudhuri,
        normalized,
        frequency_margin,
        rows[:, 4],
        required_density,
    )
    selected = select_latest_positive_frequency_root(roots)
    finite_matches = np.isfinite(N_physical)
    return {
        "nodes": nodes,
        "sampled_stored_indices": {
            "first": int(indices[0]),
            "last": int(indices[-1]),
        },
        "friedmann_density_matches": int(np.sum(finite_matches)),
        "simultaneous_sign_change_brackets": len(roots),
        "positive_frequency_brackets": sum(
            bool(root["positive_vacuum_frequency_margin"]) for root in roots
        ),
        "selected_by_preregistered_rule": selected,
        "wronskian_relative_error": float(mode_result["wronskian_relative_error"]),
        "finite": bool(
            np.all(np.isfinite(rows))
            and all(
                np.all(np.isfinite(np.asarray(quantum[key], dtype=float)))
                for key in ("rho", "pressure", "chi2")
            )
        ),
        "root_candidates": roots,
    }


def selection_refinement(levels: list[dict]) -> list[dict[str, float]]:
    output: list[dict[str, float]] = []
    for coarse, fine in zip(levels[:-1], levels[1:]):
        left = coarse["selected_by_preregistered_rule"]
        right = fine["selected_by_preregistered_rule"]
        if left is None or right is None:
            continue
        output.append(
            {
                "coarse_nodes": int(coarse["nodes"]),
                "fine_nodes": int(fine["nodes"]),
                "abs_delta_N_local": abs(
                    float(left["N_local"]) - float(right["N_local"])
                ),
                "abs_delta_N_physical": abs(
                    float(left["N_physical"]) - float(right["N_physical"])
                ),
                "relative_delta_H": abs(
                    float(left["H_Mpl"]) - float(right["H_Mpl"])
                )
                / max(abs(float(right["H_Mpl"])), 1.0e-300),
            }
        )
    return output


def build_match_pilot_report(
    apeiron_root: Path,
    config: MatchPilotConfig | None = None,
) -> dict:
    config = MatchPilotConfig() if config is None else config
    authorities = audit_authorities(apeiron_root)
    state_path = apeiron_root / "HANDOFF/CURRENT_STATE_LATEST.npz"
    frozen = decode_hard_pass_state(state_path)
    config.validate(len(frozen))
    levels = [
        evaluate_level(state_path, frozen, nodes, config) for nodes in config.levels
    ]
    refinements = selection_refinement(levels)
    gates = {
        "authority_hashes_exact": True,
        "all_levels_finite": all(level["finite"] for level in levels),
        "friedmann_density_match_exists_each_level": all(
            level["friedmann_density_matches"] > 0 for level in levels
        ),
        "simultaneous_Friedmann_Raychaudhuri_bracket_each_level": all(
            level["simultaneous_sign_change_brackets"] > 0 for level in levels
        ),
        "positive_frequency_bracket_each_level": all(
            level["positive_frequency_brackets"] > 0 for level in levels
        ),
        "deterministic_selection_each_level": all(
            level["selected_by_preregistered_rule"] is not None for level in levels
        ),
        "Wronskian_below_existing_1e_11_each_level": all(
            level["wronskian_relative_error"] < 1.0e-11 for level in levels
        ),
        "old_solver_or_physical_map_not_called": True,
        "pilot_candidate_not_released_as_seed": True,
    }
    passed = bool(all(gates.values()))
    return {
        "schema": "apeiron-ap1-m1-junction-match-pilot-v1.0",
        "updated_utc": datetime.now(timezone.utc).strftime("%Y-%m-%dT%H:%M:%SZ"),
        "classification": (
            "M1_JUNCTION_MATCH_BRACKET_PREFLIGHT_PASS_TOLERANCE_FREEZE_AND_BACKGROUND_SOLVE_OPEN"
            if passed
            else "M1_JUNCTION_MATCH_BRACKET_PREFLIGHT_NOT_PASS"
        ),
        "config": asdict(config),
        "authority_sha256": authorities,
        "matching_equations": {
            "Friedmann": "rho_std_required = 3 H^2 - rho_AP_classical - rho_q",
            "Raychaudhuri": "R_J = Hdot + 1/2[(rho+p)_AP_classical + (rho+p)_q + (rho+p)_std]",
            "selection_rule": "latest simultaneous sign-change bracket with positive inherited-vacuum frequency margin; maximize use of the frozen HARD-PASS prefix without entering its tachyonic onset",
        },
        "levels": levels,
        "selected_root_refinement": refinements,
        "gates": gates,
        "all_structural_pilot_gates_pass": passed,
        "matching_tolerances_frozen": False,
        "selected_candidate_status": "PILOT_ONLY_NOT_A_CHECKPOINT_OR_SEED",
        "new_AP1_M1_background_started": False,
        "physical_response_kernel_started": False,
        "nonpass_stored_or_used_as_seed": False,
        "equations_changed": False,
        "physics_changed": False,
        "parameters_changed": False,
        "existing_gate_thresholds_changed": False,
        "next_required": "derive and freeze junction tolerances from these pilot refinements plus an independent Decimal roundoff audit; validate on a finer held-out level before accepting any M1 seed",
        "claim_boundary": "junction-bracket pilot only; no accepted matching state, late-time background, response kernel, observable, fit or significance",
    }


def main() -> None:
    import argparse

    parser = argparse.ArgumentParser()
    parser.add_argument("apeiron_root", type=Path)
    parser.add_argument("--output", type=Path)
    parser.add_argument("--levels", type=int, nargs="+", default=[1025, 2049, 4097])
    parser.add_argument("--mode-nodes", type=int, default=8)
    parser.add_argument("--digits", type=int, default=60)
    args = parser.parse_args()
    report = build_match_pilot_report(
        args.apeiron_root,
        MatchPilotConfig(
            levels=tuple(args.levels),
            mode_nodes=args.mode_nodes,
            decimal_digits=args.digits,
        ),
    )
    rendered = json.dumps(report, indent=2) + "\n"
    if args.output is None:
        print(rendered, end="")
    else:
        args.output.write_text(rendered, encoding="utf-8")


if __name__ == "__main__":
    main()