def qualify_response_operator(
    operator: ResonantResponseOperator,
    phasors: np.ndarray,
    responses: Sequence[Sequence[complex]],
    *,
    sample_ids: Sequence[str],
    rtol: float,
) -> ResponseOperatorQualification:
    """Qualify or reject ``Phi = H conj(c)`` against held-out runs.

    Parameters
    ----------
    operator : ResonantResponseOperator
        From :func:`resonant_response_operator` [-].
    phasors : numpy.ndarray
        ``(N, K)`` excitations not among the unit runs, complex ones included [phasor].
    responses : sequence of sequence of complex
        Their direct ``Phi``; an all-NaN row marks a failed run [T].
    sample_ids : sequence of str
        One id per run [-].
    rtol : float
        Largest accepted ``||Phi_model - Phi_direct|| / ||Phi_direct||`` [-].

    Returns
    -------
    qualification : ResponseOperatorQualification
        ``rejected`` when any run exceeds ``rtol``; ``insufficient`` when any
        failed or none was given, or no run had a non-real phasor (the
        conjugate relation would be untested); otherwise ``qualified`` [-].

    Raises
    ------
    ValueError
        Shapes disagree or ``rtol`` is negative.

    Processing steps
    ----------------
    1. Set aside failed (non-finite) runs and count them.
    2. Relative vector error of ``H conj(c)`` against each direct ``Phi``.
    3. One status, requiring at least one run with a non-real phasor.

    Convention
    ----------
    The comparison is on the complex vector, phase included; matching the
    magnitudes alone would not catch the ``conj`` relation.

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

    Provenance
    ----------
    .. [1884] Issue #1884, independent combined-current validation.
    """
    if not rtol >= 0:
        raise ValueError("rtol must not be negative")
    c = np.atleast_2d(np.asarray(phasors, dtype=complex))
    direct = np.atleast_2d(np.asarray(responses, dtype=complex))
    ids = [str(s) for s in sample_ids]
    if c.shape != (len(ids), operator.H.shape[1]) or direct.shape != (len(ids), operator.H.shape[0]):
        raise ValueError("phasors, responses and sample_ids disagree with the operator")
    finite = np.all(np.isfinite(direct), axis=1)
    model = operator.response(c[finite])
    errors = {i: float(np.linalg.norm(m - d) / np.linalg.norm(d))
              for i, m, d in zip(np.asarray(ids)[finite], model, direct[finite])}
    reasons: list[str] = []
    failed = [i for i, e in errors.items() if e > rtol]
    if failed:
        status = "rejected"
        reasons.append(f"relative error above {rtol:g}: {failed}")
    elif not finite.all() or not finite.any():
        status = "insufficient"
        reasons.append(f"{int((~finite).sum())} run(s) failed; the operator is identified, not validated")
    elif not np.any(np.abs(c[finite].imag) > 0):
        status = "insufficient"
        reasons.append("no held-out run had a non-real phasor, so Phi = H conj(c) is untested")
    else:
        status = "qualified"
    return ResponseOperatorQualification(status=status, relative_errors=errors, reasons=tuple(reasons))

def resonant_response_operator(
    unit_responses: Sequence[Sequence[complex]],
    *,
    field: str,
    psi_n_rational: Sequence[float],
    m_rational: Sequence[int],
    n: int,
    groups: Sequence[str],
    amplitude_unit: float,
    provenance: Mapping[str, object] | None = None,
) -> ResonantResponseOperator:
    """The operator ``H`` from the responses to the unit excitations ``c = e_k``.

    Parameters
    ----------
    unit_responses : sequence of sequence of complex
        ``Phi`` of the run with ``c = e_k`` (group ``k`` at one phasor unit,
        phase 0, the others off), one per group, in group order [T].
    field : str
        ``Phi_res`` or ``Phi_res_v`` [-].
    psi_n_rational : sequence of float
        Surface positions, the same in every run [-].
    m_rational : sequence of int
        Poloidal mode number per surface [-].
    n : int
        Toroidal mode number [-].
    groups : sequence of str
        Coil-group names in phasor order [-].
    amplitude_unit : float
        Cosine amplitude of one phasor unit [A per turn].
    provenance : mapping, optional
        Equilibrium, code version, run identities [-].

    Returns
    -------
    operator : ResonantResponseOperator
        ``H[:, k] = Phi(e_k)``, so ``Phi = H conj(c)`` [-].

    Raises
    ------
    ValueError
        Unknown field, mismatched lengths, a non-finite response, or no group.

    Convention
    ----------
    Because ``e_k`` is real, ``Phi(e_k) = H e_k``; the conjugate shows only
    for complex excitations.  :func:`qualify_response_operator` is what
    establishes the relation for complex ``c``.

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

    Limitations
    -----------
    Built from ``K`` runs; trustworthy only once qualified on held-out complex
    runs at the same equilibrium, ``n`` and numerical settings.

    Provenance
    ----------
    .. [1884] Issue #1884, the verified linear operator.
    """
    if field not in RESONANT_FIELDS:
        raise ValueError(f"field must be one of {RESONANT_FIELDS}, got {field!r}")
    columns = [np.asarray(r, dtype=complex).ravel() for r in unit_responses]
    psi = np.asarray(psi_n_rational, dtype=float).ravel()
    m = np.asarray(m_rational, dtype=int).ravel()
    if not columns or len(columns) != len(groups):
        raise ValueError(f"{len(columns)} unit responses for {len(groups)} groups")
    if any(col.shape != psi.shape for col in columns) or m.shape != psi.shape:
        raise ValueError("every unit response must have one value per rational surface")
    h = np.column_stack(columns)
    if not np.all(np.isfinite(h)):
        raise ValueError("a unit response is not finite")
    if not amplitude_unit > 0:
        raise ValueError("amplitude_unit must be positive")
    return ResonantResponseOperator(field=field, H=h, psi_n_rational=psi, m_rational=m, n=int(n),
                                    groups=tuple(str(g) for g in groups), amplitude_unit=float(amplitude_unit),
                                    provenance=dict(provenance or {}))

def window_weights(operator: ResonantResponseOperator, window) -> np.ndarray:
    """Surface weights whose quadratic metric is the window's mean squared ``|Phi|``.

    Parameters
    ----------
    operator : ResonantResponseOperator
        Supplies the surface positions [-].
    window : vaft.process.perturbation.ResonantWindow
        From :func:`vaft.process.perturbation.resonant_windows` [-].

    Returns
    -------
    W : numpy.ndarray
        Diagonal ``N_s x N_s``, ``1 / N_in`` on the surfaces inside [-].

    Raises
    ------
    ValueError
        The window holds no surface: an empty window is not a zero response.

    Convention
    ----------
    Endpoints included, as :meth:`ResonantWindow.mask`; a surface on the
    core/edge boundary is in both.  ``sqrt(c^H Q c)`` with ``Q =
    operator.quadratic(W)`` is :func:`vaft.process.perturbation.rms_resonant_field`.

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

    Provenance
    ----------
    .. [1884] Issue #1884, Sec. 3 (core/edge metrics).
    """
    inside = window.mask(operator.psi_n_rational)
    if not inside.any():
        raise ValueError(f"the {window.name} window [{window.low:g}, {window.high:g}] holds no rational surface")
    return np.diag(inside / inside.sum()).astype(float)