ordito.edges

def faces_to_edges(
    faces: wp.array[wp.int32],
    sorted: bool = False,  # noqa: A002
) -> odt.Array2dInt32:
    """
    Directed triangle edges from a flat ``(i0, i1, i2)`` index buffer.

    For each face emits three directed edges ``(i0, i1)``, ``(i1, i2)``, ``(i2, i0)`` in
    row-major order, matching [`trimesh.geometry.faces_to_edges`][].

    Parameters
    ----------
    faces
        ``(3 * n_faces,)`` buffer of consecutive vertex-index triples.
    sorted
        If ``True``, each output row has its smaller vertex index first (undirected edges).

    Returns
    -------
    odt.Array2dInt32
        ``(3 * n_faces, 2)`` directed edges. Empty ``(0, 2)`` array when ``n_faces == 0``.

    See Also
    --------
    [`trimesh.geometry.faces_to_edges`][]

    Examples
    --------
    Three directed edges per face, in corner order; ``sorted=True`` puts the smaller index first:

    ```python
    import warp as wp

    import ordito as od

    faces = wp.array([0, 1, 2, 2, 1, 3], dtype=wp.int32)
    directed = od.edges.faces_to_edges(faces)
    print(directed.numpy().tolist())  # [[0, 1], [1, 2], [2, 0], [2, 1], [1, 3], [3, 2]]
    undirected = od.edges.faces_to_edges(faces, sorted=True)
    print(undirected.numpy().tolist())  # [[0, 1], [1, 2], [0, 2], [1, 2], [1, 3], [2, 3]]
    ```
    """
    n_faces = faces.size // 3
    edges = odt.empty_2d((n_faces * 3, 2), wp.int32, device=faces.device)
    _launch.launch(
        kernel_edges.faces_to_edges,
        dim=n_faces,
        inputs=[faces, wp.bool(sorted), edges],
        device=faces.device,
    )
    return odt.as_array2d(edges, wp.int32)

def edges_face(faces: wp.array[wp.int32]) -> wp.array[wp.int32]:
    """
    Face index for each directed edge produced by [`faces_to_edges`][ordito.edges.faces_to_edges].

    Returns an array of length ``n_faces * 3`` where entry ``3*f + k`` equals ``f``,
    matching [`trimesh.Trimesh.edges_face`][].

    Parameters
    ----------
    faces
        ``(3 * n_faces,)`` face index buffer.

    Returns
    -------
    wp.array[wp.int32]
        ``(3 * n_faces,)`` face indices on ``faces.device``.

    See Also
    --------
    [`trimesh.Trimesh.edges_face`][]

    Examples
    --------
    ```python
    import warp as wp

    import ordito as od

    faces = wp.array([0, 1, 2, 2, 1, 3], dtype=wp.int32)
    print(od.edges.edges_face(faces).numpy().tolist())  # [0, 0, 0, 1, 1, 1]
    ```
    """
    n_faces = faces.size // 3
    return arange_repeat(n_faces * 3, 3, faces.device)

def edges_unique(
    faces: wp.array[wp.int32],
    edges_sorted: odt.Array2dInt32 | None = None,
    n_vertices: int | None = None,
    *,
    validate: bool = True,
) -> tuple[odt.Array2dInt32, wp.array[wp.int32]]:
    """
    Return unique undirected edges and their inverse mapping into the sorted edge list.

    Equivalent to [`trimesh.Trimesh.edges_unique`][] and
    [`trimesh.Trimesh.edges_unique_inverse`][] computed together.

    Parameters
    ----------
    faces
        ``(3 * n_faces,)`` face index buffer.
    edges_sorted
        ``(3 * n_faces, 2)`` precomputed sorted edges (each row min-first), or ``None`` to build
        them from ``faces``.
    n_vertices
        Total number of vertices (used as the hash base, and as the radix the unique rows are
        unpacked with). When ``None`` and ``validate`` is ``True`` it is inferred from the edge
        indices with a device-host sync -- the same reduction that runs the range check.
        When ``None`` and ``validate`` is ``False`` the keys pack against
        [`constants.INDEX_RADIX_PAIR`][ordito.constants.INDEX_RADIX_PAIR] instead, which bounds
        every ``int32`` index with no reduction at all and leaves the row order unchanged.
    validate
        Whether to range-check the edge indices. With ``n_vertices`` given (and no
        ``edges_sorted``) the check is folded into the key pass and read back with the edge count,
        so it costs no extra synchronization; otherwise it is a ``ordito.reduce.minmax`` whose host
        readback serialises the device pipeline. Pass ``False``
        only where both bounds are structurally guaranteed -- a face buffer this package produced
        itself, or one an entry point has already validated. A ``False`` that is wrong does not
        raise: indices at or above ``n_vertices`` collide in the packing and silently group two
        different edges as one, and a negative index is read as a huge unsigned digit.

    Returns
    -------
    unique_edges : odt.Array2dInt32
        ``(m, 2)`` unique undirected vertex pairs, ``m <= 3 * n_faces``.
    inverse : wp.array[wp.int32]
        ``(3 * n_faces,)`` inverse indices, ``unique_edges[inverse[i]] == edges_sorted[i]``.

    Raises
    ------
    TypeError
        If ``edges_sorted`` is given and is not a rank-2 ``wp.int32`` array.
    ValueError
        If ``edges_sorted`` is given and does not have exactly two columns, if ``n_vertices`` is
        not positive, or if ``validate`` is ``True`` and an edge index is negative or reaches
        ``n_vertices``.
    RuntimeError
        If ``faces`` and ``edges_sorted`` are not all on one device.

    See Also
    --------
    [`trimesh.Trimesh.edges_unique`][]
    [`trimesh.Trimesh.edges_unique_inverse`][]

    Examples
    --------
    A cube has 12 triangles, 36 corner edges and 18 unique edges (12 sides and 6 face diagonals);
    each unique edge of a closed mesh is reached from exactly two corners:

    ```python
    import numpy as np

    import ordito as od

    vertices, faces = od.creation.box()
    unique_edges, inverse = od.edges.edges_unique(faces)
    print(unique_edges.shape)  # (18, 2)
    print(inverse.shape)  # (36,)
    print(np.bincount(inverse.numpy()).tolist() == [2] * 18)  # True
    ```
    """
    require_same_device(faces=faces, edges_sorted=edges_sorted)
    if edges_sorted is not None:
        odt.ensure_edge_pairs(edges_sorted, "edges_sorted")
    n_faces = faces.size // 3
    device = faces.device

    if n_faces == 0:
        empty_edges = odt.empty_2d((0, 2), wp.int32, device=device)
        empty_inv = _launch.empty(0, dtype=wp.int32, device=device)
        return empty_edges, empty_inv

    if edges_sorted is None:
        # The keys are built straight off ``faces``: ``face_edge_keys`` packs each corner's sorted
        # pair exactly as ``hash_indices_rows`` packs the corresponding ``faces_to_edges`` row, so
        # the ``(3 * n_faces, 2)`` table is never written or read back, and any range check reduces
        # the ``3 * n_faces`` face buffer -- the same values, half the entries.
        radix, check = _face_edge_radix(faces, n_faces, n_vertices, validate)
        return _unique_edges_from_faces(faces, n_faces, radix, check)

    if n_vertices is None:
        if validate:
            # Inferring the bound from this very buffer already reduces it, and the packing's own
            # check would reduce it again to re-test a bound derived from it -- only the negative
            # half could ever fire. One reduction answers both.
            n_vertices = od.array.index_bound(edges_sorted, require_non_negative=True)
            validate = False
        else:
            # Nothing is being checked, so the only thing the count would be used for is the
            # radix -- and the pair radix bounds every ``int32`` without reducing the rows to find
            # out. Row order is unchanged; the emit below unpacks against whichever
            # radix packed them.
            n_vertices = INDEX_RADIX_PAIR
    elif n_vertices <= 0:
        raise ValueError(f"n_vertices must be positive, got {n_vertices}")
    elif validate:
        # The packing's range check: the reduction's readback, then the same two raises as the
        # face path's key-pass check.
        low, high = od.reduce.minmax(edges_sorted)
        if low < 0:
            raise ValueError(f"edge indices must be non-negative, got a minimum of {low}")
        if high >= n_vertices:
            raise ValueError(
                f"edge indices must be less than n_vertices {n_vertices}, got a maximum of {high}"
            )
    return _unique_edges_from_sorted_rows(edges_sorted, n_vertices)

def unique_edge_count(
    faces: wp.array[wp.int32], n_vertices: int | None = None, *, validate: bool = True
) -> int:
    """
    Count the unique undirected edges, without building the edge table.

    Equals ``edges_unique(faces, n_vertices=n_vertices)[0].shape[0]``: the same sort and run scan,
    stopped before the rows and the inverse map are written.

    Parameters
    ----------
    faces
        ``(3 * n_faces,)`` flat triangle index buffer.
    n_vertices
        Vertex count bounding the indices, as in [`edges_unique`][ordito.edges.edges_unique];
        inferred from ``faces`` when ``None``.
    validate
        As in [`edges_unique`][ordito.edges.edges_unique]: check the index range. With
        ``n_vertices`` given the check rides in the count's own readback; inferring the bound
        costs one more.

    Returns
    -------
    int
        Unique undirected edge count. ``0`` for an empty mesh.

    Raises
    ------
    ValueError
        If ``n_vertices`` is not positive, or ``validate`` is ``True`` and a face index is negative
        or reaches ``n_vertices``.

    See Also
    --------
    [`edges_unique`][ordito.edges.edges_unique]
        The table this counts.

    Examples
    --------
    ```python
    import ordito as od

    vertices, faces = od.creation.box()
    print(od.edges.unique_edge_count(faces))  # 18
    ```
    """
    n_faces = faces.size // 3
    if n_faces == 0:
        return 0
    radix, check = _face_edge_radix(faces, n_faces, n_vertices, validate)
    return _sorted_edge_run_ranks(faces, n_faces, radix, check)[3]

def edges_unique_inverse(
    faces: wp.array[wp.int32],
    edges_sorted: odt.Array2dInt32 | None = None,
    n_vertices: int | None = None,
) -> wp.array[wp.int32]:
    """
    Inverse mapping from sorted edges into [`edges_unique`][ordito.edges.edges_unique].

    Maps [`faces_to_edges`][ordito.edges.faces_to_edges] (sorted) into `edges_unique`:
    ``edges_unique(faces)[edges_unique_inverse(faces)] == faces_to_edges(faces, sorted=True)``

    Parameters
    ----------
    faces
        ``(3 * n_faces,)`` face index buffer.
    edges_sorted
        ``(3 * n_faces, 2)`` precomputed sorted edges, or ``None`` to build them from ``faces``.
    n_vertices
        Total vertex count. When ``None``, inferred from ``edges_sorted`` with a device-host sync.

    Returns
    -------
    wp.array[wp.int32]
        ``(3 * n_faces,)`` inverse indices.

    Raises
    ------
    TypeError
        If ``edges_sorted`` is given and is not a rank-2 ``wp.int32`` array.
    ValueError
        If ``edges_sorted`` is given and does not have exactly two columns.
    RuntimeError
        If ``faces`` and ``edges_sorted`` are not all on one device.

    See Also
    --------
    [`trimesh.Trimesh.edges_unique_inverse`][]

    Examples
    --------
    Indexing the unique edges with the inverse reproduces the sorted corner edges:

    ```python
    import numpy as np

    import ordito as od

    vertices, faces = od.creation.box()
    unique_edges, _inverse = od.edges.edges_unique(faces)
    inverse = od.edges.edges_unique_inverse(faces)
    corner_edges = od.edges.faces_to_edges(faces, sorted=True)
    print(np.array_equal(unique_edges.numpy()[inverse.numpy()], corner_edges.numpy()))  # True
    ```
    """
    require_same_device(faces=faces, edges_sorted=edges_sorted)
    return edges_unique(faces, edges_sorted=edges_sorted, n_vertices=n_vertices)[1]

def edges_unique_length(
    vertices: wp.array[wp.vec3],
    faces: wp.array[wp.int32],
    unique_edges: odt.Array2dInt32 | None = None,
    n_vertices: int | None = None,
    *,
    validate: bool = True,
) -> wp.array[wp.float32]:
    """
    Euclidean length of each unique undirected edge.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions.
    faces
        ``(3 * n_faces,)`` face index buffer.
    unique_edges
        ``(m, 2)`` precomputed unique edges, or ``None`` to compute them from ``faces``.
    n_vertices
        Total vertex count passed to [`edges_unique`][ordito.edges.edges_unique]. Ignored when
        ``unique_edges`` is already provided. When ``None`` it is taken from ``vertices`` rather
        than inferred from ``faces``, which would cost a host readback for a count this function
        was already handed.
    validate
        Forwarded to [`edges_unique`][ordito.edges.edges_unique] when ``unique_edges`` is built
        here; ignored when ``unique_edges`` is supplied. Pass ``False`` only where the face
        indices' range is structurally guaranteed.

    Returns
    -------
    wp.array[wp.float32]
        ``(m,)`` edge lengths on ``faces.device``.

    Raises
    ------
    TypeError
        If ``unique_edges`` is not a rank-2 ``wp.int32`` array.
    ValueError
        If ``unique_edges`` does not have exactly two columns, or if ``validate`` is ``True``,
        ``unique_edges`` is omitted and a face index is negative or reaches ``n_vertices``.
    RuntimeError
        If ``vertices``, ``faces`` and ``unique_edges`` are not all on one device.

    See Also
    --------
    [`trimesh.Trimesh.edges_unique_length`][]

    Examples
    --------
    The unit cube's 18 unique edges are 12 sides of length 1 and 6 diagonals of length ``sqrt(2)``:

    ```python
    import ordito as od

    vertices, faces = od.creation.box()
    lengths = od.edges.edges_unique_length(vertices, faces)
    print(lengths.shape)  # (18,)
    print(f"{lengths.numpy().min():.4f} {lengths.numpy().max():.4f}")  # 1.0000 1.4142
    print(int((lengths.numpy() > 1.1).sum()))  # 6
    ```
    """
    require_same_device(vertices=vertices, faces=faces, unique_edges=unique_edges)
    if unique_edges is not None:
        return _edge_lengths(vertices, unique_edges, "unique_edges")
    # ``vertices.shape[0]`` when the caller gave no count: ``edges_unique`` would otherwise infer it
    # from ``faces`` with ``array.index_bound``, a device reduction plus a host readback, for a
    # number every caller of *this* function already holds in the array it passed. It is also the
    # bound ``validate`` is documented against.
    if n_vertices is None:
        n_vertices = vertices.size
    n_faces = faces.size // 3
    if n_faces == 0:
        return _launch.empty(0, dtype=wp.float32, device=vertices.device)
    # ``edges_unique``'s sort and run scan, then the lengths written at each run's rank: the edge
    # table and the corner inverse it would emit are never built.
    radix, check = _face_edge_radix(faces, n_faces, n_vertices, validate)
    keys, order, ranks, n_edges = _sorted_edge_run_ranks(faces, n_faces, radix, check)
    lengths = _launch.empty(n_edges, dtype=wp.float32, device=vertices.device)
    _launch.launch(
        kernel_edges.emit_sorted_unique_edge_lengths,
        dim=3 * n_faces,
        inputs=[vertices, faces, keys, order, ranks, lengths],
        device=vertices.device,
    )
    return lengths

def edges_length(
    vertices: wp.array[wp.vec3], faces: wp.array[wp.int32], edges_in: odt.Array2dInt32 | None = None
) -> wp.array[wp.float32]:
    """
    Euclidean length of every directed edge (one per face half-edge).

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions.
    faces
        ``(3 * n_faces,)`` face index buffer.
    edges_in
        ``(3 * n_faces, 2)`` precomputed directed edges, or ``None`` to compute them from
        ``faces``.

    Returns
    -------
    wp.array[wp.float32]
        ``(3 * n_faces,)`` edge lengths on ``faces.device``.

    Raises
    ------
    TypeError
        If ``edges_in`` is not a rank-2 ``wp.int32`` array.
    ValueError
        If ``edges_in`` does not have exactly two columns.
    RuntimeError
        If ``vertices``, ``faces`` and ``edges_in`` are not all on one device.

    Examples
    --------
    One length per corner edge, so a shared edge appears once from each of its faces:

    ```python
    import ordito as od

    vertices, faces = od.creation.box()
    lengths = od.edges.edges_length(vertices, faces)
    print(lengths.shape)  # (36,)
    print(f"{lengths.numpy().sum():.4f}")  # 40.9706
    ```
    """
    require_same_device(vertices=vertices, faces=faces, edges_in=edges_in)
    if edges_in is not None:
        return _edge_lengths(vertices, edges_in, "edges_in")
    # Without a table the rows are the halfedges, so the lengths come straight off ``faces``.
    n_halfedges = faces.size // 3 * 3
    out = _launch.empty(n_halfedges, dtype=wp.float32, device=vertices.device)
    if n_halfedges > 0:
        _launch.launch(
            kernel_edges.halfedge_lengths,
            dim=n_halfedges,
            inputs=[vertices, faces, out],
            device=vertices.device,
        )
    return out

def face_edge_lengths(vertices: wp.array[wp.vec3], faces: wp.array[wp.int32]) -> odt.Array2dFloat32:
    """
    Per-face edge lengths, in the intrinsic column order the cotangent formulas expect.

    Column ``e`` holds the length of the edge *opposite* corner ``e``, matching
    [`cotmatrix_entries_intrinsic`][ordito.laplacian.cotmatrix_entries_intrinsic] and
    ``igl::cotmatrix_entries``' intrinsic overload. Unlike
    [`edges_unique_length`][ordito.edges.edges_unique_length] this is a per-*corner* table: an
    interior edge appears twice, which is what lets the entries be perturbed per face.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` mesh vertex positions.
    faces
        ``(3 * n_faces,)`` triangle index buffer.

    Returns
    -------
    odt.Array2dFloat32
        ``(n_faces, 3)`` edge lengths on ``vertices.device``.

    Raises
    ------
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`mollify_intrinsic`][ordito.laplacian.mollify_intrinsic]
    [`edges_unique_length`][ordito.edges.edges_unique_length]

    Examples
    --------
    Column ``k`` is the edge opposite corner ``k``; on a cube every face has two sides of length 1
    and one diagonal:

    ```python
    import numpy as np

    import ordito as od

    vertices, faces = od.creation.box()
    lengths = od.edges.face_edge_lengths(vertices, faces)
    print(lengths.shape)  # (12, 3)
    print(np.sort(lengths.numpy()[0]).round(4))  # [1.     1.     1.4142]
    ```
    """
    require_same_device(vertices=vertices, faces=faces)
    device = vertices.device
    n_faces = faces.size // 3
    lengths = odt.empty_2d((n_faces, 3), wp.float32, device=device)
    if n_faces == 0:
        return odt.as_array2d(lengths, wp.float32)
    _launch.launch(
        kernel_edges.face_edge_lengths,
        dim=n_faces,
        inputs=[vertices, faces, lengths],
        device=device,
    )
    return odt.as_array2d(lengths, wp.float32)

def mean_edge_length(vertices: wp.array[wp.vec3], faces: wp.array[wp.int32]) -> float:
    """
    Mean length of the ``3 * n_faces`` per-face edges, counting a shared edge once per face.

    Every face contributes all three of its edges, so an edge shared by two faces is counted twice
    and a boundary edge once -- which is what separates this from
    [`mean_unique_edge_length`][ordito.edges.mean_unique_edge_length], the same average over each
    edge once. That weighting is deliberate, not an oversight: it is the average
    libigl's ``CurvatureCalculator::getAverageEdge`` computes, which is what
    ``igl::principal_curvature`` calls to set its sphere-search radius, and what
    [`principal_curvature`][ordito.curvature.principal_curvature] therefore uses here.

    !!! note "There are two edge averages, and they are not interchangeable"
        [`mean_unique_edge_length`][ordito.edges.mean_unique_edge_length] counts each edge once,
        matching ``igl::avg_edge_length`` and MeshLab's ``avg_edge_length``. The two agree exactly
        on a closed manifold mesh -- every edge has two incident faces there, so the doubling is
        uniform -- and diverge by a fraction of a percent on anything with a boundary or a
        non-manifold edge.

        Reach for this one when reproducing libigl's curvature; reach for the unique-edge one when
        reproducing anything else, including the heat method's timestep.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions.
    faces
        ``(3 * n_faces,)`` face index buffer.

    Returns
    -------
    float
        Mean length over the ``3 * n_faces`` per-face edges. ``0.0`` when ``n_faces == 0``.

    Raises
    ------
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`mean_unique_edge_length`][ordito.edges.mean_unique_edge_length]
        The same average taken over the unique edges instead.
    [`edges_length`][ordito.edges.edges_length]
        The per-face lengths this averages.

    Examples
    --------
    On a closed mesh this equals the unique-edge mean; on an open grid the boundary edges count once
    and the interior ones twice, so the two averages differ:

    ```python
    import ordito as od

    vertices, faces = od.creation.box()
    print(f"{od.edges.mean_edge_length(vertices, faces):.4f}")  # 1.1381
    vertices, faces = od.creation.grid(count=(2, 2))
    print(f"{od.edges.mean_edge_length(vertices, faces):.4f}")  # 1.1381
    print(f"{od.edges.mean_unique_edge_length(vertices, faces):.4f}")  # 1.0828
    ```
    """
    require_same_device(vertices=vertices, faces=faces)
    n_halfedges = faces.size // 3 * 3
    if n_halfedges == 0:
        return 0.0
    # The total is ``edge_length_sum``'s, read back once.
    return float(read_scalar(edge_length_sum(vertices, faces), 0)) / n_halfedges

def edge_length_sum(vertices: wp.array[wp.vec3], faces: wp.array[wp.int32]) -> wp.array[wp.float32]:
    """
    Total length of the ``3 * n_faces`` per-face edges, left on the device.

    The sum [`mean_edge_length`][ordito.edges.mean_edge_length] divides by ``3 * n_faces``, for a
    caller whose kernels consume it on the device and would otherwise wait on a readback to pass
    the mean in as a scalar. A shared edge counts once per face, as in the mean. Summed in
    [`sum`][ordito.reduce.sum]'s order: fixed under
    [`is_deterministic`][ordito.config.is_deterministic], arrival order otherwise.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions.
    faces
        ``(3 * n_faces,)`` face index buffer.

    Returns
    -------
    wp.array[wp.float32]
        ``(1,)`` the total on ``vertices.device``; ``0`` when ``n_faces == 0``.

    Raises
    ------
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`mean_edge_length`][ordito.edges.mean_edge_length]
        The same total over the edge count, read back.

    Examples
    --------
    ```python
    import ordito as od

    vertices, faces = od.creation.box()
    total = od.edges.edge_length_sum(vertices, faces)
    print(total.shape, f"{total.numpy()[0] / faces.size:.4f}")  # (1,) 1.1381
    ```
    """
    require_same_device(vertices=vertices, faces=faces)
    device = vertices.device
    n_halfedges = faces.size // 3 * 3
    if n_halfedges == 0:
        return _launch.zeros(1, dtype=wp.float32, device=device)
    # The lengths are formed inside the sum's own blocks: no ``(3 * n_faces,)`` length buffer.
    terms = kernel_edges.HalfedgeLengthTerms()
    terms.n = n_halfedges
    terms.vertices = vertices
    terms.faces = faces
    return od.reduce.sum_of_terms(
        kernel_edges.HALFEDGE_LENGTH_SUM, [terms], n_halfedges, device=device
    )

def mean_unique_edge_length(
    vertices: wp.array[wp.vec3], faces: wp.array[wp.int32], *, validate: bool = True
) -> float:
    """
    Mean length of the **unique** undirected edges.

    Each edge is counted once, however many faces share it. This is the definition
    ``igl::avg_edge_length`` and MeshLab's ``avg_edge_length`` both use, and the one
    ``igl::heat_geodesics`` picks its diffusion timestep from -- so it is what
    [`heat_geodesic`][ordito.heat.heat_geodesic] and
    [`vector_heat_operators`][ordito.heat.vector_heat_operators] use here.

    See [`mean_edge_length`][ordito.edges.mean_edge_length] for the per-face average, how far the
    two diverge on an open mesh, and why libigl carries both.

    Parameters
    ----------
    vertices
        ``(n_vertices,)`` vertex positions.
    faces
        ``(3 * n_faces,)`` face index buffer.
    validate
        Forwarded to [`edges_unique`][ordito.edges.edges_unique]. Pass ``False`` only where the
        face indices' range is structurally guaranteed; it removes a host readback.

    Returns
    -------
    float
        Mean length over the unique undirected edges. ``0.0`` when ``n_faces == 0``.

    Raises
    ------
    ValueError
        If ``validate`` is ``True`` and a face index is negative or reaches ``n_vertices``.
    RuntimeError
        If ``vertices`` and ``faces`` are not all on one device.

    See Also
    --------
    [`mean_edge_length`][ordito.edges.mean_edge_length]
        The same average taken over the per-face edges instead.
    [`edges_unique_length`][ordito.edges.edges_unique_length]
        The per-edge lengths this averages.

    Examples
    --------
    The cube's 18 unique edges average ``(12 + 6 * sqrt(2)) / 18``:

    ```python
    import ordito as od

    vertices, faces = od.creation.box()
    print(f"{od.edges.mean_unique_edge_length(vertices, faces):.4f}")  # 1.1381
    ```
    """
    require_same_device(vertices=vertices, faces=faces)
    n_faces = faces.size // 3
    if n_faces == 0:
        return 0.0
    # The unique edges are the sorted corner keys' run starts, so their lengths and count fold
    # straight off the sort into per-block ``(sum, count, range flag)`` partials: no emitted edge
    # table, no length buffer, and one readback carrying the total, the count and the validation.
    n_vertices = vertices.size
    n = 3 * n_faces
    keys, order = _sorted_face_edge_keys(faces, n_faces, n_vertices, n_vertices)
    partials = _launch.empty(kernel_reduce.blocks_1d(n), dtype=wp.vec3, device=vertices.device)
    _launch.launch_tiled(
        kernel_edges.unique_edge_length_partials,
        dim=[partials.size],
        inputs=[vertices, keys, n, wp.uint64(n_vertices), order, 2 * n, partials],
        block_dim=TILE_1D,
        device=vertices.device,
    )
    total, count, code = od.reduce.sum(partials)
    if validate:
        _raise_on_range_code(int(code), n_vertices)
    return float(total) / float(count)