def assess_neo_readiness(state: ResolvedTransportState,
                         surfaces: Sequence[float] = DEFAULT_SURFACES, *,
                         projection: Optional["TransportProjection"] = None) -> ReadinessReport:
    """Decide whether a resolved state can go to NEO.

    Parameters
    ----------
    state : ResolvedTransportState
        From :func:`resolve_transport_state` [-].
    surfaces : sequence of float
        The r/a NEO will solve; with an inferred-Ti gap fill, at least one must have
        an input independent of the fill [-].
    projection : TransportProjection, optional
        A ``neoclassical`` projection of this state; NEO accepts only explicit
        species, which is the state's own list [-].

    Returns
    -------
    ReadinessReport
        ``ready``/``conditional``/``insufficient``; NEO is a profile code, so no
        per-surface entries [-].

    Raises
    ------
    ValueError
        A projection of another state, for another physics, or a reduced one.

    Processing steps
    ----------------
    1. An unresolved state is ``insufficient``.
    2. ``GACODEProfile.check_neo_requirements`` (the executable-input contract) must
       report nothing missing.
    3. Declared assumptions make the state ``conditional``.

    Applicability
    -------------
    Machine-independent.
    """
    if not state.resolved:
        return ReadinessReport("neo", "insufficient", reasons=state.reasons)
    state, projected = _projected_state(state, projection, "neo")
    missing = tuple(state.profile.check_neo_requirements())
    if missing:
        return ReadinessReport("neo", "insufficient",
                               reasons=tuple(f"missing_{name}" for name in missing))
    if state.fill_probe is not None and not any(
            inferred_ti_supported(state, r, solver="neo") for r in surfaces):
        return ReadinessReport("neo", "insufficient", reasons=("no_surface_independent_of_ti_fill",))
    conditions = _state_conditions(state) + projected
    return ReadinessReport("neo", "conditional" if conditions else "ready",
                           conditions=conditions)

def assess_tglf_readiness(
    state: ResolvedTransportState,
    surfaces: Sequence[float] = DEFAULT_SURFACES,
    *,
    config: Any = None,
    projection: Optional["TransportProjection"] = None,
) -> ReadinessReport:
    """Decide whether a resolved state, and each requested surface, can go to TGLF.

    Parameters
    ----------
    state : ResolvedTransportState
        From :func:`resolve_transport_state` [-].
    surfaces : sequence of float
        Requested surfaces as r/a [-].
    config : TGLFConfig, optional
        Only ``n_species`` is consulted, through ``prepare_tglf_input`` [-].
    projection : TransportProjection, optional
        A ``turbulence`` projection of this state (:func:`project_transport_state`);
        the surfaces are then built on its profile, and a reduced projection adds
        ``species_projection_<method>`` to ``conditions``.  Default: the state's own
        species list, unchanged [-].

    Returns
    -------
    ReadinessReport
        State verdict plus one :class:`SurfaceReadiness` per surface, each carrying
        the built ``TGLFInput`` when ready [-].

    Raises
    ------
    ValueError
        A projection of another state, for another physics, or by a method TGLF
        does not accept.

    Processing steps
    ----------------
    1. An unresolved state is ``insufficient`` and no surface is built; a projection
       is checked against the state and the solver.
    2. Each surface is projected with ``prepare_tglf_input``; its documented refusal
       becomes a reason code, then ``check_tglf_requirements`` and finite gradients
       are required.
    3. A runnable state with declared assumptions is ``conditional``.

    Applicability
    -------------
    Machine-independent.
    """
    if not state.resolved:
        return ReadinessReport("tglf", "insufficient", reasons=state.reasons)
    state, projected = _projected_state(state, projection, "tglf")
    from vaft.code.gacode.tglf.inputs import LocalConversionError, prepare_tglf_input

    results = []
    for radius in surfaces:
        try:
            local = prepare_tglf_input(state.profile, float(radius), config=config)
        except LocalConversionError as error:
            results.append(SurfaceReadiness(float(radius), _surface_code(error), str(error)))
            continue
        if not inferred_ti_supported(state, float(radius), local=local):
            results.append(SurfaceReadiness(float(radius), "ti_not_inferred_here",
                                             "this surface's solver input depends on the inferred-Ti gap fill"))
            continue
        absent = local.check_tglf_requirements()
        if absent:
            code = "invalid_gradient" if set(absent) & {"rlns", "rlts"} else "local_conversion_failure"
            results.append(SurfaceReadiness(float(radius), code, ", ".join(absent), local))
            continue
        results.append(SurfaceReadiness(float(radius), "ready", "", local))
    if not any(s.ready for s in results):
        return ReadinessReport("tglf", "insufficient", reasons=("no_ready_surface",),
                               surfaces=tuple(results))
    conditions = _state_conditions(state) + projected
    if any(s.ready and s.local_input.vexb_shear is None for s in results):
        # prepare_tglf_input does not derive ExB shear yet (#553), so TGLF runs with
        # its own zero whatever rotation core_profiles carries.
        conditions = (*conditions, "no_exb_shear")
    return ReadinessReport("tglf", "conditional" if conditions else "ready",
                           conditions=conditions, surfaces=tuple(results))

def inferred_ti_supported(state: "ResolvedTransportState", r_over_a: float, *,
                          solver: str = "tglf", local: Any = None) -> bool:
    """Whether a surface's solver input is independent of the inferred-Ti gap fill.

    Parameters
    ----------
    state : ResolvedTransportState
        From :func:`resolve_transport_state` [-].
    r_over_a : float
        The surface [-].
    solver : str
        ``"tglf"`` compares the full TGLF local input; ``"neo"`` compares the ion
        temperatures and their log-gradients the way NEO reads them [-].
    local : TGLFInput, optional
        The surface's input already built from ``state.profile``, to avoid rebuilding
        it [-].

    Returns
    -------
    bool
        True when no gap was filled, or when the solver's input agrees (rtol 1e-4,
        atol 1e-4) between the profile and its ``fill_probe`` [-].

    Convention
    ----------
    GACODE differentiates on three grid points and evaluates with a spline over the
    whole profile, so a filled node can reach a surface several grid steps away.
    Rather than guess a margin, the gate builds the solver's input twice -- once from
    the profile, once from ``fill_probe``, whose gaps are filled with a different
    scale *and* offset -- and runs the surface only when nothing moves. TGLF's local
    input uses a not-a-knot spline; NEO (PROFILE_MODEL=2, ``expro_locsim``'s
    ``cub_spline1``) a natural one, so each solver is checked with its own.

    Applicability
    -------------
    Machine-independent.
    """
    probe = getattr(state, "fill_probe", None)
    if probe is None:
        return True
    if solver == "neo":
        return _neo_ion_inputs_agree(state.profile, probe, float(r_over_a))
    from vaft.code.gacode.tglf.inputs import LocalConversionError, prepare_tglf_input

    try:
        a = local if local is not None else prepare_tglf_input(state.profile, float(r_over_a))
        b = prepare_tglf_input(probe, float(r_over_a))
    except LocalConversionError:
        return False
    numeric = ("taus", "rlts", "as_", "rlns", "zs", "mass", "betae", "xnue", "zeff", "debye",
               "p_prime_loc", "q_loc", "q_prime_loc", "rmin_loc", "rmaj_loc", "kappa_loc",
               "delta_loc")
    for name in numeric:
        x = np.atleast_1d(np.asarray(getattr(a, name), dtype=float))
        y = np.atleast_1d(np.asarray(getattr(b, name), dtype=float))
        if x.shape != y.shape or not np.allclose(x, y, rtol=1e-4, atol=1e-4, equal_nan=True):
            return False
    return True

def physics_parameters(config: Any, *, exclude: Iterable[str] = ()) -> dict[str, Any]:
    """The physics half of a GACODE config, with runtime-only fields dropped.

    ``backend``, ``timeout``, ``env``, ``home``, ``executable``, ``platform`` and the
    MPI/OMP counts decide *where* a run happens, not *what* it computes, so they stay
    out of identity.

    Parameters
    ----------
    config : GACODEConfig
        A ``TGLFConfig`` or ``NEOConfig`` dataclass [-].
    exclude : iterable of str
        Further field names to leave out [-].

    Returns
    -------
    dict
        Field name to JSON-ready value, for :func:`run_identity` [-].

    Applicability
    -------------
    Machine-independent.
    """
    from dataclasses import fields as dataclass_fields

    runtime = {"backend", "timeout", "env", "home", "executable", "workdir", "args",
               "platform", "n_mpi", "n_omp", "memory_mb", *exclude}
    out = {}
    for entry in dataclass_fields(config):
        if entry.name in runtime:
            continue
        out[entry.name] = _jsonable(getattr(config, entry.name))
    return out

def project_transport_state(state: ResolvedTransportState, target_physics: str,
                            method: Optional[str] = None) -> TransportProjection:
    """One physics projection of a resolved state's species, with the profile its solver reads.

    Parameters
    ----------
    state : ResolvedTransportState
        A resolved state [-].
    target_physics : str
        A target of :data:`vaft.process.species.PROJECTION_POLICY`, e.g.
        ``"turbulence"`` (TGLF), ``"neoclassical"`` (NEO), ``"resistive"`` [-].
    method : str, optional
        A method the target allows (``"explicit"``, ``"effective_impurity"``, ...);
        default the target's own [-].

    Returns
    -------
    TransportProjection
        The species contract, ``canonical_state_id`` (the state's identity) and
        ``species_state_id``, and the projected profile for an ``explicit`` or
        ``effective_impurity`` method (``None`` for ``moments_only``) [any].

    Raises
    ------
    ValueError
        An unresolved state; a target or method the species policy refuses (an
        effective impurity for NEO); an effective impurity whose charge or mass
        varies with radius, which one GACODE species cannot carry.

    Processing steps
    ----------------
    1. :func:`transport_species_state` from the resolved profile.
    2. :func:`vaft.process.species.project_species_state` on that state.
    3. ``explicit``: the state's own profile, unchanged; ``effective_impurity``:
       the profile with the thermal impurities replaced by one pseudo-ion
       (charge, mass constant over the grid; temperature density-weighted);
       the inferred-Ti probe, if any, gets the same replacement.

    Applicability
    -------------
    Machine-independent.

    Limitations
    -----------
    The effective impurity keeps the charge density and the Z^2 moment at every
    point, not the impurity thermal pressure: its density
    ``(sum n Z)^2 / sum n Z^2`` is below ``sum n``, so the impurity ``n T`` (and
    a gyrokinetic ``p'``) is slightly lower than the explicit species give;
    ``ptot`` is left as the state carries it.

    Provenance
    ----------
    .. [issue] #1567 Sec. 6-7 and Phase C; policy from #1709.
    """
    from vaft.process.species import project_species_state

    species_state = transport_species_state(state)
    projection = project_species_state(species_state, target_physics, method)
    profile = probe = None
    if projection.method == "explicit":
        profile, probe = state.profile, state.fill_probe
    elif projection.method == "effective_impurity":
        profile = _species_profile(state.profile, species_state, projection)
        if state.fill_probe is not None:
            probe = _species_profile(state.fill_probe, species_state, projection)
    return TransportProjection(canonical_state_id=state.identity, species_state_id=species_state.state_id,
                               projection=projection, profile=profile, fill_probe=probe)

def run_identity(
    state: ResolvedTransportState,
    *,
    solver: str,
    parameters: Mapping[str, Any],
    surface: Optional[float] = None,
    solver_revision: Optional[str] = None,
    projection: Optional["TransportProjection"] = None,
) -> str:
    """Deterministic identity of one solver calculation on one resolved state.

    Parameters
    ----------
    state : ResolvedTransportState
        The upstream state; its own identity is folded in [-].
    solver : str
        ``"tglf"``, ``"neo"`` or ``"classical"`` [-].
    parameters : mapping
        The solver's effective settings (e.g. ``dataclasses.asdict(config)`` minus
        runtime-only fields) [-].
    surface : float, optional
        r/a of a local calculation [-].
    solver_revision : str, optional
        The GACODE revision the run used or will use [-].
    projection : TransportProjection, optional
        The species projection the run used (#1567).  A reduced projection's
        ``projection_id`` is folded in, so an effective-impurity run differs from
        the explicit one; an explicit projection is the state's own species list
        and gives the identity this function always gave [-].

    Returns
    -------
    str
        sha256 hex digest; equal inputs give equal identities, which is the cache key [-].

    Raises
    ------
    ValueError
        A projection of another state.

    Applicability
    -------------
    Machine-independent.
    """
    payload = {
        "state": state.identity,
        "solver": solver,
        "parameters": dict(parameters),
        "surface": None if surface is None else round(float(surface), 6),
        "solver_revision": solver_revision,
    }
    if projection is not None:
        if projection.canonical_state_id != state.identity:
            raise ValueError("the projection was made from another transport state")
        if projection.method != "explicit":
            # an explicit projection is the state's own species list: same run, same key
            payload["species_projection"] = projection.projection_id
    return _digest(payload)

def surface_toroidal_field(profile: Any, local: Any) -> float:
    """The toroidal field at a surface's centre, ``|B_centr| R_centr / R``.

    GACODE's ``B_unit`` is a flux-coordinate field (about twice B0 on VEST), not the
    field a gyro-orbit sees, so the classical baseline does not use it.

    Parameters
    ----------
    profile : GACODEProfile
        Supplies ``bcentr`` [T] and ``rcentr`` [m].
    local : TGLFInput
        Supplies ``rmaj_loc`` and the minor radius of the surface [-].

    Returns
    -------
    float
        ``|B_T|`` at ``R = rmaj_loc * a``; the poloidal field is neglected [T].

    Applicability
    -------------
    Machine-independent.

    Provenance
    ----------
    .. [1] Vacuum toroidal field ``B_T = B_0 R_0 / R``; GACODE ``input.gacode``
       ``bcentr``/``rcentr`` (https://gacode.io/input_gacode.html).
    """
    r_surface = float(local.rmaj_loc) * float(local.normalisation.minor_radius)
    return abs(float(profile.bcentr)) * float(profile.rcentr) / r_surface

def transport_partition(
    neoclassical: Mapping[str, Any],
    turbulent: Mapping[str, Any],
    *,
    turbulent_charges: Optional[Mapping[str, float]] = None,
    classical: Optional[Mapping[str, Any]] = None,
    rho_tolerance: float = 1e-3,
) -> dict[str, Any]:
    """Signed model fluxes and their magnitude partition at one surface of one state.

    Parameters
    ----------
    neoclassical : mapping
        One mapped NEO surface: ``r_over_a``, ``rho_tor_norm`` and electron/ion energy
        and particle fluxes in SI, ions keyed ``z=<charge>`` [W/m^2].
    turbulent : mapping
        One mapped TGLF surface in the same shape, ions keyed by label [W/m^2].
    turbulent_charges : mapping, optional
        Species label -> charge for the TGLF row (``{"H+": 1, "C6+": 6}``) [e].
    classical : mapping, optional
        A classical row in the same shape, added as a third component when given [-].
    rho_tolerance : float
        Largest allowed rho_tor_norm disagreement between the two rows [-].

    Returns
    -------
    dict
        ``status`` ``"available"`` with per-channel ``neo``, ``turb``, (``classical``),
        ``model`` (the modelled sum, not a measurement), ``f_neo`` =
        |neo| / (|neo| + |turb| [+ |classical|]) and ``neo_over_turb`` (``None`` below a
        floor), or ``"unavailable"`` with a reason [-].

    Applicability
    -------------
    Machine-independent. Callers join only rows with one ``state_identity``.
    """
    if neoclassical is None or turbulent is None:
        missing = "neoclassical" if neoclassical is None else "turbulent"
        return {"status": "unavailable", "reason": f"missing_{missing}_component"}
    if abs(float(neoclassical["r_over_a"]) - float(turbulent["r_over_a"])) > 1e-4:
        return {"status": "unavailable", "reason": "surfaces_differ_in_r_over_a"}
    if abs(float(neoclassical["rho_tor_norm"]) - float(turbulent["rho_tor_norm"])) > rho_tolerance:
        # Same r/a but a different flux surface means two different profiles.
        return {"status": "unavailable", "reason": "surfaces_differ_in_rho_tor_norm"}
    if classical is not None and abs(float(classical["r_over_a"]) - float(turbulent["r_over_a"])) > 1e-4:
        return {"status": "unavailable", "reason": "classical_surface_differs"}
    neo = _channels(neoclassical)
    turb = _channels(turbulent, turbulent_charges)
    cl = _channels(classical) if classical is not None else {}
    # The ion total is a channel only when both models carry the same ion species;
    # a classical row (main ion only) does not take part in it.
    for kind in ("energy", "particle"):
        ion_n = {k for k in neo if k.startswith(f"ion_{kind}_z=")}
        ion_t = {k for k in turb if k.startswith(f"ion_{kind}_z=")}
        if ion_n and ion_n == ion_t:
            neo[f"ion_{kind}_total"] = _ion_total(neo, kind)
            turb[f"ion_{kind}_total"] = _ion_total(turb, kind)
    channels = {}
    for name in sorted(set(neo) & set(turb)):
        kind = "energy" if "_energy" in name else "particle"
        n, t = neo[name], turb[name]
        parts = [abs(n), abs(t)] + ([abs(cl[name])] if name in cl else [])
        denominator = sum(parts)
        entry = {
            "neo": n, "turb": t,
            "model": n + t + (cl[name] if name in cl else 0.0),
            "f_neo": None if denominator == 0.0 else abs(n) / denominator,
            "f_turb": None if denominator == 0.0 else abs(t) / denominator,
            "neo_over_turb": None if abs(t) < _RATIO_FLOOR[kind] else n / t,
        }
        if name in cl:
            entry["classical"] = cl[name]
            entry["f_classical"] = None if denominator == 0.0 else abs(cl[name]) / denominator
        channels[name] = entry
    if not channels:
        return {"status": "unavailable", "reason": "no_common_channel"}
    return {"status": "available", "r_over_a": float(turbulent["r_over_a"]),
            "rho_tor_norm": float(turbulent["rho_tor_norm"]),
            "components": ["neo", "turb"] + (["classical"] if classical is not None else []),
            "channels": channels}

def transport_species_state(state: ResolvedTransportState):
    """The canonical species state of a resolved transport state: the species its solvers read.

    Parameters
    ----------
    state : ResolvedTransportState
        A resolved state from :func:`resolve_transport_state` [-].

    Returns
    -------
    vaft.process.species.CanonicalSpeciesState
        One component per species of ``state.profile`` (element and mass, charge,
        population from the GACODE type, density [m^-3] and temperature [eV] on the
        profile's ``rho``), with ``n_e`` and the state time [any].

    Raises
    ------
    ValueError
        An unresolved state, a profile that is itself a projection, or a species
        whose label and mass match no element.

    Convention
    ----------
    A species the converter lumped from a bundled ion (``provenance["ni"]
    ["bundled_ions"]``) or with a non-integer charge is a bundled component.  The
    element comes from the label only where its standard weight lies within 1 u
    of the stored mass and its Z_n bounds the charge, else from the mass.  The
    density origin is ``assumed`` for a policy composition closure and
    ``unlabelled`` otherwise; the state's own record says which.

    Applicability
    -------------
    Machine-independent.

    Provenance
    ----------
    .. [issue] #1567 Phase C (transport-state integration).
    """
    from vaft.process.species import CanonicalSpeciesState, SpeciesComponent

    if not state.resolved:
        raise ValueError(f"the transport state is not resolved ({', '.join(state.reasons)})")
    profile = state.profile
    if "species_projection" in (getattr(profile, "provenance", {}) or {}):
        raise ValueError("this profile is already a species projection; take the canonical state "
                         "from the resolved state it was projected from")
    ni = np.atleast_2d(np.asarray(profile.ni, dtype=float))
    ti = np.atleast_2d(np.asarray(profile.ti, dtype=float))
    charges = np.atleast_1d(np.asarray(profile.z, dtype=float))
    masses = np.atleast_1d(np.asarray(profile.mass, dtype=float))
    names = list(profile.name) or [f"ion{k}" for k in range(charges.size)]
    kinds = list(profile.type) or ["[therm]"] * charges.size
    record = dict(getattr(profile, "provenance", {}) or {}).get("ni", {})
    bundled = set((record or {}).get("bundled_ions", {}) or {})
    origin = "assumed" if state.composition.get("origin") == "policy_assumption" else "unlabelled"
    components = []
    for k, label in enumerate(names):
        charge = float(charges[k])
        nuclide, mass = _profile_nuclide(label, charge, float(masses[k]))
        is_bundled = label in bundled or abs(charge - round(charge)) > 1e-9
        components.append(SpeciesComponent(
            nuclide, charge, _POPULATION_OF_TYPE.get(str(kinds[k]).strip(), "other"),
            ni[k] * _DENSITY_SCALE, mass, temperature=ti[min(k, ti.shape[0] - 1)] * _TEMPERATURE_SCALE,
            bundled=is_bundled, origin=origin,
            source=f"transport state {state.identity[:12]}, GACODE species {label}"))
    return CanonicalSpeciesState(
        tuple(components), rho=np.asarray(profile.rho, dtype=float),
        n_e=np.asarray(profile.ne, dtype=float) * _DENSITY_SCALE, time=float(state.key.time_efit_s),
        provenance={"source": "resolved transport state profile", "state_identity": state.identity,
                    "composition": dict(state.composition)})