refactor: robustness and cleanup fixes from the 0.7.0 review

docs/api.md

@@ -74,7 +74,7 @@ the infercnvpy equivalents operate on CNV data.
 
 ### Spatial
 
-Tools for spatial data. They read a cells x cells neighbour graph from `adata.obsp`,
+Tools for spatial data. They read a cells x cells neighbor graph from `adata.obsp`,
 e.g. `obsp["spatial_connectivities"]` from {func}`squidpy.gr.spatial_neighbors_knn`.
 
 ```{eval-rst}

docs/installation.md

@@ -57,7 +57,7 @@ uv add git+https://github.com/scverse/infercnvpy.git
 
 - `infercnvpy[dask]` installs dask and the lazy readers of anndata (`anndata[dask,lazy]`) for {doc}`large and lazy datasets <large_data>`.
 - `infercnvpy[gtf]` installs [gtfparse][] for {func}`infercnvpy.io.genomic_position_from_gtf`.
-- `infercnvpy[spatial]` installs [squidpy][] for the spatial neighbour graph that {func}`infercnvpy.tl.smooth_cnv`,
+- `infercnvpy[spatial]` installs [squidpy][] for the spatial neighbor graph that {func}`infercnvpy.tl.smooth_cnv`,
   {func}`infercnvpy.tl.spatial_regions` and {func}`infercnvpy.pl.spatial_cnv` read from `adata.obsp`.
   The functions themselves work with any graph in `adata.obsp`, so squidpy is only needed to build it.
 - `infercnvpy[copykat]` installs [rpy2][] for {func}`infercnvpy.tl.copykat`.

docs/large_data.md

@@ -9,7 +9,9 @@ as a dense array. The expression matrix can be
 - a [dask](https://www.dask.org/) array, also with sparse `scipy` blocks, e.g. a matrix loaded lazily with
   {func}`anndata.experimental.read_elem_lazy` or {func}`anndata.experimental.read_lazy`, or a table from
   [spatialdata-io](https://spatialdata.scverse.org/projects/io/),
-- the matrix of a backed AnnData (`anndata.read_h5ad(..., backed="r")`).
+- the matrix of a backed AnnData (`anndata.read_h5ad(..., backed="r")`), stored as CSR or dense. Reading a block
+  of rows of a CSC matrix reads all of it, so a backed CSC matrix is read into memory once, as CSR, with a warning.
+  Convert it before writing the file (`adata.X = adata.X.tocsr()`) to read it block by block.
 
 Install the optional dependency with `pip install "infercnvpy[dask]"` to use dask arrays.
 
@@ -41,20 +43,20 @@ With `lazy=True`, the result for a dask array stays a dask array of CSR blocks,
 {meth}`dask.array.Array.map_blocks`. The noise band is computed first, with one dask computation. `lazy=True` has to
 be passed explicitly, also when a {class}`distributed.Client` is active. In this mode the dask scheduler decides how
 many blocks run at the same time, each in a single thread, and `n_jobs` is ignored. infercnvpy never starts numba
-threads inside a dask worker thread, because numba's default threading layer on macOS aborts the process when it is
-used from several threads at once.
+threads inside a dask task, with any scheduler: the scheduler's threads or processes already use the cores, and
+numba's default threading layer on macOS aborts the process when it is used from several threads at once.
 
 :::{note}
-infercnv expects log-normalised expression (see {ref}`input-data`). For a dask array, {func}`scanpy.pp.normalize_total`
-and {func}`scanpy.pp.log1p` stay lazy. A backed AnnData cannot be normalised in place by scanpy, so normalise the data
+infercnv expects log-normalized expression (see {ref}`input-data`). For a dask array, {func}`scanpy.pp.normalize_total`
+and {func}`scanpy.pp.log1p` stay lazy. A backed AnnData cannot be normalized in place by scanpy, so normalize the data
 before you write the file.
 :::
 
 ## Example: a lazily loaded table
 
 Spatial technologies often produce tables with a lazy dask `X` of CSR blocks that are indexed by Ensembl ID. Next to
 the genes, such a table can contain control probes or codewords without a genomic position, which you remove before
-normalisation.
+normalization.
 
 ```python
 import anndata as ad

docs/notebooks/reproduce_infercnv.ipynb

@@ -202,7 +202,7 @@
    "id": "artificial-animation",
    "metadata": {},
    "source": [
-    "It contains malignant cells from four tumours and two types of non-malignant cells."
+    "It contains malignant cells from four tumors and two types of non-malignant cells."
    ]
   },
   {

docs/notebooks/spatial_xenium.ipynb

@@ -9,22 +9,22 @@
     "\n",
     "Imaging-based spatial transcriptomics (Xenium, CosMx, MERSCOPE) measures a few thousand genes in hundreds of\n",
     "thousands of cells, with a few hundred counts per cell. inferCNV-style profiles of single cells are noisy at that\n",
-    "depth, but the tissue coordinates add information: neighbouring tumour cells usually share their copy number state.\n",
+    "depth, but the tissue coordinates add information: neighboring tumor cells usually share their copy number state.\n",
     "This tutorial shows the spatial layer of infercnvpy:\n",
     "\n",
-    "- {func}`infercnvpy.tl.smooth_cnv` averages each cell's CNV profile with its spatial neighbours,\n",
-    "- {func}`infercnvpy.tl.spatial_regions` makes per-cell labels (here tumour/normal calls) spatially coherent: a cell\n",
-    "  takes the label of its neighbourhood when most of its neighbours disagree with it (more than 75% at the default\n",
+    "- {func}`infercnvpy.tl.smooth_cnv` averages each cell's CNV profile with its spatial neighbors,\n",
+    "- {func}`infercnvpy.tl.spatial_regions` makes per-cell labels (here tumor/normal calls) spatially coherent: a cell\n",
+    "  takes the label of its neighborhood when most of its neighbors disagree with it (more than 75% at the default\n",
     "  `beta`),\n",
     "- {func}`infercnvpy.pl.spatial_cnv` plots the mean CNV of a chromosome or arm on the tissue.\n",
     "\n",
-    "All three read a spatial neighbour graph from `adata.obsp`, here from {func}`squidpy.gr.spatial_neighbors_knn`; any\n",
+    "All three read a spatial neighbor graph from `adata.obsp`, here from {func}`squidpy.gr.spatial_neighbors_knn`; any\n",
     "cells x cells graph works. To run this tutorial, install squidpy, gtfparse and the lazy readers of anndata with the optional extras, and s3fs to\n",
     "read the data from S3: `pip install \"infercnvpy[spatial,gtf,dask]\" s3fs`.\n",
     "\n",
-    "The ideas are not new. InSituCNV {cite}`Jensen2025` smooths infercnvpy profiles over kNN neighbours on Xenium and CosMx,\n",
-    "SlideCNA {cite}`Zhang2025` pools spatial neighbours on Slide-seq, and STARCH {cite}`Elyanow2021` and CalicoST\n",
-    "{cite}`Ma2024` make neighbouring cells prefer the same clone label. In a benchmark of CNV inference on spatial data,\n",
+    "The ideas are not new. InSituCNV {cite}`Jensen2025` smooths infercnvpy profiles over kNN neighbors on Xenium and CosMx,\n",
+    "SlideCNA {cite}`Zhang2025` pools spatial neighbors on Slide-seq, and STARCH {cite}`Elyanow2021` and CalicoST\n",
+    "{cite}`Ma2024` make neighboring cells prefer the same clone label. In a benchmark of CNV inference on spatial data,\n",
     "{cite:t}`Han2026` found that aggregating cells on imaging platforms improved CNV recovery compared to native\n",
     "single-cell resolution.\n",
     "\n",
@@ -122,8 +122,8 @@
    "id": "2f7e8f19",
    "metadata": {},
    "source": [
-    "We drop cells with fewer than 50 counts and log-normalise to 10,000 counts per cell. The CNV values are log fold\n",
-    "changes, so the normalisation target sets their scale: with the default target (the median of the remaining cells),\n",
+    "We drop cells with fewer than 50 counts and log-normalize to 10,000 counts per cell. The CNV values are log fold\n",
+    "changes, so the normalization target sets their scale: with the default target (the median of the remaining cells),\n",
     "most values are `log1p` of 0 or 1 count, which compresses the CNV values."
    ]
   },
@@ -183,7 +183,7 @@
     "expression data and annotate the clusters with marker gene scores. The Prime 5K panel leaves out many classic,\n",
     "highly expressed markers, e.g. PTPRC, KRT5 and COL1A1, which are missing from the filtered lists printed below. We\n",
     "use T, B and myeloid cell markers for immune cells, endothelial markers, fibroblast markers and squamous epithelial\n",
-    "markers, as far as they are in the panel. Immune cells, endothelial cells and fibroblasts are the reference; the epithelial cells contain the tumour."
+    "markers, as far as they are in the panel. Immune cells, endothelial cells and fibroblasts are the reference; the epithelial cells contain the tumor."
    ]
   },
   {
@@ -506,8 +506,8 @@
     "so we use running windows of 51 genes, every 5 genes, instead of the defaults\n",
     "(100 and 10). We turn off the noise filter (`dynamic_threshold=None`), because it sets small values to 0 and\n",
     "{func}`infercnvpy.tl.smooth_cnv` would then average those zeros. Then we build a spatial graph of the 10 nearest\n",
-    "neighbours of each cell with squidpy and average each cell's CNV profile with its neighbours. By default, the cell\n",
-    "itself counts as one more neighbour."
+    "neighbors of each cell with squidpy and average each cell's CNV profile with its neighbors. By default, the cell\n",
+    "itself counts as one more neighbor."
    ]
   },
   {
@@ -674,11 +674,11 @@
    "id": "e088f7f9",
    "metadata": {},
    "source": [
-    "## Tumour calls\n",
+    "## Tumor calls\n",
     "\n",
-    "We score each cell by the mean absolute value of its smoothed CNV profile and call it tumour when the score is above\n",
+    "We score each cell by the mean absolute value of its smoothed CNV profile and call it tumor when the score is above\n",
     "the 95th percentile of the scores of the reference cells, so the score and the threshold are both per cell. By\n",
-    "construction, 5% of the reference cells are called tumour, and normal cells of types outside the reference can exceed\n",
+    "construction, 5% of the reference cells are called tumor, and normal cells of types outside the reference can exceed\n",
     "the threshold as well."
    ]
   },
@@ -790,8 +790,8 @@
    "id": "8e291530",
    "metadata": {},
    "source": [
-    "Most epithelial cells are called tumour. The fibroblasts, immune and endothelial cells among the tumour calls (share\n",
-    "printed above) are the 5% of reference cells above the threshold. The tumour calls show the loss of 3p, the gain of\n",
+    "Most epithelial cells are called tumor. The fibroblasts, immune and endothelial cells among the tumor calls (share\n",
+    "printed above) are the 5% of reference cells above the threshold. The tumor calls show the loss of 3p, the gain of\n",
     "3q and losses on chromosomes 4, 5 and 9."
    ]
   },
@@ -835,8 +835,8 @@
    "source": [
     "## Spatially coherent regions\n",
     "\n",
-    "{func}`infercnvpy.tl.spatial_regions` makes labels spatially coherent: a cell takes the label of its neighbourhood\n",
-    "when enough neighbours disagree with it (for two labels and the default `beta=2`, more than 75% of them).\n",
+    "{func}`infercnvpy.tl.spatial_regions` makes labels spatially coherent: a cell takes the label of its neighborhood\n",
+    "when enough neighbors disagree with it (for two labels and the default `beta=2`, more than 75% of them).\n",
     "The calls from the smoothed profiles are already spatially coherent, and few cells change."
    ]
   },
@@ -930,9 +930,9 @@
    "source": [
     "The prior matters more for noisy labels. Calls of single cells from the unsmoothed CNV profiles, with the same rule,\n",
     "are scattered over the section (left in the figure below). With `beta=4` (a cell changes its label when more than\n",
-    "62.5% of its neighbours disagree), `tl.spatial_regions` removes about half of the tumour calls (8,788 to 4,710 in the\n",
-    "table below), including 89% of the calls of non-epithelial cells (middle). It doesn't recover tumour cells that the\n",
-    "unsmoothed profiles miss, so the regularised calls from the smoothed profiles (right) include more than twice as many\n",
+    "62.5% of its neighbors disagree), `tl.spatial_regions` removes about half of the tumor calls (8,788 to 4,710 in the\n",
+    "table below), including 89% of the calls of non-epithelial cells (middle). It doesn't recover tumor cells that the\n",
+    "unsmoothed profiles miss, so the regularized calls from the smoothed profiles (right) include more than twice as many\n",
     "epithelial cells."
    ]
   },
@@ -1064,13 +1064,13 @@
    "source": [
     "## Notes\n",
     "\n",
-    "- With 10 neighbours and the default weight of the cell itself, a cell's own profile is 1/11 of its smoothed profile, so calls from\n",
-    "  smoothed profiles are calls of neighbourhoods. Smoothing mixes the profiles of tumour and normal cells at the\n",
-    "  tumour boundary and of immune cells that infiltrate the tumour, over about `n_hops` cells. Use the unsmoothed\n",
+    "- With 10 neighbors and the default weight of the cell itself, a cell's own profile is 1/11 of its smoothed profile, so calls from\n",
+    "  smoothed profiles are calls of neighborhoods. Smoothing mixes the profiles of tumor and normal cells at the\n",
+    "  tumor boundary and of immune cells that infiltrate the tumor, over about `n_hops` cells. Use the unsmoothed\n",
     "  `X_cnv` for per-cell questions in mixed regions.\n",
-    "- The spatial graph defines the neighbourhood. A kNN graph adapts to the cell density; a radius graph\n",
+    "- The spatial graph defines the neighborhood. A kNN graph adapts to the cell density; a radius graph\n",
     "  ({func}`squidpy.gr.spatial_neighbors_radius`) keeps the physical scale fixed. `n_hops=2` in `tl.smooth_cnv`\n",
-    "  widens the neighbourhood further.\n",
+    "  widens the neighborhood further.\n",
     "- `tl.spatial_regions` changes labels only as long as a single change makes them more coherent, so the result stays\n",
     "  close to the input labels and is not a new clustering. Larger `beta` gives smoother regions."
    ]

src/infercnvpy/_numba.py

@@ -1,18 +1,20 @@
 """numba helpers: a parallel kernel with a serial twin, and `n_jobs` as the number of numba threads.
 
-A ``parallel=True`` kernel only runs on the main thread. Called from any other thread (a dask worker, a thread pool),
-the serial twin runs instead: numba's default ``workqueue`` threading layer (the only one in the macOS wheels) aborts
-the process on concurrent use, and a caller that brings its own threads should not get numba threads on top.
+A ``parallel=True`` kernel only runs on the main thread and outside dask tasks (`dask_task`). Anywhere else (a dask
+worker, a thread pool, the main thread of a dask worker process) the serial twin runs instead: numba's default
+``workqueue`` threading layer (the only one in the macOS wheels) aborts the process on concurrent use, and a caller
+that brings its own threads or processes should not get numba threads on top.
 
-All kernels use numba's disk cache, with data files named after their cache key (see `_KeyedCacheFile`).
+All kernels use numba's disk cache, with data files named after their cache key (see `_KeyedCacheFile`), and numpy's
+error model, so that both twins divide by zero alike.
 """
 
 from __future__ import annotations
 
 import hashlib
 import threading
 from contextlib import contextmanager
-from functools import update_wrapper
+from functools import update_wrapper, wraps
 from types import FunctionType
 from typing import TYPE_CHECKING, Any
 
@@ -25,7 +27,7 @@
 else:
     from numba import prange
 
-__all__ = ["n_kernel_chunks", "njit", "njit_parallel", "numba_threads", "prange"]
+__all__ = ["dask_task", "n_kernel_chunks", "njit", "njit_parallel", "numba_threads", "prange", "runs_serially"]
 
 
 class _KeyedCacheFile(IndexDataCacheFile):
@@ -58,18 +60,18 @@ def __init__(self, py_func):
 
 
 def njit(fn: Callable[..., Any], /, **options: Any) -> Any:
-    """`numba.njit` with ``nogil=True`` and the disk cache of `_KeyedCacheFile`."""
-    dispatcher = numba.njit(fn, nogil=True, **options)
-    dispatcher._cache = _KeyedFunctionCache(dispatcher.py_func)  # type: ignore[attr-defined]  # what `cache=True` sets
+    """`numba.njit` with ``nogil=True``, numpy's error model and the disk cache of `_KeyedCacheFile`."""
+    dispatcher = numba.njit(fn, nogil=True, error_model="numpy", **options)
+    if hasattr(dispatcher, "py_func"):  # with NUMBA_DISABLE_JIT=1, `njit` returns `fn` itself
+        dispatcher._cache = _KeyedFunctionCache(dispatcher.py_func)  # type: ignore[attr-defined]  # what `cache=True` sets
     return dispatcher
 
 
 class njit_parallel:
     """Compile a function with `njit` twice, ``parallel=True`` and serial.
 
-    Calling it runs the parallel version on the main thread and the serial one on any other thread, or always the
-    serial one with ``serial=True`` (for code that runs inside dask tasks). The two versions are distinct function
-    objects with distinct names, so they don't share a cache entry.
+    Calling it runs the serial version where `runs_serially` says so, otherwise the parallel one. The two versions are
+    distinct function objects with distinct names, so they don't share a cache entry.
     """
 
     def __init__(self, fn: Callable[..., Any], /):
@@ -78,12 +80,31 @@ def __init__(self, fn: Callable[..., Any], /):
         self.serial = njit(_renamed(fn, "serial"))
         update_wrapper(self, fn)
 
-    def __call__(self, *args, serial: bool = False):
-        return (self.serial if serial or not _on_main_thread() else self.parallel)(*args)
+    def __call__(self, *args):
+        return (self.serial if runs_serially() else self.parallel)(*args)
 
 
-def _on_main_thread() -> bool:
-    return threading.current_thread() is threading.main_thread()
+_scope = threading.local()
+
+
+def runs_serially() -> bool:
+    """Whether kernels called here run serially: off the main thread, or inside a dask task (see `dask_task`)."""
+    return threading.current_thread() is not threading.main_thread() or getattr(_scope, "in_dask_task", False)
+
+
+def dask_task[F: Callable[..., Any]](fn: F) -> F:
+    """Mark `fn` as a dask task: the kernels it calls run serially, whatever thread or process the scheduler uses."""
+
+    @wraps(fn)
+    def task(*args, **kwargs):
+        outer = getattr(_scope, "in_dask_task", False)
+        _scope.in_dask_task = True
+        try:
+            return fn(*args, **kwargs)
+        finally:
+            _scope.in_dask_task = outer
+
+    return task  # type: ignore[return-value]
 
 
 def _renamed(fn: FunctionType, suffix: str) -> FunctionType:
@@ -93,9 +114,9 @@ def _renamed(fn: FunctionType, suffix: str) -> FunctionType:
     return new
 
 
-def n_kernel_chunks(n_rows: int, *, serial: bool = False, per_thread: int = 4) -> int:
+def n_kernel_chunks(n_rows: int, *, per_thread: int = 4) -> int:
     """Number of row chunks for a `prange` kernel: `per_thread` per numba thread, or 1 if it runs serially."""
-    if serial or not _on_main_thread():
+    if runs_serially():
         return 1
     return max(1, min(n_rows, per_thread * numba.get_num_threads()))