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

[beabevi]

Summary

Reworks the --ceiling data ceiling (MetricsEvaluator.compute_ceiling) from the bootstrap self-split to a disjoint-split depth extrapolation. For each perturbation the real data is split without replacement into two halves at several depths (fracs=(1.0, 0.5, 0.25)), the metric suite is run on each self-split, and every metric's value-vs-depth curve is extrapolated to full depth. The result is a per-metric, unbiased estimate of the best score any model could reach given the noise in the real data.

Why

The bootstrap drew both halves from the same cells (2n draws with replacement, then split), so they weren't independent. That biases the ceiling in both directions, so it isn't a reliable upper bound:

• the shared cells make the two halves agree more than two independent samples would, which inflates the ceiling;
• the duplicate cells over-call significance, which can pull the DE recovery metrics the other way.

A disjoint (without-replacement) split gives genuinely independent halves (unbiased), at the cost of depth (each half ≤ n/2). Measuring each metric across several depths and extrapolating to full depth corrects for that depth cost.

Relationship to Spearman-Brown: for reliability-like (correlation) metrics this extrapolation is the Spearman-Brown correction: SB (2r/(1+r)) is its analytical closed form. The extrapolation generalizes the same idea to metrics that have no closed form (the DE recovery metrics), which is why it's used as the single method rather than SB.

What changed

• Remove _bootstrap_halves; add _disjoint_halves (disjoint split at a depth fraction) + a per-metric depth-curve fit (attenuation/reliability form, linear fallback), extrapolated to full depth.
• compute_ceiling returns one per-metric table (measure, with extrap_linear / resid / model diagnostics) instead of the old (results, agg_results) tuple.
• Writes a single ceiling_results.csv (was ceiling_reslts.csv).
• Docs updated (--ceiling CLI help, README, CLAUDE.md). The existing ceiling tests in tests/test_eval.py were updated for the new API, and a new tests/test_ceiling_extrap.py was added (unit tests for the extrapolation math + a disjoint-split check). pytest, ruff, and ty (on the changed files) pass.

[gemini-code-assist]

Code Review

This pull request refactors the data ceiling estimation in cell-eval to use a disjoint-split depth extrapolation method instead of a bootstrap self-split, which avoids bias in the ceiling estimates. The changes update the CLI, documentation, MetricsEvaluator class, and unit tests. The review feedback highlights three opportunities to improve robustness: validating that the fracs parameter is not empty to prevent a crash in min(), guarding against empty index lists in np.concatenate during disjoint splitting, and filtering out floating-point NaN values during metric aggregation to prevent downstream extrapolation failures.

[gemini-code-assist]

If fracs is empty, min(fracs) on line 206 will raise a ValueError: min() arg is an empty sequence. Explicitly validating that fracs is not empty prevents this crash.

Suggested change

	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}")
	fracs = tuple(sorted(set(fracs), reverse=True))
	if not fracs:
	raise ValueError("fracs must not be empty")
	if any(f <= 0.0 or f > 1.0 for f in fracs):
	raise ValueError(f"fracs must be in (0, 1]; got {fracs}")

[gemini-code-assist]

If no perturbations have enough cells to split at the given frac (e.g., if eligible is empty or cell counts are too low), a_idx and b_idx will be empty. Calling np.concatenate on an empty list raises a ValueError. Guarding against this with a descriptive error improves robustness.

        # Disjoint split has no duplicate rows, so obs names stay unique.
        if not a_idx:
            raise ValueError(
                f"No perturbations had enough cells to split at frac={frac}. "
                "Ensure your dataset has sufficient cells per perturbation."
            )
        half_real = real[np.concatenate(a_idx)].copy()
        half_pred = real[np.concatenate(b_idx)].copy()

[gemini-code-assist]

Polars' drop_nulls() drops missing values (null), but does not drop floating-point NaN values. If any metric computation returns NaN for a perturbation, np.median or np.mean will propagate the NaN, causing the subsequent curve fitting/extrapolation to fail. Filtering out non-finite values (including NaN and inf) ensures robust aggregation.

Suggested change

	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))
	arr = results[col].drop_nulls().to_numpy()
	arr = arr[np.isfinite(arr)]
	if arr.size == 0:
	continue
	out[col] = float(np.median(arr) if agg == "median" else np.mean(arr))

[abhinadduri]

can you bump the semver as this changes behavior: version = "0.8.1" -> version = "0.8.2"

[LeonHafner]

I ran it on the in-vivo spleen dataset and compared three ways of doing the depth extrapolation. Plots below: v1 = current PR, v2 = current PR with the intercept pinned to 1, v3 = fit in signal-to-noise space through the origin. Posting a quite extensive explanation here, so people (including me) still now the reasoning in a couple of days...

TL;DR: the disjoint-split + depth-extrapolation idea is sound and a real improvement over the bootstrap ceiling. But the current extrapolation produces mathematically impossible values (correlations > 1) on several metric × cell-type combinations. I'd suggest switching to the v3 formulation, which is bounded by construction, and not computing the ceiling at all for a few metrics (see the last section).

Out-of-range extrapolations across the 6 metrics × 8 cell types shown:

implementation	out of [0, 1]
v1 — current PR	9 / 48
v2 — intercept = 1	2 / 48
v3 — ours (SNR space)	0 / 48

v1 — current PR (free intercept, fit in 1/m space):

v2 — PR with intercept = 1 (fit in 1/m space, intercept pinned):

v3 — ours (intercept 0, fit in SNR space):

The model...why depth extrapolation works

The disjoint split gives us a metric $r$ (a correlation / agreement score) measured between two independent halves at reduced depth $d$ (cells per half). Under the standard measurement-error model, such a score behaves like a reliability:

$$r(d) = \frac{S}{S + N/d}$$

where $S$ is the variance of the true (biological) signal across genes and $N/d$ is the sampling-noise variance of a depth- $d$ estimate. The signal-to-noise ratio is then

$$\mathrm{SNR}(d) ;=; \frac{r}{1-r} ;=; \frac{S}{N},d \qquad\Longrightarrow\qquad \mathrm{SNR}\ \text{is linear in depth } d\ (\propto n),\ \text{through the origin,}$$

and the inverse map back to the metric is

$$r = \frac{\mathrm{SNR}}{1+\mathrm{SNR}}.$$

Because a disjoint split only reaches half depth (each half is $\le n/2$; in the code frac = 1 ↔ $n/2$ per half, frac = 2 ↔ full $n$), we measure $r$ at a few shallow depths, fit the SNR line, and evaluate it at full depth (frac = 2).

Equivalence to Spearman–Brown

Going from half depth to full depth is exactly doubling the sample size, which doubles the SNR. Substituting $\mathrm{SNR}=r/(1-r)$:

$$\mathrm{SNR}' = 2,\mathrm{SNR} \quad\Longrightarrow\quad r' = \frac{2,\mathrm{SNR}}{1+2,\mathrm{SNR}} = \frac{2r}{1+r},$$

which is the Spearman–Brown correction. So the depth extrapolation is Spearman–Brown, generalized to a least-squares fit over several depths (more robust than plugging a single point into the closed form), and it extends the same idea to metrics that have no closed form.

Why the current implementation (v1) overshoots

The PR fits the reliability in reciprocal coordinates,

$$\frac{1}{r} = a + \frac{b}{\text{frac}},$$

with a free intercept $a$. But $r \le 1$ requires $1/r \ge 1$ always, i.e. the intercept must satisfy $a \ge 1$ (equivalently: the SNR line must pass through the origin, since $\tfrac{1}{r} = 1 + \tfrac{1}{\mathrm{SNR}}$). Leaving $a$ free, a noisy / low-signal cell type fits $a < 1$ (e.g. Plasma pearson_delta: $a = -1.44$), the line drops below the $1/r = 1$ floor, and inverting gives $r > 1$. In v1 this happens for pearson_delta (Plasma → 1.50), discrimination_score_cosine (→ 1.4), clustering_agreement (→ 1.85 and −0.97), etc.

Why we favor v3 (SNR space, through the origin)

Two independent problems, and only v3 fixes both:

1. Boundedness. With the SNR fit through the origin, the extrapolated SNR is $\ge 0$, and $r = \mathrm{SNR}/(1+\mathrm{SNR})$ can never reach or exceed 1 — the ceiling of 1 is baked into the map. Pinning the PR intercept to 1 (v2) also enforces this and fixes the upper overshoots.
2. Stability near zero. v2 still fits in $1/r$ space, which explodes when the correlation is near zero ($1/r \to \infty$). This is exactly what happens for de_spearman_lfc_sig (v2 → 23.7 for B cells). v3 fits $\mathrm{SNR} = r/(1-r) \approx r$ near zero, which is well-behaved — so it stays bounded there too.

Net: v3 is the only one of the three with zero out-of-range extrapolations across every metric and cell type here, while remaining exactly Spearman–Brown on the well-behaved metrics (v1/v2/v3 all agree on roc_auc, overlap_at_N, and the larger cell types for pearson_delta). The through-origin constraint is also the physically correct one: zero cells → zero recovered signal, and infinite cells → perfect reliability.

Concretely, per cell type we fit

$$k = \frac{\sum_i \text{frac}_i \cdot \mathrm{SNR}_i}{\sum_i \text{frac}_i^{,2}}, \qquad \text{ceiling} = \frac{2k}{1+2k}.$$

Caution — some metrics should not get a ceiling like this

The whole construction assumes the metric is a reliability-like quantity in [0, 1] that improves monotonically with depth. That holds for the correlation / agreement / recovery metrics (pearson_delta, the spearman correlations, roc_auc, pr_auc, overlap@N, precision@N, de_sig_genes_recall, discrimination_score). It does not hold for:

• Count metrics — e.g. the number-of-significant-genes metric (de_nsig_counts). These are unbounded counts with best_value = NONE; there is no "maximum achievable" to extrapolate to. We should not emit a ceiling for these at all (skip them, or route on the registry's best_value).
• Quantized / non-monotone metrics — e.g. clustering_agreement. Its depth curve isn't a smooth reliability (it jumps between discrete cluster solutions), so even a bounded extrapolation is not trustworthy. Best excluded or flagged.
• Near-zero, noisy metrics — e.g. de_spearman_lfc_sig. v3 keeps it in range, but the estimate is fragile; treat with caution.
• Small cell types (few perturbations) give noisy depth curves and larger extrapolation error regardless of method.

Suggested changes

1. Fit in SNR space through the origin (v3) — or equivalently constrain the reciprocal-space intercept to $a \ge 1$.
2. Decide which metrics get a ceiling from the registry's best_value (compute for ONE; skip NONE / count metrics like de_nsig_counts).
3. Add a monotonicity / stability guard and fall back to the raw half-depth value (frac = 1) when the depth curve isn't monotone — this covers clustering_agreement-type metrics.

[LeonHafner]

Incorporating the previously discussed plots for training on 0.25, 0.5, and 0.75 while predicting at 1.0 (reaching a 1.33x extrapolation), as well as training on 0.25 and 0.5 while predicting at 1.0 (reaching a 2.0x extrapolation). Additionally, including plots that use the mean instead of the median and illustrate the behavior across datasets and methods.

---

---

---

[LeonHafner]

Adding another dataset

[LeonHafner]

Switched the ceiling estimator to the plain Spearman-Brown correction

I've pushed a commit that replaces the multi-depth reciprocal-fit extrapolation with the analytical Spearman-Brown doubling r' = 2r/(1+r).

Why. The holdout backtest — predict the held-out frac=1 from shallower depths, across 4 datasets (484 context × metric cases) — shows the plain SB correction is:

• the most accurate estimator: lowest median |error| at both reaches (SB 0.0139 @1.33× / 0.0253 @2×), best in every block.
• bounded by construction — 2r/(1+r) ∈ [0,1) for r ∈ [0,1), so 0 out-of-[0,1] predictions vs 21–31 for the free-intercept fit (no inv>0 guard, no linear fallback, no clipping needed);
• the cheapest — SB only needs the deepest split, so the depth sweep is unnecessary compute.

What changed. compute_ceiling now takes a single disjoint n/2 self-split, runs the metric pipeline once, and applies 2r/(1+r) per perturbation. The depth sweep, the 1/m = a + b/frac fit, and the linear fallback are removed.

• The correction is applied to an explicit SB_METRICS list at the top of _evaluator.py. Every other metric — error metrics, unbounded counts, and reliability metrics we don't trust the doubling for (clustering_agreement, pearson_edistance) — is emitted as NaN.
• Returns (results, agg_results) and writes ceiling_results.csv / agg_ceiling_results.csv, matching compute().
• Tests, README, and CLAUDE.md updated.

This builds directly on the disjoint-split groundwork in this PR — it just swaps the estimator that sits on top of it. @beabevi flagging since it's a change to your branch.