From 9e864182060c0da8c13376bf4d0abc50b8e3dce9 Mon Sep 17 00:00:00 2001 From: Beatrice Bevilacqua Date: Sat, 11 Jul 2026 02:21:05 +0000 Subject: [PATCH 1/4] feat(ceiling): disjoint-split depth extrapolation (replaces bootstrap) Replace the bootstrap self-split in compute_ceiling with a disjoint (without-replacement) split at multiple depths, extrapolated to full depth. The bootstrap drew both halves from the same cells, so they were not independent, which biased the ceiling in both directions (not a reliable upper bound): the shared cells make the halves over-agree (inflating it), while the duplicate cells over-call the FDR-gated DE metrics and drag the recovery metrics down. A disjoint split is unbiased but shallow (each half <= n/2); measuring each metric across depths and extrapolating to full depth corrects for it. For reliability-like (correlation) metrics this extrapolation reduces to the analytical Spearman-Brown correction (verified: extrap matches SB on those metrics); it generalizes the same idea to metrics with no closed form. - remove _bootstrap_halves; add _disjoint_halves + a per-metric depth-curve fit - compute_ceiling returns one per-metric table (was a (results, agg_results) tuple); writes a single ceiling_results.csv (was ceiling_results.csv + agg_ceiling_results.csv) - CLI --ceiling help, README, CLAUDE updated; add tests/test_ceiling_extrap.py --- CLAUDE.md | 2 +- README.md | 25 +++- src/cell_eval/_cli/_run.py | 7 +- src/cell_eval/_evaluator.py | 280 +++++++++++++++++++++++++---------- tests/test_ceiling_extrap.py | 88 +++++++++++ tests/test_eval.py | 62 +++----- 6 files changed, 333 insertions(+), 131 deletions(-) create mode 100644 tests/test_ceiling_extrap.py diff --git a/CLAUDE.md b/CLAUDE.md index 8603771..970c3b1 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -51,7 +51,7 @@ AnnData inputs (predicted + real) ### Key Abstractions -- **`MetricsEvaluator`** (`src/cell_eval/_evaluator.py`) — Main programmatic entry point. Validates input AnnData objects, computes differential expression via `pdex`, and orchestrates the metric pipeline. `compute_ceiling()` estimates a per-metric data ceiling (upper bound) from the real data alone: it bootstraps each perturbation (and the control) to 2× cells, splits into two equal halves treated as real/pred, and runs the full pipeline (DE computed in-memory). Exposed via the `run --ceiling` / `--ceiling-seed` CLI flags (additive: writes `ceiling_results.csv` / `agg_ceiling_results.csv`). +- **`MetricsEvaluator`** (`src/cell_eval/_evaluator.py`) — Main programmatic entry point. Validates input AnnData objects, computes differential expression via `pdex`, and orchestrates the metric pipeline. `compute_ceiling()` estimates a per-metric data ceiling (upper bound) from the real data alone: it splits the data into two *disjoint* halves (no cell in both) at several depths, runs the full pipeline (DE computed in-memory) on each self-split, and extrapolates each metric's value-vs-depth curve to full depth. Exposed via the `run --ceiling` / `--ceiling-seed` CLI flags (additive: writes a single per-metric `ceiling_results.csv` with the depth curve and the full-depth `extrap`). - **`MetricRegistry`** (`src/cell_eval/metrics/_registry.py`) — Global singleton `metrics_registry`. Metrics are registered with a name, type (`DE` or `ANNDATA_PAIR`), compute function, and best-value indicator. Supports both plain functions and class-based metrics requiring instantiation. diff --git a/README.md b/README.md index 61e8d42..97e3749 100644 --- a/README.md +++ b/README.md @@ -84,10 +84,21 @@ This will give you metric evaluations for each perturbation individually (`resul #### Data ceiling To estimate the *maximum* achievable score on each metric given the noise inherent in the real -data, pass `--ceiling`. This is computed from the **real data only**: per perturbation (and the -control), its cells are bootstrapped to twice their count and split into two equal halves; one -half plays "real" and the other "prediction", and the full metric suite is run on that self-split. -The result is, per metric, an upper bound on how well any model could score on this dataset. +data, pass `--ceiling`. This is computed from the **real data only**: the data is split into two +*disjoint* halves (no cell in both) at several depths, the full metric suite is run on each +self-split, and every metric's value-vs-depth curve is extrapolated to full depth. The result is, +per metric, an unbiased upper bound on how well any model could score on this dataset. + +For reliability-like (correlation) metrics this extrapolation reduces to the analytical Spearman-Brown +correction (`2r/(1+r)`), verified on real data where `extrap` matches Spearman-Brown on the +correlation metrics; it generalizes the same idea to metrics that have no closed form. + +A disjoint split is used rather than a bootstrap self-split: a bootstrap draws the two halves from +the same cells, so they are not independent, which biases the ceiling in *both* directions (so it is +not a reliable upper bound). The shared cells make the halves agree more than two independent +samples would (inflating it), while the duplicate cells over-call the FDR-gated DE metrics and drag +the recovery metrics (recall / overlap / AUC) down. The disjoint split is unbiased but shallow (each +half ≤ `n/2`), which the depth extrapolation corrects. ```bash cell-eval run \ @@ -99,11 +110,11 @@ cell-eval run \ ``` This is *additive*: it writes the normal `results.csv` / `agg_results.csv` **and** -`ceiling_results.csv` / `agg_ceiling_results.csv`. The bootstrap is reproducible via -`--ceiling-seed` (default `0`). From python, call `compute_ceiling` on the evaluator: +`ceiling_results.csv`, a per-metric table with the measured depth curve and the extrapolation +to full depth (`extrap`). The split is reproducible via `--ceiling-seed` (default `0`). From python: ```python -ceiling, ceiling_agg = evaluator.compute_ceiling(seed=0) +ceiling = evaluator.compute_ceiling(seed=0) ``` ### Score diff --git a/src/cell_eval/_cli/_run.py b/src/cell_eval/_cli/_run.py index b62a124..c42bba6 100644 --- a/src/cell_eval/_cli/_run.py +++ b/src/cell_eval/_cli/_run.py @@ -100,14 +100,15 @@ def parse_args_run(parser: ap.ArgumentParser): "--ceiling", action="store_true", help="Additionally compute a data ceiling: a real-data-only upper bound on " - "each metric, estimated by bootstrapping the real data against itself. " - "Writes ceiling_results.csv / agg_ceiling_results.csv alongside the normal results.", + "each metric, estimated by splitting the real data into two disjoint halves " + "at several depths and extrapolating each metric's depth curve to full depth. " + "Writes ceiling_results.csv alongside the normal results.", ) parser.add_argument( "--ceiling-seed", type=int, default=0, - help="Random seed for the data ceiling bootstrap [default: %(default)s]", + help="Random seed for the data ceiling disjoint split [default: %(default)s]", ) parser.add_argument( "--cpm-filter", diff --git a/src/cell_eval/_evaluator.py b/src/cell_eval/_evaluator.py index c0ac1a6..7913b93 100644 --- a/src/cell_eval/_evaluator.py +++ b/src/cell_eval/_evaluator.py @@ -1,7 +1,6 @@ import logging import multiprocessing as mp import os -import warnings from typing import Any, Literal import anndata as ad @@ -157,102 +156,145 @@ def compute_ceiling( profile: Literal["full", "vcc", "minimal", "de", "anndata", "pds"] = "full", metric_configs: dict[str, dict[str, Any]] | None = None, skip_metrics: list[str] | None = None, + fracs: tuple[float, ...] = (1.0, 0.5, 0.25), + seed: int = 0, + agg: Literal["mean", "median"] = "mean", basename: str = "ceiling_results.csv", write_csv: bool = True, break_on_error: bool = False, - seed: int = 0, - ) -> tuple[pl.DataFrame, pl.DataFrame]: + ) -> pl.DataFrame: """Estimate a data ceiling: the maximum achievable score per metric. - Uses the real data only. For each perturbation (and the control) the - cells are bootstrapped to twice their count and split into two equal - halves; one half is treated as "real" and the other as "prediction". - Running the normal metric pipeline on that self-split yields, per metric, - an upper bound on how well any model could score given the noise inherent - in the real data. - - Outputs mirror :meth:`compute` (``ceiling_results.csv`` / - ``agg_ceiling_results.csv``). The bootstrap DE is computed in-memory and - not written to disk. The same ``pdex_kwargs`` and ``allow_discrete`` used - for the main evaluation are reused so the ceiling is directly comparable. + Uses the real data only. The real data is split into two *disjoint* halves + (no cell in both) at several depths, each metric is measured on that + self-split at each depth, and the metric-vs-depth curve is extrapolated to + full depth - an unbiased estimate of the ceiling any model could reach + given the noise inherent in the real data. For reliability-like + (correlation) metrics this reduces to the analytical Spearman-Brown + correction; it generalizes the same idea to metrics with no closed form. + + A disjoint split is used (rather than a bootstrap self-split) because a + bootstrap draws the two halves from the same cells, so they are not + independent - which biases the ceiling in both directions (not a reliable + upper bound): the shared cells make the halves over-agree (inflating it), + while the duplicate cells over-call the FDR-gated DE metrics and drag the + recovery metrics down. The cost of a disjoint split is depth (each half is + at most ``n/2``), which the depth extrapolation corrects for. + + For each ``frac`` in ``fracs`` every perturbation's cells (and the + control's) are shuffled and split without replacement into two halves of + ``floor(frac * n/2)`` cells. ``frac=1`` uses all cells (each half ``n/2``); + the full-depth target (each half ``n``) is ``frac=2``, where the curve is + extrapolated to. The same ``pdex_kwargs`` / ``allow_discrete`` / ``skip_de`` + as the main evaluation are reused so the ceiling is directly comparable, + and the sweep DE is computed in-memory (never written to disk). + + Returns a per-metric table with the measured depth curve (``m@``) + and the extrapolation to full depth (``extrap``, alongside an + ``extrap_linear`` companion, the fit residual and the model used). """ - logger.info(f"Computing data ceiling (seed={seed})") - half_real, half_pred = self._bootstrap_halves(seed) - - ceiling_pair = PerturbationAnndataPair( - real=half_real, - pred=half_pred, - control_pert=self.anndata_pair.control_pert, - pert_col=self.anndata_pair.pert_col, - embed_key=self.anndata_pair.embed_key, - ) + fracs = tuple(sorted(set(fracs), reverse=True)) + if any(f <= 0.0 or f > 1.0 for f in fracs): + raise ValueError(f"fracs must be in (0, 1]; got {fracs}") - if self._skip_de: - ceiling_de = None - else: - ceiling_de = _build_de_comparison( - anndata_pair=ceiling_pair, - num_threads=self._num_threads, - allow_discrete=self._allow_discrete, - outdir=None, # keep the bootstrap DE in-memory; never persisted - prefix=None, - pdex_kwargs=dict(self._pdex_kwargs), + # Fix the perturbation set across depths so the curve isn't confounded by + # small perts dropping out at shallow depths: keep only perts (incl. the + # control) with enough cells to split at the shallowest requested depth. + pert_col = self.anndata_pair.pert_col + control = self.anndata_pair.control_pert + counts = self.anndata_pair.real.obs[pert_col].value_counts() + min_cells = int(np.ceil(2.0 / min(fracs))) + eligible = {str(p) for p, c in counts.items() if c >= min_cells} + if control not in eligible: + raise ValueError( + f"Control '{control}' has too few cells to split at " + f"frac={min(fracs)} (needs >= {min_cells})." + ) + dropped = {str(p) for p in counts.index} - eligible + if dropped: + logger.warning( + f"Depth-extrapolation ceiling: dropping {len(dropped)} perturbation(s) " + f"with < {min_cells} cells (too few to split at frac={min(fracs)})." ) - pipeline = MetricPipeline( - profile=profile, - metric_configs=metric_configs, - break_on_error=break_on_error, - ) - if skip_metrics is not None: - pipeline.skip_metrics(skip_metrics) - pipeline.compute_de_metrics(ceiling_de) - pipeline.compute_anndata_metrics(ceiling_pair) - results = pipeline.get_results() - agg_results = pipeline.get_agg_results() - - if write_csv: - self._write_results(results, agg_results, basename) - - return results, agg_results + curve: dict[float, dict[str, float]] = {} + for frac in fracs: + logger.info(f"Ceiling depth sweep: frac={frac:g} (seed={seed})") + half_real, half_pred = self._disjoint_halves(seed, frac, eligible) + pair = PerturbationAnndataPair( + real=half_real, + pred=half_pred, + control_pert=control, + pert_col=pert_col, + embed_key=self.anndata_pair.embed_key, + ) + de = None + if not self._skip_de: + de = _build_de_comparison( + anndata_pair=pair, + num_threads=self._num_threads, + allow_discrete=self._allow_discrete, + outdir=None, # keep the sweep DE in-memory; never persisted + prefix=None, + pdex_kwargs=dict(self._pdex_kwargs), + ) + pipeline = MetricPipeline( + profile=profile, + metric_configs=metric_configs, + break_on_error=break_on_error, + ) + if skip_metrics is not None: + pipeline.skip_metrics(skip_metrics) + pipeline.compute_de_metrics(de) + pipeline.compute_anndata_metrics(pair) + curve[frac] = _aggregate_metric_values(pipeline.get_results(), agg) - def _bootstrap_halves(self, seed: int) -> tuple[ad.AnnData, ad.AnnData]: - """Build two same-size bootstrap halves of the real data. + table = _extrapolate_ceiling_curve(curve, target_frac=2.0) - Resampling is stratified per perturbation (including the control): each - group's ``n`` cells are drawn ``2n`` times with replacement and split - into two halves of ``n`` cells. This guarantees both halves carry every - perturbation plus the control with the same per-perturbation membership - as the real data, so the resulting ``PerturbationAnndataPair`` validates - and the bootstrap DE keeps the same statistical power. + if write_csv: + prefix = self.prefix.replace("/", "-") if self.prefix is not None else None + outname = basename.replace("/", "-") + outpath = os.path.join( + self.outdir, f"{prefix}_{outname}" if prefix else outname + ) + logger.info(f"Writing depth-extrapolated ceiling to {outpath}") + table.write_csv(outpath) + + return table + + def _disjoint_halves( + self, seed: int, frac: float, eligible: set[str] | None = None + ) -> tuple[ad.AnnData, ad.AnnData]: + """Split the real data into two *disjoint* halves at depth ``frac``. + + Each perturbation's cells are shuffled and split without replacement into + two halves of ``floor(frac * n/2)`` cells each, so no cell appears in both + halves - the independence a bootstrap self-split lacks. Because the shuffle + is seeded per call and the group order is stable, shallower depths are + nested prefixes of deeper ones (monotone subsampling). ``eligible`` (if + given) restricts to a fixed perturbation set so every depth uses the same + perts. """ real = self.anndata_pair.real pert_col = self.anndata_pair.pert_col - rng = np.random.default_rng(seed) - # Group row positions per perturbation in a single pass (`observed=True` - # keeps only perturbations actually present). `.indices` yields positional - # indices, so they can index the AnnData directly. - half_real_idx: list[np.ndarray] = [] - half_pred_idx: list[np.ndarray] = [] - for _pert, idx in real.obs.groupby(pert_col, observed=True).indices.items(): - draws = rng.choice(idx, size=2 * idx.size, replace=True) - half_real_idx.append(draws[: idx.size]) - half_pred_idx.append(draws[idx.size :]) - - # Sampling with replacement duplicates obs names; anndata warns about the - # non-unique index on slice, so silence that one known-benign warning and - # make the names unique immediately afterwards. - with warnings.catch_warnings(): - warnings.filterwarnings( - "ignore", message="Observation names are not unique" - ) - half_real = real[np.concatenate(half_real_idx)].copy() - half_pred = real[np.concatenate(half_pred_idx)].copy() - half_real.obs_names_make_unique() - half_pred.obs_names_make_unique() - + a_idx: list[np.ndarray] = [] + b_idx: list[np.ndarray] = [] + for pert, idx in real.obs.groupby(pert_col, observed=True).indices.items(): + if eligible is not None and str(pert) not in eligible: + continue + idx = np.array(idx) + rng.shuffle(idx) + h = int(frac * (idx.size // 2)) + if h < 1: + continue + a_idx.append(idx[:h]) + b_idx.append(idx[h : 2 * h]) + + # Disjoint split has no duplicate rows, so obs names stay unique. + half_real = real[np.concatenate(a_idx)].copy() + half_pred = real[np.concatenate(b_idx)].copy() return half_real, half_pred def _write_results( @@ -281,6 +323,82 @@ def _write_results( agg_results.write_csv(agg_outpath) +def _aggregate_metric_values(results: pl.DataFrame, agg: str) -> dict[str, float]: + """Collapse the per-perturbation results to one scalar per metric.""" + out: dict[str, float] = {} + if results.is_empty(): + return out + for col in results.columns: + if col == "perturbation" or not results[col].dtype.is_numeric(): + continue + arr = results[col].drop_nulls().to_numpy() + if arr.size == 0: + continue + out[col] = float(np.median(arr) if agg == "median" else np.mean(arr)) + return out + + +def _extrapolate_metric( + fracs: np.ndarray, values: np.ndarray, target: float +) -> dict[str, float | str | None]: + """Extrapolate one metric's depth curve to ``target`` (full depth = 2). + + Primary model is the reliability/attenuation form ``1/m = a + b/frac`` (the + standard measurement-error model, in which reliability grows with depth); it is + only well-defined for reliability-like metrics (``0 < m < 1``). Otherwise we + fall back to a linear fit (the downsampling curves were observed to be ~linear). + """ + mask = np.isfinite(values) + f, m = fracs[mask], values[mask] + out: dict[str, float | str | None] = { + "extrap": None, + "extrap_linear": None, + "resid": None, + "model": None, + } + if f.size < 2: + return out + # Linear fit m = a + b*frac (matches the ~linear downsampling behaviour). + a_lin, b_lin = np.linalg.lstsq(np.vstack([np.ones_like(f), f]).T, m, rcond=None)[0] + out["extrap_linear"] = float(a_lin + b_lin * target) + # Attenuation fit 1/m = a + b/frac; reduces to SB for a single deepest point. + if np.all((m > 0.0) & (m < 1.0)): + design = np.vstack([np.ones_like(f), 1.0 / f]).T + a_att, b_att = np.linalg.lstsq(design, 1.0 / m, rcond=None)[0] + inv = a_att + b_att / target + out["extrap"] = float(1.0 / inv) if inv > 0 else None + out["resid"] = float(np.sqrt(np.mean((design @ [a_att, b_att] - 1.0 / m) ** 2))) + out["model"] = "attenuation" + else: + out["extrap"] = out["extrap_linear"] + out["model"] = "linear" + return out + + +def _extrapolate_ceiling_curve( + curve: dict[float, dict[str, float]], target_frac: float = 2.0 +) -> pl.DataFrame: + """Build the per-metric ceiling table from the measured depth curve. + + Columns: ``metric``, the measured value at each depth (``m@``), the + extrapolation to full depth (``extrap``) plus its ``extrap_linear`` companion, + the fit residual (``resid``) and the model used (``model``). + """ + fracs = sorted(curve.keys(), reverse=True) + metrics = sorted({m for depth in curve.values() for m in depth}) + f_arr = np.array(fracs, dtype=float) + + rows: list[dict[str, Any]] = [] + for metric in metrics: + vals = np.array([curve[fr].get(metric, np.nan) for fr in fracs], dtype=float) + row: dict[str, Any] = {"metric": metric} + for fr in fracs: + row[f"m@{fr:g}"] = curve[fr].get(metric) + row.update(_extrapolate_metric(f_arr, vals, target_frac)) + rows.append(row) + return pl.DataFrame(rows) + + def _build_anndata_pair( real: ad.AnnData | str, pred: ad.AnnData | str, diff --git a/tests/test_ceiling_extrap.py b/tests/test_ceiling_extrap.py new file mode 100644 index 0000000..f0f6dc5 --- /dev/null +++ b/tests/test_ceiling_extrap.py @@ -0,0 +1,88 @@ +import shutil + +import numpy as np +import pytest + +from cell_eval import MetricsEvaluator +from cell_eval._evaluator import _extrapolate_ceiling_curve, _extrapolate_metric +from cell_eval.data import CONTROL_VAR, PERT_COL, build_random_anndata + +OUTDIR = "TEST_OUTPUT_EXTRAP" + + +def test_extrapolation_on_attenuation_curve(): + """A true reliability curve m(f)=f/(f+k) must extrapolate to its full-depth + value 2/(2+k) (frac=2) with a clean attenuation fit.""" + k = 1.0 + fracs = np.array([1.0, 0.5, 0.25]) + values = fracs / (fracs + k) # attenuation / reliability form + + out = _extrapolate_metric(fracs, values, target=2.0) + expected = 2.0 / (2.0 + k) # value at frac=2 (full depth) + + assert out["model"] == "attenuation" + assert out["extrap"] == pytest.approx(expected, abs=1e-6) + assert out["resid"] == pytest.approx(0.0, abs=1e-9) + + +def test_extrapolation_flat_curve_stays_flat(): + """A depth-independent (flat) metric must extrapolate to the same value.""" + fracs = np.array([1.0, 0.5, 0.25]) + values = np.array([0.5, 0.5, 0.5]) + + out = _extrapolate_metric(fracs, values, target=2.0) + extrap = out["extrap"] + assert isinstance(extrap, float) + assert extrap == pytest.approx(0.5, abs=1e-6) + + +def test_extrapolation_linear_fallback_outside_unit_range(): + """Metrics not in (0,1) (e.g. counts) can't use the attenuation form and + fall back to a linear fit.""" + fracs = np.array([1.0, 0.5, 0.25]) + values = np.array([1000.0, 500.0, 250.0]) # count-like, linear in depth + + out = _extrapolate_metric(fracs, values, target=2.0) + assert out["model"] == "linear" + # linear through (0.25,250),(0.5,500),(1,1000) -> slope 1000 -> f=2 : 2000 + assert out["extrap"] == pytest.approx(2000.0, rel=1e-3) + + +def test_extrapolate_ceiling_curve_table_shape(): + curve = { + 1.0: {"pert_r": 0.5, "de_spearman_sig": 0.5}, + 0.5: {"pert_r": 1 / 3, "de_spearman_sig": 0.5}, + 0.25: {"pert_r": 0.2, "de_spearman_sig": 0.5}, + } + table = _extrapolate_ceiling_curve(curve, target_frac=2.0) + + assert set(table["metric"]) == {"pert_r", "de_spearman_sig"} + for col in ("metric", "m@1", "m@0.5", "m@0.25", "extrap", "resid", "model"): + assert col in table.columns + assert "sb" not in table.columns # SB is a validation tool, not shipped output + + +def test_disjoint_halves_share_no_cells(): + adata_real = build_random_anndata() + evaluator = MetricsEvaluator( + adata_pred=adata_real.copy(), + adata_real=adata_real, + control_pert=CONTROL_VAR, + pert_col=PERT_COL, + outdir=OUTDIR, + skip_de=True, + ) + counts = evaluator.anndata_pair.real.obs[PERT_COL].value_counts() + eligible = {str(p) for p in counts.index} + half_real, half_pred = evaluator._disjoint_halves( + seed=0, frac=1.0, eligible=eligible + ) + + # disjoint: no original cell (by name) appears in both halves + assert set(half_real.obs_names).isdisjoint(set(half_pred.obs_names)) + # both halves carry the control + every eligible perturbation + assert CONTROL_VAR in set(half_real.obs[PERT_COL].astype(str)) + assert set(half_real.obs[PERT_COL].astype(str)) == set( + half_pred.obs[PERT_COL].astype(str) + ) + shutil.rmtree(OUTDIR) diff --git a/tests/test_eval.py b/tests/test_eval.py index 9e094aa..c8887df 100644 --- a/tests/test_eval.py +++ b/tests/test_eval.py @@ -405,32 +405,6 @@ def _assert_results_close(a, b) -> None: assert np.allclose(num_a, num_b, equal_nan=True) -def test_ceiling_bootstrap_halves_membership(): - """Each half must carry every perturbation + control with the real counts.""" - adata_real = build_random_anndata() - evaluator = MetricsEvaluator( - adata_pred=adata_real.copy(), - adata_real=adata_real, - control_pert=CONTROL_VAR, - pert_col=PERT_COL, - outdir=OUTDIR, - skip_de=True, # membership only depends on the bootstrap, skip pdex - ) - - half_real, half_pred = evaluator._bootstrap_halves(seed=0) - - real_counts = evaluator.anndata_pair.real.obs[PERT_COL].value_counts().to_dict() - assert CONTROL_VAR in real_counts - for half in (half_real, half_pred): - half_counts = half.obs[PERT_COL].value_counts().to_dict() - # same set of perturbations (incl. control) and same per-pert membership - assert half_counts == real_counts - # sampling with replacement must yield unique obs names - assert half.obs_names.is_unique - - shutil.rmtree(OUTDIR) - - def test_eval_ceiling(): adata_real = build_random_anndata() adata_pred = adata_real.copy() @@ -441,11 +415,11 @@ def test_eval_ceiling(): pert_col=PERT_COL, outdir=OUTDIR, ) - results, agg_results = evaluator.compute_ceiling(break_on_error=True) - assert results.height > 0 - assert agg_results.height > 0 + table = evaluator.compute_ceiling(break_on_error=True) + assert table.height > 0 + for col in ("metric", "extrap"): + assert col in table.columns assert os.path.exists(f"{OUTDIR}/ceiling_results.csv") - assert os.path.exists(f"{OUTDIR}/agg_ceiling_results.csv") shutil.rmtree(OUTDIR) @@ -462,7 +436,6 @@ def test_eval_ceiling_prefix(): ) evaluator.compute_ceiling(break_on_error=True) assert os.path.exists(f"{OUTDIR}/arbitrary_ceiling_results.csv") - assert os.path.exists(f"{OUTDIR}/arbitrary_agg_ceiling_results.csv") shutil.rmtree(OUTDIR) @@ -479,6 +452,7 @@ def test_eval_ceiling_profiles(): for profile in KNOWN_PROFILES: evaluator.compute_ceiling( profile=profile, + fracs=(1.0, 0.5), break_on_error=True, write_csv=False, ) @@ -498,12 +472,13 @@ def test_eval_ceiling_pds_skips_de(): skip_de=True, ) assert evaluator.de_comparison is None - results, _ = evaluator.compute_ceiling( + table = evaluator.compute_ceiling( profile="pds", + fracs=(1.0, 0.5), break_on_error=True, write_csv=False, ) - assert results.height > 0 + assert table.height > 0 shutil.rmtree(OUTDIR) @@ -518,9 +493,19 @@ def test_eval_ceiling_reproducible(): outdir=OUTDIR, num_threads=1, # deterministic reductions ) - r1, _ = evaluator.compute_ceiling(seed=7, write_csv=False, break_on_error=True) - r2, _ = evaluator.compute_ceiling(seed=7, write_csv=False, break_on_error=True) - _assert_results_close(r1, r2) + r1 = evaluator.compute_ceiling( + seed=7, fracs=(1.0, 0.5), write_csv=False, break_on_error=True + ).sort("metric") + r2 = evaluator.compute_ceiling( + seed=7, fracs=(1.0, 0.5), write_csv=False, break_on_error=True + ).sort("metric") + assert r1["metric"].to_list() == r2["metric"].to_list() + for col in ("extrap", "extrap_linear"): + assert np.allclose( + np.array(r1[col].to_list(), dtype=float), + np.array(r2[col].to_list(), dtype=float), + equal_nan=True, + ) shutil.rmtree(OUTDIR) @@ -544,16 +529,15 @@ def test_eval_ceiling_does_not_clobber_de(): with open(pred_de_path, "rb") as fh: before_pred = fh.read() - evaluator.compute_ceiling(seed=0, break_on_error=True) + evaluator.compute_ceiling(seed=0, fracs=(1.0, 0.5), break_on_error=True) with open(real_de_path, "rb") as fh: assert fh.read() == before_real with open(pred_de_path, "rb") as fh: assert fh.read() == before_pred - # ceiling outputs exist, but no ceiling DE artifacts are written + # ceiling output exists, but no ceiling DE artifacts are written assert os.path.exists(f"{OUTDIR}/ceiling_results.csv") - assert os.path.exists(f"{OUTDIR}/agg_ceiling_results.csv") assert not os.path.exists(f"{OUTDIR}/ceiling_real_de.csv") assert not os.path.exists(f"{OUTDIR}/ceiling_pred_de.csv") shutil.rmtree(OUTDIR) From 245d2da15faaa85891b3d683236e15dde971860e Mon Sep 17 00:00:00 2001 From: Leon Hafner Date: Tue, 21 Jul 2026 21:58:02 +0000 Subject: [PATCH 2/4] feat(ceiling): single disjoint split + Spearman-Brown (drop depth extrapolation) Replace the multi-depth reciprocal-fit extrapolation with a single disjoint self-split at n/2 corrected by the analytical Spearman-Brown doubling r' = 2r/(1+r). A holdout backtest across four datasets showed the plain Spearman-Brown correction is both the most accurate estimator and the cheapest (it needs only the deepest split), so the depth sweep and the free-intercept reciprocal fit are removed. - compute_ceiling: one _disjoint_halves(seed) split -> metrics -> Spearman-Brown - correction applied to an explicit SB_METRICS list of reliability metrics; error metrics, unbounded counts, and reliability metrics off the list (clustering_agreement, pearson_edistance) are emitted as NaN - returns (results, agg_results); writes ceiling_results.csv / agg_ceiling_results.csv, matching compute() - tests, README, and CLAUDE.md updated --- CLAUDE.md | 2 +- README.md | 26 +-- src/cell_eval/_cli/_run.py | 5 +- src/cell_eval/_evaluator.py | 315 ++++++++++++++--------------------- tests/test_ceiling.py | 113 +++++++++++++ tests/test_ceiling_extrap.py | 88 ---------- tests/test_eval.py | 40 ++--- 7 files changed, 279 insertions(+), 310 deletions(-) create mode 100644 tests/test_ceiling.py delete mode 100644 tests/test_ceiling_extrap.py diff --git a/CLAUDE.md b/CLAUDE.md index 970c3b1..45e43f1 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -51,7 +51,7 @@ AnnData inputs (predicted + real) ### Key Abstractions -- **`MetricsEvaluator`** (`src/cell_eval/_evaluator.py`) — Main programmatic entry point. Validates input AnnData objects, computes differential expression via `pdex`, and orchestrates the metric pipeline. `compute_ceiling()` estimates a per-metric data ceiling (upper bound) from the real data alone: it splits the data into two *disjoint* halves (no cell in both) at several depths, runs the full pipeline (DE computed in-memory) on each self-split, and extrapolates each metric's value-vs-depth curve to full depth. Exposed via the `run --ceiling` / `--ceiling-seed` CLI flags (additive: writes a single per-metric `ceiling_results.csv` with the depth curve and the full-depth `extrap`). +- **`MetricsEvaluator`** (`src/cell_eval/_evaluator.py`) — Main programmatic entry point. Validates input AnnData objects, computes differential expression via `pdex`, and orchestrates the metric pipeline. `compute_ceiling()` estimates a per-metric data ceiling (upper bound) from the real data alone: it splits the data into two *disjoint* halves of `n/2` cells (no cell in both), runs the full pipeline (DE computed in-memory) on that self-split, and maps the reliability metrics in the explicit `SB_METRICS` list from half depth to full depth with the Spearman-Brown correction `r'=2r/(1+r)` (all other metrics — error, counts, `clustering_agreement`, `pearson_edistance` — emitted as NaN). Returns `(results, agg_results)`. Exposed via the `run --ceiling` / `--ceiling-seed` CLI flags (additive: writes `ceiling_results.csv` / `agg_ceiling_results.csv`). - **`MetricRegistry`** (`src/cell_eval/metrics/_registry.py`) — Global singleton `metrics_registry`. Metrics are registered with a name, type (`DE` or `ANNDATA_PAIR`), compute function, and best-value indicator. Supports both plain functions and class-based metrics requiring instantiation. diff --git a/README.md b/README.md index 97e3749..43a9d92 100644 --- a/README.md +++ b/README.md @@ -84,21 +84,23 @@ This will give you metric evaluations for each perturbation individually (`resul #### Data ceiling To estimate the *maximum* achievable score on each metric given the noise inherent in the real -data, pass `--ceiling`. This is computed from the **real data only**: the data is split into two -*disjoint* halves (no cell in both) at several depths, the full metric suite is run on each -self-split, and every metric's value-vs-depth curve is extrapolated to full depth. The result is, -per metric, an unbiased upper bound on how well any model could score on this dataset. - -For reliability-like (correlation) metrics this extrapolation reduces to the analytical Spearman-Brown -correction (`2r/(1+r)`), verified on real data where `extrap` matches Spearman-Brown on the -correlation metrics; it generalizes the same idea to metrics that have no closed form. +data, pass `--ceiling`. This is computed from the **real data only**: each perturbation's cells +(and the control's) are split into two *disjoint* halves of `n/2` cells (no cell in both), one half +plays "real" and the other "prediction", and the full metric suite is run on that self-split. Each +reliability metric is then mapped from half depth back to full depth by the analytical +Spearman-Brown correction `r' = 2r/(1+r)`. The result is, per metric, an unbiased upper bound on how +well any model could score on this dataset. A disjoint split is used rather than a bootstrap self-split: a bootstrap draws the two halves from the same cells, so they are not independent, which biases the ceiling in *both* directions (so it is not a reliable upper bound). The shared cells make the halves agree more than two independent samples would (inflating it), while the duplicate cells over-call the FDR-gated DE metrics and drag the recovery metrics (recall / overlap / AUC) down. The disjoint split is unbiased but shallow (each -half ≤ `n/2`), which the depth extrapolation corrects. +half `n/2`), which the Spearman-Brown doubling corrects. + +The correction is applied only to a fixed set of reliability metrics (the `SB_METRICS` list in +`_evaluator.py`); every other metric — error metrics, unbounded counts, and reliability metrics +left off that list (`clustering_agreement`, `pearson_edistance`) — is reported as `NaN`. ```bash cell-eval run \ @@ -110,11 +112,11 @@ cell-eval run \ ``` This is *additive*: it writes the normal `results.csv` / `agg_results.csv` **and** -`ceiling_results.csv`, a per-metric table with the measured depth curve and the extrapolation -to full depth (`extrap`). The split is reproducible via `--ceiling-seed` (default `0`). From python: +`ceiling_results.csv` / `agg_ceiling_results.csv`. The split is reproducible via `--ceiling-seed` +(default `0`). From python, call `compute_ceiling` on the evaluator: ```python -ceiling = evaluator.compute_ceiling(seed=0) +ceiling, ceiling_agg = evaluator.compute_ceiling(seed=0) ``` ### Score diff --git a/src/cell_eval/_cli/_run.py b/src/cell_eval/_cli/_run.py index c42bba6..f73b1d8 100644 --- a/src/cell_eval/_cli/_run.py +++ b/src/cell_eval/_cli/_run.py @@ -101,8 +101,9 @@ def parse_args_run(parser: ap.ArgumentParser): action="store_true", help="Additionally compute a data ceiling: a real-data-only upper bound on " "each metric, estimated by splitting the real data into two disjoint halves " - "at several depths and extrapolating each metric's depth curve to full depth. " - "Writes ceiling_results.csv alongside the normal results.", + "of n/2 cells and applying the Spearman-Brown correction (2r/(1+r)) to map " + "each reliability metric back to full depth. Writes ceiling_results.csv / " + "agg_ceiling_results.csv alongside the normal results.", ) parser.add_argument( "--ceiling-seed", diff --git a/src/cell_eval/_evaluator.py b/src/cell_eval/_evaluator.py index 7913b93..f8b1b36 100644 --- a/src/cell_eval/_evaluator.py +++ b/src/cell_eval/_evaluator.py @@ -18,6 +18,38 @@ logger = logging.getLogger(__name__) +# Metrics that receive the Spearman-Brown ceiling correction (r' = 2r/(1+r)): +# the bounded, higher-is-better reliability metrics for which doubling the depth +# is meaningful and empirically accurate. Every OTHER metric in the ceiling output +# - error metrics, unbounded counts, and reliability metrics where the doubling is +# not trustworthy (e.g. clustering_agreement, pearson_edistance) - is emitted as +# NaN. Edit this set to change which metrics are corrected (names must match the +# metric column names produced by the pipeline). +SB_METRICS = frozenset( + { + "pearson_delta", + "discrimination_score_l1", + "discrimination_score_l2", + "discrimination_score_cosine", + "overlap_at_N", + "overlap_at_50", + "overlap_at_100", + "overlap_at_200", + "overlap_at_500", + "precision_at_N", + "precision_at_50", + "precision_at_100", + "precision_at_200", + "precision_at_500", + "de_spearman_sig", + "de_spearman_lfc_sig", + "de_direction_match", + "de_sig_genes_recall", + "pr_auc", + "roc_auc", + } +) + def _available_cpus() -> int: """Return CPUs the current process is allowed to use. @@ -156,124 +188,91 @@ def compute_ceiling( profile: Literal["full", "vcc", "minimal", "de", "anndata", "pds"] = "full", metric_configs: dict[str, dict[str, Any]] | None = None, skip_metrics: list[str] | None = None, - fracs: tuple[float, ...] = (1.0, 0.5, 0.25), - seed: int = 0, - agg: Literal["mean", "median"] = "mean", basename: str = "ceiling_results.csv", write_csv: bool = True, break_on_error: bool = False, - ) -> pl.DataFrame: + seed: int = 0, + ) -> tuple[pl.DataFrame, pl.DataFrame]: """Estimate a data ceiling: the maximum achievable score per metric. - Uses the real data only. The real data is split into two *disjoint* halves - (no cell in both) at several depths, each metric is measured on that - self-split at each depth, and the metric-vs-depth curve is extrapolated to - full depth - an unbiased estimate of the ceiling any model could reach - given the noise inherent in the real data. For reliability-like - (correlation) metrics this reduces to the analytical Spearman-Brown - correction; it generalizes the same idea to metrics with no closed form. - - A disjoint split is used (rather than a bootstrap self-split) because a - bootstrap draws the two halves from the same cells, so they are not - independent - which biases the ceiling in both directions (not a reliable - upper bound): the shared cells make the halves over-agree (inflating it), - while the duplicate cells over-call the FDR-gated DE metrics and drag the - recovery metrics down. The cost of a disjoint split is depth (each half is - at most ``n/2``), which the depth extrapolation corrects for. - - For each ``frac`` in ``fracs`` every perturbation's cells (and the - control's) are shuffled and split without replacement into two halves of - ``floor(frac * n/2)`` cells. ``frac=1`` uses all cells (each half ``n/2``); - the full-depth target (each half ``n``) is ``frac=2``, where the curve is - extrapolated to. The same ``pdex_kwargs`` / ``allow_discrete`` / ``skip_de`` - as the main evaluation are reused so the ceiling is directly comparable, - and the sweep DE is computed in-memory (never written to disk). - - Returns a per-metric table with the measured depth curve (``m@``) - and the extrapolation to full depth (``extrap``, alongside an - ``extrap_linear`` companion, the fit residual and the model used). + Uses the real data only. Each perturbation's cells (and the control's) are + split into two *disjoint* halves of ``n/2`` cells - no cell in both - and + one half is treated as "real", the other as "prediction". Running the + normal metric pipeline on that self-split measures each metric's + reliability at half depth; the Spearman-Brown correction ``r' = 2r/(1+r)`` + then maps it to full depth (``n``), an unbiased upper bound on how well any + model could score given the noise inherent in the real data. + + A *disjoint* split is used (rather than a bootstrap self-split) because a + bootstrap draws both halves from the same cells, so they are not + independent - which biases the ceiling in both directions: shared cells + make the halves over-agree (inflating it), while duplicate cells over-call + the FDR-gated DE metrics and drag the recovery metrics down. The cost of a + disjoint split is depth (each half is ``n/2``), which the Spearman-Brown + doubling corrects for. + + The correction is applied only to the reliability metrics listed in the + module-level ``SB_METRICS`` set (bounded, higher-is-better, and empirically + well-behaved under doubling). Every other metric - error metrics, unbounded + counts, and reliability metrics left off that list (``clustering_agreement``, + ``pearson_edistance``) - is emitted as ``NaN`` (no defensible ceiling). Outputs + mirror :meth:`compute` (``ceiling_results.csv`` / + ``agg_ceiling_results.csv``); the self-split DE is computed in-memory and + never written. The same ``pdex_kwargs`` and ``allow_discrete`` as the main + evaluation are reused so the ceiling is directly comparable. """ - fracs = tuple(sorted(set(fracs), reverse=True)) - if any(f <= 0.0 or f > 1.0 for f in fracs): - raise ValueError(f"fracs must be in (0, 1]; got {fracs}") + logger.info(f"Computing data ceiling (seed={seed})") + half_real, half_pred = self._disjoint_halves(seed) + + ceiling_pair = PerturbationAnndataPair( + real=half_real, + pred=half_pred, + control_pert=self.anndata_pair.control_pert, + pert_col=self.anndata_pair.pert_col, + embed_key=self.anndata_pair.embed_key, + ) - # Fix the perturbation set across depths so the curve isn't confounded by - # small perts dropping out at shallow depths: keep only perts (incl. the - # control) with enough cells to split at the shallowest requested depth. - pert_col = self.anndata_pair.pert_col - control = self.anndata_pair.control_pert - counts = self.anndata_pair.real.obs[pert_col].value_counts() - min_cells = int(np.ceil(2.0 / min(fracs))) - eligible = {str(p) for p, c in counts.items() if c >= min_cells} - if control not in eligible: - raise ValueError( - f"Control '{control}' has too few cells to split at " - f"frac={min(fracs)} (needs >= {min_cells})." - ) - dropped = {str(p) for p in counts.index} - eligible - if dropped: - logger.warning( - f"Depth-extrapolation ceiling: dropping {len(dropped)} perturbation(s) " - f"with < {min_cells} cells (too few to split at frac={min(fracs)})." + if self._skip_de: + ceiling_de = None + else: + ceiling_de = _build_de_comparison( + anndata_pair=ceiling_pair, + num_threads=self._num_threads, + allow_discrete=self._allow_discrete, + outdir=None, # keep the self-split DE in-memory; never persisted + prefix=None, + pdex_kwargs=dict(self._pdex_kwargs), ) - curve: dict[float, dict[str, float]] = {} - for frac in fracs: - logger.info(f"Ceiling depth sweep: frac={frac:g} (seed={seed})") - half_real, half_pred = self._disjoint_halves(seed, frac, eligible) - pair = PerturbationAnndataPair( - real=half_real, - pred=half_pred, - control_pert=control, - pert_col=pert_col, - embed_key=self.anndata_pair.embed_key, - ) - de = None - if not self._skip_de: - de = _build_de_comparison( - anndata_pair=pair, - num_threads=self._num_threads, - allow_discrete=self._allow_discrete, - outdir=None, # keep the sweep DE in-memory; never persisted - prefix=None, - pdex_kwargs=dict(self._pdex_kwargs), - ) - pipeline = MetricPipeline( - profile=profile, - metric_configs=metric_configs, - break_on_error=break_on_error, - ) - if skip_metrics is not None: - pipeline.skip_metrics(skip_metrics) - pipeline.compute_de_metrics(de) - pipeline.compute_anndata_metrics(pair) - curve[frac] = _aggregate_metric_values(pipeline.get_results(), agg) + pipeline = MetricPipeline( + profile=profile, + metric_configs=metric_configs, + break_on_error=break_on_error, + ) + if skip_metrics is not None: + pipeline.skip_metrics(skip_metrics) + pipeline.compute_de_metrics(ceiling_de) + pipeline.compute_anndata_metrics(ceiling_pair) - table = _extrapolate_ceiling_curve(curve, target_frac=2.0) + # Half-depth self-split scores -> full-depth ceiling via Spearman-Brown. + results = _spearman_brown_correct(pipeline.get_results()) + agg_results = results.drop("perturbation").describe() if write_csv: - prefix = self.prefix.replace("/", "-") if self.prefix is not None else None - outname = basename.replace("/", "-") - outpath = os.path.join( - self.outdir, f"{prefix}_{outname}" if prefix else outname - ) - logger.info(f"Writing depth-extrapolated ceiling to {outpath}") - table.write_csv(outpath) - - return table - - def _disjoint_halves( - self, seed: int, frac: float, eligible: set[str] | None = None - ) -> tuple[ad.AnnData, ad.AnnData]: - """Split the real data into two *disjoint* halves at depth ``frac``. - - Each perturbation's cells are shuffled and split without replacement into - two halves of ``floor(frac * n/2)`` cells each, so no cell appears in both - halves - the independence a bootstrap self-split lacks. Because the shuffle - is seeded per call and the group order is stable, shallower depths are - nested prefixes of deeper ones (monotone subsampling). ``eligible`` (if - given) restricts to a fixed perturbation set so every depth uses the same - perts. + self._write_results(results, agg_results, basename) + + return results, agg_results + + def _disjoint_halves(self, seed: int) -> tuple[ad.AnnData, ad.AnnData]: + """Split the real data into two *disjoint* halves of ``n/2`` cells each. + + Each perturbation's cells (including the control's) are shuffled and split + without replacement into two halves of ``floor(n/2)`` cells - so no cell + appears in both halves, giving the independence a bootstrap self-split + lacks. Perturbations with fewer than 2 cells cannot be split and are + dropped from both halves. The resulting half depth (``n/2``) is corrected + back to full depth by the Spearman-Brown doubling in + :meth:`compute_ceiling`. """ real = self.anndata_pair.real pert_col = self.anndata_pair.pert_col @@ -281,16 +280,13 @@ def _disjoint_halves( a_idx: list[np.ndarray] = [] b_idx: list[np.ndarray] = [] - for pert, idx in real.obs.groupby(pert_col, observed=True).indices.items(): - if eligible is not None and str(pert) not in eligible: - continue - idx = np.array(idx) - rng.shuffle(idx) - h = int(frac * (idx.size // 2)) + for _pert, idx in real.obs.groupby(pert_col, observed=True).indices.items(): + perm = rng.permutation(np.asarray(idx)) + h = perm.size // 2 if h < 1: - continue - a_idx.append(idx[:h]) - b_idx.append(idx[h : 2 * h]) + continue # < 2 cells: cannot form two disjoint halves + a_idx.append(perm[:h]) + b_idx.append(perm[h : 2 * h]) # Disjoint split has no duplicate rows, so obs names stay unique. half_real = real[np.concatenate(a_idx)].copy() @@ -323,80 +319,25 @@ def _write_results( agg_results.write_csv(agg_outpath) -def _aggregate_metric_values(results: pl.DataFrame, agg: str) -> dict[str, float]: - """Collapse the per-perturbation results to one scalar per metric.""" - out: dict[str, float] = {} - if results.is_empty(): - return out - for col in results.columns: - if col == "perturbation" or not results[col].dtype.is_numeric(): - continue - arr = results[col].drop_nulls().to_numpy() - if arr.size == 0: - continue - out[col] = float(np.median(arr) if agg == "median" else np.mean(arr)) - return out - - -def _extrapolate_metric( - fracs: np.ndarray, values: np.ndarray, target: float -) -> dict[str, float | str | None]: - """Extrapolate one metric's depth curve to ``target`` (full depth = 2). +def _spearman_brown_correct(results: pl.DataFrame) -> pl.DataFrame: + """Map half-depth self-split scores to the full-depth ceiling. - Primary model is the reliability/attenuation form ``1/m = a + b/frac`` (the - standard measurement-error model, in which reliability grows with depth); it is - only well-defined for reliability-like metrics (``0 < m < 1``). Otherwise we - fall back to a linear fit (the downsampling curves were observed to be ~linear). + Applies the Spearman-Brown prophecy ``r' = 2r/(1+r)`` - the reliability of a + test of doubled length - to the reliability metrics listed in ``SB_METRICS``. + Every other column (error metrics, unbounded counts, and reliability metrics + not in that set) is emitted as ``NaN``, since a Spearman-Brown ceiling has no + defensible meaning there. """ - mask = np.isfinite(values) - f, m = fracs[mask], values[mask] - out: dict[str, float | str | None] = { - "extrap": None, - "extrap_linear": None, - "resid": None, - "model": None, - } - if f.size < 2: - return out - # Linear fit m = a + b*frac (matches the ~linear downsampling behaviour). - a_lin, b_lin = np.linalg.lstsq(np.vstack([np.ones_like(f), f]).T, m, rcond=None)[0] - out["extrap_linear"] = float(a_lin + b_lin * target) - # Attenuation fit 1/m = a + b/frac; reduces to SB for a single deepest point. - if np.all((m > 0.0) & (m < 1.0)): - design = np.vstack([np.ones_like(f), 1.0 / f]).T - a_att, b_att = np.linalg.lstsq(design, 1.0 / m, rcond=None)[0] - inv = a_att + b_att / target - out["extrap"] = float(1.0 / inv) if inv > 0 else None - out["resid"] = float(np.sqrt(np.mean((design @ [a_att, b_att] - 1.0 / m) ** 2))) - out["model"] = "attenuation" - else: - out["extrap"] = out["extrap_linear"] - out["model"] = "linear" - return out - - -def _extrapolate_ceiling_curve( - curve: dict[float, dict[str, float]], target_frac: float = 2.0 -) -> pl.DataFrame: - """Build the per-metric ceiling table from the measured depth curve. - - Columns: ``metric``, the measured value at each depth (``m@``), the - extrapolation to full depth (``extrap``) plus its ``extrap_linear`` companion, - the fit residual (``resid``) and the model used (``model``). - """ - fracs = sorted(curve.keys(), reverse=True) - metrics = sorted({m for depth in curve.values() for m in depth}) - f_arr = np.array(fracs, dtype=float) - - rows: list[dict[str, Any]] = [] - for metric in metrics: - vals = np.array([curve[fr].get(metric, np.nan) for fr in fracs], dtype=float) - row: dict[str, Any] = {"metric": metric} - for fr in fracs: - row[f"m@{fr:g}"] = curve[fr].get(metric) - row.update(_extrapolate_metric(f_arr, vals, target_frac)) - rows.append(row) - return pl.DataFrame(rows) + nan = float("nan") + exprs: list[pl.Expr] = [] + for col in results.columns: + if col == "perturbation": + continue + if col in SB_METRICS: + exprs.append((2.0 * pl.col(col) / (1.0 + pl.col(col))).alias(col)) + else: + exprs.append(pl.lit(nan).alias(col)) + return results.with_columns(exprs) if exprs else results def _build_anndata_pair( diff --git a/tests/test_ceiling.py b/tests/test_ceiling.py new file mode 100644 index 0000000..f83d68d --- /dev/null +++ b/tests/test_ceiling.py @@ -0,0 +1,113 @@ +import shutil + +import numpy as np +import polars as pl +import pytest + +from cell_eval import MetricsEvaluator +from cell_eval._evaluator import SB_METRICS, _spearman_brown_correct +from cell_eval.data import CONTROL_VAR, PERT_COL, build_random_anndata + +OUTDIR = "TEST_OUTPUT_CEILING" + + +def test_spearman_brown_doubling_on_reliability_metrics(): + """SB doubling r' = 2r/(1+r) is applied to reliability (best_value == ONE) + metric columns.""" + df = pl.DataFrame( + { + "perturbation": ["a", "b"], + "pearson_delta": [0.5, 1.0 / 3.0], + "overlap_at_N": [0.6, 0.2], + } + ) + out = _spearman_brown_correct(df) + # 2*0.5/(1+0.5)=2/3 ; 2*(1/3)/(1+1/3)=0.5 + assert out["pearson_delta"].to_list() == pytest.approx([2 / 3, 0.5], abs=1e-6) + # 2*0.6/1.6=0.75 ; 2*0.2/1.2=1/3 + assert out["overlap_at_N"].to_list() == pytest.approx([0.75, 1 / 3], abs=1e-6) + + +def test_non_reliability_and_excluded_metrics_are_nan(): + """Error metrics, unbounded counts, and the excluded reliabilities + (clustering_agreement, pearson_edistance) are emitted as NaN.""" + df = pl.DataFrame( + { + "perturbation": ["a"], + "pearson_delta": [0.5], # reliability -> SB + "mse": [0.01], # error metric -> NaN + "mae": [0.02], # error metric -> NaN + "de_nsig_counts_real": [10.0], # unbounded count -> NaN + "clustering_agreement": [0.4], # excluded reliability -> NaN + "pearson_edistance": [0.9], # excluded reliability -> NaN + } + ) + out = _spearman_brown_correct(df) + assert out["pearson_delta"][0] == pytest.approx(2 / 3, abs=1e-6) + for col in ( + "mse", + "mae", + "de_nsig_counts_real", + "clustering_agreement", + "pearson_edistance", + ): + assert np.isnan(out[col][0]), col + # SB_METRICS is the explicit inclusion list; the excluded ones are absent + assert "pearson_delta" in SB_METRICS + assert "clustering_agreement" not in SB_METRICS + assert "pearson_edistance" not in SB_METRICS + + +def test_disjoint_halves_share_no_cells(): + adata_real = build_random_anndata() + evaluator = MetricsEvaluator( + adata_pred=adata_real.copy(), + adata_real=adata_real, + control_pert=CONTROL_VAR, + pert_col=PERT_COL, + outdir=OUTDIR, + skip_de=True, + ) + half_real, half_pred = evaluator._disjoint_halves(seed=0) + + # disjoint: no original cell (by name) appears in both halves + assert set(half_real.obs_names).isdisjoint(set(half_pred.obs_names)) + # both halves carry the control + every (splittable) perturbation + assert CONTROL_VAR in set(half_real.obs[PERT_COL].astype(str)) + assert set(half_real.obs[PERT_COL].astype(str)) == set( + half_pred.obs[PERT_COL].astype(str) + ) + shutil.rmtree(OUTDIR) + + +def test_compute_ceiling_end_to_end(): + """compute_ceiling returns (results, agg): reliability metrics are SB-corrected + and bounded in [0, 1]; error metrics come back as NaN.""" + adata_real = build_random_anndata() + evaluator = MetricsEvaluator( + adata_pred=adata_real.copy(), + adata_real=adata_real, + control_pert=CONTROL_VAR, + pert_col=PERT_COL, + outdir=OUTDIR, + skip_de=True, + ) + results, agg = evaluator.compute_ceiling( + profile="anndata", write_csv=False, break_on_error=True + ) + assert results.height > 0 + assert "perturbation" in results.columns + + assert "pearson_delta" in results.columns + pv = results["pearson_delta"].drop_nulls().to_numpy() + assert np.all(pv <= 1.0 + 1e-9) # SB doubling can never exceed 1 + + # error metrics are emitted as NaN (no defensible SB ceiling) + for col in ("mse", "mae"): + if col in results.columns: + assert np.all(np.isnan(results[col].to_numpy())) + # excluded reliability metrics are NaN too + for col in ("clustering_agreement", "pearson_edistance"): + if col in results.columns: + assert np.all(np.isnan(results[col].to_numpy())) + shutil.rmtree(OUTDIR) diff --git a/tests/test_ceiling_extrap.py b/tests/test_ceiling_extrap.py deleted file mode 100644 index f0f6dc5..0000000 --- a/tests/test_ceiling_extrap.py +++ /dev/null @@ -1,88 +0,0 @@ -import shutil - -import numpy as np -import pytest - -from cell_eval import MetricsEvaluator -from cell_eval._evaluator import _extrapolate_ceiling_curve, _extrapolate_metric -from cell_eval.data import CONTROL_VAR, PERT_COL, build_random_anndata - -OUTDIR = "TEST_OUTPUT_EXTRAP" - - -def test_extrapolation_on_attenuation_curve(): - """A true reliability curve m(f)=f/(f+k) must extrapolate to its full-depth - value 2/(2+k) (frac=2) with a clean attenuation fit.""" - k = 1.0 - fracs = np.array([1.0, 0.5, 0.25]) - values = fracs / (fracs + k) # attenuation / reliability form - - out = _extrapolate_metric(fracs, values, target=2.0) - expected = 2.0 / (2.0 + k) # value at frac=2 (full depth) - - assert out["model"] == "attenuation" - assert out["extrap"] == pytest.approx(expected, abs=1e-6) - assert out["resid"] == pytest.approx(0.0, abs=1e-9) - - -def test_extrapolation_flat_curve_stays_flat(): - """A depth-independent (flat) metric must extrapolate to the same value.""" - fracs = np.array([1.0, 0.5, 0.25]) - values = np.array([0.5, 0.5, 0.5]) - - out = _extrapolate_metric(fracs, values, target=2.0) - extrap = out["extrap"] - assert isinstance(extrap, float) - assert extrap == pytest.approx(0.5, abs=1e-6) - - -def test_extrapolation_linear_fallback_outside_unit_range(): - """Metrics not in (0,1) (e.g. counts) can't use the attenuation form and - fall back to a linear fit.""" - fracs = np.array([1.0, 0.5, 0.25]) - values = np.array([1000.0, 500.0, 250.0]) # count-like, linear in depth - - out = _extrapolate_metric(fracs, values, target=2.0) - assert out["model"] == "linear" - # linear through (0.25,250),(0.5,500),(1,1000) -> slope 1000 -> f=2 : 2000 - assert out["extrap"] == pytest.approx(2000.0, rel=1e-3) - - -def test_extrapolate_ceiling_curve_table_shape(): - curve = { - 1.0: {"pert_r": 0.5, "de_spearman_sig": 0.5}, - 0.5: {"pert_r": 1 / 3, "de_spearman_sig": 0.5}, - 0.25: {"pert_r": 0.2, "de_spearman_sig": 0.5}, - } - table = _extrapolate_ceiling_curve(curve, target_frac=2.0) - - assert set(table["metric"]) == {"pert_r", "de_spearman_sig"} - for col in ("metric", "m@1", "m@0.5", "m@0.25", "extrap", "resid", "model"): - assert col in table.columns - assert "sb" not in table.columns # SB is a validation tool, not shipped output - - -def test_disjoint_halves_share_no_cells(): - adata_real = build_random_anndata() - evaluator = MetricsEvaluator( - adata_pred=adata_real.copy(), - adata_real=adata_real, - control_pert=CONTROL_VAR, - pert_col=PERT_COL, - outdir=OUTDIR, - skip_de=True, - ) - counts = evaluator.anndata_pair.real.obs[PERT_COL].value_counts() - eligible = {str(p) for p in counts.index} - half_real, half_pred = evaluator._disjoint_halves( - seed=0, frac=1.0, eligible=eligible - ) - - # disjoint: no original cell (by name) appears in both halves - assert set(half_real.obs_names).isdisjoint(set(half_pred.obs_names)) - # both halves carry the control + every eligible perturbation - assert CONTROL_VAR in set(half_real.obs[PERT_COL].astype(str)) - assert set(half_real.obs[PERT_COL].astype(str)) == set( - half_pred.obs[PERT_COL].astype(str) - ) - shutil.rmtree(OUTDIR) diff --git a/tests/test_eval.py b/tests/test_eval.py index c8887df..2684031 100644 --- a/tests/test_eval.py +++ b/tests/test_eval.py @@ -415,11 +415,15 @@ def test_eval_ceiling(): pert_col=PERT_COL, outdir=OUTDIR, ) - table = evaluator.compute_ceiling(break_on_error=True) - assert table.height > 0 - for col in ("metric", "extrap"): - assert col in table.columns + results, agg = evaluator.compute_ceiling(break_on_error=True) + assert results.height > 0 + assert "perturbation" in results.columns + # a reliability metric is SB-corrected; the doubling can never exceed 1 + assert "pearson_delta" in results.columns + pv = results["pearson_delta"].drop_nulls().to_numpy() + assert np.all(pv <= 1.0 + 1e-9) assert os.path.exists(f"{OUTDIR}/ceiling_results.csv") + assert os.path.exists(f"{OUTDIR}/agg_ceiling_results.csv") shutil.rmtree(OUTDIR) @@ -452,7 +456,6 @@ def test_eval_ceiling_profiles(): for profile in KNOWN_PROFILES: evaluator.compute_ceiling( profile=profile, - fracs=(1.0, 0.5), break_on_error=True, write_csv=False, ) @@ -472,13 +475,12 @@ def test_eval_ceiling_pds_skips_de(): skip_de=True, ) assert evaluator.de_comparison is None - table = evaluator.compute_ceiling( + results, _agg = evaluator.compute_ceiling( profile="pds", - fracs=(1.0, 0.5), break_on_error=True, write_csv=False, ) - assert table.height > 0 + assert results.height > 0 shutil.rmtree(OUTDIR) @@ -493,17 +495,15 @@ def test_eval_ceiling_reproducible(): outdir=OUTDIR, num_threads=1, # deterministic reductions ) - r1 = evaluator.compute_ceiling( - seed=7, fracs=(1.0, 0.5), write_csv=False, break_on_error=True - ).sort("metric") - r2 = evaluator.compute_ceiling( - seed=7, fracs=(1.0, 0.5), write_csv=False, break_on_error=True - ).sort("metric") - assert r1["metric"].to_list() == r2["metric"].to_list() - for col in ("extrap", "extrap_linear"): - assert np.allclose( - np.array(r1[col].to_list(), dtype=float), - np.array(r2[col].to_list(), dtype=float), + r1, _ = evaluator.compute_ceiling(seed=7, write_csv=False, break_on_error=True) + r2, _ = evaluator.compute_ceiling(seed=7, write_csv=False, break_on_error=True) + r1 = r1.sort("perturbation") + r2 = r2.sort("perturbation") + assert r1["perturbation"].to_list() == r2["perturbation"].to_list() + for col in (c for c in r1.columns if c != "perturbation"): + np.testing.assert_allclose( + r1[col].to_numpy().astype(float), + r2[col].to_numpy().astype(float), equal_nan=True, ) shutil.rmtree(OUTDIR) @@ -529,7 +529,7 @@ def test_eval_ceiling_does_not_clobber_de(): with open(pred_de_path, "rb") as fh: before_pred = fh.read() - evaluator.compute_ceiling(seed=0, fracs=(1.0, 0.5), break_on_error=True) + evaluator.compute_ceiling(seed=0, break_on_error=True) with open(real_de_path, "rb") as fh: assert fh.read() == before_real From 81e4af0f739c030ffeb3e40829bde4bd45ab603b Mon Sep 17 00:00:00 2001 From: Leon Hafner Date: Tue, 21 Jul 2026 22:21:09 +0000 Subject: [PATCH 3/4] chore: bump version to 0.8.2 --- pyproject.toml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/pyproject.toml b/pyproject.toml index c256a3c..fb653db 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "cell-eval" -version = "0.8.1" +version = "0.8.2" description = "Evaluation metrics for single-cell perturbation predictions" readme = "README.md" authors = [ From 459b52a82633a51598598605f41060a0654e80b4 Mon Sep 17 00:00:00 2001 From: Leon Hafner Date: Thu, 23 Jul 2026 17:15:54 +0000 Subject: [PATCH 4/4] refactor(ceiling): apply Spearman-Brown to the per-context aggregate Average each metric over perturbations first, then apply r' = 2r/(1+r) to that mean (SB-of-mean) - matching the validated backtest - instead of correcting each perturbation individually and then aggregating (mean-of-SB); the two differ by a Jensen gap for the concave SB transform. ceiling_results.csv now holds the raw per-perturbation self-split; agg_ceiling_results.csv holds the SB-corrected per-metric ceiling. --- CLAUDE.md | 2 +- README.md | 13 +++++++------ src/cell_eval/_evaluator.py | 28 ++++++++++++++++------------ tests/test_ceiling.py | 26 +++++++++++++------------- tests/test_eval.py | 9 +++++---- 5 files changed, 42 insertions(+), 36 deletions(-) diff --git a/CLAUDE.md b/CLAUDE.md index 45e43f1..568001d 100644 --- a/CLAUDE.md +++ b/CLAUDE.md @@ -51,7 +51,7 @@ AnnData inputs (predicted + real) ### Key Abstractions -- **`MetricsEvaluator`** (`src/cell_eval/_evaluator.py`) — Main programmatic entry point. Validates input AnnData objects, computes differential expression via `pdex`, and orchestrates the metric pipeline. `compute_ceiling()` estimates a per-metric data ceiling (upper bound) from the real data alone: it splits the data into two *disjoint* halves of `n/2` cells (no cell in both), runs the full pipeline (DE computed in-memory) on that self-split, and maps the reliability metrics in the explicit `SB_METRICS` list from half depth to full depth with the Spearman-Brown correction `r'=2r/(1+r)` (all other metrics — error, counts, `clustering_agreement`, `pearson_edistance` — emitted as NaN). Returns `(results, agg_results)`. Exposed via the `run --ceiling` / `--ceiling-seed` CLI flags (additive: writes `ceiling_results.csv` / `agg_ceiling_results.csv`). +- **`MetricsEvaluator`** (`src/cell_eval/_evaluator.py`) — Main programmatic entry point. Validates input AnnData objects, computes differential expression via `pdex`, and orchestrates the metric pipeline. `compute_ceiling()` estimates a per-metric data ceiling (upper bound) from the real data alone: it splits the data into two *disjoint* halves of `n/2` cells (no cell in both), runs the full pipeline (DE computed in-memory) on that self-split, and averages each metric over perturbations and maps the reliability metrics in the explicit `SB_METRICS` list from half depth to full depth with the Spearman-Brown correction `r'=2r/(1+r)` (all other metrics — error, counts, `clustering_agreement`, `pearson_edistance` — emitted as NaN). Returns `(results, agg_results)`: the raw per-perturbation self-split and the SB-corrected per-metric ceiling. Exposed via the `run --ceiling` / `--ceiling-seed` CLI flags (additive: writes `ceiling_results.csv` / `agg_ceiling_results.csv`). - **`MetricRegistry`** (`src/cell_eval/metrics/_registry.py`) — Global singleton `metrics_registry`. Metrics are registered with a name, type (`DE` or `ANNDATA_PAIR`), compute function, and best-value indicator. Supports both plain functions and class-based metrics requiring instantiation. diff --git a/README.md b/README.md index 43a9d92..b9836fe 100644 --- a/README.md +++ b/README.md @@ -86,10 +86,10 @@ This will give you metric evaluations for each perturbation individually (`resul To estimate the *maximum* achievable score on each metric given the noise inherent in the real data, pass `--ceiling`. This is computed from the **real data only**: each perturbation's cells (and the control's) are split into two *disjoint* halves of `n/2` cells (no cell in both), one half -plays "real" and the other "prediction", and the full metric suite is run on that self-split. Each -reliability metric is then mapped from half depth back to full depth by the analytical -Spearman-Brown correction `r' = 2r/(1+r)`. The result is, per metric, an unbiased upper bound on how -well any model could score on this dataset. +plays "real" and the other "prediction", and the full metric suite is run on that self-split. +Averaging each metric over perturbations and applying the analytical Spearman-Brown correction +`r' = 2r/(1+r)` maps that per-context mean from half depth back to full depth. The result is, per +metric, an unbiased upper bound on how well any model could score on this dataset. A disjoint split is used rather than a bootstrap self-split: a bootstrap draws the two halves from the same cells, so they are not independent, which biases the ceiling in *both* directions (so it is @@ -112,8 +112,9 @@ cell-eval run \ ``` This is *additive*: it writes the normal `results.csv` / `agg_results.csv` **and** -`ceiling_results.csv` / `agg_ceiling_results.csv`. The split is reproducible via `--ceiling-seed` -(default `0`). From python, call `compute_ceiling` on the evaluator: +`ceiling_results.csv` (the raw per-perturbation self-split) / `agg_ceiling_results.csv` (the +SB-corrected per-metric ceiling). The split is reproducible via `--ceiling-seed` (default `0`). From +python, call `compute_ceiling` on the evaluator: ```python ceiling, ceiling_agg = evaluator.compute_ceiling(seed=0) diff --git a/src/cell_eval/_evaluator.py b/src/cell_eval/_evaluator.py index f8b1b36..4bcf711 100644 --- a/src/cell_eval/_evaluator.py +++ b/src/cell_eval/_evaluator.py @@ -198,10 +198,11 @@ def compute_ceiling( Uses the real data only. Each perturbation's cells (and the control's) are split into two *disjoint* halves of ``n/2`` cells - no cell in both - and one half is treated as "real", the other as "prediction". Running the - normal metric pipeline on that self-split measures each metric's - reliability at half depth; the Spearman-Brown correction ``r' = 2r/(1+r)`` - then maps it to full depth (``n``), an unbiased upper bound on how well any - model could score given the noise inherent in the real data. + normal metric pipeline on that self-split measures each metric per + perturbation at half depth; averaging over perturbations and applying the + Spearman-Brown correction ``r' = 2r/(1+r)`` maps that per-context mean to + full depth (``n``), an unbiased upper bound on how well any model could + score given the noise inherent in the real data. A *disjoint* split is used (rather than a bootstrap self-split) because a bootstrap draws both halves from the same cells, so they are not @@ -215,11 +216,12 @@ def compute_ceiling( module-level ``SB_METRICS`` set (bounded, higher-is-better, and empirically well-behaved under doubling). Every other metric - error metrics, unbounded counts, and reliability metrics left off that list (``clustering_agreement``, - ``pearson_edistance``) - is emitted as ``NaN`` (no defensible ceiling). Outputs - mirror :meth:`compute` (``ceiling_results.csv`` / - ``agg_ceiling_results.csv``); the self-split DE is computed in-memory and - never written. The same ``pdex_kwargs`` and ``allow_discrete`` as the main - evaluation are reused so the ceiling is directly comparable. + ``pearson_edistance``) - is emitted as ``NaN`` (no defensible ceiling). + ``ceiling_results.csv`` holds the raw per-perturbation self-split scores; + ``agg_ceiling_results.csv`` holds the SB-corrected per-metric ceiling. The + self-split DE is computed in-memory and never written. The same + ``pdex_kwargs`` and ``allow_discrete`` as the main evaluation are reused so + the ceiling is directly comparable. """ logger.info(f"Computing data ceiling (seed={seed})") half_real, half_pred = self._disjoint_halves(seed) @@ -254,9 +256,11 @@ def compute_ceiling( pipeline.compute_de_metrics(ceiling_de) pipeline.compute_anndata_metrics(ceiling_pair) - # Half-depth self-split scores -> full-depth ceiling via Spearman-Brown. - results = _spearman_brown_correct(pipeline.get_results()) - agg_results = results.drop("perturbation").describe() + # Spearman-Brown ceiling on the per-context AGGREGATE: average each metric + # over perturbations, then map that mean from half depth to full depth with + # r' = 2r/(1+r). results keeps the raw per-perturbation self-split scores. + results = pipeline.get_results() + agg_results = _spearman_brown_correct(results.drop("perturbation").mean()) if write_csv: self._write_results(results, agg_results, basename) diff --git a/tests/test_ceiling.py b/tests/test_ceiling.py index f83d68d..f8b0703 100644 --- a/tests/test_ceiling.py +++ b/tests/test_ceiling.py @@ -81,8 +81,9 @@ def test_disjoint_halves_share_no_cells(): def test_compute_ceiling_end_to_end(): - """compute_ceiling returns (results, agg): reliability metrics are SB-corrected - and bounded in [0, 1]; error metrics come back as NaN.""" + """compute_ceiling returns (results, agg): `results` is the raw per-perturbation + self-split; `agg` is the SB-corrected per-context ceiling (one row) - reliability + metrics bounded in [0, 1], error / excluded metrics NaN.""" adata_real = build_random_anndata() evaluator = MetricsEvaluator( adata_pred=adata_real.copy(), @@ -95,19 +96,18 @@ def test_compute_ceiling_end_to_end(): results, agg = evaluator.compute_ceiling( profile="anndata", write_csv=False, break_on_error=True ) + # per-perturbation self-split measurements assert results.height > 0 assert "perturbation" in results.columns - assert "pearson_delta" in results.columns - pv = results["pearson_delta"].drop_nulls().to_numpy() - assert np.all(pv <= 1.0 + 1e-9) # SB doubling can never exceed 1 + # the ceiling is the SB-corrected aggregate (mean over perturbations) - one row + assert agg.height == 1 + assert "pearson_delta" in agg.columns + cv = agg["pearson_delta"].to_numpy() + assert np.all(cv <= 1.0 + 1e-9) # SB doubling can never exceed 1 - # error metrics are emitted as NaN (no defensible SB ceiling) - for col in ("mse", "mae"): - if col in results.columns: - assert np.all(np.isnan(results[col].to_numpy())) - # excluded reliability metrics are NaN too - for col in ("clustering_agreement", "pearson_edistance"): - if col in results.columns: - assert np.all(np.isnan(results[col].to_numpy())) + # error metrics and excluded reliabilities have no ceiling -> NaN in the aggregate + for col in ("mse", "mae", "clustering_agreement", "pearson_edistance"): + if col in agg.columns: + assert np.all(np.isnan(agg[col].to_numpy())) shutil.rmtree(OUTDIR) diff --git a/tests/test_eval.py b/tests/test_eval.py index 2684031..47c88a1 100644 --- a/tests/test_eval.py +++ b/tests/test_eval.py @@ -418,10 +418,11 @@ def test_eval_ceiling(): results, agg = evaluator.compute_ceiling(break_on_error=True) assert results.height > 0 assert "perturbation" in results.columns - # a reliability metric is SB-corrected; the doubling can never exceed 1 - assert "pearson_delta" in results.columns - pv = results["pearson_delta"].drop_nulls().to_numpy() - assert np.all(pv <= 1.0 + 1e-9) + # the ceiling (SB of the per-context mean) lives in the aggregate and, being a + # doubling of a reliability, can never exceed 1 + assert "pearson_delta" in agg.columns + cv = agg["pearson_delta"].drop_nulls().to_numpy() + assert np.all(cv <= 1.0 + 1e-9) assert os.path.exists(f"{OUTDIR}/ceiling_results.csv") assert os.path.exists(f"{OUTDIR}/agg_ceiling_results.csv") shutil.rmtree(OUTDIR)