ENH: add diag_indices, tril_indices, triu_indices

Resolves #686. Adds the three index-generating functions that numpy,
jax, and cupy all have but that are missing from the array-api
standard and (so far) from this library.

Signatures follow array-api conventions: parameter `offset` (matching
`xp.linalg.diagonal`) instead of numpy's `k`; keyword-only arguments
for everything except `n`; `xp` is required (these functions have no
input array to infer from, following the `default_dtype` precedent).

Delegation:
- numpy/cupy/jax: forward directly (signatures match verbatim).
- dask: has tril/triu_indices but no diag_indices.
- torch: has tril/triu_indices but with (row, col, *, offset) signature
  returning a 2xN tensor rather than a tuple; delegation translates.
  No torch.diag_indices exists; falls through to generic.
- sparse, array-api-strictest: fall through to generic; marked xfail
  on those backends (no nonzero / data-dependent shapes).

Generic implementation uses `xp.arange` + broadcasting + `xp.nonzero`
for the triangle variants. Validation (n >= 0, ndim >= 1, m >= 0)
happens in the delegation layer so all backends produce consistent
ValueErrors.

Also fixes a pre-existing bug in tests/conftest.py's NumPyReadOnly
wrapper: `type(o)(*gen)` worked for namedtuples but failed for plain
tuples of length >= 2. Exposed here because these are the first
functions in the library that return a tuple of arrays.

Author: bruAristimunha

src/array_api_extra/__init__.py

@@ -5,6 +5,7 @@
     atleast_nd,
     cov,
     create_diagonal,
+    diag_indices,
     expand_dims,
     isclose,
     isin,
@@ -15,6 +16,8 @@
     searchsorted,
     setdiff1d,
     sinc,
+    tril_indices,
+    triu_indices,
     union1d,
 )
 from ._lib._at import at
@@ -40,6 +43,7 @@
     "cov",
     "create_diagonal",
     "default_dtype",
+    "diag_indices",
     "expand_dims",
     "isclose",
     "isin",
@@ -53,5 +57,7 @@
     "searchsorted",
     "setdiff1d",
     "sinc",
+    "tril_indices",
+    "triu_indices",
     "union1d",
 ]

src/array_api_extra/_delegation.py

@@ -22,16 +22,27 @@
     "atleast_nd",
     "cov",
     "create_diagonal",
+    "diag_indices",
     "expand_dims",
     "isclose",
     "nan_to_num",
     "one_hot",
     "pad",
     "searchsorted",
     "sinc",
+    "tril_indices",
+    "triu_indices",
 ]
 
 
+def _check_nonneg(**kwargs: int | None) -> None:  # numpydoc ignore=PR01
+    """Raise ``ValueError`` if any keyword argument is a negative int."""
+    for name, value in kwargs.items():
+        if value is not None and value < 0:
+            msg = f"{name} must be non-negative, got {value}"
+            raise ValueError(msg)
+
+
 def atleast_nd(x: Array, /, *, ndim: int, xp: ModuleType | None = None) -> Array:
     """
     Recursively expand the dimension of an array to at least `ndim`.
@@ -238,6 +249,47 @@ def create_diagonal(
     return _funcs.create_diagonal(x, offset=offset, xp=xp)
 
 
+def diag_indices(n: int, /, *, ndim: int = 2, xp: ModuleType) -> tuple[Array, ...]:
+    """
+    Return the indices to access the main diagonal of an array.
+
+    Equivalent to ``numpy.diag_indices``.
+
+    Parameters
+    ----------
+    n : int
+        The size of each dimension of the (hyper-)cube ``(n, n, ..., n)``
+        that the returned indices index into.
+    ndim : int, optional
+        The number of dimensions. Default: ``2``.
+    xp : array_namespace
+        The standard-compatible namespace to create the indices in.
+
+    Returns
+    -------
+    tuple of array
+        ``ndim`` 1-D integer arrays of length ``n`` that together index
+        the main diagonal of an array of shape ``(n,) * ndim``.
+
+    Examples
+    --------
+    >>> import array_api_strict as xp
+    >>> import array_api_extra as xpx
+    >>> rows, cols = xpx.diag_indices(3, xp=xp)
+    >>> rows
+    Array([0, 1, 2], dtype=array_api_strict.int64)
+    >>> cols
+    Array([0, 1, 2], dtype=array_api_strict.int64)
+    """
+    _check_nonneg(n=n)
+    if ndim < 1:
+        msg = f"ndim must be >= 1, got {ndim}"
+        raise ValueError(msg)
+    if is_numpy_namespace(xp) or is_cupy_namespace(xp) or is_jax_namespace(xp):
+        return xp.diag_indices(n, ndim=ndim)
+    return _funcs.diag_indices(n, ndim=ndim, xp=xp)
+
+
 def expand_dims(
     a: Array, /, *, axis: int | tuple[int, ...] = (0,), xp: ModuleType | None = None
 ) -> Array:
@@ -1150,3 +1202,109 @@ def union1d(a: Array, b: Array, /, *, xp: ModuleType | None = None) -> Array:
         return xp.union1d(a, b)
 
     return _funcs.union1d(a, b, xp=xp)
+
+
+def tril_indices(
+    n: int, /, *, offset: int = 0, m: int | None = None, xp: ModuleType
+) -> tuple[Array, Array]:
+    """
+    Return the indices of the lower triangle of an ``(n, m)`` array.
+
+    Equivalent to ``numpy.tril_indices`` with parameter ``k`` renamed to
+    ``offset`` to match ``xp.linalg.diagonal``'s naming.
+
+    Parameters
+    ----------
+    n : int
+        The row dimension of the array.
+    offset : int, optional
+        Diagonal offset; ``0`` (default) is the main diagonal. Corresponds
+        to ``k`` in ``numpy.tril_indices``.
+    m : int, optional
+        The column dimension. If ``None`` (default), assumed equal to `n`.
+    xp : array_namespace
+        The standard-compatible namespace to create the indices in.
+
+    Returns
+    -------
+    tuple of array
+        Row and column indices ``(rows, cols)`` of the lower triangle of
+        the ``(n, m)`` matrix, shifted by `offset`.
+
+    Examples
+    --------
+    >>> import array_api_strict as xp
+    >>> import array_api_extra as xpx
+    >>> rows, cols = xpx.tril_indices(3, xp=xp)
+    >>> rows
+    Array([0, 1, 1, 2, 2, 2], dtype=array_api_strict.int64)
+    >>> cols
+    Array([0, 0, 1, 0, 1, 2], dtype=array_api_strict.int64)
+    """
+    _check_nonneg(n=n, m=m)
+    if (
+        is_numpy_namespace(xp)
+        or is_cupy_namespace(xp)
+        or is_jax_namespace(xp)
+        or is_dask_namespace(xp)
+    ):
+        return xp.tril_indices(n, k=offset, m=m)
+    if is_torch_namespace(xp):
+        # `torch.tril_indices` returns a 2xN tensor, not a tuple, and
+        # takes (row, col) rather than (n, *, m=None).
+        cols = n if m is None else m
+        idx = xp.tril_indices(n, cols, offset=offset)
+        return (idx[0], idx[1])
+    return _funcs.tril_indices(n, offset=offset, m=m, xp=xp)
+
+
+def triu_indices(
+    n: int, /, *, offset: int = 0, m: int | None = None, xp: ModuleType
+) -> tuple[Array, Array]:
+    """
+    Return the indices of the upper triangle of an ``(n, m)`` array.
+
+    Equivalent to ``numpy.triu_indices`` with parameter ``k`` renamed to
+    ``offset`` to match ``xp.linalg.diagonal``'s naming.
+
+    Parameters
+    ----------
+    n : int
+        The row dimension of the array.
+    offset : int, optional
+        Diagonal offset; ``0`` (default) is the main diagonal. Corresponds
+        to ``k`` in ``numpy.triu_indices``.
+    m : int, optional
+        The column dimension. If ``None`` (default), assumed equal to `n`.
+    xp : array_namespace
+        The standard-compatible namespace to create the indices in.
+
+    Returns
+    -------
+    tuple of array
+        Row and column indices ``(rows, cols)`` of the upper triangle of
+        the ``(n, m)`` matrix, shifted by `offset`.
+
+    Examples
+    --------
+    >>> import array_api_strict as xp
+    >>> import array_api_extra as xpx
+    >>> rows, cols = xpx.triu_indices(3, xp=xp)
+    >>> rows
+    Array([0, 0, 0, 1, 1, 2], dtype=array_api_strict.int64)
+    >>> cols
+    Array([0, 1, 2, 1, 2, 2], dtype=array_api_strict.int64)
+    """
+    _check_nonneg(n=n, m=m)
+    if (
+        is_numpy_namespace(xp)
+        or is_cupy_namespace(xp)
+        or is_jax_namespace(xp)
+        or is_dask_namespace(xp)
+    ):
+        return xp.triu_indices(n, k=offset, m=m)
+    if is_torch_namespace(xp):
+        cols = n if m is None else m
+        idx = xp.triu_indices(n, cols, offset=offset)
+        return (idx[0], idx[1])
+    return _funcs.triu_indices(n, offset=offset, m=m, xp=xp)

src/array_api_extra/_lib/_funcs.py

@@ -28,13 +28,16 @@
     "broadcast_shapes",
     "cov",
     "create_diagonal",
+    "diag_indices",
     "expand_dims",
     "kron",
     "nunique",
     "pad",
     "searchsorted",
     "setdiff1d",
     "sinc",
+    "tril_indices",
+    "triu_indices",
 ]
 
 
@@ -346,6 +349,41 @@ def create_diagonal(
     return xp.reshape(diag, (*batch_dims, n, n))
 
 
+def diag_indices(
+    n: int, /, *, ndim: int = 2, xp: ModuleType
+) -> tuple[Array, ...]:  # numpydoc ignore=PR01,RT01
+    """See docstring in array_api_extra._delegation."""
+    idx = xp.arange(n)
+    return (idx,) * ndim
+
+
+def _tri_indices(
+    n: int, *, offset: int, m: int | None, upper: bool, xp: ModuleType
+) -> tuple[Array, Array]:  # numpydoc ignore=PR01,RT01
+    """Shared implementation for `tril_indices` and `triu_indices`."""
+    cols = n if m is None else m
+    rows = xp.arange(n)[:, None]
+    cols_a = xp.arange(cols)[None, :]
+    delta = cols_a - rows
+    mask = delta >= offset if upper else delta <= offset
+    r, c = xp.nonzero(mask)
+    return (r, c)
+
+
+def tril_indices(
+    n: int, /, *, offset: int = 0, m: int | None = None, xp: ModuleType
+) -> tuple[Array, Array]:  # numpydoc ignore=PR01,RT01
+    """See docstring in array_api_extra._delegation."""
+    return _tri_indices(n, offset=offset, m=m, upper=False, xp=xp)
+
+
+def triu_indices(
+    n: int, /, *, offset: int = 0, m: int | None = None, xp: ModuleType
+) -> tuple[Array, Array]:  # numpydoc ignore=PR01,RT01
+    """See docstring in array_api_extra._delegation."""
+    return _tri_indices(n, offset=offset, m=m, upper=True, xp=xp)
+
+
 def default_dtype(
     xp: ModuleType,
     kind: Literal[

tests/conftest.py

@@ -96,9 +96,12 @@ def as_readonly(o: T) -> T:  # numpydoc ignore=PR01,RT01
                 # Cannot interpret as a data type
                 return o
 
-            # This works with namedtuples too
             if isinstance(o, tuple | list):
-                return type(o)(*(as_readonly(i) for i in o))  # type: ignore[arg-type,return-value] # pyright: ignore[reportArgumentType]
+                # namedtuple wants positional args; plain tuple/list wants an iterable.
+                items = (as_readonly(i) for i in o)
+                if hasattr(o, "_fields"):
+                    return type(o)(*items)  # type: ignore[arg-type,return-value] # pyright: ignore[reportArgumentType]
+                return type(o)(items)  # type: ignore[return-value]
 
             return o