feat(ceiling): disjoint-split depth extrapolation (replaces bootstrap)

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.
 

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

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",

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@<frac>``)
+        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@<frac>``), 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,