move rng to be scanpy like

.gitignore

@@ -7,6 +7,7 @@ __pycache__/
 /benchmarks/
 test-data/
 .vscode/
+PanGPA.log
 
 # Distribution / packaging
 /dist/

docs/release-notes/0.17.0.md

@@ -8,6 +8,7 @@
 * Split {func}`~rapids_singlecell.gr.calculate_niche` into {func}`~rapids_singlecell.gr.calculate_niche_neighborhood`, {func}`~rapids_singlecell.gr.calculate_niche_utag` and {func}`~rapids_singlecell.gr.calculate_niche_cellcharter`, with ``mask``, ``library_key`` and cross-flavor ``min_niche_size``, following {mod}`squidpy` {pr}`758` {smaller}`S Dicks`
 * Speed up {func}`~rapids_singlecell.pp.harmony_integrate` and make it reproducible by seeding k-means from a deterministic `float64` fit on a bounded, batch-stratified random subsample instead of a non-reproducible `float32` fit over all cells. `dtype` now defaults to `numpy.float32` {pr}`756` {smaller}`S Dicks`
 * Derive unset {func}`~rapids_singlecell.pp.harmony_integrate` stopping rules from `flavor`: `harmony2` follows Harmony2 defaults, `harmony1` still follows harmony-pytorch {pr}`756` {smaller}`S Dicks`
+* Add support for {class}`numpy.random.Generator` to all functions previously accepting a ``random_state`` parameter, which is renamed to ``rng`` following `scanpy`. Passing ``random_state`` still works and keeps producing the same results; passing both raises {smaller}`S Dicks`
 
 ```{rubric} Misc
 ```
@@ -17,3 +18,7 @@
 ```{rubric} Deprecations
 ```
 * Deprecate {func}`~rapids_singlecell.gr.calculate_niche`; the new functions replace ``copy`` with ``inplace`` and drop ``gmm_init`` and the leiden ``random_state`` {pr}`758` {smaller}`S Dicks`
+
+```{rubric} Removals
+```
+* Remove the ``random_seed`` parameter of {func}`~rapids_singlecell.pp.scrublet_simulate_doublets`, following `scanpy`. Use ``rng`` (or ``random_state``) instead {smaller}`S Dicks`

src/rapids_singlecell/_utils/_random.py

@@ -0,0 +1,96 @@
+from __future__ import annotations
+
+from collections.abc import Sequence
+from functools import wraps
+from typing import TYPE_CHECKING
+
+import numpy as np
+
+if TYPE_CHECKING:
+    from collections.abc import Callable
+
+__all__ = [
+    "RNGLike",
+    "SeedLike",
+    "_LegacyRandom",
+    "_LegacyRng",
+    "_accepts_legacy_random_state",
+    "_legacy_random_state",
+    "_seed_from_rng",
+]
+
+type SeedLike = int | np.integer | Sequence[int] | np.random.SeedSequence
+type RNGLike = np.random.Generator | np.random.BitGenerator
+type _LegacyRandom = int | np.random.RandomState | None
+
+_SEED_BOUND = 2**32
+"""cuML, cuGraph and CuPy all take a 32-bit unsigned seed."""
+
+
+class _LegacyRng:
+    """Marks a seed that arrived through the superseded `random_state` argument.
+
+    Unlike scanpy's class of the same name, this is only a marker. Integer-only
+    GPU consumers use it to recover the legacy seed, while host-side consumers
+    can recover the original :class:`~numpy.random.RandomState` object and
+    continue its exact stream.
+    """
+
+    __slots__ = ("arg",)
+
+    def __init__(self, arg: _LegacyRandom) -> None:
+        self.arg = arg
+
+    def __repr__(self) -> str:
+        return f"_LegacyRng({self.arg!r})"
+
+
+def _accepts_legacy_random_state[**P, R](
+    default: _LegacyRandom, /
+) -> Callable[[Callable[P, R]], Callable[P, R]]:
+    """Let a function taking `rng` still be called with `random_state`.
+
+    A `random_state` argument is wrapped in a :class:`_LegacyRng` and passed as
+    `rng`. Passing both is an error. If neither is given, `default` is used.
+    """
+
+    def decorator(func: Callable[P, R]) -> Callable[P, R]:
+        @wraps(func)
+        def wrapper(*args: P.args, **kwargs: P.kwargs) -> R:
+            match "random_state" in kwargs, "rng" in kwargs:
+                case True, True:
+                    raise TypeError("Specify at most one of `rng` and `random_state`.")
+                case True, False:
+                    kwargs["rng"] = _LegacyRng(kwargs.pop("random_state"))
+                case False, False:
+                    kwargs["rng"] = _LegacyRng(default)
+            return func(*args, **kwargs)
+
+        return wrapper
+
+    return decorator
+
+
+def _seed_from_rng(rng: SeedLike | RNGLike | _LegacyRng | None, /) -> int | None:
+    """The integer seed to hand to cuML, cuGraph, CuPy or a kernel.
+
+    A `random_state` integer is forwarded untouched so that existing calls keep
+    producing the same results; anything else is drawn from the generator.
+    """
+    if isinstance(rng, _LegacyRng):
+        if rng.arg is None:
+            return None
+        if isinstance(rng.arg, (int, np.integer)):
+            return int(rng.arg)
+        return int(rng.arg.randint(0, _SEED_BOUND))
+    return int(np.random.default_rng(rng).integers(0, _SEED_BOUND))
+
+
+def _legacy_random_state(
+    rng: SeedLike | RNGLike | _LegacyRng | None, /
+) -> _LegacyRandom:
+    """The value to hand to host-side APIs that still take a `random_state`."""
+    if isinstance(rng, _LegacyRng):
+        return rng.arg
+    [bit_generator] = np.random.default_rng(rng).bit_generator.spawn(1)
+    return np.random.RandomState(bit_generator)

src/rapids_singlecell/preprocessing/_harmony_integrate.py

@@ -6,12 +6,20 @@
 import cupy as cp
 import numpy as np
 
+from rapids_singlecell._utils._random import (
+    RNGLike,
+    SeedLike,
+    _accepts_legacy_random_state,
+    _seed_from_rng,
+)
+
 if TYPE_CHECKING:
     from anndata import AnnData
 
     from ._harmony import COLSUM_ALGO
 
 
+@_accepts_legacy_random_state(0)
 def harmony_integrate(
     adata: AnnData,
     key: str | list[str],
@@ -34,7 +42,7 @@ def harmony_integrate(
     correction_method: Literal["fast", "batched"] | None = None,
     colsum_algo: COLSUM_ALGO | None = None,
     block_proportion: float = 0.05,
-    random_state: int = 0,
+    rng: SeedLike | RNGLike | None = None,
     verbose: bool = False,
 ) -> None:
     """Integrate different experiments using the Harmony algorithm :cite:p:`Korsunsky2019,Patikas2026`.
@@ -158,8 +166,9 @@ def harmony_integrate(
         Proportion of cells updated per clustering sub-iteration.
         Smaller values produce more stochastic updates.
         Larger values are faster but may converge to different solutions.
-    random_state
-        Random seed for reproducibility.
+    rng
+        Random seed or :class:`~numpy.random.Generator` for reproducibility.
+        The superseded `random_state` argument is still accepted.
     verbose
         Whether to print benchmarking and convergence information.
 
@@ -169,6 +178,8 @@ def harmony_integrate(
     containing principal components adjusted by Harmony
     such that different experiments are integrated.
     """
+    random_state = _seed_from_rng(rng)
+
     from ._harmony import harmonize
 
     # Resolve flavor into internal flags

src/rapids_singlecell/preprocessing/_neighbors/__init__.py

@@ -6,6 +6,12 @@
 import cupy as cp
 import numpy as np
 
+from rapids_singlecell._utils._random import (
+    RNGLike,
+    SeedLike,
+    _accepts_legacy_random_state,
+    _seed_from_rng,
+)
 from rapids_singlecell.preprocessing._neighbors._helper import (
     _check_metrics,
     _check_neighbors_X,
@@ -32,13 +38,14 @@
 ]
 
 
+@_accepts_legacy_random_state(0)
 def neighbors(
     adata: AnnData,
     n_neighbors: int = 15,
     n_pcs: int | None = None,
     *,
     use_rep: str | None = None,
-    random_state: AnyRandom = 0,
+    rng: SeedLike | RNGLike | None = None,
     algorithm: _Algorithms = "brute",
     metric: _Metrics = "euclidean",
     metric_kwds: Mapping[str, Any] = MappingProxyType({}),
@@ -70,8 +77,9 @@ def neighbors(
         If None, the representation is chosen automatically: For .n_vars < 50, .X
         is used, otherwise `'X_pca'` is used. If `'X_pca'` is not present, it's
         computed with default parameters or `n_pcs` if present.
-    random_state
-        A numpy random seed.
+    rng
+        Random seed or :class:`~numpy.random.Generator` for reproducibility.
+        The superseded `random_state` argument is still accepted.
     algorithm
         The query algorithm to use. Valid options are:
             `'brute'`
@@ -174,6 +182,8 @@ def neighbors(
             neighbors.
 
     """
+    random_state = _seed_from_rng(rng)
+
     adata = adata.copy() if copy else adata
 
     if adata.is_view:
@@ -246,14 +256,15 @@ def neighbors(
     return adata if copy else None
 
 
+@_accepts_legacy_random_state(0)
 def bbknn(
     adata: AnnData,
     neighbors_within_batch: int = 3,
     n_pcs: int | None = None,
     *,
     batch_key: str | None = None,
     use_rep: str | None = None,
-    random_state: AnyRandom = 0,
+    rng: SeedLike | RNGLike | None = None,
     algorithm: _Algorithms_bbknn = "brute",
     metric: _Metrics = "euclidean",
     metric_kwds: Mapping[str, Any] = MappingProxyType({}),
@@ -284,8 +295,9 @@ def bbknn(
         If `None`, the representation is chosen automatically: For `.n_vars < 50`, `.X`
         is used, otherwise `'X_pca'` is used. If `'X_pca'` is not present, it's
         computed with default parameters or `n_pcs` if present.
-    random_state
-        A numpy random seed.
+    rng
+        Random seed or :class:`~numpy.random.Generator` for reproducibility.
+        The superseded `random_state` argument is still accepted.
     algorithm
         The query algorithm to use. Valid options are:
 
@@ -359,6 +371,8 @@ def bbknn(
     connectivities and distances.
     """
 
+    random_state = _seed_from_rng(rng)
+
     if batch_key is None:
         raise ValueError("Please provide a batch key to perform batch-balanced KNN.")
 

src/rapids_singlecell/preprocessing/_pca.py

@@ -12,6 +12,12 @@
 from scipy.sparse import issparse
 
 from rapids_singlecell._compat import DaskArray
+from rapids_singlecell._utils._random import (
+    RNGLike,
+    SeedLike,
+    _accepts_legacy_random_state,
+    _seed_from_rng,
+)
 from rapids_singlecell.get import X_to_GPU, _check_mask, _get_obs_rep
 
 from ._utils import _check_gpu_X
@@ -62,14 +68,15 @@ def _resolve_mask_var(
     return mask_var_param, _check_mask(adata, mask_var, "var")
 
 
+@_accepts_legacy_random_state(0)
 def pca(
     data: AnnData | ArrayTypesDask,
     n_comps: int | None = None,
     *,
     layer: str = None,
     zero_center: bool = True,
     svd_solver: str | None = None,
-    random_state: int | None = 0,
+    rng: SeedLike | RNGLike | None = None,
     mask_var: NDArray[np.bool] | str | None = _empty,
     use_highly_variable: bool | None = None,
     dtype: str = "float32",
@@ -145,8 +152,9 @@ def pca(
         `'jacobi'`
             cuML: Jacobi iterative solver. Faster but less accurate. For dense arrays only.
 
-    random_state
-        Random state for initialization.
+    rng
+        Random seed or :class:`~numpy.random.Generator` for initialization.
+        The superseded `random_state` argument is still accepted.
 
     mask_var
         Mask to use for the PCA computation.
@@ -216,6 +224,8 @@ def pca(
                 Explained variance, equivalent to the eigenvalues of the \
                 covariance matrix.
     """
+    random_state = _seed_from_rng(rng)
+
     if not isinstance(data, AnnData):
         if layer is not None:
             raise ValueError("`layer` can only be used with an AnnData object.")