"""Structural preflight for the AP1-R2c implicit Lobatto DAE."""
from __future__ import annotations

import json
from pathlib import Path

import numpy as np

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, match_seed_N, standard_stress
from chi_background_closure import ChiParameters, physical_shells
from planck2018_neutrino_closure import Planck2018Pilot


def lobatto_grid_and_D(nodes: int, delta_N: float) -> tuple[np.ndarray, np.ndarray]:
    if nodes < 5 or not delta_N > 0.0:
        raise ValueError("valid Lobatto grid required")
    j = np.arange(nodes)
    x = -np.cos(np.pi * j / (nodes - 1))
    weights = (-1.0) ** j
    weights[[0, -1]] *= 0.5
    differences = x[:, None] - x[None, :]
    D_x = np.zeros((nodes, nodes))
    mask = ~np.eye(nodes, dtype=bool)
    ratios = weights[None, :] / weights[:, None]
    np.divide(ratios, differences, out=D_x, where=mask)
    D_x[np.diag_indices(nodes)] = -np.sum(D_x, axis=1)
    N = 0.5 * delta_N * (x + 1.0)
    return N, (2.0 / delta_N) * D_x


def positive_friedmann_root(total_density: np.ndarray | float) -> np.ndarray:
    total = np.asarray(total_density, dtype=float)
    if np.any(~np.isfinite(total)) or np.any(total <= 0.0):
        raise ValueError("positive finite Friedmann density required")
    return np.sqrt(total / 3.0)


def dae_dimensions(nodes: int) -> dict[str, int | bool]:
    unknowns = 7 * nodes + 1
    equations = 4 * nodes + 3 * nodes + 1
    return {"unknowns": unknowns, "equations": equations,
            "square": unknowns == equations}


def build_preflight(state_path: Path, nodes: int = 17,
                    delta_N: float = 1.0e-7) -> dict:
    N, D = lobatto_grid_and_D(nodes, delta_N)
    x = 2.0 * N / delta_N - 1.0
    polynomial_errors = {}
    for degree in range(5):
        values = x**degree
        # Compare in the dimensionless Lobatto coordinate to avoid attaching
        # an artificial absolute tolerance to the 2/delta_N scale.
        numerical = (0.5 * delta_N) * (D @ values)
        exact = np.zeros_like(x) if degree == 0 else degree * x**(degree - 1)
        scale = max(float(np.max(np.abs(exact))), 1.0)
        polynomial_errors[str(degree)] = float(np.max(np.abs(numerical - exact)) / scale)

    frozen = decode_hard_pass_state(state_path)
    seed = frozen[-1, :5]
    prefix = inherited_trajectory(state_path, 1025)
    p = ChiParameters()
    k, weights = physical_shells(p.Lambda, nodes=8)
    inherited = completed_history_decimal(prefix, k, weights, 60)
    rho_q = float(inherited["completed_history"]["rho"][-1])
    pilot = Planck2018Pilot()
    N_seed = match_seed_N(seed, rho_q, pilot)
    rho_std = float(standard_stress(np.array([N_seed]), pilot)[0][0])
    total = classical_terms(seed)["rho"] + rho_std + rho_q
    H_root = float(positive_friedmann_root(total))
    H_relative_error = abs(H_root / float(seed[4]) - 1.0)
    dims = dae_dimensions(nodes)
    gates = {
        "lobatto_polynomials_0_to_4_below_1e_10": max(polynomial_errors.values()) < 1.0e-10,
        "square_residual_system": bool(dims["square"]),
        "positive_seed_friedmann_root": H_root > 0.0,
        "seed_friedmann_root_relative_error_below_1e_10": H_relative_error < 1.0e-10,
        "old_solver_or_physical_map_called": False,
    }
    passed = (
        all(value for key, value in gates.items()
            if key != "old_solver_or_physical_map_called")
        and not gates["old_solver_or_physical_map_called"]
    )
    return {
        "schema": "apeiron-ap1-r2c-implicit-dae-preflight-v1.0",
        "classification": "IMPLICIT_DAE_FORMULATION_PREFLIGHT_PASS_SOLVE_NOT_RUN" if passed else "IMPLICIT_DAE_PREFLIGHT_NOT_PASS",
        "grid": {"nodes": nodes, "delta_N": delta_N,
                 "min_spacing": float(np.min(np.diff(N))),
                 "max_spacing": float(np.max(np.diff(N)))},
        "polynomial_derivative_relative_errors": polynomial_errors,
        "dimensions": dims,
        "seed": {"N_physical": N_seed, "rho_q": rho_q,
                 "H_frozen_Mpl": float(seed[4]), "H_DAE_root_Mpl": H_root,
                 "H_root_relative_error": H_relative_error},
        "gates": gates,
        "next_required": "assemble full 7n+1 residual and safeguarded Newton-Krylov solve",
        "claim_boundary": "structural DAE preflight only; no candidate PASS or observable released",
    }


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(build_preflight(args.state), indent=2) + "\n"
    if args.output:
        args.output.write_text(rendered, encoding="utf-8")
    else:
        print(rendered, end="")


if __name__ == "__main__":
    main()