(feat): read_elem_as_dask method

docs/api.md

@@ -121,6 +121,7 @@ Low level methods for reading and writing elements of an `AnnData` object to a s
 
    experimental.read_elem
    experimental.write_elem
+   experimental.read_elem_as_dask
 ```
 
 Utilities for customizing the IO process:

docs/release-notes/0.11.0.md

@@ -8,6 +8,7 @@
 * Add `should_remove_unused_categories` option to `anndata.settings` to override current behavior.  Default is `True` (i.e., previous behavior).  Please refer to the [documentation](https://anndata.readthedocs.io/en/latest/generated/anndata.settings.html) for usage.  {pr}`1340` {user}`ilan-gold`
 * `scipy.sparse.csr_array` and `scipy.sparse.csc_array` are now supported when constructing `AnnData` objects {pr}`1028` {user}`ilan-gold` {user}`isaac-virshup`
 * Add `should_check_uniqueness` option to `anndata.settings` to override current behavior.  Default is `True` (i.e., previous behavior).  Please refer to the [documentation](https://anndata.readthedocs.io/en/latest/generated/anndata.settings.html) for usage.  {pr}`1507` {user}`ilan-gold`
+* Add :func:`~anndata.experimental.read_elem_as_dask` function to handle i/o with sparse and dense arrays {pr}`1469` {user}`ilan-gold`
 * Add functionality to write from GPU {class}`dask.array.Array` to disk {pr}`1550` {user}`ilan-gold`
 
 #### Bugfix

pyproject.toml

@@ -66,7 +66,7 @@ dev = [
     "pytest-xdist",
 ]
 doc = [
-    "sphinx>=4.4",
+    "sphinx>=7.4.6",
     "sphinx-book-theme>=1.1.0",
     "sphinx-autodoc-typehints>=2.2.0",
     "sphinx-issues",

src/anndata/_core/file_backing.py

@@ -3,7 +3,7 @@
 import weakref
 from collections.abc import Mapping
 from functools import singledispatch
-from pathlib import Path
+from pathlib import Path, PurePosixPath
 from typing import TYPE_CHECKING
 
 import h5py
@@ -175,3 +175,18 @@ def _(x):
 @filename.register(ZarrGroup)
 def _(x):
     return x.store.path
+
+
+@singledispatch
+def get_elem_name(x):
+    raise NotImplementedError(f"Not implemented for {type(x)}")
+
+
+@get_elem_name.register(h5py.Group)
+def _(x):
+    return x.name
+
+
+@get_elem_name.register(ZarrGroup)
+def _(x):
+    return PurePosixPath(x.path).name

src/anndata/_io/specs/__init__.py

@@ -1,21 +1,25 @@
 from __future__ import annotations
 
-from . import methods
+from . import lazy_methods, methods
 from .registry import (
+    _LAZY_REGISTRY,  # noqa: F401
     _REGISTRY,  # noqa: F401
     IOSpec,
     Reader,
     Writer,
     get_spec,
     read_elem,
+    read_elem_as_dask,
     write_elem,
 )
 
 __all__ = [
     "methods",
+    "lazy_methods",
     "write_elem",
     "get_spec",
     "read_elem",
+    "read_elem_as_dask",
     "Reader",
     "Writer",
     "IOSpec",

src/anndata/_io/specs/lazy_methods.py

@@ -0,0 +1,173 @@
+from __future__ import annotations
+
+from contextlib import contextmanager
+from functools import wraps
+from pathlib import Path
+from typing import TYPE_CHECKING
+
+import h5py
+import numpy as np
+from scipy import sparse
+
+import anndata as ad
+from anndata._core.file_backing import filename, get_elem_name
+from anndata.compat import H5Array, H5Group, ZarrArray, ZarrGroup
+
+from .registry import _LAZY_REGISTRY, IOSpec
+
+if TYPE_CHECKING:
+    from collections.abc import Callable, Mapping, Sequence
+    from typing import Concatenate, Literal, ParamSpec, TypeVar
+
+    from anndata.compat import DaskArray
+
+    from .registry import DaskReader
+
+    BlockInfo2D = Mapping[
+        Literal[None],
+        dict[str, tuple[tuple[int, int], tuple[int, int]]],
+    ]
+
+    BlockInfoND = Mapping[
+        Literal[None],
+        dict[str, Sequence[tuple[int, int]]],
+    ]
+
+    B = TypeVar("B", BlockInfo2D, BlockInfoND)
+    P = ParamSpec("P")
+    R = TypeVar("R")
+
+
+@contextmanager
+def maybe_open_h5(path_or_group: Path | ZarrGroup, elem_name: str):
+    if not isinstance(path_or_group, Path):
+        yield path_or_group
+        return
+    file = h5py.File(path_or_group, "r")
+    try:
+        yield file[elem_name]
+    finally:
+        file.close()
+
+
+_DEFAULT_STRIDE = 1000
+
+
+def compute_chunk_layout_for_axis_shape(
+    chunk_axis_shape: int, full_axis_shape: int
+) -> tuple[int, ...]:
+    n_strides, rest = np.divmod(full_axis_shape, chunk_axis_shape)
+    chunk = (chunk_axis_shape,) * n_strides
+    if rest > 0:
+        chunk += (rest,)
+    return chunk
+
+
+def require_block_info(
+    f: Callable[Concatenate[B, P], R],
+) -> Callable[Concatenate[B | None, P], R]:
+    @wraps(f)
+    def wrapper(block_info: B | None = None, *args: P.args, **kwargs: P.kwargs) -> R:
+        if block_info is None:
+            msg = "Block info is required"
+            raise ValueError(msg)
+        return f(block_info, *args, **kwargs)
+
+    return wrapper
+
+
+@_LAZY_REGISTRY.register_read(H5Group, IOSpec("csc_matrix", "0.1.0"))
+@_LAZY_REGISTRY.register_read(H5Group, IOSpec("csr_matrix", "0.1.0"))
+@_LAZY_REGISTRY.register_read(ZarrGroup, IOSpec("csc_matrix", "0.1.0"))
+@_LAZY_REGISTRY.register_read(ZarrGroup, IOSpec("csr_matrix", "0.1.0"))
+def read_sparse_as_dask(
+    elem: H5Group | ZarrGroup,
+    *,
+    _reader: DaskReader,
+    chunks: tuple[int, ...] | None = None,  # only tuple[int, int] is supported here
+) -> DaskArray:
+    import dask.array as da
+
+    path_or_group = Path(filename(elem)) if isinstance(elem, H5Group) else elem
+    elem_name = get_elem_name(elem)
+    shape: tuple[int, int] = tuple(elem.attrs["shape"])
+    dtype = elem["data"].dtype
+    is_csc: bool = elem.attrs["encoding-type"] == "csc_matrix"
+
+    stride: int = _DEFAULT_STRIDE
+    major_dim, minor_dim = (1, 0) if is_csc else (0, 1)
+    if chunks is not None:
+        if len(chunks) != 2:
+            raise ValueError("`chunks` must be a tuple of two integers")
+        if chunks[minor_dim] != shape[minor_dim]:
+            raise ValueError(
+                "Only the major axis can be chunked. "
+                f"Try setting chunks to {((-1, _DEFAULT_STRIDE) if is_csc else (_DEFAULT_STRIDE, -1))}"
+            )
+        stride = chunks[major_dim]
+
+    @require_block_info
+    def make_dask_chunk(block_info: BlockInfo2D):
+        # We need to open the file in each task since `dask` cannot share h5py objects when using `dask.distributed`
+        # https://github.com/scverse/anndata/issues/1105
+        with maybe_open_h5(path_or_group, elem_name) as f:
+            mtx = ad.experimental.sparse_dataset(f)
+            range_i, range_j = block_info[None]["array-location"]
+            chunk = mtx[slice(*range_i), slice(*range_j)]
+        return chunk
+
+    shape_minor, shape_major = shape if is_csc else shape[::-1]
+    chunks_major = compute_chunk_layout_for_axis_shape(stride, shape_major)
+    chunks_minor = (shape_minor,)
+    chunk_layout = (
+        (chunks_minor, chunks_major) if is_csc else (chunks_major, chunks_minor)
+    )
+    memory_format = sparse.csc_matrix if is_csc else sparse.csr_matrix
+    da_mtx = da.map_blocks(
+        make_dask_chunk,
+        dtype=dtype,
+        chunks=chunk_layout,
+        meta=memory_format((0, 0), dtype=dtype),
+    )
+    return da_mtx
+
+
+@_LAZY_REGISTRY.register_read(H5Array, IOSpec("array", "0.2.0"))
+def read_h5_array(
+    elem: H5Array, *, _reader: DaskReader, chunks: tuple[int, ...] | None = None
+) -> DaskArray:
+    import dask.array as da
+
+    path = Path(elem.file.filename)
+    elem_name: str = elem.name
+    shape = tuple(elem.shape)
+    dtype = elem.dtype
+    chunks: tuple[int, ...] = (
+        chunks if chunks is not None else (_DEFAULT_STRIDE,) * len(shape)
+    )
+
+    @require_block_info
+    def make_dask_chunk(block_info: BlockInfoND):
+        with maybe_open_h5(path, elem_name) as f:
+            idx = ()
+            for i in range(len(shape)):
+                (start, stop) = block_info[None]["array-location"][i]
+                idx += (slice(start, stop),)
+            return f[idx]
+
+    chunk_layout = tuple(
+        compute_chunk_layout_for_axis_shape(chunks[i], shape[i])
+        for i in range(len(shape))
+    )
+
+    return da.map_blocks(make_dask_chunk, dtype=dtype, chunks=chunk_layout)
+
+
+@_LAZY_REGISTRY.register_read(ZarrArray, IOSpec("array", "0.2.0"))
+def read_zarr_array(
+    elem: ZarrArray, *, _reader: DaskReader, chunks: tuple[int, ...] | None = None
+) -> DaskArray:
+    chunks: tuple[int, ...] = chunks if chunks is not None else elem.chunks
+    import dask.array as da
+
+    return da.from_zarr(elem, chunks=chunks)