From 4cc54a411da5de368bb666c3296c69982b70b51a Mon Sep 17 00:00:00 2001
From: noamteyssier <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 14 Jul 2026 10:04:41 -0700
Subject: [PATCH 01/12] fix: better handling of sched_getaffinity with
 conditional availability

---
 src/pdex/_utils.py | 10 ++++++----
 1 file changed, 6 insertions(+), 4 deletions(-)

diff --git a/src/pdex/_utils.py b/src/pdex/_utils.py
index dbb33e4..6a185bc 100644
--- a/src/pdex/_utils.py
+++ b/src/pdex/_utils.py
@@ -20,10 +20,12 @@ def _available_cpus() -> int:
     CPU even when affinity or a cgroup restricts the process to far fewer, causing
     ``numba.set_num_threads`` to raise.
     """
-    try:
-        return len(os.sched_getaffinity(0))
-    except AttributeError:
-        return mp.cpu_count()
+    if hasattr(os, "sched_getaffinity"):
+        try:
+            return len(os.sched_getaffinity(0))
+        except AttributeError:
+            pass
+    return mp.cpu_count()
 
 
 def set_numba_threadpool(threads: int = 0):

From 9ad0e85886d558f1eb20213a55f92dac648791f8 Mon Sep 17 00:00:00 2001
From: noamteyssier <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 14 Jul 2026 10:56:57 -0700
Subject: [PATCH 02/12] refactor: improve path for one-v-rest

---
 src/pdex/__init__.py | 68 ++++++++++++++++++++++---------
 src/pdex/_math.py    | 95 ++++++++++++++++++++++++++++++++++++++------
 2 files changed, 132 insertions(+), 31 deletions(-)

diff --git a/src/pdex/__init__.py b/src/pdex/__init__.py
index e9a7e56..fec9514 100644
--- a/src/pdex/__init__.py
+++ b/src/pdex/__init__.py
@@ -14,11 +14,15 @@
 
 from pdex._math import (
     bulk_matrix_arithmetic,
+    bulk_matrix_pre_transform_mean,
     cpm_bulk,
+    cpm_from_gene_means,
     log2_fold_change,
     mwu,
+    mwu_one_vs_rest,
     percent_change,
     pseudobulk,
+    pseudobulk_from_pre_mean,
 )
 
 from ._utils import _detect_is_log1p, set_numba_threadpool
@@ -594,43 +598,71 @@ def _pdex_all(
     cpm_filter: float | None = None,
 ) -> pl.DataFrame:
     unique_groups, unique_group_indices = _unique_groups(adata.obs, groupby)
-    log.info("Found %d groups for 1-vs-rest comparison", len(unique_groups))
+    n_groups = len(unique_groups)
+    log.info("Found %d groups for 1-vs-rest comparison", n_groups)
+
+    if n_groups < 2:
+        raise ValueError(f"mode='all' requires at least 2 groups, found {n_groups}")
 
     feature_names = adata.var_names
 
+    # group ∪ rest is always the full (non-filtered) dataset in "all" mode, so the
+    # matrix is materialized once here rather than once per group (see CLAUDE.md).
+    valid_mask = np.flatnonzero(unique_group_indices >= 0)
+    codes_valid = unique_group_indices[valid_mask]
+    n_valid = valid_mask.size
+
+    global_matrix = _isolate_matrix(adata, valid_mask)
+
+    # One-shot 1-vs-rest MWU: each gene is ranked once, not once per group.
+    mwu_result = mwu_one_vs_rest(global_matrix, codes_valid, n_groups)
+    all_statistic = mwu_result.statistic
+    all_pvalue = np.asarray(mwu_result.pvalue).clip(0, 1)
+
+    # Pseudobulk pre-transform mean, computed once; "rest" is derived per group via
+    # (global_sum - group_sum) / n_rest rather than pseudobulk() on a fresh rest slice.
+    global_pre_mean = bulk_matrix_pre_transform_mean(
+        global_matrix, geometric_mean=geometric_mean, is_log1p=is_log1p
+    )
+    global_pre_sum = global_pre_mean * n_valid
+
+    if cpm_filter is not None:
+        global_arith_mean = bulk_matrix_arithmetic(global_matrix, is_log1p)
+        global_arith_sum = global_arith_mean * n_valid
+
     results = []
     for group_idx in tqdm(
-        range(len(unique_groups)),
+        range(n_groups),
         desc="Running parallel differential expression (1 vs Rest)",
     ):
         group_name = unique_groups[group_idx]
 
-        group_mask = np.flatnonzero(unique_group_indices == group_idx)
-        rest_mask = np.flatnonzero(
-            (unique_group_indices != group_idx) & (unique_group_indices >= 0)
-        )
+        local_group_mask = np.flatnonzero(codes_valid == group_idx)
+        n_group = local_group_mask.size
+        n_rest = n_valid - n_group
 
-        group_matrix = _isolate_matrix(adata, group_mask)
-        rest_matrix = _isolate_matrix(adata, rest_mask)
+        group_matrix = global_matrix[local_group_mask]
 
-        group_bulk = pseudobulk(
+        group_pre_mean = bulk_matrix_pre_transform_mean(
             group_matrix, geometric_mean=geometric_mean, is_log1p=is_log1p
         )
-        rest_bulk = pseudobulk(
-            rest_matrix, geometric_mean=geometric_mean, is_log1p=is_log1p
-        )
+        rest_pre_mean = (global_pre_sum - group_pre_mean * n_group) / n_rest
+
+        group_bulk = pseudobulk_from_pre_mean(group_pre_mean, geometric_mean)
+        rest_bulk = pseudobulk_from_pre_mean(rest_pre_mean, geometric_mean)
 
         lfc = log2_fold_change(group_bulk, rest_bulk, epsilon)
         pc = percent_change(group_bulk, rest_bulk, epsilon)
-        mwu_result = mwu(group_matrix, rest_matrix)
 
-        mwu_statistic = mwu_result.statistic
-        mwu_pvalue = np.asarray(mwu_result.pvalue).clip(0, 1)
+        mwu_statistic = all_statistic[group_idx]
+        mwu_pvalue = all_pvalue[group_idx]
 
         if cpm_filter is None:
             keep = None
         else:
-            rest_cpm = cpm_bulk(rest_matrix, is_log1p)
+            group_arith_mean = bulk_matrix_arithmetic(group_matrix, is_log1p)
+            rest_arith_mean = (global_arith_sum - group_arith_mean * n_group) / n_rest
+            rest_cpm = cpm_from_gene_means(rest_arith_mean)
             keep = _cpm_keep_mask(group_matrix, rest_cpm, is_log1p, cpm_filter)
 
         results.append(
@@ -639,8 +671,8 @@ def _pdex_all(
                 feature_names=feature_names,
                 target_bulk=group_bulk,
                 ref_bulk=rest_bulk,
-                target_membership=group_mask.size,
-                ref_membership=rest_mask.size,
+                target_membership=n_group,
+                ref_membership=n_rest,
                 lfc=lfc,
                 pc=pc,
                 pvalue=mwu_pvalue,
diff --git a/src/pdex/_math.py b/src/pdex/_math.py
index c67e321..a007a7a 100644
--- a/src/pdex/_math.py
+++ b/src/pdex/_math.py
@@ -4,6 +4,8 @@
     MannWhitneyUResult,
     SparseColumnIndex,
     mannwhitneyu_columns,
+    mannwhitneyu_one_vs_rest,
+    mannwhitneyu_one_vs_rest_sparse,
     mannwhitneyu_sparse,
 )
 from scipy.sparse import csr_matrix
@@ -84,6 +86,24 @@ def pseudobulk(
     return bulk_matrix_arithmetic(matrix, is_log1p=is_log1p)
 
 
+def _log1p_transformed_mean(
+    matrix: np.ndarray | csr_matrix, is_log1p: bool, axis=0
+) -> np.ndarray:
+    """The pre-``expm1`` intermediate that :func:`bulk_matrix_geometric` back-transforms.
+
+    This quantity is a plain mean, hence linear across row partitions: it can be
+    recovered for an unmaterialized subset of rows (e.g. "all rows except group G")
+    via ``(global_sum - subset_sum) / n_subset`` without re-reading those rows.
+    """
+    if is_log1p:
+        return np.array(matrix.mean(axis=axis)).flatten()
+    if isinstance(matrix, csr_matrix):
+        m = matrix.copy()
+        np.log1p(m.data, out=m.data)
+        return np.array(m.mean(axis=axis)).flatten()
+    return _log1p_col_mean(np.asarray(matrix, dtype=np.float64))
+
+
 def bulk_matrix_geometric(
     matrix: np.ndarray | csr_matrix, is_log1p: bool, axis=0
 ) -> np.ndarray:
@@ -94,15 +114,42 @@ def bulk_matrix_geometric(
     Both paths return values in count space.  For sparse matrices only the
     stored values are transformed (log1p(0) = 0, so sparsity is preserved).
     """
-    if is_log1p:
-        log_mean = np.array(matrix.mean(axis=axis)).flatten()
-    elif isinstance(matrix, csr_matrix):
-        m = matrix.copy()
-        np.log1p(m.data, out=m.data)
-        log_mean = np.array(m.mean(axis=axis)).flatten()
-    else:
-        log_mean = _log1p_col_mean(np.asarray(matrix, dtype=np.float64))
-    return _expm1_vec(log_mean)
+    return _expm1_vec(_log1p_transformed_mean(matrix, is_log1p, axis=axis))
+
+
+def bulk_matrix_pre_transform_mean(
+    matrix: np.ndarray | csr_matrix, geometric_mean: bool, is_log1p: bool
+) -> np.ndarray:
+    """The linear (pre-back-transform) quantity :func:`pseudobulk` is built from.
+
+    Being a plain mean, it is safe to scale by row count and subtract across row
+    partitions (``global_sum - group_sum = rest_sum``), letting a "rest of the
+    dataset" pseudobulk be derived without ever materializing those rows.
+    Pair with :func:`pseudobulk_from_pre_mean` to recover the exact value
+    :func:`pseudobulk` would have returned.
+    """
+    if geometric_mean:
+        return _log1p_transformed_mean(matrix, is_log1p)
+    return bulk_matrix_arithmetic(matrix, is_log1p)
+
+
+def pseudobulk_from_pre_mean(pre_mean: np.ndarray, geometric_mean: bool) -> np.ndarray:
+    """Back-transform :func:`bulk_matrix_pre_transform_mean`'s output to match ``pseudobulk()``."""
+    if geometric_mean:
+        return _expm1_vec(pre_mean)
+    return pre_mean
+
+
+def cpm_from_gene_means(gene_means: np.ndarray) -> np.ndarray:
+    """Normalize per-gene pooled means to counts-per-million.
+
+    Split out of :func:`cpm_bulk` so callers that already have the pooled
+    arithmetic-mean sum (e.g. derived via a global-sum-minus-subset-sum trick)
+    can skip recomputing it. See :func:`cpm_bulk` for the full semantics.
+    """
+    total = float(gene_means.sum())
+    denom = total if total != 0.0 else 1.0
+    return gene_means / denom * 1e6
 
 
 def cpm_bulk(matrix: np.ndarray | csr_matrix, is_log1p: bool) -> np.ndarray:
@@ -123,10 +170,7 @@ def cpm_bulk(matrix: np.ndarray | csr_matrix, is_log1p: bool) -> np.ndarray:
     so every gene's CPM is ``0.0`` (and is dropped by any positive threshold) rather
     than ``inf``/``nan``. Returns a flat float64 array of length ``n_genes``.
     """
-    gene_means = bulk_matrix_arithmetic(matrix, is_log1p)
-    total = float(gene_means.sum())
-    denom = total if total != 0.0 else 1.0
-    return gene_means / denom * 1e6
+    return cpm_from_gene_means(bulk_matrix_arithmetic(matrix, is_log1p))
 
 
 @nb.njit(parallel=True)
@@ -180,3 +224,28 @@ def mwu(
         return mannwhitneyu_sparse(x, y)
     else:
         return mannwhitneyu_columns(x, y)
+
+
+def mwu_one_vs_rest(
+    matrix: np.ndarray | csr_matrix, codes: np.ndarray, n_groups: int
+) -> MannWhitneyUResult:
+    """One-shot 1-vs-rest Mann-Whitney U test: each gene is ranked once against the
+    full ``matrix``, and every group's statistic/p-value (vs. all other rows) is
+    derived from that single ranking rather than by re-ranking per group.
+
+    Thin wrapper over ``numba_mwu.mannwhitneyu_one_vs_rest``/``_sparse`` (the
+    one-vs-rest generalization of ``mwu()``'s pairwise kernels, upstreamed from
+    this module — see numba_mwu's CLAUDE.md for the algorithm).
+
+    ``matrix`` must contain only the rows to be compared (e.g. filtered/-1 codes
+    already excluded); ``codes[i]`` gives the group id of row ``i`` in
+    ``[0, n_groups)``. Returns a :class:`MannWhitneyUResult` whose ``statistic``
+    and ``pvalue`` are each ``(n_groups, n_genes)`` arrays — row ``g`` is group
+    ``g``'s test against every other row in ``matrix``.
+    """
+    codes = np.ascontiguousarray(codes, dtype=np.int64)
+    if isinstance(matrix, csr_matrix):
+        return mannwhitneyu_one_vs_rest_sparse(matrix, codes, n_groups=n_groups)
+    return mannwhitneyu_one_vs_rest(
+        np.asarray(matrix, dtype=np.float64), codes, n_groups=n_groups
+    )

From aecc697137234048000724799284796317b4309a Mon Sep 17 00:00:00 2001
From: noamteyssier <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 14 Jul 2026 10:57:43 -0700
Subject: [PATCH 03/12] tests: testing one-v-rest with multiple groups and
 edge-cases

---
 tests/conftest.py  |  45 ++++++++++++++
 tests/test_math.py | 135 ++++++++++++++++++++++++++++++++++++++++
 tests/test_pdex.py | 149 +++++++++++++++++++++++++++++++++++++++++++++
 3 files changed, 329 insertions(+)

diff --git a/tests/conftest.py b/tests/conftest.py
index cb225a3..bae2eef 100644
--- a/tests/conftest.py
+++ b/tests/conftest.py
@@ -118,3 +118,48 @@ def small_adata_backed(small_adata, tmp_path):
     path = tmp_path / "test.h5ad"
     small_adata.write_h5ad(path)
     return ad.read_h5ad(path, backed="r")
+
+
+@pytest.fixture
+def multi_group_adata(rng):
+    """Synthetic AnnData for stress-testing 1-vs-rest across many uneven groups.
+
+    6 groups (sizes 1, 3, 7, 5, 2, 6 -> 24 cells), 8 genes. Most genes have
+    group-specific shifts; gene_6 is heavily tied (two distinct values) and
+    gene_7 is all-zero, exercising the degenerate (s_sq <= 0) MWU case.
+    """
+    group_sizes = {"g0": 1, "g1": 3, "g2": 7, "g3": 5, "g4": 2, "g5": 6}
+    offsets = {"g0": 0, "g1": 2, "g2": 0, "g3": 5, "g4": 0, "g5": -2}
+    n_genes = 8
+
+    obs_groups = np.concatenate(
+        [np.repeat(name, size) for name, size in group_sizes.items()]
+    )
+    n_cells = len(obs_groups)
+
+    X = rng.poisson(lam=5, size=(n_cells, n_genes)).astype(np.float64)
+    idx = 0
+    for name, size in group_sizes.items():
+        X[idx : idx + size] += offsets[name]
+        idx += size
+    X = np.clip(X, 0, None)
+
+    # gene_6: heavily tied (only two distinct values across all cells)
+    X[:, 6] = (rng.random(n_cells) > 0.5).astype(np.float64) * 3.0
+    # gene_7: zero everywhere (degenerate MWU case)
+    X[:, 7] = 0.0
+
+    obs = pd.DataFrame(
+        {"guide": obs_groups},
+        index=np.array([f"cell_{i}" for i in range(n_cells)]),
+    )
+    var = pd.DataFrame(index=np.array([f"gene_{i}" for i in range(n_genes)]))
+    return ad.AnnData(X=X, obs=obs, var=var)
+
+
+@pytest.fixture
+def multi_group_adata_sparse(multi_group_adata):
+    """multi_group_adata with sparse CSR X matrix."""
+    adata = multi_group_adata.copy()
+    adata.X = csr_matrix(adata.X)
+    return adata
diff --git a/tests/test_math.py b/tests/test_math.py
index 171ecf6..6bff5a4 100644
--- a/tests/test_math.py
+++ b/tests/test_math.py
@@ -1,13 +1,21 @@
 """Tests for pdex._math (log2_fold_change, percent_change, bulk_matrix_geometric)."""
 
 import numpy as np
+import pytest
+from scipy import stats
 from scipy.sparse import csr_matrix
 
 from pdex._math import (
+    bulk_matrix_arithmetic,
     bulk_matrix_geometric,
+    bulk_matrix_pre_transform_mean,
     cpm_bulk,
+    cpm_from_gene_means,
     log2_fold_change,
+    mwu_one_vs_rest,
     percent_change,
+    pseudobulk,
+    pseudobulk_from_pre_mean,
 )
 
 
@@ -233,3 +241,130 @@ def test_scale_invariant(self):
         """Uniformly rescaling counts does not change the CPM (ratio cancels)."""
         x = np.array([[1.0, 3.0, 7.0], [3.0, 5.0, 2.0], [5.0, 1.0, 4.0]])
         np.testing.assert_allclose(cpm_bulk(x, False), cpm_bulk(100.0 * x, False))
+
+
+class TestCpmFromGeneMeans:
+    """Tests for cpm_from_gene_means, the normalize-only half of cpm_bulk."""
+
+    def test_matches_cpm_bulk(self):
+        x = np.array([[1.0, 3.0], [3.0, 5.0], [5.0, 1.0]])
+        gene_means = bulk_matrix_arithmetic(x, is_log1p=False)
+        np.testing.assert_allclose(
+            cpm_from_gene_means(gene_means), cpm_bulk(x, is_log1p=False)
+        )
+
+
+class TestBulkMatrixPreTransformMean:
+    """Tests for bulk_matrix_pre_transform_mean / pseudobulk_from_pre_mean.
+
+    These power the "all" mode's global-sum-minus-group-sum optimization: the
+    round trip must exactly reproduce ``pseudobulk()``, and the linear quantity
+    must be safely additive across row partitions.
+    """
+
+    @pytest.mark.parametrize("geometric_mean", [True, False])
+    @pytest.mark.parametrize("is_log1p", [True, False])
+    def test_round_trip_matches_pseudobulk(self, geometric_mean, is_log1p):
+        rng = np.random.default_rng(0)
+        counts = rng.poisson(5, size=(6, 4)).astype(np.float64)
+        x = np.log1p(counts) if is_log1p else counts
+
+        pre_mean = bulk_matrix_pre_transform_mean(
+            x, geometric_mean=geometric_mean, is_log1p=is_log1p
+        )
+        result = pseudobulk_from_pre_mean(pre_mean, geometric_mean)
+        expected = pseudobulk(x, geometric_mean=geometric_mean, is_log1p=is_log1p)
+        np.testing.assert_allclose(result, expected)
+
+    def test_global_minus_group_equals_rest(self):
+        """global_sum - group_sum = rest_sum — the core "all" mode algebraic trick."""
+        rng = np.random.default_rng(1)
+        x = rng.poisson(5, size=(10, 3)).astype(np.float64)
+        group_mask = np.zeros(10, dtype=bool)
+        group_mask[:4] = True
+        n_total, n_group = x.shape[0], int(group_mask.sum())
+        n_rest = n_total - n_group
+
+        pre_mean_global = bulk_matrix_pre_transform_mean(
+            x, geometric_mean=True, is_log1p=False
+        )
+        pre_mean_group = bulk_matrix_pre_transform_mean(
+            x[group_mask], geometric_mean=True, is_log1p=False
+        )
+        rest_pre_mean = (pre_mean_global * n_total - pre_mean_group * n_group) / n_rest
+        expected = bulk_matrix_pre_transform_mean(
+            x[~group_mask], geometric_mean=True, is_log1p=False
+        )
+        np.testing.assert_allclose(rest_pre_mean, expected, rtol=1e-10)
+
+
+def _assert_mwu_one_vs_rest_matches_scipy(matrix, dense_reference, codes, n_groups):
+    result = mwu_one_vs_rest(matrix, codes, n_groups)
+    for g in range(n_groups):
+        group_mask = codes == g
+        for j in range(dense_reference.shape[1]):
+            col = dense_reference[:, j]
+            gv = col[group_mask]
+            rv = col[~group_mask]
+            scipy_result = stats.mannwhitneyu(
+                gv, rv, alternative="two-sided", method="asymptotic"
+            )
+            np.testing.assert_allclose(
+                result.statistic[g, j], scipy_result.statistic, rtol=1e-8, atol=1e-8
+            )
+            np.testing.assert_allclose(
+                result.pvalue[g, j], scipy_result.pvalue, rtol=1e-8, atol=1e-8
+            )
+    return result
+
+
+class TestMwuOneVsRest:
+    """Tests for the one-shot global-rank 1-vs-rest MWU kernel used by mode='all'."""
+
+    def test_dense_matches_scipy_multi_group(self):
+        rng = np.random.default_rng(2)
+        x = rng.poisson(5, size=(15, 3)).astype(np.float64)
+        x[:5] += 4
+        x[10:] = np.clip(x[10:] - 2, 0, None)
+        codes = np.array([0] * 5 + [1] * 5 + [2] * 5)
+        _assert_mwu_one_vs_rest_matches_scipy(x, x, codes, 3)
+
+    def test_dense_with_ties_matches_scipy(self):
+        """Ties spanning multiple groups: tie correction is shared across all groups."""
+        x = np.array([[1.0], [1.0], [2.0], [2.0], [3.0], [1.0]])
+        codes = np.array([0, 0, 1, 1, 2, 2])
+        _assert_mwu_one_vs_rest_matches_scipy(x, x, codes, 3)
+
+    def test_group_of_size_one_matches_scipy(self):
+        x = np.array([[1.0], [2.0], [3.0], [4.0], [5.0]])
+        codes = np.array([0, 1, 1, 1, 1])
+        _assert_mwu_one_vs_rest_matches_scipy(x, x, codes, 2)
+
+    def test_all_zero_column_gives_pvalue_one(self):
+        """A fully-tied (all-zero) gene needs no special-case: s_sq == 0 -> p = 1.0."""
+        x = np.zeros((6, 1))
+        codes = np.array([0, 0, 1, 1, 2, 2])
+        result = mwu_one_vs_rest(x, codes, 3)
+        np.testing.assert_allclose(result.pvalue, 1.0)
+        assert np.isfinite(result.statistic).all()
+
+    def test_sparse_matches_dense(self):
+        rng = np.random.default_rng(3)
+        x = rng.poisson(2, size=(12, 4)).astype(np.float64)
+        x[x < 1.5] = 0.0
+        codes = np.array([0] * 4 + [1] * 3 + [2] * 5)
+        dense_result = mwu_one_vs_rest(x, codes, 3)
+        sparse_result = mwu_one_vs_rest(csr_matrix(x), codes, 3)
+        np.testing.assert_allclose(
+            dense_result.statistic, sparse_result.statistic, rtol=1e-8
+        )
+        np.testing.assert_allclose(
+            dense_result.pvalue, sparse_result.pvalue, rtol=1e-8, atol=1e-10
+        )
+
+    def test_sparse_matches_scipy(self):
+        rng = np.random.default_rng(4)
+        x = rng.poisson(2, size=(12, 3)).astype(np.float64)
+        x[x < 1.5] = 0.0
+        codes = np.array([0] * 4 + [1] * 3 + [2] * 5)
+        _assert_mwu_one_vs_rest_matches_scipy(csr_matrix(x), x, codes, 3)
diff --git a/tests/test_pdex.py b/tests/test_pdex.py
index c0eb5f7..4e0465e 100644
--- a/tests/test_pdex.py
+++ b/tests/test_pdex.py
@@ -2,6 +2,7 @@
 
 import anndata as ad
 import numpy as np
+import pandas as pd
 import polars as pl
 import pytest
 from scipy import stats
@@ -285,6 +286,154 @@ def test_sparse_dense_agreement(self, small_adata, small_adata_sparse):
             )
 
 
+def _assert_all_mode_matches_scipy_per_group(result, adata_dense, skip_genes=()):
+    """Cross-checks every (group, gene) pair against an independently-computed
+    scipy.stats.mannwhitneyu — a stronger oracle than diffing against pdex's own
+    prior implementation, since it validates the multi-group rank-sum reduction
+    from first principles.
+
+    ``skip_genes`` excludes genes where pdex deliberately diverges from raw scipy:
+    a fully-tied column (all values equal) makes scipy's variance term 0/0 -> NaN,
+    while pdex (matching numba_mwu's convention) defines this degenerate case as
+    p=1.0 (see test_all_zero_gene_pvalue_is_one).
+    """
+    X_dense = np.asarray(adata_dense.X)
+    obs = adata_dense.obs
+    for group_name in obs["guide"].unique():
+        group_mask = (obs["guide"] == group_name).to_numpy()
+        rest_mask = ~group_mask
+        group_rows = result.filter(pl.col("target") == group_name)
+        for gene_idx, gene in enumerate(adata_dense.var_names):
+            if gene in skip_genes:
+                continue
+            gv = X_dense[group_mask, gene_idx]
+            rv = X_dense[rest_mask, gene_idx]
+            scipy_result = stats.mannwhitneyu(
+                gv, rv, alternative="two-sided", method="asymptotic"
+            )
+            row = group_rows.filter(pl.col("feature") == gene)
+            np.testing.assert_allclose(
+                row["statistic"][0], scipy_result.statistic, rtol=1e-6, atol=1e-8
+            )
+            np.testing.assert_allclose(
+                row["p_value"][0], scipy_result.pvalue, rtol=1e-6, atol=1e-8
+            )
+
+
+class TestPdexAllModeMultiGroup:
+    """Stress tests for the 1-vs-rest global-rank optimization across many groups,
+    uneven group sizes (including size 1), ties, and an all-zero gene."""
+
+    def test_matches_scipy_per_group_dense(self, multi_group_adata):
+        result = pdex(multi_group_adata, groupby="guide", mode="all", is_log1p=False)
+        _assert_all_mode_matches_scipy_per_group(
+            result, multi_group_adata, skip_genes={"gene_7"}
+        )
+
+    def test_matches_scipy_per_group_sparse(
+        self, multi_group_adata, multi_group_adata_sparse
+    ):
+        result = pdex(
+            multi_group_adata_sparse, groupby="guide", mode="all", is_log1p=False
+        )
+        _assert_all_mode_matches_scipy_per_group(
+            result, multi_group_adata, skip_genes={"gene_7"}
+        )
+
+    def test_sparse_dense_agreement(self, multi_group_adata, multi_group_adata_sparse):
+        dense_result = pdex(
+            multi_group_adata, groupby="guide", mode="all", is_log1p=False
+        ).sort(["target", "feature"])
+        sparse_result = pdex(
+            multi_group_adata_sparse, groupby="guide", mode="all", is_log1p=False
+        ).sort(["target", "feature"])
+
+        for col in [
+            "p_value",
+            "statistic",
+            "fold_change",
+            "percent_change",
+            "target_mean",
+            "ref_mean",
+        ]:
+            np.testing.assert_allclose(
+                dense_result[col].to_numpy(),
+                sparse_result[col].to_numpy(),
+                rtol=1e-6,
+                atol=1e-8,
+                err_msg=f"Mismatch in column {col}",
+            )
+
+    def test_all_zero_gene_pvalue_is_one(self, multi_group_adata):
+        """gene_7 is zero everywhere: the degenerate s_sq <= 0 case needs no
+        special-casing but must still yield p=1.0 and a finite statistic."""
+        result = pdex(multi_group_adata, groupby="guide", mode="all", is_log1p=False)
+        zero_gene_rows = result.filter(pl.col("feature") == "gene_7")
+        assert (zero_gene_rows["p_value"] == 1.0).all()
+        assert np.isfinite(zero_gene_rows["statistic"].to_numpy()).all()
+
+    def test_membership_counts(self, multi_group_adata):
+        result = pdex(multi_group_adata, groupby="guide", mode="all", is_log1p=False)
+        n_total = multi_group_adata.n_obs
+        for group_name in multi_group_adata.obs["guide"].unique():
+            expected_group = (multi_group_adata.obs["guide"] == group_name).sum()
+            expected_rest = n_total - expected_group
+            group_rows = result.filter(pl.col("target") == group_name)
+            assert group_rows["target_membership"].unique().to_list() == [
+                expected_group
+            ]
+            assert group_rows["ref_membership"].unique().to_list() == [expected_rest]
+
+    def test_filtered_group_excluded_from_both_sides(self, rng):
+        """Cells with a NaN/empty groupby value are excluded from every group and
+        from 'rest', for every comparison (no test previously exercised this at
+        the pdex() level for mode='all')."""
+        n_per = 5
+        n_genes = 3
+        groups = np.array(["A"] * n_per + ["B"] * n_per + [""] * n_per)
+        X = rng.poisson(lam=5, size=(len(groups), n_genes)).astype(np.float64)
+        obs = pd.DataFrame(
+            {"guide": groups}, index=[f"cell_{i}" for i in range(len(groups))]
+        )
+        var = pd.DataFrame(index=[f"gene_{i}" for i in range(n_genes)])
+        adata = ad.AnnData(X=X, obs=obs, var=var)
+
+        result = pdex(adata, groupby="guide", mode="all", is_log1p=False)
+        assert set(result["target"].unique().to_list()) == {"A", "B"}
+        # rest for "A" is only "B" (n_per cells), not "B" + filtered cells
+        a_rows = result.filter(pl.col("target") == "A")
+        assert a_rows["target_membership"].unique().to_list() == [n_per]
+        assert a_rows["ref_membership"].unique().to_list() == [n_per]
+
+    def test_single_group_raises(self, rng):
+        n_genes = 3
+        X = rng.poisson(lam=5, size=(6, n_genes)).astype(np.float64)
+        obs = pd.DataFrame(
+            {"guide": ["only"] * 6}, index=[f"cell_{i}" for i in range(6)]
+        )
+        var = pd.DataFrame(index=[f"gene_{i}" for i in range(n_genes)])
+        adata = ad.AnnData(X=X, obs=obs, var=var)
+
+        with pytest.raises(ValueError, match="at least 2 groups"):
+            pdex(adata, groupby="guide", mode="all", is_log1p=False)
+
+    def test_cpm_filter_drops_floor_across_many_groups(self, multi_group_adata):
+        """Extends TestCpmFilter::test_all_mode_drops_floor to 6 uneven groups:
+        the all-zero gene_7 is dropped for every group under the new rest_cpm
+        (global-sum-minus-group-sum) derivation."""
+        result = pdex(
+            multi_group_adata,
+            groupby="guide",
+            mode="all",
+            is_log1p=False,
+            cpm_filter=5,
+        )
+        assert "gene_7" not in result["feature"].to_list()
+        assert set(result["target"].to_list()) == set(
+            multi_group_adata.obs["guide"].unique()
+        )
+
+
 class TestPdexOnTargetMode:
     """Tests for pdex(..., mode='on_target')."""
 

From e67d08ee070371133ea1a0e97541733a4611aed7 Mon Sep 17 00:00:00 2001
From: noamteyssier <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 14 Jul 2026 10:57:55 -0700
Subject: [PATCH 04/12] dep: update numba-mwu

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 9824f54..6c01551 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -8,7 +8,7 @@ requires-python = ">=3.11"
 dependencies = [
     "anndata>=0.12.10",
     "numba>=0.64.0",
-    "numba-mwu>=0.1.1",
+    "numba-mwu>=0.2.0",
     "numpy>=2.4.2",
     "pandas>=2.3.3",
     "polars>=1.38.1",

From 02f1413df319d2d3a8e715e8ed14bfa31d63d33f Mon Sep 17 00:00:00 2001
From: noamteyssier <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 14 Jul 2026 10:58:03 -0700
Subject: [PATCH 05/12] doc: update claude

---
 CLAUDE.md | 8 ++++++--
 1 file changed, 6 insertions(+), 2 deletions(-)

diff --git a/CLAUDE.md b/CLAUDE.md
index 89ac38c..47aa6ab 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -73,15 +73,19 @@ checked empirically (inspect the per-gene CPM distribution).
 | File                   | Role                                                                                                    |
 | ---------------------- | ------------------------------------------------------------------------------------------------------- |
 | `src/pdex/__init__.py` | `pdex()` entry point and full pipeline logic                                                            |
-| `src/pdex/_math.py`    | Numba JIT-compiled `fold_change()`, `percent_change()`, and `mwu()` wrappers; `pseudobulk()` dispatcher; `cpm_bulk()` pooled-CPM view for the filter |
+| `src/pdex/_math.py`    | Numba JIT-compiled `fold_change()`, `percent_change()`, and `mwu()`/`mwu_one_vs_rest()` wrappers over `numba-mwu`; `pseudobulk()` dispatcher; `cpm_bulk()` pooled-CPM view for the filter; the `bulk_matrix_pre_transform_mean()`/`pseudobulk_from_pre_mean()`/`cpm_from_gene_means()` trio powering `"all"` mode's global-sum optimization |
 | `src/pdex/_utils.py`   | `set_numba_threadpool()` — sets Numba thread count before JIT warmup; `_available_cpus()` — affinity-aware CPU count (respects cgroup/SLURM limits); `_detect_is_log1p()` heuristic |
 
 ### Performance Design
 
 - Numba JIT compilation accelerates per-cell/per-gene math (`fold_change`, `percent_change`, `_log1p_col_mean`, `_expm1_vec`)
-- `numba-mwu` (external dep) provides a Numba-accelerated Mann-Whitney U implementation
+- `numba-mwu` (external dep, `>=0.2.0`) provides Numba-accelerated Mann-Whitney U kernels for **both** the pairwise case (`mannwhitneyu_columns`/`mannwhitneyu_sparse`, used by `"ref"` and `"on_target"` modes) and the one-vs-rest case (`mannwhitneyu_one_vs_rest`/`_sparse`, used by `"all"` mode) — the one-vs-rest kernels originated in pdex and were upstreamed into `numba-mwu` since the optimization is domain-agnostic (see that package's CLAUDE.md for the algorithm).
 - Sparse CSR matrices are handled by reusing pre-computed non-targeting column indices to avoid redundant dense conversion
 - Parallelism is controlled via `threads` passed to `set_numba_threadpool()`
+- **`"all"` mode (1-vs-rest) is a one-shot computation, not a per-group loop over `numba-mwu`.** Because `group ∪ rest` is always the full (non-filtered) dataset regardless of which group is being tested, `_pdex_all` materializes the expression matrix exactly once (instead of once per group) and:
+  - Ranks each gene once across all cells via `mwu_one_vs_rest()` (`_math.py`, a thin wrapper over `numba_mwu.mannwhitneyu_one_vs_rest`/`_sparse`) and reduces to every group's rank-sum in the same pass, rather than re-ranking group+rest from scratch per group.
+  - Derives each group's "rest" pseudobulk and CPM algebraically as `(global_sum - group_sum) / n_rest` (`bulk_matrix_pre_transform_mean()`/`pseudobulk_from_pre_mean()`/`cpm_from_gene_means()`) instead of recomputing over a freshly sliced "rest" matrix.
+  - This turns `_pdex_all` from `O(n_groups × n_obs)` matrix I/O and ranking into ~`O(n_obs)` total, which matters most for screens with many groups (e.g. guides) and/or large cell counts.
 
 ### Output Schema
 

From 8c82028cca940e6a06888e95cbbe37909a1f1acc Mon Sep 17 00:00:00 2001
From: noamteyssier <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 14 Jul 2026 11:18:08 -0700
Subject: [PATCH 06/12] fix: potential floating point noise and regression test

---
 src/pdex/__init__.py | 12 ++++++++--
 tests/test_pdex.py   | 55 ++++++++++++++++++++++++++++++++++++++++++++
 2 files changed, 65 insertions(+), 2 deletions(-)

diff --git a/src/pdex/__init__.py b/src/pdex/__init__.py
index fec9514..290cf4f 100644
--- a/src/pdex/__init__.py
+++ b/src/pdex/__init__.py
@@ -646,7 +646,12 @@ def _pdex_all(
         group_pre_mean = bulk_matrix_pre_transform_mean(
             group_matrix, geometric_mean=geometric_mean, is_log1p=is_log1p
         )
-        rest_pre_mean = (global_pre_sum - group_pre_mean * n_group) / n_rest
+        # Clip: rest_pre_mean is a mean of non-negative data, so it can never be
+        # legitimately negative. This guards against floating-point noise which can
+        # otherwise corrupt log2_fold_change/percent_change.
+        rest_pre_mean = np.clip(
+            (global_pre_sum - group_pre_mean * n_group) / n_rest, 0, None
+        )
 
         group_bulk = pseudobulk_from_pre_mean(group_pre_mean, geometric_mean)
         rest_bulk = pseudobulk_from_pre_mean(rest_pre_mean, geometric_mean)
@@ -661,7 +666,10 @@ def _pdex_all(
             keep = None
         else:
             group_arith_mean = bulk_matrix_arithmetic(group_matrix, is_log1p)
-            rest_arith_mean = (global_arith_sum - group_arith_mean * n_group) / n_rest
+            # Same floating-point-cancellation guard as rest_pre_mean above.
+            rest_arith_mean = np.clip(
+                (global_arith_sum - group_arith_mean * n_group) / n_rest, 0, None
+            )
             rest_cpm = cpm_from_gene_means(rest_arith_mean)
             keep = _cpm_keep_mask(group_matrix, rest_cpm, is_log1p, cpm_filter)
 
diff --git a/tests/test_pdex.py b/tests/test_pdex.py
index 4e0465e..177d974 100644
--- a/tests/test_pdex.py
+++ b/tests/test_pdex.py
@@ -889,6 +889,61 @@ def test_one_sided_zero_still_infinite(self, small_adata):
         np.testing.assert_allclose(gene0["percent_change"].to_numpy(), -1.0)
 
 
+class TestAllModeRestMeanFloatingPointCancellation:
+    """Regression test for `_pdex_all`'s "rest = global - group" derivation.
+
+    ``rest_pre_mean``/``rest_arith_mean`` are each computed as the difference
+    of two independently-rounded means (see __init__.py::_pdex_all). For a
+    gene expressed only in the target group (exactly zero in "rest"), this
+    difference is mathematically exactly zero but can land a hair below zero
+    in floating point (e.g. -3e-15) — impossible for a mean of non-negative
+    data, but enough to corrupt log2_fold_change (silently masked to 0.0, as
+    if "unexpressed in both groups") and percent_change (a large *negative*
+    finite value) when epsilon=0.0. `_pdex_all` clips both to ``[0, None]``.
+    """
+
+    def test_one_sided_zero_via_subtraction_is_not_masked(self):
+        n_group, n_rest = 3857, 555
+        rng = np.random.default_rng(1)
+        group_vals = np.log1p(rng.poisson(10, size=n_group).astype(np.float64) + 1)
+
+        X = np.zeros((n_group + n_rest, 1))
+        X[:n_group, 0] = group_vals
+        obs = pd.DataFrame(
+            {"guide": ["A"] * n_group + ["B"] * n_rest},
+            index=[f"c{i}" for i in range(n_group + n_rest)],
+        )
+        var = pd.DataFrame(index=["gene_0"])
+        adata = ad.AnnData(X=X, obs=obs, var=var)
+
+        result = pdex(
+            adata,
+            groupby="guide",
+            mode="all",
+            is_log1p=True,
+            epsilon=0.0,
+            geometric_mean=True,
+        )
+        row = result.filter(pl.col("target") == "A")
+
+        ref_mean = row["ref_mean"][0]
+        lfc = row["log2_fold_change"][0]
+        pc = row["percent_change"][0]
+
+        # The clip's direct invariant: a mean of non-negative data can never
+        # be negative, regardless of which way the floating-point noise breaks.
+        assert ref_mean >= 0.0
+
+        # gene_0 is expressed in the target and unexpressed in rest, so this
+        # must read as a large *positive* change -- never NaN, and never
+        # masked to 0.0 (the distinct "unexpressed in both groups" case,
+        # which does not apply here since target_mean > 0).
+        assert not np.isnan(lfc)
+        assert np.isposinf(lfc) or lfc > 10
+        assert not np.isnan(pc)
+        assert pc > 0
+
+
 def _pairs(df) -> set:
     """Set of (target, feature) tuples in a result frame."""
     return set(zip(df["target"].to_list(), df["feature"].to_list()))

From afe08ebcfe2ca291972b786e9cc2fd664c1363df Mon Sep 17 00:00:00 2001
From: noamteyssier <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 14 Jul 2026 11:21:17 -0700
Subject: [PATCH 07/12] fix: only run ci once on PR

---
 .github/workflows/ci.yml | 5 ++++-
 1 file changed, 4 insertions(+), 1 deletion(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 958bbb7..41212ab 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -1,6 +1,9 @@
 name: "CI"
 
-on: [push, pull_request]
+on:
+  push:
+    branches: [main]
+  pull_request:
 
 jobs:
   all_jobs:

From 38cec6e4fa3213d7c056963203c4a9faf985d2a3 Mon Sep 17 00:00:00 2001
From: noamteyssier <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 14 Jul 2026 11:31:28 -0700
Subject: [PATCH 08/12] ci: only trigger semver-check when PR to main

---
 .github/workflows/ci.yml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/.github/workflows/ci.yml b/.github/workflows/ci.yml
index 41212ab..d8957ca 100644
--- a/.github/workflows/ci.yml
+++ b/.github/workflows/ci.yml
@@ -21,7 +21,7 @@ jobs:
 
   semver-check:
     runs-on: ubuntu-latest
-    if: github.event_name == 'pull_request'
+    if: github.event_name == 'pull_request' && github.base_ref == 'main'
 
     steps:
       - uses: actions/checkout@v4

From 64386ff372126ae8dc7468b5628813640e04dbf6 Mon Sep 17 00:00:00 2001
From: noamteyssier <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 14 Jul 2026 11:38:46 -0700
Subject: [PATCH 09/12] refactor: completely deprecated bare fold_change column

---
 CLAUDE.md            |  5 ++---
 src/pdex/__init__.py | 23 ++++-------------------
 tests/test_pdex.py   | 27 +++++++++++++--------------
 3 files changed, 19 insertions(+), 36 deletions(-)

diff --git a/CLAUDE.md b/CLAUDE.md
index 47aa6ab..b1004f6 100644
--- a/CLAUDE.md
+++ b/CLAUDE.md
@@ -73,12 +73,12 @@ checked empirically (inspect the per-gene CPM distribution).
 | File                   | Role                                                                                                    |
 | ---------------------- | ------------------------------------------------------------------------------------------------------- |
 | `src/pdex/__init__.py` | `pdex()` entry point and full pipeline logic                                                            |
-| `src/pdex/_math.py`    | Numba JIT-compiled `fold_change()`, `percent_change()`, and `mwu()`/`mwu_one_vs_rest()` wrappers over `numba-mwu`; `pseudobulk()` dispatcher; `cpm_bulk()` pooled-CPM view for the filter; the `bulk_matrix_pre_transform_mean()`/`pseudobulk_from_pre_mean()`/`cpm_from_gene_means()` trio powering `"all"` mode's global-sum optimization |
+| `src/pdex/_math.py`    | Numba JIT-compiled `log2_fold_change()`, `percent_change()`, and `mwu()`/`mwu_one_vs_rest()` wrappers over `numba-mwu`; `pseudobulk()` dispatcher; `cpm_bulk()` pooled-CPM view for the filter; the `bulk_matrix_pre_transform_mean()`/`pseudobulk_from_pre_mean()`/`cpm_from_gene_means()` trio powering `"all"` mode's global-sum optimization |
 | `src/pdex/_utils.py`   | `set_numba_threadpool()` — sets Numba thread count before JIT warmup; `_available_cpus()` — affinity-aware CPU count (respects cgroup/SLURM limits); `_detect_is_log1p()` heuristic |
 
 ### Performance Design
 
-- Numba JIT compilation accelerates per-cell/per-gene math (`fold_change`, `percent_change`, `_log1p_col_mean`, `_expm1_vec`)
+- Numba JIT compilation accelerates per-cell/per-gene math (`log2_fold_change`, `percent_change`, `_log1p_col_mean`, `_expm1_vec`)
 - `numba-mwu` (external dep, `>=0.2.0`) provides Numba-accelerated Mann-Whitney U kernels for **both** the pairwise case (`mannwhitneyu_columns`/`mannwhitneyu_sparse`, used by `"ref"` and `"on_target"` modes) and the one-vs-rest case (`mannwhitneyu_one_vs_rest`/`_sparse`, used by `"all"` mode) — the one-vs-rest kernels originated in pdex and were upstreamed into `numba-mwu` since the optimization is domain-agnostic (see that package's CLAUDE.md for the algorithm).
 - Sparse CSR matrices are handled by reusing pre-computed non-targeting column indices to avoid redundant dense conversion
 - Parallelism is controlled via `threads` passed to `set_numba_threadpool()`
@@ -99,7 +99,6 @@ The returned Polars DataFrame (or pandas DataFrame when `as_pandas=True`) has co
 | `ref_mean`          | float | Pseudobulk mean for the reference, always in natural (count) space    |
 | `target_membership` | int   | Number of cells in the target group                                   |
 | `ref_membership`    | int   | Number of cells in the reference                                      |
-| `fold_change`       | float | **Deprecated** alias for `log2_fold_change` (identical values). Retained for one release; emits a `FutureWarning` on every `pdex(...)` call and will be removed in pdex 0.3.0. |
 | `log2_fold_change`  | float | log2((target_mean + epsilon) / (ref_mean + epsilon)) — computed from pseudobulk means. `epsilon` defaults to `1e-9` (finite-guard), so by default there are no `±inf`/`NaN`: one-sided zeros become large-but-finite and `0/0` is `0.0`. With `epsilon == 0`, `0/0` is still defined as `0.0` (not `NaN`) but one-sided zeros yield `±inf`. |
 | `percent_change`    | float | (target_mean - ref_mean) / (ref_mean + epsilon) — computed from pseudobulk means. With the default `epsilon=1e-9` there are no non-finite values; with `epsilon == 0`, `0/0` is `0.0` (not `NaN`) and a zero reference with nonzero target yields `+inf`. |
 | `p_value`           | float | Mann-Whitney U p-value (per-cell vectors)                             |
diff --git a/src/pdex/__init__.py b/src/pdex/__init__.py
index 290cf4f..e8c0d8e 100644
--- a/src/pdex/__init__.py
+++ b/src/pdex/__init__.py
@@ -256,7 +256,7 @@ def pdex(
     epsilon:
         Pseudocount added to the **native (count-space) means** — the denominator
         (and, for ``log2_fold_change``, the numerator) — before computing
-        ``fold_change`` and ``percent_change``. It is never applied to the CPM view
+        ``log2_fold_change`` and ``percent_change``. It is never applied to the CPM view
         used by ``cpm_filter``. When ``epsilon > 0``, extreme values from near-zero
         reference means (scRNA-seq sparsity artifact) are dampened toward zero, and
         one-sided zeros become large-but-finite instead of ``±inf``. Has no effect on
@@ -307,7 +307,7 @@ def pdex(
     pl.DataFrame | pd.DataFrame
         One row per (group, feature) pair with columns: ``target``, ``feature``,
         ``target_mean``, ``ref_mean``, ``target_membership``, ``ref_membership``,
-        ``fold_change``, ``log2_fold_change``, ``percent_change``, ``p_value``,
+        ``log2_fold_change``, ``percent_change``, ``p_value``,
         ``statistic``, ``fdr``.
 
         ``target_mean`` and ``ref_mean`` are always in **natural (count) space**.
@@ -322,12 +322,8 @@ def pdex(
         and ``percent_change`` define this as ``0.0`` (no change) rather than
         ``NaN``. One-sided zeros still produce ``±inf``.
 
-        ``fold_change`` is a **deprecated** alias for ``log2_fold_change``
-        (identical values). It is retained for one release to ease migration
-        and will be removed in pdex 0.3.0. New code should read
-        ``log2_fold_change`` directly. A :class:`FutureWarning` is emitted
-        on every ``pdex(...)`` call.  The MWU ``p_value`` and
-        ``statistic`` are computed directly on the per-cell expression vectors.
+        The MWU ``p_value`` and ``statistic`` are computed directly on the
+        per-cell expression vectors.
 
         For ``mode="ref"``, the reference group itself is excluded from the output.
 
@@ -345,14 +341,6 @@ def pdex(
     if epsilon < 0:
         raise ValueError(f"epsilon must be non-negative, got {epsilon}")
 
-    warnings.warn(
-        "The `fold_change` column in pdex output is deprecated and will be "
-        "removed in pdex 0.3.0. Use `log2_fold_change` instead — it contains "
-        "the same values (`log2(target_mean / ref_mean)`).",
-        FutureWarning,
-        stacklevel=2,
-    )
-
     # Set the global threadpool for numba
     set_numba_threadpool(threads)
 
@@ -502,7 +490,6 @@ def _assemble_group_frame(
             "ref_mean": ref_mean,
             "target_membership": np.full(feature.shape[0], target_membership),
             "ref_membership": np.full(feature.shape[0], ref_membership),
-            "fold_change": lfc,
             "log2_fold_change": lfc,
             "percent_change": pc,
             "p_value": pvalue,
@@ -790,7 +777,6 @@ def _pdex_on_target(
                 "ref_mean": ref_mean,
                 "target_membership": group_mask.size,
                 "ref_membership": ref_membership,
-                "fold_change": lfc,
                 "log2_fold_change": lfc,
                 "percent_change": pc,
                 "p_value": p_value,
@@ -809,7 +795,6 @@ def _pdex_on_target(
                 "ref_mean": pl.Float64,
                 "target_membership": pl.Int64,
                 "ref_membership": pl.Int64,
-                "fold_change": pl.Float64,
                 "log2_fold_change": pl.Float64,
                 "percent_change": pl.Float64,
                 "p_value": pl.Float64,
diff --git a/tests/test_pdex.py b/tests/test_pdex.py
index 177d974..a5f38b5 100644
--- a/tests/test_pdex.py
+++ b/tests/test_pdex.py
@@ -17,7 +17,6 @@
     "ref_mean",
     "target_membership",
     "ref_membership",
-    "fold_change",
     "log2_fold_change",
     "percent_change",
     "p_value",
@@ -87,7 +86,7 @@ def test_fold_change_sign(self, small_adata):
         for group_name in ["A", "B"]:
             group_rows = result.filter(pl.col("target") == group_name)
             # Mean fold change should be positive since we boosted these groups
-            mean_fc = group_rows["fold_change"].mean()
+            mean_fc = group_rows["log2_fold_change"].mean()
             assert mean_fc > 0, f"Expected positive fold change for group {group_name}"  # type: ignore
 
     def test_statistics_against_scipy(self, small_adata):
@@ -177,7 +176,7 @@ def test_sparse_dense_agreement(self, small_adata, small_adata_sparse):
 
         assert dense_result.shape == sparse_result.shape
 
-        for col in ["p_value", "statistic", "fold_change", "percent_change"]:
+        for col in ["p_value", "statistic", "log2_fold_change", "percent_change"]:
             np.testing.assert_allclose(
                 dense_result[col].to_numpy(),
                 sparse_result[col].to_numpy(),
@@ -237,7 +236,7 @@ def test_fold_change_sign(self, small_adata):
         """Group B was boosted the most, so its fold change vs rest should be positive."""
         result = pdex(small_adata, groupby="guide", mode="all", is_log1p=False)
         group_b_rows = result.filter(pl.col("target") == "B")
-        mean_fc = group_b_rows["fold_change"].mean()
+        mean_fc = group_b_rows["log2_fold_change"].mean()
         assert mean_fc > 0  # type: ignore
 
     def test_statistics_against_scipy(self, small_adata):
@@ -277,7 +276,7 @@ def test_sparse_dense_agreement(self, small_adata, small_adata_sparse):
 
         assert dense_result.shape == sparse_result.shape
 
-        for col in ["p_value", "statistic", "fold_change", "percent_change"]:
+        for col in ["p_value", "statistic", "log2_fold_change", "percent_change"]:
             np.testing.assert_allclose(
                 dense_result[col].to_numpy(),
                 sparse_result[col].to_numpy(),
@@ -351,7 +350,7 @@ def test_sparse_dense_agreement(self, multi_group_adata, multi_group_adata_spars
         for col in [
             "p_value",
             "statistic",
-            "fold_change",
+            "log2_fold_change",
             "percent_change",
             "target_mean",
             "ref_mean",
@@ -572,7 +571,7 @@ def test_sparse_dense_agreement(self, on_target_adata, on_target_adata_sparse):
         )
 
         assert dense_result.shape == sparse_result.shape
-        for col in ["p_value", "statistic", "fold_change", "percent_change"]:
+        for col in ["p_value", "statistic", "log2_fold_change", "percent_change"]:
             np.testing.assert_allclose(
                 dense_result[col].to_numpy(),
                 sparse_result[col].to_numpy(),
@@ -744,7 +743,7 @@ def test_both_log1p_paths_agree(self, small_adata, small_adata_log1p):
         """pdex on raw counts with is_log1p=False and on log1p counts with is_log1p=True
         should yield identical results across all output columns.
 
-        Pseudobulk means back-transform to the same count space, so fold_change and
+        Pseudobulk means back-transform to the same count space, so log2_fold_change and
         percent_change must match. The MWU statistic and p_value operate on the raw
         cell-level values (which differ between the two inputs), so they are NOT
         expected to match — only the pseudobulk-derived columns are tested here.
@@ -763,7 +762,7 @@ def test_both_log1p_paths_agree(self, small_adata, small_adata_log1p):
             is_log1p=True,
             geometric_mean=True,
         )
-        for col in ["target_mean", "ref_mean", "fold_change", "percent_change"]:
+        for col in ["target_mean", "ref_mean", "log2_fold_change", "percent_change"]:
             np.testing.assert_allclose(
                 raw_result[col].to_numpy(),
                 log_result[col].to_numpy(),
@@ -797,7 +796,7 @@ def test_ref_mode_backed_matches_inmemory(self, small_adata, small_adata_backed)
         inmem = pdex(small_adata, groupby="guide", mode="ref", is_log1p=False)
         backed = pdex(small_adata_backed, groupby="guide", mode="ref", is_log1p=False)
         assert inmem.shape == backed.shape
-        for col in ["p_value", "statistic", "fold_change", "percent_change"]:
+        for col in ["p_value", "statistic", "log2_fold_change", "percent_change"]:
             np.testing.assert_allclose(
                 inmem[col].to_numpy(),
                 backed[col].to_numpy(),
@@ -809,7 +808,7 @@ def test_all_mode_backed_matches_inmemory(self, small_adata, small_adata_backed)
         inmem = pdex(small_adata, groupby="guide", mode="all", is_log1p=False)
         backed = pdex(small_adata_backed, groupby="guide", mode="all", is_log1p=False)
         assert inmem.shape == backed.shape
-        for col in ["p_value", "statistic", "fold_change", "percent_change"]:
+        for col in ["p_value", "statistic", "log2_fold_change", "percent_change"]:
             np.testing.assert_allclose(
                 inmem[col].to_numpy(),
                 backed[col].to_numpy(),
@@ -849,7 +848,7 @@ def test_zero_in_both_is_zero_not_nan(self, small_adata, mode):
 
         assert (gene0["target_mean"].to_numpy() == 0).all()
         assert (gene0["ref_mean"].to_numpy() == 0).all()
-        for col in ["log2_fold_change", "fold_change", "percent_change"]:
+        for col in ["log2_fold_change", "percent_change"]:
             values = gene0[col].to_numpy()
             assert not np.isnan(values).any(), f"{col} contains NaN"
             np.testing.assert_array_equal(values, 0.0)
@@ -870,7 +869,7 @@ def test_on_target_zero_in_both_is_zero_not_nan(self, on_target_adata):
         row = result.filter(pl.col("target") == "A")
         assert row["target_mean"].to_numpy()[0] == 0
         assert row["ref_mean"].to_numpy()[0] == 0
-        for col in ["log2_fold_change", "fold_change", "percent_change"]:
+        for col in ["log2_fold_change", "percent_change"]:
             value = row[col].to_numpy()[0]
             assert not np.isnan(value), f"{col} is NaN"
             assert value == 0.0
@@ -1129,7 +1128,7 @@ def test_backed_sparse_matches_inmemory_sparse(
         for col in [
             "target_mean",
             "ref_mean",
-            "fold_change",
+            "log2_fold_change",
             "percent_change",
             "p_value",
             "statistic",

From b09fbf84400dc85dd001ea92c2029e7d3e3d01b5 Mon Sep 17 00:00:00 2001
From: noamteyssier <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 14 Jul 2026 11:44:30 -0700
Subject: [PATCH 10/12] chore(semver): bump

---
 pyproject.toml | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/pyproject.toml b/pyproject.toml
index 6c01551..ca2a3a1 100644
--- a/pyproject.toml
+++ b/pyproject.toml
@@ -1,6 +1,6 @@
 [project]
 name = "pdex"
-version = "0.2.5"
+version = "0.3.0"
 description = "Parallel differential expression for single-cell perturbation sequencing"
 readme = "README.md"
 authors = [{ name = "noam teyssier", email = "noam.teyssier@arcinstitute.org" }]

From 585c5e338a8a0ca3f5fc01067d151917c33667b8 Mon Sep 17 00:00:00 2001
From: noamteyssier <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 14 Jul 2026 11:46:18 -0700
Subject: [PATCH 11/12] docs: remove deprecated fold_change from docs

---
 README.md | 1 -
 1 file changed, 1 deletion(-)

diff --git a/README.md b/README.md
index 3f54929..747e070 100644
--- a/README.md
+++ b/README.md
@@ -103,7 +103,6 @@ Returns a Polars DataFrame (or pandas if `as_pandas=True`) with one row per (gro
 | `ref_mean`          | Pseudobulk mean for the reference (count space)    |
 | `target_membership` | Number of cells in the target group                |
 | `ref_membership`    | Number of cells in the reference                   |
-| `fold_change`       | **Deprecated alias** for `log2_fold_change` (identical values). Will be removed in pdex 0.3.0. |
 | `log2_fold_change`  | log2((target_mean + epsilon) / (ref_mean + epsilon)). With default `epsilon=1e-9`, one-sided zeros are large finite values; with `epsilon=0.0`, one-sided zeros yield `±inf`. Genes unexpressed in both groups (0/0) report `0.0`, not `NaN`. |
 | `percent_change`    | (target_mean - ref_mean) / (ref_mean + epsilon). With default `epsilon=1e-9`, zero-reference cases are finite; with `epsilon=0.0`, a zero reference with nonzero target yields `+inf`. Genes unexpressed in both groups (0/0) report `0.0`, not `NaN`. |
 | `p_value`           | Mann-Whitney U p-value                             |

From 5e37804ecf3399bff79da6d93c15578b72aaa19d Mon Sep 17 00:00:00 2001
From: noamteyssier <22600644+noamteyssier@users.noreply.github.com>
Date: Tue, 14 Jul 2026 11:46:58 -0700
Subject: [PATCH 12/12] docs: added changelog

---
 CHANGELOG.md | 48 ++++++++++++++++++++++++++++++++++++++++++++++++
 1 file changed, 48 insertions(+)
 create mode 100644 CHANGELOG.md

diff --git a/CHANGELOG.md b/CHANGELOG.md
new file mode 100644
index 0000000..067c46f
--- /dev/null
+++ b/CHANGELOG.md
@@ -0,0 +1,48 @@
+# Changelog
+
+All notable changes to this project are documented in this file.
+
+## [0.3.0] 
+
+### Removed
+
+- **Breaking:** the deprecated `fold_change` output column has been removed. It was an
+  alias for `log2_fold_change` (identical values); use `log2_fold_change` directly. The
+  `FutureWarning` previously emitted on every `pdex(...)` call is also gone.
+
+### Changed
+
+- `"all"` mode (1-vs-rest) is now a genuine one-shot computation instead of a per-group
+  loop: the expression matrix is materialized once, each gene is ranked once via
+  `mwu_one_vs_rest()`, and each group's "rest" pseudobulk/CPM is derived algebraically
+  from the global sum rather than by re-slicing and re-ranking a fresh "rest" matrix.
+  This turns `_pdex_all` from `O(n_groups × n_obs)` into ~`O(n_obs)`, which matters most
+  for screens with many groups (e.g. guides) and/or large cell counts.
+- Bumped the `numba-mwu` dependency floor to `>=0.2.0` (required by the `"all"`-mode
+  optimization above).
+
+### Fixed
+
+- Guarded against floating-point cancellation noise in the algebraic "rest" pseudobulk
+  and CPM derivation (`_pdex_all`): the rest mean of non-negative data can never be
+  legitimately negative, so it is now clipped to `>= 0` before feeding
+  `log2_fold_change`/`percent_change`, preventing spurious negative values from
+  floating-point noise.
+- `_available_cpus()` now checks `hasattr(os, "sched_getaffinity")` before calling it,
+  instead of relying on `AttributeError` from the call itself, fixing a case where the
+  attribute exists but raises for an unrelated reason on some platforms.
+
+### Tests
+
+- Added regression coverage for `"all"` mode with multiple groups and edge cases
+  (`tests/test_math.py`, `tests/test_pdex.py`, `tests/conftest.py`).
+
+### CI
+
+- `semver-check` now only runs on pull requests targeting `main` (previously it could
+  run — and fail spuriously — on unrelated PR bases).
+- CI workflows now avoid redundant duplicate runs on a single PR push.
+
+## [0.2.5] - previous release
+
+See git history prior to this file's introduction for earlier changes.