diff --git a/src/squidpy/gr/_ppatterns.py b/src/squidpy/gr/_ppatterns.py index a4416a1b2..c242a26d3 100644 --- a/src/squidpy/gr/_ppatterns.py +++ b/src/squidpy/gr/_ppatterns.py @@ -13,11 +13,12 @@ from scanpy import logging as logg from scanpy.metrics import gearys_c, morans_i from scipy import stats -from scipy.sparse import spmatrix +from scipy.sparse import issparse, spmatrix from sklearn.metrics import pairwise_distances from sklearn.preprocessing import normalize from spatialdata import SpatialData from statsmodels.stats.multitest import multipletests +from tqdm.auto import tqdm from squidpy._compat import old_positionals from squidpy._constants._constants import SpatialAutocorr @@ -27,12 +28,11 @@ NDArrayA, RNGLike, SeedLike, - Signal, - SigQueue, deprecated_params, deprecated_randomness_param, - get_n_processes, - parallelize, + get_n_numba_threads, + numba_threads, + thread_map, ) from squidpy._validators import assert_key_in_adata, assert_positive from squidpy.gr._utils import ( @@ -53,6 +53,14 @@ fp = np.float32 bl = nt.boolean +# Permutation entries held at once. Permutations are reused across features, so they are +# materialized rather than redrawn per feature; the cap keeps that buffer off `n_perms`, which +# users raise for FDR resolution. Each block re-pays one scatter per feature, so the cost is set +# by the block size this leaves: measured within a few percent up to 10M cells, where the block is +# still tens of permutations, and ~1.3-1.6x at 100M, where it falls to two. Raise the cap if that +# scale matters. It changes no result, see `test_spatial_autocorr_perm_blocks`. +_PERM_BLOCK_SIZE = 2**28 + @d.dedent @inject_docs(key=Key.obsp.spatial_conn(), sp=SpatialAutocorr) @@ -74,6 +82,7 @@ "show_progress_bar", ) @deprecated_randomness_param +@deprecated_params({"backend": "1.10.0"}) def spatial_autocorr( adata: AnnData | SpatialData, *, @@ -90,7 +99,6 @@ def spatial_autocorr( use_raw: bool = False, copy: bool = False, n_jobs: int | None = None, - backend: str = "loky", show_progress_bar: bool = True, table_key: str | None = None, ) -> pd.DataFrame | None: @@ -110,6 +118,10 @@ def spatial_autocorr( %(rng_versionchanged)s + .. versionchanged:: 1.8.4 + Permutations run on a thread pool, and ``n_jobs = None`` now uses all ``NUMBA_NUM_THREADS`` + threads instead of one process. Pass ``n_jobs = 1`` for the old serial default. + Parameters ---------- %(adata)s @@ -149,7 +161,8 @@ def spatial_autocorr( Which attribute of :class:`~anndata.AnnData` to access. See ``genes`` parameter for more information. %(rng)s %(copy)s - %(parallelize)s + %(n_jobs_threads)s + %(show_progress_bar)s Returns ------- @@ -184,6 +197,10 @@ def extract_X(adata: AnnData, genes: str | Sequence[str] | None) -> tuple[NDArra genes = [genes] if not use_raw: + # full var_names in order: `adata[:, genes]` would copy X to hand back the same + # columns it was given. + if len(genes) == adata.n_vars and np.array_equal(np.asarray(genes), adata.var_names.values): + return (adata.X if layer is None else adata.layers[layer]).T, genes subset = adata[:, genes] return (subset.X if layer is None else subset.layers[layer]).T, genes if adata.raw is None: @@ -232,31 +249,24 @@ def extract_obsm(adata: AnnData, ixs: int | Sequence[int] | None) -> tuple[NDArr else: raise NotImplementedError(f"Mode `{mode}` is not yet implemented.") - g = adata.obsp[connectivity_key].copy() + g = adata.obsp[connectivity_key].tocsr(copy=True) if transformation: # row-normalize normalize(g, norm="l1", axis=1, copy=False) - score = params["func"](g, vals) # type: ignore - - n_jobs = get_n_processes(n_jobs) - start = logg.info(f"Calculating {mode}'s statistic for `{n_perms}` permutations using `{n_jobs}` core(s)") + n_jobs = get_n_numba_threads(n_jobs) + start = logg.info(f"Calculating {mode}'s statistic for `{n_perms}` permutations using `{n_jobs}` thread(s)") if n_perms is not None: assert_positive(n_perms, name="n_perms") - perms = list(np.arange(n_perms)) - rngs = np.random.default_rng(rng).spawn(n_perms) - - score_perms = parallelize( - _score_helper, - collection=perms, - extractor=np.concatenate, - n_jobs=n_jobs, - backend=backend, - show_progress_bar=show_progress_bar, - )(mode=mode, g=g, vals=vals, rngs=rngs) + # the observed score comes from the same kernel as the permuted ones, so exact ties + # in the tally below stay exact + score, score_perms = _score_perms( + g, vals, mode=mode, n_perms=n_perms, rng=rng, n_jobs=n_jobs, show_progress_bar=show_progress_bar + ) else: - score_perms = None + score, score_perms = params["func"](g, vals), None # type: ignore - with np.errstate(divide="ignore"): + # a degenerate feature has `var_sim == 0`, so its z score is an invalid 0/0, not a divide + with np.errstate(divide="ignore", invalid="ignore"): pval_results = _p_value_calc(score, score_perms, g, params) data_dict: dict[str, Any] = {str(params["stat"]): score, **pval_results} @@ -264,8 +274,14 @@ def extract_obsm(adata: AnnData, ixs: int | Sequence[int] | None) -> tuple[NDArr if corr_method is not None: for pv in filter(lambda x: "pval" in x, df.columns): - _, pvals_adj, _, _ = multipletests(df[pv].values, alpha=0.05, method=corr_method) - df[f"{pv}_{corr_method}"] = pvals_adj + pvals = df[pv].to_numpy(dtype=np.float64) + # a degenerate feature has no p-value to correct, and `multipletests` would spread its + # NaN over every other gene, so correct the defined ones among themselves + defined = np.isfinite(pvals) + adj = np.full(pvals.shape, np.nan) + if defined.any(): + _, adj[defined], _, _ = multipletests(pvals[defined], alpha=0.05, method=corr_method) + df[f"{pv}_{corr_method}"] = adj df.sort_values(by=params["stat"], ascending=params["ascending"], inplace=True) @@ -278,30 +294,234 @@ def extract_obsm(adata: AnnData, ixs: int | Sequence[int] | None) -> tuple[NDArr _save_data(adata, attr="uns", key=mode_str + stat_str, data=df, time=start) -def _score_helper( - perms: Sequence[int], - mode: SpatialAutocorr, - g: spmatrix, - vals: NDArrayA, - rngs: Sequence[np.random.Generator], - *, - queue: SigQueue | None = None, -) -> pd.DataFrame: - score_perms = np.empty((len(perms), vals.shape[0])) - func = morans_i if mode == SpatialAutocorr.MORAN else gearys_c +@njit(parallel=False, nogil=True, cache=True) +def _autocorr_perms( # noqa: PLR0917, numba requires positional arguments + tptr: NDArrayA, + tind: NDArrayA, + tdat: NDArrayA, + w_sum: NDArrayA, + col_sum: NDArrayA, + nz: NDArrayA, + xv: NDArrayA, + w: float, + perms: NDArrayA, + moran: bool, +) -> NDArrayA: + """Permutation scores of one feature over the row-permuted graph, one per row of ``perms``. - for i, p in enumerate(perms): - rng = rngs[p] - idx_shuffle = rng.permutation(g.shape[0]) - score_perms[i, :] = func(g[idx_shuffle, :], vals) + Both statistics sum, over the rows of ``g[perm, :]``, a per-row quantity that depends on the + permutation only through *which* row is visited, so those quantities are accumulated once per + row of the unpermuted graph rather than re-derived for every permutation. - if queue is not None: - queue.put(Signal.UPDATE) + ``nz`` lists the cells where this feature is non-zero. Both statistics split into a term that + only touches those cells and a term that is permutation-invariant (``sum(perm)`` over a bijection + is ``sum``), so the per-permutation loop costs ``O(len(nz))`` rather than ``O(n_cells)``. Pass + ``arange(n)`` for a dense feature. - if queue is not None: - queue.put(Signal.FINISH) + ``xv`` must already be :class:`numpy.float64`, and ``tdat`` the float64 graph weights, + matching the casts :mod:`scanpy` applies before its own kernels. - return score_perms + Compiled serial on purpose: `_score_perms` parallelizes across features instead, which measured + faster than an in-kernel ``prange`` even with fewer features than workers. ``nogil`` is what + lets that pool enter it concurrently, which a ``parallel=True`` kernel cannot survive under + numba's default threading layer. + """ + n_perms, n = perms.shape + out = np.empty(n_perms, dtype=np.float64) + + # Independence from ``n_jobs`` comes from `_score_perms`: permutation ``p`` is always drawn + # from ``rngs[p]``, and each feature is owned by one worker. + n_nz = nz.shape[0] + x_bar = 0.0 + for t in range(n_nz): + x_bar += xv[t] + x_bar /= n + + # sum_i (x_i - x_bar)^2: Moran's denominator, and Geary's. + css = 0.0 + for t in range(n_nz): # looped rather than ((xv - x_bar) ** 2).sum() to skip the temporary + d = xv[t] - x_bar + css += d * d + css += (n - n_nz) * x_bar * x_bar # every zero entry contributes x_bar^2 + + if moran: + # I = n / W * sum_ij w_ij z_i z_j / sum_i z_i^2. The inner row sum is the spatial lag, + # which the permutation only reindexes. + lag = np.zeros(n, dtype=np.float64) + # lag = g @ z = (g @ x) - x_bar * rowsum(g). Only columns where x != 0 contribute, so + # this scatters over the feature's non-zeros rather than sweeping the whole graph. + for t in range(n_nz): + j = nz[t] + xj = xv[t] + for e in range(tptr[j], tptr[j + 1]): + lag[tind[e]] += tdat[e] * xj + for k in range(n): + lag[k] -= x_bar * w_sum[k] + # sum_i lag[perm[i]] * z[i] = sum_{x[i]!=0} lag[perm[i]] * x[i] - x_bar * sum_k lag[k]. + # The second term is permutation-invariant, so it is hoisted out of the loop. + # sum_k lag[k] = sum_j x[j] * colsum(g)[j] - x_bar * W: O(nnz), not O(n_cells). + lag_tot = 0.0 + for t in range(n_nz): + lag_tot += xv[t] * col_sum[nz[t]] + lag_tot -= x_bar * w + for p in range(n_perms): + inum = 0.0 + for t in range(n_nz): + inum += lag[perms[p, nz[t]]] * xv[t] + out[p] = n / w * (inum - x_bar * lag_tot) / css + else: + # C = (n - 1) * sum_ij w_ij (x_i - x_j)^2 / (2 W sum_i (x_i - x_bar)^2). Expanding the + # square splits each row into (sum w, sum w x_j, sum w x_j^2), all permutation-independent. + wx = np.zeros(n, dtype=np.float64) + wx2 = np.zeros(n, dtype=np.float64) + for t in range(n_nz): + j = nz[t] + xj = xv[t] + xj2 = xj * xj + for e in range(tptr[j], tptr[j + 1]): + k = tind[e] + weight = tdat[e] + wx[k] += weight * xj + wx2[k] += weight * xj2 + denom = 2.0 * w * css + # sum_i wx2[perm[i]] is permutation-invariant; the other two terms vanish where x[i] == 0. + wx2_tot = 0.0 + for t in range(n_nz): + wx2_tot += xv[t] * xv[t] * col_sum[nz[t]] + for p in range(n_perms): + total = 0.0 + for t in range(n_nz): + k = perms[p, nz[t]] + xi = xv[t] + total += xi * xi * w_sum[k] - 2.0 * xi * wx[k] + out[p] = (n - 1) * (total + wx2_tot) / denom + return out + + +def _chunks(n: int, workers: int) -> list[range]: + """Contiguous index runs over ``n``, ~4 per worker so the pool can balance uneven items.""" + step = max(1, -(-n // (workers * 4))) + return [range(s, min(s + step, n)) for s in range(0, n, step)] + + +def _score_perms( + g: spmatrix, + vals: NDArrayA | spmatrix, + *, + mode: SpatialAutocorr, + n_perms: int, + rng: SeedLike | RNGLike | None, + n_jobs: int, + show_progress_bar: bool, +) -> tuple[NDArrayA, dict[str, Any]]: + """Observed scores ``(n_features,)`` and the permutation accumulators `_p_value_calc` needs. + + The permutation scores are only ever reduced along the permutation axis, so the + ``(n_perms, n_features)`` matrix is never materialized: a tally, plus sums of the scores + shifted by the observed one, carry the same information in ``O(n_features)``. + """ + n_cells = g.shape[0] + # Match the casts scanpy applies to its own inputs, so the kernel sees the same numbers. + g = g.astype(np.float64, copy=False) + w = g.data.sum() + rngs = np.random.default_rng(rng).spawn(n_perms) + # Blocks of at most `_PERM_BLOCK_SIZE` entries; every feature is revisited once per block, so + # the per-block scatter is re-paid. More than one block only once n_perms * n_cells exceeds it. + block = int(np.clip(_PERM_BLOCK_SIZE // max(n_cells, 1), 1, n_perms)) + buffer = np.empty((block, n_cells), dtype=np.int32) + identity = np.arange(n_cells, dtype=np.int32)[None] + all_cells = np.arange(n_cells, dtype=np.int32) # `nz` for a dense feature + w_sum = np.asarray(g.sum(axis=1)).ravel() # feature-independent: compute once, not per call + col_sum = np.asarray(g.sum(axis=0)).ravel() # ditto; turns the O(n) totals into O(nnz) + + moran = mode == SpatialAutocorr.MORAN + sparse_vals = issparse(vals) + if sparse_vals: + # ``vals`` arrives as ``X.T``, i.e. CSC, whose row slicing is O(nnz) rather than O(nnz_row); + # one conversion here makes the per-feature extraction below ~180x cheaper. + vals = vals.tocsr() + n_features = vals.shape[0] + # The kernel walks g by column, so it needs the transpose. One O(nnz) pass, and it removes + # the per-feature dense vector the row-wise form needed. + gt = g.T.tocsr() + tptr, tind, tdat = gt.indptr, gt.indices, gt.data + # Constant features have a zero denominator; scanpy drops them and reports `nan`, so seed with it. + score = np.full(n_features, np.nan, dtype=np.float64) + count_ge = np.zeros(n_features, dtype=np.int64) # integer tally of `sims >= score`, so exact + # Mean and variance accumulate as sums of `sims - score`. Both are plain sums, so blocks merge + # with `+=` and no Welford state is carried. `s1/n` is the observed score's distance from the + # permutation mean, so `s2/n - (s1/n)**2` loses about eps*z**2: ~1e-13 at z=27, against ~1e-16 + # for a two-pass variance. Worth it for the simpler merge, but it is a real difference, and it + # makes the variance exactly 0 when every permutation scores alike (see `_p_value_calc`). + s1 = np.zeros(n_features, dtype=np.float64) + s2 = np.zeros(n_features, dtype=np.float64) + n_blocks = -(-n_perms // block) + # Features are independent and the kernel is nogil, so the parallelism is here rather than + # inside the kernel. That held even with fewer features than workers, where an in-kernel + # `prange` over permutations still measured slower. + pool_workers = max(1, min(n_jobs, n_features)) + + def run_feature(m: int, perms: NDArrayA, first: bool) -> None: + if sparse_vals: + # `vals` is CSR, so this feature's non-zero cells and values are already contiguous. + # The kernel reads them directly and never densifies the feature. + lo_m, hi_m = vals.indptr[m], vals.indptr[m + 1] + nz = vals.indices[lo_m:hi_m].astype(np.int32, copy=False) + xv = np.ascontiguousarray(vals.data[lo_m:hi_m], dtype=np.float64) + else: + nz = all_cells + # both statistics are invariant to a shift, and centring here keeps the kernel's + # `g @ x - x_bar * w_sum` from cancelling on a feature whose mean dwarfs its spread. + # A sparse feature cannot reach that ratio, so it pays nothing for this. + xv = np.asarray(vals[m], dtype=np.float64) + xv = xv - xv.mean() + if len(nz) == 0 or (xv.min() == xv.max() and (len(nz) == n_cells or xv[0] == 0.0)): + return # constant feature: scanpy drops it and reports `nan` + args = (tptr, tind, tdat, w_sum, col_sum, nz, xv, w) + if first: + score[m] = _autocorr_perms(*args, identity, moran)[0] + sims = _autocorr_perms(*args, perms, moran) + count_ge[m] += int((sims >= score[m]).sum()) + shifted = sims - score[m] + s1[m] += shifted.sum() + s2[m] += (shifted * shifted).sum() + + with ( + numba_threads(1), # the kernel is serial; the pool below owns the parallelism + tqdm(total=n_blocks * n_features, unit="feature", disable=not show_progress_bar) as pbar, + ): + for lo in range(0, n_perms, block): + perms = buffer[: min(block, n_perms - lo)] + first = lo == 0 + + def fill(chunk: range, perms: NDArrayA = perms, lo: int = lo) -> None: + for i in chunk: + perms[i] = rngs[lo + i].permutation(n_cells) + + # Drained before the feature map below starts, so no worker reads `perms` while + # another is still filling it. + thread_map(fill, _chunks(len(perms), pool_workers), n_jobs=pool_workers) + + # Each worker takes a contiguous run of features rather than one at a time; a single + # feature is short enough that per-item dispatch cost about as much as the work. Each + # `m` belongs to one worker, so the writes in `run_feature` never collide. The defaults + # bind this block's values rather than the loop's last ones. + def run_chunk(chunk: range, perms: NDArrayA = perms, first: bool = first) -> int: + for m in chunk: + run_feature(m, perms, first) + return len(chunk) + + for done in thread_map(run_chunk, _chunks(n_features, pool_workers), n_jobs=pool_workers): + pbar.update(done) + # Constant features never ran, so they keep scanpy's `nan` rather than a zero mean/variance. + ran = ~np.isnan(score) + mean_shift = np.where(ran, s1 / n_perms, np.nan) + return score, { + "n_perms": n_perms, + "count_ge": count_ge, + "mean": score + mean_shift, + "var": np.where(ran, s2 / n_perms - mean_shift**2, np.nan), + } @njit(parallel=True, fastmath=True, cache=True) @@ -469,7 +689,7 @@ def _find_min_max(spatial: NDArrayA) -> tuple[float, float]: def _p_value_calc( score: NDArrayA, - sims: NDArrayA | None, + sims: dict[str, Any] | None, weights: spmatrix | NDArrayA, params: dict[str, Any], ) -> dict[str, Any]: @@ -481,7 +701,10 @@ def _p_value_calc( score (n_features,). sims - (n_simulations, n_features). + Permutation accumulators from `_score_perms`: ``n_perms`` and, per feature, the + ``count_ge`` tally plus the ``mean``/``var`` of the permuted scores. + weights + The spatial connectivity graph, used for the analytic (normality) p-value. params Object to store relevant function parameters. @@ -500,24 +723,25 @@ def _p_value_calc( if sims is None: return results - n_perms = sims.shape[0] - large_perm = (sims >= score).sum(axis=0) + n_perms = sims["n_perms"] + large_perm = sims["count_ge"].copy() # copy: the fold below writes in place # subtract total perm for negative values large_perm[(n_perms - large_perm) < large_perm] = n_perms - large_perm[(n_perms - large_perm) < large_perm] # get p-value based on permutation p_sim: NDArrayA = (large_perm + 1) / (n_perms + 1) # get p-value based on standard normal approximation from permutations - e_score_sim = sims.sum(axis=0) / n_perms - se_score_sim = sims.std(axis=0) + e_score_sim = sims["mean"] + var_sim = sims["var"] + se_score_sim = np.sqrt(var_sim) z_sim = (score - e_score_sim) / se_score_sim - p_z_sim = np.empty(z_sim.shape) + # NaN where a feature's permutations are all identical, so `var_sim` is exactly 0; the masks + # below skip those and `np.empty` would leave them uninitialized. + p_z_sim = np.full(z_sim.shape, np.nan) p_z_sim[z_sim > 0] = 1 - stats.norm.cdf(z_sim[z_sim > 0]) p_z_sim[z_sim <= 0] = stats.norm.cdf(z_sim[z_sim <= 0]) - var_sim = np.var(sims, axis=0) - results["pval_z_sim"] = p_z_sim results["pval_sim"] = p_sim results["var_sim"] = var_sim @@ -555,7 +779,7 @@ def _analytic_pval(score: NDArrayA, g: spmatrix | NDArrayA, params: dict[str, An seScore_norm = Vscore_norm ** (1 / 2.0) z_norm = (score - params["expected"]) / seScore_norm - p_norm = np.empty(score.shape) + p_norm = np.full(score.shape, np.nan) # constant features have a NaN score and match neither mask p_norm[z_norm > 0] = 1 - stats.norm.cdf(z_norm[z_norm > 0]) p_norm[z_norm <= 0] = stats.norm.cdf(z_norm[z_norm <= 0]) diff --git a/tests/graph/test_ppatterns.py b/tests/graph/test_ppatterns.py index aac5f1b85..d7c8caeef 100644 --- a/tests/graph/test_ppatterns.py +++ b/tests/graph/test_ppatterns.py @@ -9,7 +9,7 @@ from squidpy._constants._pkg_constants import Key from squidpy.gr import co_occurrence, spatial_autocorr -from squidpy.gr._ppatterns import _find_min_max +from squidpy.gr._ppatterns import _autocorr_perms, _find_min_max, _score_perms MORAN_K = "moranI" GEARY_C = "gearyC" @@ -37,22 +37,17 @@ def test_spatial_autocorr_seq_par(dummy_adata: AnnData, mode: str): assert "pval_norm_fdr_bh" in dummy_adata.uns[UNS_KEY] assert dummy_adata.uns[UNS_KEY].columns.shape == (4,) assert df.columns.shape == (9,) - # test pval_norm same - # will need to increase the tolerance because numba parallel computations might not be exactly the same - # these pval_norms don't use the seed anyway so the difference is not due to the seed - # see https://github.com/scverse/squidpy/issues/1030 for more details np.testing.assert_allclose(df["pval_norm"].values, df_parallel["pval_norm"].values, atol=1e-12) # test highly variable assert dummy_adata.uns[UNS_KEY].shape != df.shape # assert idx are sorted and contain same elements assert not np.array_equal(idx_df, idx_adata) np.testing.assert_array_equal(sorted(idx_df), sorted(idx_adata)) - # check parallel gives same results - # each permutation now gets its own seed (spawned from a SeedSequence), so the - # simulated p-values no longer depend on how the permutations are split across jobs - np.testing.assert_allclose(df["pval_sim"].values, df_parallel["pval_sim"].values, atol=1e-12) - np.testing.assert_allclose(df["pval_z_sim"].values, df_parallel["pval_z_sim"].values, atol=1e-12) - np.testing.assert_allclose(df["var_sim"].values, df_parallel["var_sim"].values, atol=1e-12) + # each permutation draws from its own spawned generator, so the simulated columns do not + # depend on how the features were split across workers + df_parallel = df_parallel.loc[df.index] # align in case the stat-based sort ties differently + for col in ("pval_sim", "pval_z_sim", "var_sim"): + np.testing.assert_allclose(df[col].values, df_parallel[col].values, atol=1e-12) @pytest.mark.parametrize("mode", ["moran", "geary"]) @@ -77,12 +72,6 @@ def test_spatial_autocorr_reproducibility(dummy_adata: AnnData, n_jobs: int, mod # assert fdr correction in adata.uns assert "pval_sim_fdr_bh" in df_1 assert "pval_norm_fdr_bh" in dummy_adata.uns[UNS_KEY] - # test pval_norm same - # will need to increase the tolerance because numba parallel computations might not be exactly the same - # see https://github.com/scverse/squidpy/issues/1030 for more details about the tolerance - # these pval_norms don't use the seed anyway so the difference is not due to the seed - np.testing.assert_allclose(df_1["pval_norm"].values, df_2["pval_norm"].values, atol=1e-12) - np.testing.assert_allclose(df_1["var_norm"].values, df_2["var_norm"].values, atol=1e-12) assert dummy_adata.uns[UNS_KEY].columns.shape == (4,) assert df_2.columns.shape == (9,) # test highly variable @@ -95,16 +84,168 @@ def test_spatial_autocorr_reproducibility(dummy_adata: AnnData, n_jobs: int, mod @pytest.mark.parametrize("mode", ["moran", "geary"]) -def test_spatial_autocorr_n_jobs_invariance(dummy_adata: AnnData, mode: str): - """The number of workers must not change the permutation-based results (seed spawned per permutation).""" +def test_spatial_autocorr_degenerate_feature_is_nan(mode: str): + """A gene whose permutations all score alike has no z-test, and must not poison the FDR column. + + Its variance is exactly zero, so the z-score is undefined. Filling `pval_z_sim` with + `np.empty` left those entries at whatever memory held, which reads as a significant + p-value, and a NaN handed to `multipletests` spreads over every other gene. + """ + import scipy.sparse as sps + + from squidpy.gr import spatial_neighbors_knn + + rng = np.random.default_rng(0) + n = 400 + X = np.zeros((n, 3), dtype=np.float32) + X[5, 0] = 4.0 # expressed in a single cell: every permutation gives the same score + X[:, 1] = rng.poisson(2, n) + X[7, 2] = 1.0 + adata = AnnData(sps.csr_matrix(X)) + adata.var_names = ["solo", "normal", "solo2"] + adata.obsm["spatial"] = rng.random((n, 2)) + spatial_neighbors_knn(adata, n_neighs=6) + + df = spatial_autocorr(adata, mode=mode, n_perms=13, rng=3, copy=True, show_progress_bar=False) + degenerate, real = ["solo", "solo2"], "normal" + + assert (df.loc[degenerate, "var_sim"] == 0.0).all() + assert df.loc[degenerate, "pval_z_sim"].isna().all() + assert df.loc[degenerate, "pval_z_sim_fdr_bh"].isna().all() + # the gene with a defined z-test keeps one, and its correction is unaffected + assert np.isfinite(df.loc[real, "pval_z_sim"]) + assert np.isfinite(df.loc[real, "pval_z_sim_fdr_bh"]) + # the permutation p-value is a tally, so it stays defined for every gene + assert df["pval_sim"].notna().all() + + +@pytest.mark.parametrize(("mode", "stat"), [("moran", "I"), ("geary", "C")]) +def test_spatial_autocorr_offset_feature_precision(mode: str, stat: str): + """Both statistics ignore a constant shift, so a large offset must not move the score. + + The kernel rebuilds the centred quantities by subtraction, e.g. `g @ x - x_bar * w_sum`, + which cancels away roughly `(mean / sd) ** 2` digits. Dense features are centred before + they reach it; without that, an offset of 1e6 moved Moran's I by several percent. + """ + from squidpy.gr import spatial_neighbors_knn + + rng = np.random.default_rng(1) + n = 400 + base = rng.standard_normal(n) + adata = AnnData(np.stack([base, base + 1e6], axis=1).astype(np.float64)) + adata.var_names = ["plain", "offset"] + adata.obsm["spatial"] = rng.random((n, 2)) + spatial_neighbors_knn(adata, n_neighs=6) + + # `n_perms` routes the observed score through the kernel rather than scanpy + df = spatial_autocorr(adata, mode=mode, n_perms=1, rng=0, copy=True, show_progress_bar=False) + # 1e-8 is the floor for recovering a unit-scale value from one offset by 1e6, not slack: + # without the centring the difference was 5e-2 + np.testing.assert_allclose(df.loc["offset", stat], df.loc["plain", stat], rtol=1e-8) + + +def test_spatial_autocorr_ties_match_scanpy(): + """A gene with one non-zero count ties most permutations exactly; ties must count as in scanpy.""" + import scipy.sparse as sps + from scanpy.metrics import gearys_c + from sklearn.preprocessing import normalize + + from squidpy.gr import spatial_neighbors_knn + + rng = np.random.default_rng(0) + X = np.zeros((300, 2), dtype=np.float32) + X[0, 0] = 5 + X[:, 1] = rng.poisson(2, 300) + adata = AnnData(sps.csr_matrix(X)) + adata.obsm["spatial"] = rng.random((300, 2)) + spatial_neighbors_knn(adata, n_neighs=6) + df = spatial_autocorr(adata, mode="geary", n_perms=50, rng=0, copy=True).loc[adata.var_names] + + # the pre-numba algorithm: scanpy on the row-permuted graph, one spawned generator per permutation + g = normalize(adata.obsp["spatial_connectivities"], norm="l1", axis=1) + vals = adata.X.T + obs = gearys_c(g, vals) + sims = np.stack([gearys_c(g[gen.permutation(300), :], vals) for gen in np.random.default_rng(0).spawn(50)]) + large = (sims >= obs).sum(axis=0) + large = np.minimum(large, 50 - large) + assert large[0] > 0 # the scenario really has ties + np.testing.assert_array_equal(df["pval_sim"].values, (large + 1) / 51) + + +@pytest.mark.parametrize("mode", ["moran", "geary"]) +def test_spatial_autocorr_perm_blocks(dummy_adata: AnnData, mode: str, monkeypatch): + """Drawing the permutations block by block must not change the result.""" + import squidpy.gr._ppatterns as ppatterns + + kw = {"mode": mode, "copy": True, "rng": 42, "n_perms": 50} + expected = spatial_autocorr(dummy_adata, **kw) + monkeypatch.setattr(ppatterns, "_PERM_BLOCK_SIZE", 7 * dummy_adata.n_obs) # 8 blocks, the last one short + assert_frame_equal(spatial_autocorr(dummy_adata, **kw), expected) + + +def test_spatial_autocorr_full_gene_list_reordered(dummy_adata: AnnData): + """A full-length but reordered `genes` must not take the identity fast path in `extract_X`. + + `extract_X` skips `adata[:, genes]` when the selection is every gene in `var_names` order. + Dropping the order check from that guard leaves a list of the right length taking the fast + path, which returns `X` in var order while labelling the rows in the caller's order. + """ + genes = list(dummy_adata.var_names) + kw = {"mode": "geary", "n_perms": 20, "rng": 0, "copy": True, "show_progress_bar": False} + fwd = spatial_autocorr(dummy_adata, genes=genes, **kw) + rev = spatial_autocorr(dummy_adata, genes=genes[::-1], **kw) + + # a gene's statistic cannot depend on the order the caller listed the genes in + assert set(fwd.index) == set(rev.index) + np.testing.assert_allclose(fwd["C"], rev.loc[fwd.index, "C"], rtol=1e-12) + + +@pytest.mark.parametrize("mode", ["moran", "geary"]) +def test_spatial_autocorr_csc_connectivities(dummy_adata: AnnData, mode: str): + """A CSC graph must give the CSR result: `normalize(axis=1, copy=False)` leaves CSC untouched.""" + key = Key.obsp.spatial_conn() kw = {"mode": mode, "copy": True, "rng": 42, "n_perms": 50} - df_serial = spatial_autocorr(dummy_adata, n_jobs=1, **kw) - df_parallel = spatial_autocorr(dummy_adata, n_jobs=2, **kw) + csc = dummy_adata.copy() + csc.obsp[key] = csc.obsp[key].tocsc() + before = csc.obsp[key].copy() - # align on the gene index in case the stat-based sort order ties differently - df_parallel = df_parallel.loc[df_serial.index] - for col in ["pval_sim", "pval_z_sim", "var_sim"]: - np.testing.assert_allclose(df_serial[col].values, df_parallel[col].values, atol=1e-12) + assert_frame_equal(spatial_autocorr(csc, **kw), spatial_autocorr(dummy_adata, **kw)) + # row-normalization must not leak back into the caller's graph + np.testing.assert_array_equal(csc.obsp[key].toarray(), before.toarray()) + + +def test_spatial_autocorr_v183_positional_backend(dummy_adata: AnnData): + """A v1.8.3 positional call through ``backend`` binds every value and warns about ``backend``. + + Delete this together with the ``backend`` shim. The positional signature it pins only exists + because `@deprecated_params` keeps accepting `backend` in its v1.8.3 slot; once that is dropped + for 1.10.0 there is no positional form left to protect and the call below starts raising. + """ + kw = {"mode": "moran", "n_perms": 20, "rng": 0, "copy": True, "n_jobs": 1, "show_progress_bar": False} + expected = spatial_autocorr(dummy_adata, **kw) + args = ( + "spatial_connectivities", + None, + "moran", + True, + 20, + False, + "fdr_bh", + "X", + None, + 0, + False, + True, + 1, + "loky", + False, + ) + with pytest.warns(FutureWarning) as record: + got = spatial_autocorr(dummy_adata, *args) + messages = [str(w.message) for w in record] + assert any("`backend`" in m for m in messages) + assert any("`seed`" in m for m in messages) + assert_frame_equal(got, expected) @pytest.mark.parametrize("mode", ["moran", "geary"]) @@ -226,3 +367,98 @@ def test_use_raw(dummy_adata: AnnData): df = spatial_autocorr(dummy_adata, use_raw=True, copy=True) np.testing.assert_equal(sorted(df.index), sorted(var_names)) + + +@pytest.mark.parametrize("mode", ["moran", "geary"]) +def test_score_perms_matches_scanpy_per_permutation(mode: str): + """The kernel reuses per-row sums across permutations; pin it to the naive scanpy reference. + + ``_score_perms`` accumulates each row of the graph once and reindexes those sums per + permutation. The reference below is the definition it replaced: build ``g[perm, :]`` and hand it + to :mod:`scanpy` for every permutation. Agreement has to hold to floating-point reordering only. + """ + from scanpy.metrics import gearys_c, morans_i + from scipy.sparse import csr_matrix + from sklearn.neighbors import kneighbors_graph + from sklearn.preprocessing import normalize + + from squidpy._constants._constants import SpatialAutocorr + + n, n_genes, n_perms = 300, 4, 6 + rng = np.random.default_rng(0) + g = csr_matrix(kneighbors_graph(rng.random((n, 2)), 5, mode="connectivity")) + normalize(g, norm="l1", axis=1, copy=False) + vals = rng.random((n_genes, n), dtype=np.float32) + + autocorr = SpatialAutocorr(mode) + moran = autocorr == SpatialAutocorr.MORAN + observed, got = _score_perms(g, vals, mode=autocorr, n_perms=n_perms, rng=0, n_jobs=1, show_progress_bar=False) + + func = morans_i if autocorr == SpatialAutocorr.MORAN else gearys_c + expected = np.stack([func(g[gen.permutation(n), :], vals) for gen in np.random.default_rng(0).spawn(n_perms)]) + np.testing.assert_allclose(observed, func(g, vals), rtol=1e-9) + + # `_score_perms` only keeps reductions of the permutation scores, so the per-permutation + # comparison happens one level down, against the kernel that still returns them. + gg = g.astype(np.float64, copy=False) + gt = gg.T.tocsr() + perms = np.stack([gen.permutation(n).astype(np.int32) for gen in np.random.default_rng(0).spawn(n_perms)]) + sims = np.stack( + [ + _autocorr_perms( + gt.indptr, + gt.indices, + gt.data, + np.asarray(gg.sum(axis=1)).ravel(), + np.asarray(gg.sum(axis=0)).ravel(), + np.arange(n, dtype=np.int32), + np.ascontiguousarray(vals[m], np.float64), + gg.data.sum(), + perms, + moran, + ) + for m in range(n_genes) + ] + ).T + assert sims.shape == (n_perms, n_genes) + np.testing.assert_allclose(sims, expected, rtol=1e-9) + + # and the accumulators must be those same scores, reduced + assert got["n_perms"] == n_perms + np.testing.assert_array_equal(got["count_ge"], (sims >= observed).sum(axis=0)) + np.testing.assert_allclose(got["mean"], expected.mean(axis=0), rtol=1e-9) + np.testing.assert_allclose(got["var"], expected.var(axis=0), rtol=1e-9) + + +@pytest.mark.parametrize("mode", ["moran", "geary"]) +def test_score_perms_thread_invariant(mode: str): + """Permutation scores must not depend on `n_jobs`, which sizes the thread pool.""" + from scipy.sparse import csr_matrix + from sklearn.neighbors import kneighbors_graph + from sklearn.preprocessing import normalize + + from squidpy._constants._constants import SpatialAutocorr + + n, n_perms = 300, 8 + rng = np.random.default_rng(0) + g = csr_matrix(kneighbors_graph(rng.random((n, 2)), 5, mode="connectivity")) + normalize(g, norm="l1", axis=1, copy=False) + vals = rng.random((3, n), dtype=np.float32) + + autocorr = SpatialAutocorr(mode) + serial = _score_perms(g, vals, mode=autocorr, n_perms=n_perms, rng=0, n_jobs=1, show_progress_bar=False) + threaded = _score_perms(g, vals, mode=autocorr, n_perms=n_perms, rng=0, n_jobs=4, show_progress_bar=False) + np.testing.assert_array_equal(serial[0], threaded[0]) + assert serial[1].keys() == threaded[1].keys() + for key in serial[1]: + np.testing.assert_array_equal(serial[1][key], threaded[1][key], err_msg=key) + + +def test_spatial_autocorr_backend_deprecated(dummy_adata: AnnData): + """``backend`` no longer selects a process pool; it warns and is ignored.""" + with pytest.warns(FutureWarning, match=r"`backend`.*deprecated"): + with_backend = spatial_autocorr( + dummy_adata, copy=True, n_perms=10, rng=0, backend="loky", show_progress_bar=False + ) + without = spatial_autocorr(dummy_adata, copy=True, n_perms=10, rng=0, show_progress_bar=False) + assert_frame_equal(with_backend, without)