def interchange_comparison(profile: InterchangeProfile, psi_n: Sequence[float], di: Sequence[float],
                           dr: Sequence[float], *, label: str) -> dict[str, np.ndarray]:
    """Equilibrium-native ``D_I``, ``D_R`` against a stability code's ``di``, ``dr``.

    Parameters
    ----------
    profile : InterchangeProfile
        From :func:`mercier_ggj_profile` [-].
    psi_n : sequence of float
        The code's normalized flux grid [-].
    di : sequence of float
        The code's Mercier index on that grid, DCON's sign [-].
    dr : sequence of float
        The code's resistive-interchange index [-].
    label : str
        Which code and run the comparison is against, kept on the result [-].

    Returns
    -------
    comparison : dict
        On the profile's levels: ``psi_n``, ``D_I``, ``di``, ``D_R``, ``dr``,
        the absolute differences ``dD_I``, ``dD_R``, the code's implied
        ``H`` (``0.5 - sqrt(dr - di)``, the root nearer the equilibrium's), and
        ``label`` [-].

    Raises
    ------
    ValueError
        The code's grid is not increasing or the arrays disagree in length.

    Processing steps
    ----------------
    1. Interpolate the code's ``di``, ``dr`` linearly onto the profile's levels
       inside its own range; NaN outside it.
    2. Differences and the code's implied ``H``, choosing the root of
       ``(H - 1/2)^2 = dr - di`` nearer the equilibrium's ``H``.

    Convention
    ----------
    Both in DCON's sign convention; no QA cut, no verdict.

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

    Provenance
    ----------
    .. [1902] Issue #1902, comparison with DCON/RDCON on a reference slice.
    """
    x = np.asarray(psi_n, dtype=float)
    di, dr = np.asarray(di, dtype=float), np.asarray(dr, dtype=float)
    if di.shape != x.shape or dr.shape != x.shape or np.any(np.diff(x) <= 0):
        raise ValueError("psi_n must be increasing and di, dr on it")
    lv = profile.psi_n
    inside = (lv >= x[0]) & (lv <= x[-1])
    # anti-alias: spatial interpolation over normalized flux, not time -- no sample rate to reduce
    di_i = np.where(inside, np.interp(lv, x, di), np.nan)
    dr_i = np.where(inside, np.interp(lv, x, dr), np.nan)
    root = np.sqrt(np.clip(dr_i - di_i, 0.0, None))
    h_code = np.where(np.abs(0.5 - root - profile.H) <= np.abs(0.5 + root - profile.H), 0.5 - root, 0.5 + root)
    return {"psi_n": lv, "D_I": profile.D_I, "di": di_i, "D_R": profile.D_R, "dr": dr_i,
            "dD_I": profile.D_I - di_i, "dD_R": profile.D_R - dr_i, "H": profile.H, "H_code": h_code,
            "label": label}

def mercier_ggj_from_ods(ods: Any, time_slice: int = 0, *, levels_norm: Sequence[float] | None = None,
                         min_points: int = 40, oversample: int = 2) -> InterchangeProfile:
    """Mercier / GGJ indices of one ODS equilibrium slice.

    Parameters
    ----------
    ods : ODS
        Holding ``equilibrium.time_slice[time_slice]`` with ``profiles_2d[0]``
        ``psi`` on a rectangular grid, ``profiles_1d`` ``psi``, ``f``,
        ``pressure``, ``q``, and the axis and boundary flux and axis position [-].
    time_slice : int, optional
        The slice to read [-].
    levels_norm : sequence of float, optional
        Surfaces to evaluate; by default 0.05 to 0.98 in steps of 0.01, the
        range the reference comparison covers [-].
    min_points : int, optional
        Passed to :func:`mercier_ggj_profile`; a numerical-resolution control [-].
    oversample : int, optional
        Passed to :func:`mercier_ggj_profile`; a numerical-resolution control [-].

    Returns
    -------
    profile : InterchangeProfile
        As :func:`mercier_ggj_profile`, with the slice index and time in its
        provenance [-].

    Raises
    ------
    KeyError
        A required path is missing.

    Processing steps
    ----------------
    1. Read the slice's flux map, profiles and axis, refusing a missing path.
    2. Convert the flux to per radian and the profile grid to normalized flux.
    3. Evaluate with :func:`mercier_ggj_profile` and record the slice and time.

    Input semantics
    ---------------
    The ODS flux is converted to per radian with
    :func:`vaft.data.eqdsk.ods_psi_to_wb_per_radian_factor`, which reads the
    record's own convention.

    Output semantics
    ----------------
    Raw equilibrium-derived indices with no QA cut; comparing with a
    stability code is :func:`interchange_comparison`'s job.

    Convention
    ----------
    As :func:`mercier_ggj_profile`: DCON's sign and orientation.

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

    Provenance
    ----------
    .. [1902] Issue #1902, evaluation from one reconstructed equilibrium ODS.
    """
    inputs = _ods_slice_inputs(ods, time_slice)
    levels = np.arange(0.05, 0.985, 0.01) if levels_norm is None else np.asarray(levels_norm, dtype=float)
    profile = mercier_ggj_profile(levels_norm=levels, min_points=min_points, oversample=oversample, **inputs)
    return InterchangeProfile(**{**profile.__dict__, "provenance": {**profile.provenance, "time_slice": int(time_slice),
                                                                     "time": _slice_time(ods, time_slice)}})

def s_alpha_from_ods(ods: Any, time_slice: int = 0, *, levels_norm: Sequence[float] | None = None,
                     R0: float | None = None) -> SAlphaProfile:
    """Per-surface s-alpha coordinates of one ODS equilibrium slice.

    Parameters
    ----------
    ods : ODS
        As :func:`mercier_ggj_from_ods` [-].
    time_slice : int, optional
        The slice to read [-].
    levels_norm : sequence of float, optional
        Surfaces; by default 0.05 to 0.98 in steps of 0.01 [-].
    R0 : float, optional
        Reference major radius; the magnetic axis by default (a convention) [m].

    Returns
    -------
    profile : SAlphaProfile
        As :func:`s_alpha_profile`, with the slice and time in provenance [-].

    Raises
    ------
    KeyError
        A required path is missing.

    Processing steps
    ----------------
    1. Read the slice (per-radian flux) as :func:`mercier_ggj_from_ods`.
    2. Evaluate with :func:`s_alpha_profile` and record the slice and time.

    Input semantics
    ---------------
    One reconstructed equilibrium slice; nothing else of the ODS is read.

    Output semantics
    ----------------
    Raw ``s``, ``alpha`` and the model's boundaries; no verdict beyond the
    model's own criterion.

    Convention
    ----------
    As :func:`s_alpha_profile`.

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

    Provenance
    ----------
    .. [1902] Issue #1902, ODS-facing per-surface s-alpha evaluation.
    """
    inputs = _ods_slice_inputs(ods, time_slice)
    levels = np.arange(0.05, 0.985, 0.01) if levels_norm is None else np.asarray(levels_norm, dtype=float)
    profile = s_alpha_profile(levels_norm=levels, R0=R0, **inputs)
    return SAlphaProfile(**{**profile.__dict__, "provenance": {**profile.provenance, "time_slice": int(time_slice),
                                                               "time": _slice_time(ods, time_slice)}})