diff --git a/src/cell_eval/_cli/_run.py b/src/cell_eval/_cli/_run.py index 8cb4fdc..d91ee56 100644 --- a/src/cell_eval/_cli/_run.py +++ b/src/cell_eval/_cli/_run.py @@ -2,6 +2,7 @@ import importlib.metadata import logging import os +import sys from typing import Any from .. import KNOWN_PROFILES @@ -18,8 +19,11 @@ def parse_args_run(parser: ap.ArgumentParser): "-ap", "--adata-pred", type=str, - help="Path to the predicted adata object to evaluate", - required=True, + help="Path to the predicted adata object to evaluate. Optional in " + "ceiling-only mode (omit when passing --ceiling to compute just the " + "real-data ceiling without a prediction). In that mode only the ceiling " + "outputs are written - no results.csv and no DE tables.", + required=False, ) parser.add_argument( "-ar", @@ -167,6 +171,19 @@ def run_evaluation(args: ap.Namespace): skip_metrics = args.skip_metrics.split(",") if args.skip_metrics else None + # Ceiling-only mode: no prediction supplied, so only the real-data ceiling + # can be computed. Requires --ceiling to make the intent explicit. + ceiling_only = args.adata_pred is None + if ceiling_only and not args.ceiling: + # Usage error, so exit like argparse would (message on stderr, status 2) + # rather than surfacing a traceback for a wrong invocation. + print( + "cell-eval run: error: --adata-pred is required unless --ceiling is " + "passed (ceiling-only mode)", + file=sys.stderr, + ) + sys.exit(2) + # DE knobs forwarded to pdex (cpm_filter off by default; epsilon default 0.0). pdex_kwargs: dict[str, Any] = { "cpm_filter": args.cpm_filter, @@ -175,14 +192,16 @@ def run_evaluation(args: ap.Namespace): if args.celltype_col is not None: real = ad.read_h5ad(args.adata_real) - pred = ad.read_h5ad(args.adata_pred) - real_split = split_anndata_on_celltype(real, args.celltype_col) - pred_split = split_anndata_on_celltype(pred, args.celltype_col) - assert len(real_split) == len(pred_split), ( - f"Number of celltypes in real and pred anndata must match: {len(real_split)} != {len(pred_split)}" - ) + if ceiling_only: + pred_split = {ct: None for ct in real_split} + else: + pred = ad.read_h5ad(args.adata_pred) + pred_split = split_anndata_on_celltype(pred, args.celltype_col) + assert len(real_split) == len(pred_split), ( + f"Number of celltypes in real and pred anndata must match: {len(real_split)} != {len(pred_split)}" + ) for ct in real_split.keys(): real_ct = real_split[ct] @@ -202,12 +221,13 @@ def run_evaluation(args: ap.Namespace): skip_de=args.profile == "pds", pdex_kwargs=pdex_kwargs, ) - evaluator.compute( - profile=args.profile, - metric_configs=metric_kwargs, - skip_metrics=skip_metrics, - basename="results.csv", - ) + if not ceiling_only: + evaluator.compute( + profile=args.profile, + metric_configs=metric_kwargs, + skip_metrics=skip_metrics, + basename="results.csv", + ) if args.ceiling: evaluator.compute_ceiling( profile=args.profile, @@ -231,12 +251,13 @@ def run_evaluation(args: ap.Namespace): skip_de=args.profile == "pds", pdex_kwargs=pdex_kwargs, ) - evaluator.compute( - profile=args.profile, - metric_configs=metric_kwargs, - skip_metrics=skip_metrics, - basename="results.csv", - ) + if not ceiling_only: + evaluator.compute( + profile=args.profile, + metric_configs=metric_kwargs, + skip_metrics=skip_metrics, + basename="results.csv", + ) if args.ceiling: evaluator.compute_ceiling( profile=args.profile, diff --git a/src/cell_eval/_evaluator.py b/src/cell_eval/_evaluator.py index 5881684..17cc1b1 100644 --- a/src/cell_eval/_evaluator.py +++ b/src/cell_eval/_evaluator.py @@ -71,8 +71,10 @@ class MetricsEvaluator: Arguments ========= - adata_pred: ad.AnnData | str - Predicted anndata object or path to anndata object. + adata_pred: ad.AnnData | str | None + Predicted anndata object or path to anndata object. May be ``None`` to run + in ceiling-only mode (the data ceiling is estimated from the real data + alone); in that case only :meth:`compute_ceiling` is available. adata_real: ad.AnnData | str Real anndata object or path to anndata object. de_pred: pl.DataFrame | str | None = None @@ -100,7 +102,7 @@ class MetricsEvaluator: def __init__( self, - adata_pred: ad.AnnData | str, + adata_pred: ad.AnnData | str | None, adata_real: ad.AnnData | str, de_pred: pl.DataFrame | str | None = None, de_real: pl.DataFrame | str | None = None, @@ -132,6 +134,25 @@ def __init__( self._skip_de = skip_de self._pdex_kwargs = pdex_kwargs or {} + # Ceiling-only mode: the data ceiling is estimated from the real data + # alone, so no prediction is required. When adata_pred is None we skip + # building the main de_comparison (compute() is unavailable) - only + # compute_ceiling() may be called. + self.ceiling_only = adata_pred is None + + # Precomputed DE cannot be reused in ceiling-only mode: the main comparison + # is skipped, and the ceiling computes DE on its own disjoint halves. Warn + # (rather than fail) so a stray argument does not break an otherwise valid + # run, but the user is not left believing their table was used. + if self.ceiling_only: + for name, value in (("de_pred", de_pred), ("de_real", de_real)): + if value is not None: + logger.warning( + f"{name} is ignored in ceiling-only mode (adata_pred=None): " + f"the ceiling computes differential expression on its own " + f"disjoint halves of the real data." + ) + self.anndata_pair = _build_anndata_pair( real=adata_real, pred=adata_pred, @@ -140,7 +161,7 @@ def __init__( allow_discrete=allow_discrete, ) - if skip_de: + if skip_de or self.ceiling_only: self.de_comparison = None else: self.de_comparison = _build_de_comparison( @@ -166,6 +187,11 @@ def compute( write_csv: bool = True, break_on_error: bool = False, ) -> tuple[pl.DataFrame, pl.DataFrame]: + if self.ceiling_only: + raise ValueError( + "compute() requires a prediction (adata_pred). This evaluator was " + "created without one (ceiling-only mode); call compute_ceiling() instead." + ) pipeline = MetricPipeline( profile=profile, metric_configs=metric_configs, @@ -399,7 +425,7 @@ def _spearman_brown_correct(results: pl.DataFrame) -> pl.DataFrame: def _build_anndata_pair( real: ad.AnnData | str, - pred: ad.AnnData | str, + pred: ad.AnnData | str | None, control_pert: str, pert_col: str, allow_discrete: bool = False, @@ -413,11 +439,24 @@ def _build_anndata_pair( # Cast float16 to float32 since NUMBA (used by pdex) does not support float16 _cast_float16_to_float32(real, which="real") - _cast_float16_to_float32(pred, which="pred") # Validate that the input is normalized and log-transformed _convert_to_normlog(real, which="real", allow_discrete=allow_discrete) - _convert_to_normlog(pred, which="pred", allow_discrete=allow_discrete) + + # Ceiling-only mode: no prediction supplied. The data ceiling reads only + # `.real`, so mirror real into pred to satisfy the pair (it is never scored). + # + # INVARIANT: this aliases the SAME object - `pair.real is pair.pred`. It is not + # a copy, because copying a matrix that is never scored would double peak memory + # for nothing. Safe only because ceiling-only mode blocks `compute()` and + # `compute_ceiling()` derives fresh copies of both halves from `.real`. Anything + # that mutates `.pred` in place would therefore corrupt `.real`: take a copy + # first, or gate the write on `adata_pred is not None`. + if pred is None: + pred = real + else: + _cast_float16_to_float32(pred, which="pred") + _convert_to_normlog(pred, which="pred", allow_discrete=allow_discrete) # Build the anndata pair return PerturbationAnndataPair( diff --git a/tests/test_ceiling.py b/tests/test_ceiling.py index 419b322..cd4b7dc 100644 --- a/tests/test_ceiling.py +++ b/tests/test_ceiling.py @@ -175,3 +175,99 @@ def test_compute_ceiling_end_to_end(): if col in agg.columns: assert np.all(np.isnan(agg[col].to_numpy())) shutil.rmtree(OUTDIR) + + +def _ceiling_only_evaluator(adata_real, **kwargs): + return MetricsEvaluator( + adata_pred=None, + adata_real=adata_real, + control_pert=CONTROL_VAR, + pert_col=PERT_COL, + outdir=OUTDIR, + **kwargs, + ) + + +def test_ceiling_only_mode_shape(): + """adata_pred=None is ceiling-only: the main DE comparison is skipped, the pair's + pred side is the real object itself (a placeholder that is never scored), and + compute() refuses rather than silently scoring real against itself.""" + adata_real = build_random_anndata() + evaluator = _ceiling_only_evaluator(adata_real) + + assert evaluator.ceiling_only + assert evaluator.de_comparison is None # main comparison skipped + # documented invariant: the placeholder aliases real, it is not a copy + assert evaluator.anndata_pair.real is evaluator.anndata_pair.pred + + with pytest.raises(ValueError, match="ceiling-only mode"): + evaluator.compute(profile="anndata", write_csv=False) + + shutil.rmtree(OUTDIR, ignore_errors=True) + + +def test_ceiling_only_matches_ceiling_with_a_prediction(): + """The ceiling is a property of the real data alone, so supplying a prediction + must not change it: the same seed yields the same ceiling either way.""" + adata_real = build_random_anndata() + + with_pred = MetricsEvaluator( + adata_pred=adata_real.copy(), + adata_real=adata_real, + control_pert=CONTROL_VAR, + pert_col=PERT_COL, + outdir=OUTDIR, + skip_de=True, + ) + _, agg_pred = with_pred.compute_ceiling( + profile="anndata", write_csv=False, break_on_error=True, seed=0 + ) + + ceiling_only = _ceiling_only_evaluator(adata_real, skip_de=True) + _, agg_only = ceiling_only.compute_ceiling( + profile="anndata", write_csv=False, break_on_error=True, seed=0 + ) + + assert agg_only.columns == agg_pred.columns + for col in agg_pred.columns: + a, b = agg_pred[col].to_numpy(), agg_only[col].to_numpy() + assert np.all(np.isnan(a) == np.isnan(b)), col + mask = ~np.isnan(a) + assert np.allclose(a[mask], b[mask], rtol=1e-12, atol=0.0), col + + shutil.rmtree(OUTDIR, ignore_errors=True) + + +def test_cli_rejects_missing_prediction_without_ceiling(): + """Omitting --adata-pred without --ceiling leaves nothing to compute. That is a + usage error, so it exits 2 with a message rather than raising a traceback.""" + import argparse + + from cell_eval._cli._run import run_evaluation + + args = argparse.Namespace( + adata_pred=None, + ceiling=False, + embed_key=None, + skip_metrics=None, + num_threads=1, + ) + with pytest.raises(SystemExit) as exc: + run_evaluation(args) + assert exc.value.code == 2 + + +def test_ceiling_only_warns_that_precomputed_de_is_unused(caplog): + """Precomputed DE cannot be reused in ceiling-only mode (the ceiling runs DE on + its own halves), so it warns for both sides rather than failing or going quiet.""" + adata_real = build_random_anndata() + de = pl.DataFrame({"target": ["a"], "feature": ["g"], "p_value": [0.5]}) + + with caplog.at_level("WARNING"): + _ceiling_only_evaluator(adata_real, de_pred=de, de_real=de, skip_de=True) + + warnings = [r.message for r in caplog.records if r.levelname == "WARNING"] + assert any("de_pred is ignored in ceiling-only mode" in m for m in warnings) + assert any("de_real is ignored in ceiling-only mode" in m for m in warnings) + + shutil.rmtree(OUTDIR, ignore_errors=True)