def canonical_sign(V: np.ndarray) -> np.ndarray:
    """Flip each column so its largest-magnitude entry is positive.

    An eigenvector is defined up to sign; without a rule two runs of the same
    build can disagree, and a stored basis cannot be compared with a fresh one.

    Parameters
    ----------
    V : numpy.ndarray
        Modes as columns, ``(n, m)`` [any].

    Returns
    -------
    numpy.ndarray
        The same modes with each column's sign fixed [any].

    Convention
    ----------
    The sign rule of this module (``SIGN_RULE = "max_abs_positive"``): the
    largest-magnitude entry of every column is made positive, ties going to the
    lowest index (``argmax`` on ``|v|``).  A pure function of the column, so
    repeated runs agree bitwise.

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

    Limitations
    -----------
    Two entries of exactly equal magnitude and opposite sign make the rule
    index-dependent rather than ill-defined; a near-degenerate *eigenspace* is
    the real ambiguity, and :func:`build_wall_mode_basis` refuses that case
    instead of signing it.

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    """
    V = np.array(V, dtype=float, copy=True)
    if V.ndim != 2:
        raise WallModeError("modes must be a 2-D array of columns")
    pivots = np.argmax(np.abs(V), axis=0)
    signs = np.sign(V[pivots, np.arange(V.shape[1])])
    signs[signs == 0.0] = 1.0
    return V * signs[None, :]

def check_wall_mode_basis(basis: WallModeBasis, R_mat: np.ndarray, M_mat: np.ndarray) -> dict[str, Any]:
    """Metrics of a basis against the matrices it claims to diagonalize; no verdict.

    Numbers a caller or a report can threshold as it sees fit; nothing here
    decides whether a basis is good enough.

    Parameters
    ----------
    basis : WallModeBasis
        The eigenbasis [-].
    R_mat : numpy.ndarray
        Diagonal loop resistance, ``(N, N)`` or its diagonal [Ohm].
    M_mat : numpy.ndarray
        Passive-passive inductance in code naming, ``(N, N)`` [H].

    Returns
    -------
    dict
        ``n_elements``, ``n_modes``, ``r_r_identity_error``
        (``max|V^T R V - I|``), ``l_r_symmetry_error``,
        ``max_segment_residual``, ``min_relative_gap`` and ``coupling`` [-].

        ``coupling`` holds, per segment pair, the Frobenius norm of that
        off-diagonal block of ``L_r`` relative to the geometric mean of the two
        diagonal blocks, so a reader can see which segments actually talk to
        each other [-].

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    """
    ops = reduced_operators(basis, R_mat, M_mat)
    M_tot = ops.L_r.shape[0]
    identity_error = float(np.max(np.abs(ops.R_r - np.eye(M_tot))))
    offsets = np.cumsum([0] + [seg.size for seg in basis.segments])
    coupling: dict[str, float] = {}
    for i, a in enumerate(basis.segments):
        for j, b in enumerate(basis.segments):
            if j <= i:
                continue
            block = ops.L_r[offsets[i]:offsets[i + 1], offsets[j]:offsets[j + 1]]
            scale = np.sqrt(
                np.linalg.norm(ops.L_r[offsets[i]:offsets[i + 1], offsets[i]:offsets[i + 1]])
                * np.linalg.norm(ops.L_r[offsets[j]:offsets[j + 1], offsets[j]:offsets[j + 1]])
            )
            coupling[f"{a.id}-{b.id}"] = float(np.linalg.norm(block) / scale) if scale > 0 else float("nan")
    return {
        "n_elements": basis.n_elements,
        "n_modes": basis.n_modes(),
        "r_r_identity_error": identity_error,
        "l_r_symmetry_error": float(np.max(np.abs(ops.L_r - ops.L_r.T)) / np.max(np.abs(ops.L_r))),
        "max_segment_residual": max(seg.residual for seg in basis.segments),
        "min_relative_gap": min(seg.min_relative_gap for seg in basis.segments),
        "coupling": coupling,
    }

def combined_operators(
    V: np.ndarray,
    R_mat: np.ndarray,
    M_mat: np.ndarray,
    L_mat: np.ndarray | None = None,
    *,
    label: str = "pattern",
) -> ReducedWall:
    """Reduced operators for an arbitrary R-orthonormal basis.

    The same projections as :func:`reduced_operators`, for a basis that is not
    (only) segment eigenmodes -- an enrichment by :func:`moment_patterns`, say,
    or a POD basis.

    Input semantics
    ---------------
    Element space: the assembled wall matrices and a basis whose columns are wall
    patterns.

    Output semantics
    ----------------
    Mode space: operators of the basis's own coordinates, labelled ``(label, k)``
    rather than by segment.

    Parameters
    ----------
    V : numpy.ndarray
        R-orthonormal basis as columns, ``(N, m)`` [Ohm**-0.5].
    R_mat : numpy.ndarray
        Diagonal loop resistance, ``(N, N)`` or its diagonal [Ohm].
    M_mat : numpy.ndarray
        Passive-passive inductance in code naming, ``(N, N)`` [H].
    L_mat : numpy.ndarray or None, optional
        Passive-to-source coupling in code naming, ``(N, n_src)`` [H].
    label : str, optional
        Label component of every coefficient's ``(label, k)`` name [-].

    Returns
    -------
    ReducedWall
        ``L_r`` [H], ``R_r`` [Ohm], ``M_r`` [H] when a source coupling was given,
        and the labels [-].

    Convention
    ----------
    ``R_r`` is computed, not assumed: a basis that is not R-orthonormal in the
    matrices passed here shows up as a non-identity rather than silently changing
    the dynamics.  Use :func:`orthonormalize_r` to make one that is.

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    """
    V = np.asarray(V, dtype=float)
    r = _diagonal_resistance(R_mat, where="R_mat")
    if V.ndim != 2 or V.shape[0] != r.size:
        raise WallModeError(f"basis has shape {V.shape}; expected ({r.size}, m)")
    L = np.asarray(M_mat, dtype=float)
    return ReducedWall(
        L_r=V.T @ L @ V,
        R_r=V.T @ (r[:, None] * V),
        M_r=None if L_mat is None else V.T @ np.asarray(L_mat, dtype=float),
        labels=tuple((str(label), k) for k in range(V.shape[1])),
        keep=(np.arange(V.shape[1], dtype=np.int64),),
    )

def global_time_constants(basis: WallModeBasis, M_mat: np.ndarray) -> np.ndarray:
    """The whole wall's decay times, from the full-rank reduced inductance.

    With ``R_r = I`` the reduced pencil is an ordinary symmetric eigenproblem, so
    ``eigvalsh(L_r)`` are the global L/R times.  They are not the segments' local
    times: the difference between the two is exactly the inter-segment coupling
    the basis keeps.

    Parameters
    ----------
    basis : WallModeBasis
        The eigenbasis [-].
    M_mat : numpy.ndarray
        Passive-passive inductance in code naming, ``(N, N)`` [H].

    Returns
    -------
    numpy.ndarray
        Global decay times, descending, ``(N,)`` [s].

    Convention
    ----------
    Full rank by construction: every mode of every segment takes part, because a
    truncated basis would report the decay times of a different wall.

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    """
    V = basis.V(None)
    L_r = V.T @ np.asarray(M_mat, dtype=float) @ V
    L_r = 0.5 * (L_r + L_r.T)
    return np.sort(np.linalg.eigvalsh(L_r))[::-1]

def moment_patterns(R_mat: np.ndarray, M_mat: np.ndarray, L_mat: np.ndarray, order: int = 1) -> np.ndarray:
    """Drive-independent wall patterns from the source coupling.

    An R-orthonormal basis of the block Krylov space
    ``span{R^{-1} M, (R^{-1} L) R^{-1} M, ...}`` up to ``order`` blocks.  The
    first block is the resistive limit -- the wall current a constant source ramp
    settles into, ``I_w = -R^{-1} M dI_src/dt`` -- and each further block is the
    next inductive correction of the slowly driven response: the Laplace-domain
    moments of the wall's transfer function, matched at zero frequency.

    Built as a block Arnoldi iteration in the R inner product, so the blocks stay
    independent where the raw powers would collapse onto the slowest mode; a
    block that adds no new direction ends the iteration early and the basis is
    narrower than ``order * n_src``.

    Parameters
    ----------
    R_mat : numpy.ndarray
        Diagonal loop resistance, ``(N, N)`` or its diagonal [Ohm].
    M_mat : numpy.ndarray
        Passive-passive inductance in code naming; the physics ``L``, ``(N, N)``
        [H].
    L_mat : numpy.ndarray
        Passive-to-source coupling in code naming; the physics ``M``,
        ``(N, n_src)`` [H].
    order : int, optional
        How many Krylov blocks to build [-].

    Returns
    -------
    numpy.ndarray
        R-orthonormal patterns as columns, ``(N, m)`` with
        ``m <= order * n_src`` [Ohm**-0.5].

    Convention
    ----------
    Argument naming follows :mod:`vaft.process.electromagnetics`, not the
    physics: ``M_mat`` is the passive-passive inductance (physics ``L``) and
    ``L_mat`` the passive-to-source coupling (physics ``M``).  ``R_mat`` is the
    diagonal loop resistance either way.

    Applicability
    -------------
    Machine-independent.  The comparison quoted below is a VEST observation, not
    a property of the method.

    Limitations
    -----------
    A zero-frequency expansion: on the packaged VEST wall ten resistive patterns
    reproduce the probe response of a real PF drive to about 1 %, where 150
    eigenmodes are needed for the same, but fast transients remain the
    eigenmodes' territory.

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    .. [vfit10] The reduced-order validation study (VEST-Tokamak/vfit#10, vaft
       #494), which asks which modes to keep; this module answers no part of that
       question by itself.
    """
    if int(order) < 1:
        raise WallModeError("order must be at least 1")
    r = _diagonal_resistance(R_mat, where="R_mat")
    L = np.asarray(M_mat, dtype=float)
    Q = orthonormalize_r(np.asarray(L_mat, dtype=float) / r[:, None], R_mat)
    blocks = [Q]
    for _ in range(int(order) - 1):
        X = (L @ blocks[-1]) / r[:, None]
        V = np.hstack(blocks)
        X = X - V @ (V.T @ (r[:, None] * X))          # R-orthogonal to everything so far
        X = X - V @ (V.T @ (r[:, None] * X))          # twice, for the usual reason
        Q = orthonormalize_r(X, R_mat, rtol=1e-8)
        if Q.shape[1] == 0:
            break
        blocks.append(Q)
    return canonical_sign(np.hstack(blocks))

def orthonormalize_r(X: np.ndarray, R_mat: np.ndarray, *, rtol: float = 1e-10) -> np.ndarray:
    """An R-orthonormal basis of ``span(X)``, canonical in sign.

    QR of ``R^{1/2} X`` with the dependent columns dropped, so a set of patterns
    that repeats a direction does not return a singular basis.

    Parameters
    ----------
    X : numpy.ndarray
        Patterns as columns, ``(N, m)`` [any].
    R_mat : numpy.ndarray
        Diagonal loop resistance, ``(N, N)`` or its diagonal [Ohm].
    rtol : float, optional
        A column whose triangular-factor diagonal falls below this times the
        largest is dependent and is dropped [-].

    Returns
    -------
    numpy.ndarray
        R-orthonormal, sign-canonical basis, ``(N, m')`` with ``m' <= m``
        [Ohm**-0.5].

    Defaults
    --------
    ``rtol = 1e-10`` is a numerical convenience: the rank cut of the QR in double
    precision.

    Convention
    ----------
    Modes are R-orthonormal (``v^T R v = 1``) and sign-canonical
    (:func:`canonical_sign`), so amplitudes carry units of sqrt(W) and ``a^T a``
    is the wall current's ohmic dissipation.  This is the same normalization the eigenbasis carries, which
    is what lets a pattern basis and an eigenbasis be combined in one reduction
    (:func:`combined_operators`).

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    """
    r = _diagonal_resistance(R_mat, where="R_mat")
    X = np.asarray(X, dtype=float)
    if X.ndim != 2 or X.shape[0] != r.size:
        raise WallModeError(f"patterns have shape {X.shape}; expected ({r.size}, m)")
    w = np.sqrt(r)[:, None]
    Q, T = np.linalg.qr(w * X)
    diag = np.abs(np.diag(T))
    independent = diag > rtol * max(diag.max() if diag.size else 0.0, 1e-300)
    return canonical_sign(Q[:, independent] / w)

def project(
    basis: WallModeBasis, I_w: np.ndarray, R_mat: np.ndarray, keep: Sequence[np.ndarray] | None = None
) -> np.ndarray:
    """Modal amplitudes ``a = V^T R I_w``.

    Exact for a current that lies inside the retained subspace; for a truncated
    basis it is the R-orthogonal projection, which is the least-dissipation one.

    Input semantics
    ---------------
    Element space: a wall current, one entry per element, measured or solved.

    Output semantics
    ----------------
    Mode space: amplitudes in the R inner product, so ``a^T a`` is the current's
    ohmic dissipation.

    Parameters
    ----------
    basis : WallModeBasis
        The eigenbasis [-].
    I_w : numpy.ndarray
        Wall current, ``(N,)`` or ``(n_times, N)`` [A].
    R_mat : numpy.ndarray
        Diagonal loop resistance, ``(N, N)`` or its diagonal [Ohm].
    keep : sequence of numpy.ndarray or None, optional
        Retained mode indices per segment; every mode when ``None`` [-].

    Returns
    -------
    numpy.ndarray
        Amplitudes, ``(M_tot,)`` or ``(n_times, M_tot)`` [W**0.5].

    Convention
    ----------
    Modes are R-orthonormal (``v^T R v = 1``) and sign-canonical
    (:func:`canonical_sign`), so amplitudes carry units of sqrt(W) and ``a^T a``
    is the wall current's ohmic dissipation.  Two other normalizations are in use elsewhere and
    convert from this one: Euclidean-normalized modes have ``a_E = a ||v||_2``,
    inductance-normalized ones ``a_L = a sqrt(tau)``.

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    """
    V = basis.V(keep)
    r = _diagonal_resistance(R_mat, where="R_mat")
    I_w = np.asarray(I_w, dtype=float)
    if I_w.ndim == 1:
        return V.T @ (r * I_w)
    return (I_w * r[None, :]) @ V

def reconstruct(basis: WallModeBasis, a: np.ndarray, keep: Sequence[np.ndarray] | None = None) -> np.ndarray:
    """Wall currents ``I_w = V a`` from modal amplitudes.

    The inverse direction of :func:`project`; exact only for the retained
    subspace, which is what :func:`reconstruction_error` measures.

    Input semantics
    ---------------
    Mode space: amplitudes of the retained modes.

    Output semantics
    ----------------
    Element space: a wall current on every element, whatever the retained order.

    Parameters
    ----------
    basis : WallModeBasis
        The eigenbasis [-].
    a : numpy.ndarray
        Amplitudes, ``(M_tot,)`` or ``(n_times, M_tot)`` [W**0.5].
    keep : sequence of numpy.ndarray or None, optional
        Retained mode indices per segment; every mode when ``None`` [-].

    Returns
    -------
    numpy.ndarray
        Wall current, ``(N,)`` or ``(n_times, N)`` [A].

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    """
    V = basis.V(keep)
    a = np.asarray(a, dtype=float)
    return V @ a if a.ndim == 1 else a @ V.T

def reconstruction_error(
    I_w: np.ndarray, I_rec: np.ndarray, R_mat: np.ndarray, basis: WallModeBasis | None = None
) -> dict[str, Any]:
    """How far a reconstruction is from the full current, in the norms that matter.

    Parameters
    ----------
    I_w : numpy.ndarray
        The full wall current, ``(N,)`` or ``(n_times, N)`` [A].
    I_rec : numpy.ndarray
        The reconstruction, same shape [A].
    R_mat : numpy.ndarray
        Diagonal loop resistance, ``(N, N)`` or its diagonal [Ohm].
    basis : WallModeBasis or None, optional
        When given, the same measures are reported per segment [-].

    Returns
    -------
    dict
        ``relative_l2``, ``relative_dissipation`` and ``max_abs``, plus one
        such entry per segment id when a basis was given [-].

        ``max_abs`` is the largest absolute element error, in amperes; the
        other two are relative.

    Convention
    ----------
    ``relative_dissipation`` weighs each element by its resistance -- the natural
    norm of this basis, where the error's ohmic power is compared with the
    current's -- and ``relative_l2`` is the plain Euclidean measure.  They differ
    by however unevenly the wall's resistance is distributed, so a reported
    error means nothing without saying which of the two it is.

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    """
    r = _diagonal_resistance(R_mat, where="R_mat")
    I_w = np.asarray(I_w, dtype=float).reshape(-1, r.size)
    I_rec = np.asarray(I_rec, dtype=float).reshape(-1, r.size)
    err = I_rec - I_w

    def _measures(e: np.ndarray, ref: np.ndarray, w: np.ndarray) -> dict[str, float]:
        l2 = float(np.linalg.norm(e) / max(np.linalg.norm(ref), 1e-300))
        diss = float(np.sqrt(np.sum(e**2 * w) / max(np.sum(ref**2 * w), 1e-300)))
        return {"relative_l2": l2, "relative_dissipation": diss, "max_abs": float(np.max(np.abs(e)))}

    out: dict[str, Any] = _measures(err, I_w, r[None, :])
    if basis is not None:
        out["segments"] = {
            seg.id: _measures(err[:, seg.index], I_w[:, seg.index], r[None, seg.index])
            for seg in basis.segments
        }
    return out

def reduce_response(
    G_full: np.ndarray, basis: WallModeBasis, keep: Sequence[np.ndarray] | None = None
) -> np.ndarray:
    """Reduce any linear response with the wall elements as its columns.

    Probe field, flux-loop flux, grid psi or B, an EFIT table -- one matmul, no
    per-class logic, so every response class is reduced by the same basis.

    Input semantics
    ---------------
    Element space: a response matrix mapping every wall element to observations.

    Output semantics
    ----------------
    Mode space: the same map from the retained modal amplitudes.

    Parameters
    ----------
    G_full : numpy.ndarray
        Response with wall elements as columns, ``(n_obs, N)`` [any/A].
    basis : WallModeBasis
        The eigenbasis [-].
    keep : sequence of numpy.ndarray or None, optional
        Retained mode indices per segment; every mode when ``None`` [-].

    Returns
    -------
    numpy.ndarray
        ``G_red = G_full V_seg``, ``(n_obs, M_tot)`` [any/W**0.5].

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    """
    G_full = np.asarray(G_full, dtype=float)
    if G_full.shape[-1] != basis.n_elements:
        raise WallModeError(
            f"response has {G_full.shape[-1]} wall columns, basis has {basis.n_elements} elements"
        )
    return G_full @ basis.V(keep)

def reduced_operators(
    basis: WallModeBasis,
    R_mat: np.ndarray,
    M_mat: np.ndarray,
    L_mat: np.ndarray | None = None,
    keep: Sequence[np.ndarray] | None = None,
) -> ReducedWall:
    """Project the full wall matrices onto the retained modes.

    The global half of the contract.  Because the *full* matrices are projected,
    the mutual inductance between segments survives in the off-diagonal blocks of
    ``L_r``; only the diagonal blocks are the local ``diag(tau_g)``.

    Input semantics
    ---------------
    Element space: the assembled wall matrices, one row and column per wall
    element, as :mod:`vaft.process.electromagnetics` builds them.

    Output semantics
    ----------------
    Mode space: operators of the retained modal coordinates, whose amplitudes
    carry units of sqrt(W).  A solve on them is the reduced circuit
    (:func:`solve_reduced_eddy`); to come back to element space, use
    :func:`reconstruct`.

    Parameters
    ----------
    basis : WallModeBasis
        The eigenbasis [-].
    R_mat : numpy.ndarray
        Diagonal loop resistance, ``(N, N)`` or its diagonal [Ohm].
    M_mat : numpy.ndarray
        Passive-passive inductance in code naming, ``(N, N)`` [H].
    L_mat : numpy.ndarray or None, optional
        Passive-to-source coupling in code naming, ``(N, n_src)`` [H].
    keep : sequence of numpy.ndarray or None, optional
        Retained mode indices per segment; every mode when ``None`` [-].

    Returns
    -------
    ReducedWall
        ``L_r`` ``(M_tot, M_tot)`` [H], ``R_r`` ``(M_tot, M_tot)`` [Ohm],
        ``M_r`` ``(M_tot, n_src)`` [H] when a source coupling was given, and the
        ``(segment_id, k)`` label of every coefficient [-].

    Convention
    ----------
    Argument naming follows :mod:`vaft.process.electromagnetics`, not the
    physics: ``M_mat`` is the passive-passive inductance (physics ``L``) and
    ``L_mat`` the passive-to-source coupling (physics ``M``).  ``R_mat`` is the
    diagonal loop resistance either way.

    ``R_r`` is computed rather than assumed to be the identity, so a basis built
    on other matrices than the ones passed here shows up as a non-identity
    instead of silently changing the dynamics.

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    """
    V = basis.V(keep)
    r = _diagonal_resistance(R_mat, where="R_mat")
    L = np.asarray(M_mat, dtype=float)
    L_r = V.T @ L @ V
    R_r = V.T @ (r[:, None] * V)
    M_r = None if L_mat is None else V.T @ np.asarray(L_mat, dtype=float)
    return ReducedWall(L_r=L_r, R_r=R_r, M_r=M_r, labels=basis.labels(keep), keep=basis._keep(keep))

def select_all(basis: WallModeBasis) -> tuple[np.ndarray, ...]:
    """Keep every mode of every segment: the full-rank selection.

    Parameters
    ----------
    basis : WallModeBasis
        The eigenbasis [-].

    Returns
    -------
    tuple of numpy.ndarray
        One index array per segment, holding every mode index [-].

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    """
    return basis._keep(None)

def select_by_score(
    basis: WallModeBasis, score: np.ndarray, M: int, keep: Sequence[np.ndarray] | None = None
) -> tuple[np.ndarray, ...]:
    """Keep the highest-scoring modes across segments.

    The global counterpart of :func:`select_slowest`, for any ranking -- usually
    one of :func:`mode_scores`.

    Parameters
    ----------
    basis : WallModeBasis
        The eigenbasis [-].
    score : numpy.ndarray
        One value per coefficient of ``basis.labels(keep)`` [-].
    M : int
        How many modes to keep in total [-].
    keep : sequence of numpy.ndarray or None, optional
        Candidate modes per segment; every mode when ``None`` [-].

    Returns
    -------
    tuple of numpy.ndarray
        One index array per segment [-].

    Convention
    ----------
    Ties resolve toward the earlier label, so the selection is deterministic for
    a given score vector.

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    .. [vfit10] The reduced-order validation study (VEST-Tokamak/vfit#10, vaft
       #494), which asks which modes to keep; this module answers no part of that
       question by itself.
    """
    labels = basis.labels(keep)
    score = np.asarray(score, dtype=float).reshape(-1)
    if score.size != len(labels):
        raise WallModeError(f"score has {score.size} entries for {len(labels)} candidate modes")
    order = np.argsort(-score, kind="stable")[: max(int(M), 0)]
    chosen = set(labels[i] for i in order)
    return tuple(
        np.array(sorted(k for (seg_id, k) in chosen if seg_id == seg.id), dtype=np.int64)
        for seg in basis.segments
    )

def select_slowest(basis: WallModeBasis, M: int | Sequence[int]) -> tuple[np.ndarray, ...]:
    """Keep the slowest modes, globally or per segment.

    The selection the spectrum alone suggests; whether it is the right one for a
    given drive is what :func:`mode_scores` exists to question.

    Parameters
    ----------
    basis : WallModeBasis
        The eigenbasis [-].
    M : int or sequence of int
        A count taken globally across segments by descending decay time, or one
        count per segment [-].

    Returns
    -------
    tuple of numpy.ndarray
        One index array per segment [-].

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    .. [vfit10] The reduced-order validation study (VEST-Tokamak/vfit#10, vaft
       #494), which asks which modes to keep; this module answers no part of that
       question by itself.
    """
    if isinstance(M, (int, np.integer)):
        tau = basis.tau()
        labels = basis.labels()
        chosen = set(labels[i] for i in np.argsort(-tau, kind="stable")[: int(M)])
        return tuple(
            np.array([k for k in range(seg.size) if (seg.id, k) in chosen], dtype=np.int64)
            for seg in basis.segments
        )
    counts = list(M)
    if len(counts) != len(basis.segments):
        raise WallModeError(f"per-segment M has {len(counts)} entries for {len(basis.segments)} segments")
    return tuple(np.arange(min(int(m), seg.size), dtype=np.int64) for seg, m in zip(basis.segments, counts))

def select_tau_range(basis: WallModeBasis, tau_min: float, tau_max: float = np.inf) -> tuple[np.ndarray, ...]:
    """Keep the modes whose decay time lies in a range, per segment.

    Parameters
    ----------
    basis : WallModeBasis
        The eigenbasis [-].
    tau_min : float
        Lower bound, inclusive [s].
    tau_max : float, optional
        Upper bound, inclusive [s].

    Returns
    -------
    tuple of numpy.ndarray
        One index array per segment [-].

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    """
    return tuple(
        np.flatnonzero((seg.tau >= tau_min) & (seg.tau <= tau_max)).astype(np.int64)
        for seg in basis.segments
    )

def solve_reduced_eddy(
    reduced: ReducedWall,
    drive: np.ndarray,
    time: np.ndarray,
    *,
    V: np.ndarray | None = None,
    dt_sub: float = 5.0e-5,
    method: str = "auto",
) -> tuple[np.ndarray, np.ndarray | None]:
    """Integrate the reduced circuit ``L_r da/dt + R_r a = -M_r dI_src/dt``.

    The same integrator as the full wall
    (:func:`vaft.process.electromagnetics.solve_eddy_currents`), fed the projected
    operators, so the only difference between the two solutions is the retained
    subspace: with every mode kept they agree to rounding.

    Input semantics
    ---------------
    Mode space: reduced operators and a source drive in element-independent
    coordinates.

    Output semantics
    ----------------
    Mode space amplitudes, and -- when the retained ``V`` is given -- the
    reconstructed element-space wall current alongside them.

    Parameters
    ----------
    reduced : ReducedWall
        The projected operators; ``M_r`` must be present [-].
    drive : numpy.ndarray
        Source currents, ``(n_times, n_src)`` on ``time`` [A].
    time : numpy.ndarray
        Time grid of the drive [s].
    V : numpy.ndarray or None, optional
        The retained modes, ``(N, M_tot)``; when given, the wall current is
        reconstructed too [Ohm**-0.5].
    dt_sub : float, optional
        Sub-step of the integrator [s].
    method : str, optional
        Integration method, forwarded to the full-wall solver [-].

    Returns
    -------
    a : numpy.ndarray
        Modal amplitudes, ``(n_times, M_tot)`` [W**0.5].
    I_w : numpy.ndarray or None
        ``V a``, ``(n_times, N)``, or ``None`` when no ``V`` was given [A].

    Defaults
    --------
    ``dt_sub = 50 us`` is a validated-workflow default: the sub-step the full-wall
    solver is used with in the VEST eddy-current workflow.  ``method="auto"``
    leaves the choice to that solver.

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    .. [vfit10] The reduced-order validation study (VEST-Tokamak/vfit#10, vaft
       #494), which asks which modes to keep; this module answers no part of that
       question by itself.
    """
    from vaft.process.electromagnetics import solve_eddy_currents

    if reduced.M_r is None:
        raise WallModeError("solve_reduced_eddy needs the projected source coupling M_r")
    a = solve_eddy_currents(
        reduced.R_r, reduced.M_r, reduced.L_r,
        np.asarray(drive, dtype=float), np.asarray(time, dtype=float),
        dt_sub=dt_sub, method=method,
    )
    I_w = None if V is None else a @ np.asarray(V, dtype=float).T
    return a, I_w

def subspace_angles_r(V_a: np.ndarray, V_b: np.ndarray, R_mat: np.ndarray) -> np.ndarray:
    """Principal angles between two mode subspaces in the R inner product.

    Two bases of a near-degenerate eigenspace differ by a rotation, so the right
    question is whether they span the same space, not whether their columns
    match.

    Parameters
    ----------
    V_a : numpy.ndarray
        First set of modes as columns, ``(N, m_a)`` [Ohm**-0.5].
    V_b : numpy.ndarray
        Second set, ``(N, m_b)`` [Ohm**-0.5].
    R_mat : numpy.ndarray
        Diagonal loop resistance, ``(N, N)`` or its diagonal [Ohm].

    Returns
    -------
    numpy.ndarray
        Principal angles, ascending, ``(min(m_a, m_b),)`` [rad].

    Convention
    ----------
    The angles are taken in the R metric, not the Euclidean one: both sets are
    mapped by ``R^{1/2}`` first, so the ordinary angles of the mapped columns are
    the R-metric angles of the originals.  A zero angle means a shared direction;
    ``pi/2`` means an R-orthogonal one.

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

    Provenance
    ----------
    .. [473] Issue #473, which built the segment-wise eigenbasis, and #500, which
       added the selection, scoring and moment tools.
    .. [vfit8] The reduced vessel-wall contract, VEST-Tokamak/vfit#8: local
       eigenbasis, global electromagnetic dynamics.
    """
    from scipy.linalg import subspace_angles

    r = _diagonal_resistance(R_mat, where="R_mat")
    w = np.sqrt(r)[:, None]
    return subspace_angles(w * np.asarray(V_a, dtype=float), w * np.asarray(V_b, dtype=float))