def check_stability(algorithm: Algorithm) -> StabilityResult:
    """Check the stability of an algorithm's matching.

    Args:
        algorithm: An Algorithm instance that has been executed.

    Returns:
        StabilityResult with is_stable, is_individually_rational, and blocking_pairs.
    """
    ir = is_individually_rational(algorithm.proposers, algorithm.responders)
    bp = find_blocking_pairs(algorithm.proposers, algorithm.responders)
    return StabilityResult(
        is_stable=ir and len(bp) == 0,
        is_individually_rational=ir,
        blocking_pairs=bp,
    )

def create_matching(
    proposer_preferences: dict[str, list[str]],
    responder_preferences: dict[str, list[str]],
) -> MatchingResult:
    """Create a matching from preference dictionaries.

    This is the simplest way to use the library. Provide preference rankings
    for each side and get back a MatchingResult.

    Persons not listed in a preference list are considered unacceptable.
    Each person is automatically added to their own preference list at the end
    (making everyone they listed acceptable).

    Args:
        proposer_preferences: Mapping of proposer names to ordered list of responder names.
        responder_preferences: Mapping of responder names to ordered list of proposer names.

    Returns:
        MatchingResult with the matching outcome.

    Example:
        >>> result = create_matching(
        ...     proposer_preferences={
        ...         "alice": ["bob", "charlie"],
        ...         "dave": ["charlie", "bob"],
        ...     },
        ...     responder_preferences={
        ...         "bob": ["alice", "dave"],
        ...         "charlie": ["dave", "alice"],
        ...     },
        ... )
        >>> result.matches
        {'alice': 'bob', 'dave': 'charlie'}
    """
    algorithm = _build_algorithm(proposer_preferences, responder_preferences)
    return algorithm.execute()

def find_blocking_pairs(proposers: list[Proposer], responders: list[Responder]) -> list[tuple[str, str]]:  # noqa: ARG001
    """Find all blocking pairs in a matching.

    A blocking pair (p, r) exists when proposer p and responder r both prefer
    each other over their current matches.

    Args:
        proposers: List of proposers in the matching.
        responders: List of responders in the matching.

    Returns:
        List of (proposer_name, responder_name) blocking pairs.
    """
    blocking: list[tuple[str, str]] = []
    for proposer in proposers:
        if not (bool(proposer.preferences) and proposer.is_matched):
            continue

        better_than_match = proposer.preferences[: proposer.preferences.index(proposer.match)]
        for responder in better_than_match:
            if not isinstance(responder, Responder):
                continue

            match responder.is_matched:
                case False:
                    blocking.append((proposer.name, responder.name))
                case True if (
                    bool(responder.preferences)
                    and all(item in responder.preferences for item in [proposer, responder.match])
                    and responder.preferences.index(proposer) < responder.preferences.index(responder.match)
                ):
                    blocking.append((proposer.name, responder.name))

    return blocking

def is_individually_rational(proposers: list[Proposer], responders: list[Responder]) -> bool:
    """Check if the matching is individually rational.

    A matching is individually rational if every matched person finds their
    match acceptable (ranked at or above themselves).

    Args:
        proposers: List of proposers in the matching.
        responders: List of responders in the matching.

    Returns:
        True if individually rational, False otherwise.
    """
    persons = proposers + responders
    return all(person.match is None or person.is_acceptable(person.match) for person in persons)

def fifo_selector(free: Sequence[int]) -> int:
    """Return the first index — pops in arrival order (FIFO)."""
    del free
    return 0

def gale_shapley(proposer_rank: NDArray[np.integer], responder_rank: NDArray[np.integer]) -> NDArray[np.int16]:
    """Run proposer-optimal deferred acceptance.

    Args:
        proposer_rank: ``(n, n)`` array where ``proposer_rank[i, j]`` is
            the 1-indexed position of responder ``j`` in proposer ``i``'s
            preference list. Each row must be a permutation of ``1..n``.
        responder_rank: ``(n, n)`` array with the symmetric convention.

    Returns:
        ``match`` of shape ``(n,)``, dtype ``int16``, where ``match[i]`` is
        the responder paired with proposer ``i``.

    Raises:
        ValueError: if the two arrays don't have the same square shape, or
            any row is not a permutation of ``1..n``.
    """
    return gale_shapley_traced(proposer_rank, responder_rank).match

def gale_shapley_traced(
    proposer_rank: NDArray[np.integer],
    responder_rank: NDArray[np.integer],
    *,
    selector: Selector = lifo_selector,
) -> GSStats:
    """Run proposer-optimal deferred acceptance and return per-proposer statistics.

    Mechanically identical to :func:`gale_shapley` (sequential McVitie-Wilson
    on the free-proposer pool), but exposes the proposal count and a
    pluggable proposer-selection rule. The final matching is invariant under
    selector choice (Knuth); the proposal count is not.

    Args:
        proposer_rank: ``(n, n)`` array where ``proposer_rank[i, j]`` is
            the 1-indexed position of responder ``j`` in proposer ``i``'s
            preference list. Each row must be a permutation of ``1..n``.
        responder_rank: ``(n, n)`` array with the symmetric convention.
        selector: picks which free proposer acts next; defaults to
            :func:`lifo_selector`. See :data:`Selector`.

    Returns:
        :class:`GSStats` with the matching, total proposal count, and
        per-proposer counts.

    Raises:
        ValueError: if the two arrays don't have the same square shape, or
            any row is not a permutation of ``1..n``.
        IndexError: if ``selector`` returns an index outside the free list.
    """
    _validate_rank_matrices(proposer_rank, responder_rank)
    n = proposer_rank.shape[0]
    next_proposal = np.zeros(n, dtype=np.int16)
    responder_match = np.full(n, -1, dtype=np.int16)
    proposer_pref = np.argsort(proposer_rank, axis=1).astype(np.int16)
    free: list[int] = list(range(n))
    while free:
        p = free.pop(selector(free))
        r = int(proposer_pref[p, next_proposal[p]])
        next_proposal[p] += 1
        current = int(responder_match[r])
        if current == -1:
            responder_match[r] = p
        elif responder_rank[r, p] < responder_rank[r, current]:
            responder_match[r] = p
            free.append(current)
        else:
            free.append(p)
    return GSStats(
        match=np.argsort(responder_match).astype(np.int16),
        proposals=int(next_proposal.sum()),
        proposals_per_proposer=next_proposal.copy(),
    )

def lifo_selector(free: Sequence[int]) -> int:
    """Return the last index — equivalent to ``list.pop()`` (LIFO)."""
    del free
    return -1

def men_optimal_gs(men_rank: NDArray[np.integer], women_rank: NDArray[np.integer]) -> NDArray[np.int16]:
    """Return the men-optimal stable matching (conventionally men propose, women dispose)."""
    return gale_shapley(men_rank, women_rank)

def men_optimal_traced(
    men_rank: NDArray[np.integer],
    women_rank: NDArray[np.integer],
    *,
    selector: Selector = lifo_selector,
) -> GSStats:
    """Men-optimal stable matching with per-man proposal stats.

    Equivalent to ``gale_shapley_traced(men_rank, women_rank, selector=...)``;
    provided for symmetry with :func:`men_optimal_gs`.

    Args:
        men_rank: ``(n, n)`` rank matrix for men (proposers).
        women_rank: ``(n, n)`` rank matrix for women (responders).
        selector: see :data:`Selector`.

    Returns:
        :class:`GSStats` where ``match`` and ``proposals_per_proposer`` are
        both men-indexed.
    """
    return gale_shapley_traced(men_rank, women_rank, selector=selector)

def random_selector(rng: np.random.Generator) -> Selector:
    """Build a uniform-random selector backed by a numpy ``Generator``.

    Roth-Vande Vate-style: pick a free proposer uniformly at random each step.
    Final matching is invariant to selector choice (Knuth's order independence)
    but the proposal count is not — random order is a useful baseline against
    LIFO/FIFO for amortized analysis.

    Args:
        rng: a numpy ``Generator`` (use ``np.random.default_rng(seed)`` for
            determinism).

    Returns:
        A :data:`Selector` closure that draws ``rng.integers(len(free))`` per call.
    """

    def _selector(free: Sequence[int]) -> int:
        return int(rng.integers(len(free)))

    return _selector

def women_optimal_gs(men_rank: NDArray[np.integer], women_rank: NDArray[np.integer]) -> NDArray[np.int16]:
    """Return the women-optimal stable matching, still in men-indexed form ``match[m] = w``."""
    return np.argsort(gale_shapley(women_rank, men_rank)).astype(np.int16)

def women_optimal_traced(
    men_rank: NDArray[np.integer],
    women_rank: NDArray[np.integer],
    *,
    selector: Selector = lifo_selector,
) -> GSStats:
    """Women-optimal stable matching with per-woman proposal stats.

    Runs :func:`gale_shapley_traced` with women as the proposing side. The
    returned ``match`` is re-indexed to men-indexed form for consistency with
    :func:`women_optimal_gs`. ``proposals_per_proposer`` remains
    **women-indexed**, since women are the proposers in this run.

    Args:
        men_rank: ``(n, n)`` rank matrix for men.
        women_rank: ``(n, n)`` rank matrix for women.
        selector: see :data:`Selector`.

    Returns:
        :class:`GSStats` with ``match[m] = w`` (men-indexed) and
        ``proposals_per_proposer[w] = number of proposals woman w made``
        (women-indexed). ``proposals`` is the unambiguous total.
    """
    stats = gale_shapley_traced(women_rank, men_rank, selector=selector)
    return GSStats(
        match=np.argsort(stats.match).astype(np.int16),
        proposals=stats.proposals,
        proposals_per_proposer=stats.proposals_per_proposer,
    )

def find_blocking_pairs(
    proposer_rank: NDArray[np.integer],
    responder_rank: NDArray[np.integer],
    match: NDArray[np.integer],
) -> list[tuple[int, int]]:
    """Return every ``(proposer, responder)`` pair that blocks ``match``.

    A blocking pair ``(p, r)`` is one where ``p`` strictly prefers ``r`` to
    his current partner and ``r`` strictly prefers ``p`` to her current
    partner.

    Args:
        proposer_rank: ``(n, n)`` rank matrix for the proposing side.
        responder_rank: ``(n, n)`` rank matrix for the responding side.
        match: ``(n,)`` men-indexed matching.

    Returns:
        List of 0-indexed ``(proposer, responder)`` blocking pairs.
    """
    n = proposer_rank.shape[0]
    partner_of_responder = np.argsort(match).astype(np.int16)
    blocking: list[tuple[int, int]] = []
    for p in range(n):
        current_r = int(match[p])
        for r in range(n):
            if r == current_r:
                continue
            if (
                proposer_rank[p, r] < proposer_rank[p, current_r]
                and responder_rank[r, p] < responder_rank[r, partner_of_responder[r]]
            ):
                blocking.append((p, r))
    return blocking

def is_stable(
    proposer_rank: NDArray[np.integer],
    responder_rank: NDArray[np.integer],
    match: NDArray[np.integer],
) -> bool:
    """Return ``True`` iff ``match`` is a perfect matching with no blocking pair."""
    n = proposer_rank.shape[0]
    if match.shape[0] != n:
        return False
    if np.unique(match).shape[0] != n:
        return False
    return len(find_blocking_pairs(proposer_rank, responder_rank, match)) == 0

def is_stable_batch(
    proposer_rank: NDArray[np.integer],
    responder_rank: NDArray[np.integer],
    matchings: NDArray[np.integer],
) -> NDArray[np.bool_]:
    """Vectorized stability check across a batch of matchings.

    Memory use is O(K * n**2) where K is the batch size. For K = n! (full
    permutation enumeration) the tensor grows as O(n**2 * n!), ~12 MB at
    n=8 and ~360 MB at n=10 — caller's responsibility to keep K bounded.

    Args:
        proposer_rank: ``(n, n)`` rank matrix.
        responder_rank: ``(n, n)`` rank matrix.
        matchings: ``(K, n)`` array of candidate men-indexed matchings.

    Returns:
        ``(K,)`` bool array.
    """
    arange_n = np.arange(proposer_rank.shape[0], dtype=np.int16)
    partner_rank_proposer = proposer_rank[arange_n[np.newaxis, :], matchings]
    partner_of_responder = np.argsort(matchings, axis=1)
    partner_rank_responder = responder_rank[arange_n[np.newaxis, :], partner_of_responder]

    block_p = proposer_rank[np.newaxis, :, :] < partner_rank_proposer[:, :, np.newaxis]
    block_r = responder_rank.T[np.newaxis, :, :] < partner_rank_responder[:, np.newaxis, :]
    has_block = np.any(block_p & block_r, axis=(1, 2))
    return ~has_block

def apply_rotation(
    matching: NDArray[np.integer],
    rotation: NDArray[np.integer],
) -> NDArray[np.int16]:
    """Apply a rotation to a stable matching, producing the next matching.

    Args:
        matching: ``(n,)`` men-indexed stable matching.
        rotation: ``(r, 2)`` array from :func:`exposed_rotations`.

    Returns:
        A fresh ``(n,)`` int16 matching. Input arrays are not mutated.
    """
    new_matching = matching.astype(np.int16, copy=True)
    for i in range(int(rotation.shape[0])):
        new_matching[int(rotation[i, 0])] = int(rotation[i, 1])
    return new_matching

def enumerate_stable_matchings(
    men_rank: NDArray[np.integer],
    women_rank: NDArray[np.integer],
    method: Literal["rotation", "brute"] = "rotation",
    max_n: int = _DEFAULT_MAX_N_BRUTE,
    batch_size: int = _DEFAULT_BATCH_SIZE,
) -> NDArray[np.int16]:
    """Return every stable matching of the instance.

    Args:
        men_rank: ``(n, n)`` men's rank matrix.
        women_rank: ``(n, n)`` women's rank matrix.
        method: ``"rotation"`` (default) uses the Gusfield-Irving rotation
            algorithm and scales to large n; ``"brute"`` uses vectorized
            permutation enumeration and is capped at ``max_n``.
        max_n: only consulted when ``method="brute"``.
        batch_size: only consulted when ``method="brute"``.

    Returns:
        ``(|L|, n)`` int16 array of men-indexed stable matchings.

    Raises:
        ValueError: if the rank matrices are not same-shape square permutation
            matrices, or ``method`` is unknown, or brute-force is requested
            with ``n > max_n``.
    """
    _validate_rank_matrices(men_rank, women_rank)
    if method == "rotation":
        return _enumerate_via_rotations(men_rank, women_rank)
    if method == "brute":
        return _enumerate_via_brute_force(men_rank, women_rank, max_n=max_n, batch_size=batch_size)
    raise ValueError(f"Unknown method {method!r}. Expected 'rotation' or 'brute'.")

def exposed_rotations(
    men_rank: NDArray[np.integer],
    women_rank: NDArray[np.integer],
    matching: NDArray[np.integer],
) -> list[NDArray[np.int16]]:
    """Return the list of rotations currently exposed at ``matching``.

    Each rotation is an ``(r, 2)`` int16 array; row ``i`` is ``[m_i, w_{i+1}]``,
    where ``w_{i+1}`` is the partner ``m_i`` acquires if the rotation is eliminated.
    Applying the rotation amounts to assigning ``new_matching[m_i] = w_{i+1}``
    cyclically.

    At the men-optimal matching every rotation exposed here corresponds to a
    distinct one-step move down the lattice; at the women-optimal matching no
    rotations are exposed (you have reached the infimum).

    Args:
        men_rank: ``(n, n)`` men's rank matrix.
        women_rank: ``(n, n)`` women's rank matrix.
        matching: ``(n,)`` men-indexed stable matching.

    Returns:
        List of rotations; may be empty if no rotation is exposed.

    Raises:
        ValueError: if the rank matrices are not same-shape square permutation
            matrices.
    """
    _validate_rank_matrices(men_rank, women_rank)
    next_target, next_man = _compute_rotation_pointers(matching, men_rank, women_rank)
    cycles = _find_cycles(next_man)
    rotations: list[NDArray[np.int16]] = []
    for cycle in cycles:
        rot = np.empty((len(cycle), 2), dtype=np.int16)
        for i, m in enumerate(cycle):
            rot[i, 0] = m
            rot[i, 1] = next_target[m]
        rotations.append(rot)
    return rotations