diff --git a/docs/index.md b/docs/index.md
index 53ff8f2..170df23 100644
--- a/docs/index.md
+++ b/docs/index.md
@@ -6,6 +6,7 @@
 :hidden: true
 :maxdepth: 1
 
+tutorials.md
 api.md
 changelog.md
 contributing.md
diff --git a/docs/notebooks/lung_example.ipynb b/docs/notebooks/lung_example.ipynb
new file mode 100644
index 0000000..6c8a8f3
--- /dev/null
+++ b/docs/notebooks/lung_example.ipynb
@@ -0,0 +1,159 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "wei6g2uemfc",
+   "source": "# Benchmarking lung integration\n\nHere we walkthrough applying integration benchmarking metrics on the lung atlas example from the [scIB paper](https://www.nature.com/articles/s41592-021-01336-8). This mirrors the [scib-metrics tutorial](https://scib-metrics.readthedocs.io/en/stable/notebooks/lung_example.html), but uses scib-rapids for GPU-accelerated computation.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "k1ufrfd0xe",
+   "source": "import time\n\nimport numpy as np\nimport pandas as pd\nimport scanpy as sc\n\nimport scib_rapids\nfrom scib_rapids.nearest_neighbors import NeighborsResults, pynndescent",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "u7l9ngw927b",
+   "source": "## Load and preprocess data",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "y9bsu2y7k9k",
+   "source": "adata = sc.read(\n    \"data/lung_atlas.h5ad\",\n    backup_url=\"https://figshare.com/ndownloader/files/24539942\",\n)\nadata",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "code",
+   "id": "6wp4m7qcn7t",
+   "source": "sc.pp.highly_variable_genes(adata, n_top_genes=2000, flavor=\"cell_ranger\", batch_key=\"batch\")\nsc.tl.pca(adata, n_comps=30, use_highly_variable=True)",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "q7ppcgda4vk",
+   "source": "We subset to the highly variable genes so that each method has the same input.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "36kjgh5fyb4",
+   "source": "adata = adata[:, adata.var.highly_variable].copy()\nadata.obsm[\"Unintegrated\"] = adata.obsm[\"X_pca\"]",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "eah0tzhtdzg",
+   "source": "## Run integration methods\n\nHere we run a few embedding-based methods. By focusing on embedding-based methods, we can substantially reduce the runtime of the benchmarking metrics.\n\n### Harmony",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "sfvgtfdgrnk",
+   "source": "from harmony import harmonize\n\nadata.obsm[\"Harmony\"] = harmonize(adata.obsm[\"X_pca\"], adata.obs, batch_key=\"batch\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "42bnd5dmqke",
+   "source": "### scVI",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "stz57vrnvmd",
+   "source": "import scvi\n\nscvi.model.SCVI.setup_anndata(adata, layer=\"counts\", batch_key=\"batch\")\nvae = scvi.model.SCVI(adata, gene_likelihood=\"nb\", n_layers=2, n_latent=30)\nvae.train()\nadata.obsm[\"scVI\"] = vae.get_latent_representation()",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7ko3cf0gmx4",
+   "source": "### scANVI",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "te2a3x7vzkr",
+   "source": "lvae = scvi.model.SCANVI.from_scvi_model(\n    vae,\n    adata=adata,\n    labels_key=\"cell_type\",\n    unlabeled_category=\"Unknown\",\n)\nlvae.train(max_epochs=20, n_samples_per_label=100)\nadata.obsm[\"scANVI\"] = lvae.get_latent_representation()",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0dyu1vtgmkns",
+   "source": "## Compute all metrics\n\nWe define a helper function that computes all metrics for a given embedding and returns a results dictionary.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "wlmt2ehwmo",
+   "source": "def compute_all_metrics(\n    adata,\n    embedding_key,\n    batch_key=\"batch\",\n    label_key=\"cell_type\",\n    pre_integrated_key=\"X_pca\",\n    n_neighbors=90,\n):\n    \"\"\"Compute all scib-rapids metrics for a single embedding.\"\"\"\n    X = np.asarray(adata.obsm[embedding_key], dtype=np.float32)\n    labels = np.asarray(pd.Categorical(adata.obs[label_key]).codes)\n    batch = np.asarray(pd.Categorical(adata.obs[batch_key]).codes)\n\n    # Compute nearest neighbors\n    nn = pynndescent(X, n_neighbors=n_neighbors)\n\n    results = {}\n\n    # Bio conservation\n    results[\"Isolated labels\"] = scib_rapids.isolated_labels(X, labels, batch)\n    kmeans = scib_rapids.nmi_ari_cluster_labels_kmeans(X, labels)\n    results[\"KMeans NMI\"] = kmeans[\"nmi\"]\n    results[\"KMeans ARI\"] = kmeans[\"ari\"]\n    results[\"Silhouette label\"] = scib_rapids.silhouette_label(X, labels)\n    results[\"cLISI\"] = scib_rapids.clisi_knn(nn, labels)\n\n    # Batch correction\n    results[\"Silhouette batch\"] = scib_rapids.silhouette_batch(X, labels, batch)\n    results[\"iLISI\"] = scib_rapids.ilisi_knn(nn, batch)\n    results[\"KBET\"] = scib_rapids.kbet_per_label(nn, batch, labels)\n    results[\"Graph connectivity\"] = scib_rapids.graph_connectivity(nn, labels)\n\n    # PCR comparison (requires pre-integrated embedding)\n    if pre_integrated_key is not None:\n        X_pre = np.asarray(adata.obsm[pre_integrated_key], dtype=np.float32)\n        results[\"PCR comparison\"] = scib_rapids.pcr_comparison(X_pre, X, batch, categorical=True)\n    else:\n        results[\"PCR comparison\"] = 0.0\n\n    return results",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "cprkjbq718h",
+   "source": "## Perform the benchmark\n\nRun all metrics on each embedding.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "l56sgynqc7",
+   "source": "embedding_keys = [\"Unintegrated\", \"Harmony\", \"scVI\", \"scANVI\"]\n\nall_results = {}\nfor key in embedding_keys:\n    print(f\"Computing metrics for {key}...\")\n    t0 = time.time()\n    all_results[key] = compute_all_metrics(adata, key)\n    elapsed = time.time() - t0\n    print(f\"  Done in {elapsed:.1f}s\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "4imqdpx10pt",
+   "source": "## Visualize the results",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "nbx89tlqvq9",
+   "source": "df = pd.DataFrame(all_results).T\ndf.index.name = \"Embedding\"\n\n# Add metric type annotations\nbio_metrics = [\"Isolated labels\", \"KMeans NMI\", \"KMeans ARI\", \"Silhouette label\", \"cLISI\"]\nbatch_metrics = [\"Silhouette batch\", \"iLISI\", \"KBET\", \"Graph connectivity\", \"PCR comparison\"]\n\n# Compute aggregate scores\ndf[\"Bio conservation\"] = df[bio_metrics].mean(axis=1)\ndf[\"Batch correction\"] = df[batch_metrics].mean(axis=1)\ndf[\"Total\"] = 0.6 * df[\"Bio conservation\"] + 0.4 * df[\"Batch correction\"]\n\ndf.round(4)",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "code",
+   "id": "adwoc8otaih",
+   "source": "df[[\"Bio conservation\", \"Batch correction\", \"Total\"]].sort_values(\"Total\", ascending=False)",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
\ No newline at end of file
diff --git a/docs/notebooks/quickstart.ipynb b/docs/notebooks/quickstart.ipynb
new file mode 100644
index 0000000..94b2344
--- /dev/null
+++ b/docs/notebooks/quickstart.ipynb
@@ -0,0 +1,201 @@
+{
+ "cells": [
+  {
+   "cell_type": "markdown",
+   "id": "hfcdfmavktq",
+   "source": "# Quickstart: scib-rapids metrics\n\nThis tutorial demonstrates how to use scib-rapids to compute single-cell integration benchmarking metrics. scib-rapids provides the same metrics as [scib-metrics](https://github.com/YosefLab/scib-metrics), but uses CuPy/RAPIDS for GPU acceleration instead of JAX.\n\nThe API is identical — if you're familiar with scib-metrics, you can use scib-rapids as a drop-in replacement.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "b3qo8n0m9ym",
+   "source": "import numpy as np\nimport scib_rapids\nfrom scib_rapids.nearest_neighbors import NeighborsResults, pynndescent",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "0rcwtblzizq",
+   "source": "## Generate example data\n\nWe create a simple toy dataset with 1000 cells, 50 features, 5 cell types, and 3 batches. In practice, `X` would be a PCA embedding from an AnnData object.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "9zqcbsj3ins",
+   "source": "rng = np.random.default_rng(42)\nn_cells, n_features, n_labels, n_batches = 1000, 50, 5, 3\n\nX = rng.normal(size=(n_cells, n_features)).astype(np.float32)\nlabels = rng.integers(0, n_labels, size=n_cells)\nbatch = rng.integers(0, n_batches, size=n_cells)\n\nprint(f\"X: {X.shape}, labels: {np.unique(labels)}, batches: {np.unique(batch)}\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "ysiihlzae5i",
+   "source": "## Compute nearest neighbors\n\nSeveral metrics (LISI, kBET, graph connectivity, Leiden clustering) require a precomputed k-nearest neighbor graph. scib-rapids provides a `pynndescent` wrapper that returns a `NeighborsResults` dataclass.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "rwf2yqookf",
+   "source": "nn = pynndescent(X, n_neighbors=30)\nprint(f\"NeighborsResults: indices {nn.indices.shape}, distances {nn.distances.shape}\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "wxytsr0hga",
+   "source": "## Bio-conservation metrics\n\nThese metrics measure how well biological signal (cell type identity) is preserved after integration.\n\n### Silhouette scores",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "0f3wpmxrfpfe",
+   "source": "# Average Silhouette Width — measures cell type separation\n# Returns a score in [0, 1] (rescaled) where higher = better separation\nasw = scib_rapids.silhouette_label(X, labels)\nprint(f\"Silhouette label (ASW): {asw:.4f}\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "code",
+   "id": "db5lpk6t86l",
+   "source": "# Isolated labels — score for cell types present in few batches\niso = scib_rapids.isolated_labels(X, labels, batch)\nprint(f\"Isolated labels: {iso:.4f}\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "bppaoyqjq2p",
+   "source": "### Clustering-based metrics (NMI / ARI)",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "2daq4xiyto6",
+   "source": "# K-means clustering agreement with true labels\nkmeans_result = scib_rapids.nmi_ari_cluster_labels_kmeans(X, labels)\nprint(f\"KMeans NMI: {kmeans_result['nmi']:.4f}\")\nprint(f\"KMeans ARI: {kmeans_result['ari']:.4f}\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "code",
+   "id": "gk87gnmbmgj",
+   "source": "# Leiden clustering agreement (uses the kNN graph)\nleiden_result = scib_rapids.nmi_ari_cluster_labels_leiden(nn, labels)\nprint(f\"Leiden NMI: {leiden_result['nmi']:.4f}\")\nprint(f\"Leiden ARI: {leiden_result['ari']:.4f}\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "3jzyghql9hp",
+   "source": "### cLISI — cell-type Local Inverse Simpson Index",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "kgmn3rzv6un",
+   "source": "# cLISI — how well cell types are preserved in neighborhoods\n# Higher = better cell type preservation (scaled to [0, 1])\nclisi = scib_rapids.clisi_knn(nn, labels)\nprint(f\"cLISI: {clisi:.4f}\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "7lhgb3o3ylo",
+   "source": "## Batch-correction metrics\n\nThese metrics measure how well batch effects have been removed.\n\n### Silhouette batch & BRAS",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "5c1ahtgwqzv",
+   "source": "# Silhouette batch — measures batch mixing within cell types\nsil_batch = scib_rapids.silhouette_batch(X, labels, batch)\nprint(f\"Silhouette batch: {sil_batch:.4f}\")\n\n# BRAS — Batch Removal Adapted Silhouette (uses cosine distance, mean_other)\nbras_score = scib_rapids.bras(X, labels, batch)\nprint(f\"BRAS: {bras_score:.4f}\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "997uerovcmo",
+   "source": "### iLISI — integration LISI",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "g7mq1ji756g",
+   "source": "# iLISI — how well batches are mixed in neighborhoods\n# Higher = better batch integration (scaled to [0, 1])\nilisi = scib_rapids.ilisi_knn(nn, batch)\nprint(f\"iLISI: {ilisi:.4f}\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "qs7qfkjvhri",
+   "source": "### kBET",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "7ia6pfojmbi",
+   "source": "# kBET — batch effect test via chi-square statistics\n# acceptance_rate close to 1 means batches are well mixed\nacceptance_rate, test_stats, p_values = scib_rapids.kbet(nn, batch)\nprint(f\"kBET acceptance rate: {acceptance_rate:.4f}\")\n\n# Per-label kBET (uses diffusion distances within each cell type)\nkbet_score = scib_rapids.kbet_per_label(nn, batch, labels)\nprint(f\"kBET per label: {kbet_score:.4f}\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "uc1tbd9hl3i",
+   "source": "### Graph connectivity",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "q6fqb045uer",
+   "source": "# Graph connectivity — fraction of cells in the largest connected component\n# per cell type subgraph. Higher = better.\ngc = scib_rapids.graph_connectivity(nn, labels)\nprint(f\"Graph connectivity: {gc:.4f}\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "hmm14u2au94",
+   "source": "### PCR comparison",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "iz8vintxvyn",
+   "source": "# PCR comparison — reduction in batch variance explained by PCA\n# Compares pre- vs post-integration embeddings\nX_post = rng.normal(size=X.shape).astype(np.float32)\npcr = scib_rapids.pcr_comparison(X, X_post, batch, categorical=True)\nprint(f\"PCR comparison: {pcr:.4f}\")",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  },
+  {
+   "cell_type": "markdown",
+   "id": "vro1fk08j3j",
+   "source": "## Summary table\n\nCollect all metrics into a single table.",
+   "metadata": {}
+  },
+  {
+   "cell_type": "code",
+   "id": "a67mgvo76v",
+   "source": "import pandas as pd\n\nresults = {\n    \"Metric\": [\n        \"Silhouette label\", \"Isolated labels\",\n        \"KMeans NMI\", \"KMeans ARI\",\n        \"Leiden NMI\", \"Leiden ARI\",\n        \"cLISI\",\n        \"Silhouette batch\", \"BRAS\",\n        \"iLISI\", \"kBET\", \"kBET per label\",\n        \"Graph connectivity\", \"PCR comparison\",\n    ],\n    \"Type\": [\n        \"Bio conservation\", \"Bio conservation\",\n        \"Bio conservation\", \"Bio conservation\",\n        \"Bio conservation\", \"Bio conservation\",\n        \"Bio conservation\",\n        \"Batch correction\", \"Batch correction\",\n        \"Batch correction\", \"Batch correction\", \"Batch correction\",\n        \"Batch correction\", \"Batch correction\",\n    ],\n    \"Score\": [\n        asw, iso,\n        kmeans_result[\"nmi\"], kmeans_result[\"ari\"],\n        leiden_result[\"nmi\"], leiden_result[\"ari\"],\n        clisi,\n        sil_batch, bras_score,\n        ilisi, acceptance_rate, kbet_score,\n        gc, pcr,\n    ],\n}\n\ndf = pd.DataFrame(results).set_index(\"Metric\")\ndf[\"Score\"] = df[\"Score\"].map(\"{:.4f}\".format)\ndf",
+   "metadata": {},
+   "execution_count": null,
+   "outputs": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 3",
+   "language": "python",
+   "name": "python3"
+  },
+  "language_info": {
+   "name": "python",
+   "version": "3.11.0"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 5
+}
\ No newline at end of file
diff --git a/docs/tutorials.md b/docs/tutorials.md
new file mode 100644
index 0000000..e05eeb2
--- /dev/null
+++ b/docs/tutorials.md
@@ -0,0 +1,9 @@
+# Tutorials
+
+```{toctree}
+:hidden: true
+:maxdepth: 1
+
+notebooks/quickstart
+notebooks/lung_example
+```
diff --git a/scripts/benchmark_report.py b/scripts/benchmark_report.py
new file mode 100644
index 0000000..0175752
--- /dev/null
+++ b/scripts/benchmark_report.py
@@ -0,0 +1,301 @@
+#!/usr/bin/env python3
+"""Benchmark report: compare scib-rapids vs scib-metrics speed and numerical agreement.
+
+Runs each metric in separate subprocesses (one per venv) so that JAX and CuPy
+don't interfere.  Outputs a formatted table with timings, speedup, and the
+numerical difference between the two backends.
+
+Usage:
+    python scripts/benchmark_report.py                 # default 500 cells
+    python scripts/benchmark_report.py --n-cells 5000  # larger dataset
+"""
+
+from __future__ import annotations
+
+import argparse
+import json
+import subprocess
+import sys
+import tempfile
+import textwrap
+import time
+from dataclasses import dataclass, field
+from pathlib import Path
+
+import numpy as np
+from sklearn.neighbors import NearestNeighbors
+
+# ---------------------------------------------------------------------------
+# Configuration
+# ---------------------------------------------------------------------------
+_ENVS = Path("/home/inference/environments")
+RAPIDS_PYTHON = str(_ENVS / "scib-rapids" / "bin" / "python")
+METRICS_PYTHON = str(_ENVS / "scib-metrics" / "bin" / "python")
+
+WARMUP = 1
+REPEATS = 3
+
+
+# ---------------------------------------------------------------------------
+# Data classes
+# ---------------------------------------------------------------------------
+@dataclass
+class MetricResult:
+    name: str
+    rapids_time: float
+    metrics_time: float
+    rapids_value: float | np.ndarray
+    metrics_value: float | np.ndarray
+    max_abs_diff: float = 0.0
+    max_rel_diff: float = 0.0
+
+    @property
+    def speedup(self) -> float:
+        if self.rapids_time == 0:
+            return float("inf")
+        return self.metrics_time / self.rapids_time
+
+
+# ---------------------------------------------------------------------------
+# Subprocess helpers
+# ---------------------------------------------------------------------------
+def _run_in_venv(python: str, script: str, timeout: int = 300) -> str:
+    result = subprocess.run(
+        [python, "-c", script],
+        capture_output=True, text=True, timeout=timeout,
+    )
+    if result.returncode != 0:
+        raise RuntimeError(
+            f"Script failed (rc={result.returncode}).\n"
+            f"--- stderr ---\n{result.stderr[:2000]}"
+        )
+    return result.stdout.strip()
+
+
+def _run_metric_timed(python: str, pkg: str, data_dir: str, code: str) -> tuple[float, dict]:
+    """Run metric code in a subprocess and return (min_time, results_dict)."""
+    wrapper = textwrap.dedent(f"""\
+        import json, time, warnings
+        import numpy as np
+        warnings.filterwarnings("ignore")
+
+        data_dir = "{data_dir}"
+        pkg = "{pkg}"
+
+        X = np.load(data_dir + "/X.npy")
+        labels = np.load(data_dir + "/labels.npy")
+        batch = np.load(data_dir + "/batch.npy")
+        dists = np.load(data_dir + "/dists.npy")
+        inds = np.load(data_dir + "/inds.npy")
+        X_post = np.load(data_dir + "/X_post.npy")
+
+        if pkg == "scib_metrics":
+            import scib_metrics as lib
+            from scib_metrics.nearest_neighbors import NeighborsResults
+        else:
+            import scib_rapids as lib
+            from scib_rapids.nearest_neighbors import NeighborsResults
+
+        nn = NeighborsResults(indices=inds, distances=dists)
+
+        def _run():
+            {textwrap.indent(code, "            ").strip()}
+            return results
+
+        # warmup
+        for _ in range({WARMUP}):
+            _run()
+
+        # timed repeats
+        times = []
+        for _ in range({REPEATS}):
+            t0 = time.perf_counter()
+            results = _run()
+            times.append(time.perf_counter() - t0)
+
+        best = min(times)
+
+        out = {{"_time": best}}
+        for k, v in results.items():
+            if isinstance(v, np.ndarray):
+                path = data_dir + f"/{{pkg}}_{{k}}.npy"
+                np.save(path, v)
+                out[k] = {{"type": "array", "path": path}}
+            else:
+                out[k] = {{"type": "scalar", "value": float(v)}}
+        print(json.dumps(out))
+    """)
+    raw = _run_in_venv(python, wrapper)
+    meta = json.loads(raw)
+    t = meta.pop("_time")
+    loaded = {}
+    for k, info in meta.items():
+        if info["type"] == "array":
+            loaded[k] = np.load(info["path"])
+        else:
+            loaded[k] = info["value"]
+    return t, loaded
+
+
+# ---------------------------------------------------------------------------
+# Metric definitions
+# ---------------------------------------------------------------------------
+METRICS: list[tuple[str, str]] = [
+    ("silhouette_label", 'results = {"score": lib.silhouette_label(X, labels)}'),
+    ("silhouette_batch", 'results = {"score": lib.silhouette_batch(X, labels, batch)}'),
+    ("bras", 'results = {"score": lib.bras(X, labels, batch)}'),
+    ("isolated_labels", 'results = {"score": lib.isolated_labels(X, labels, batch)}'),
+    ("lisi_knn", 'results = {"lisi": np.asarray(lib.lisi_knn(nn, labels, perplexity=10))}'),
+    ("ilisi_knn", 'results = {"score": lib.ilisi_knn(nn, batch, perplexity=10)}'),
+    ("clisi_knn", 'results = {"score": lib.clisi_knn(nn, labels, perplexity=10)}'),
+    ("kbet", textwrap.dedent("""\
+        acc, stats, pvals = lib.kbet(nn, batch)
+        results = {"acc": acc, "stats": np.asarray(stats), "pvals": np.asarray(pvals)}""")),
+    ("kbet_per_label", 'results = {"score": lib.kbet_per_label(nn, batch, labels)}'),
+    ("graph_connectivity", 'results = {"score": lib.graph_connectivity(nn, labels)}'),
+    ("pcr_comparison", 'results = {"score": lib.pcr_comparison(X, X_post, batch, categorical=True)}'),
+    ("nmi_ari_leiden", textwrap.dedent("""\
+        r = lib.nmi_ari_cluster_labels_leiden(nn, labels, optimize_resolution=False, resolution=1.0)
+        results = {"nmi": r["nmi"], "ari": r["ari"]}""")),
+]
+
+
+# ---------------------------------------------------------------------------
+# Main
+# ---------------------------------------------------------------------------
+def _compute_diffs(r: dict, m: dict) -> tuple[float, float]:
+    """Return (max_abs_diff, max_rel_diff) across all keys."""
+    max_abs = 0.0
+    max_rel = 0.0
+    for k in r:
+        rv, mv = np.asarray(r[k], dtype=np.float64), np.asarray(m[k], dtype=np.float64)
+        abs_diff = float(np.max(np.abs(rv - mv)))
+        with np.errstate(divide="ignore", invalid="ignore"):
+            rel_diff = np.where(mv != 0, np.abs((rv - mv) / mv), 0.0)
+        rel_max = float(np.max(rel_diff))
+        max_abs = max(max_abs, abs_diff)
+        max_rel = max(max_rel, rel_max)
+    return max_abs, max_rel
+
+
+def _representative_value(d: dict) -> str:
+    """Pick a single representative value string from a results dict."""
+    for k in ("score", "acc", "nmi"):
+        if k in d:
+            return f"{d[k]:.6f}"
+    # first key
+    k0 = next(iter(d))
+    v = d[k0]
+    if isinstance(v, np.ndarray):
+        return f"array(mean={np.mean(v):.6f})"
+    return f"{v:.6f}"
+
+
+def main():
+    parser = argparse.ArgumentParser(description="Benchmark scib-rapids vs scib-metrics")
+    parser.add_argument("--n-cells", type=int, default=500, help="Number of cells")
+    parser.add_argument("--n-features", type=int, default=50, help="Number of features")
+    parser.add_argument("--n-neighbors", type=int, default=30, help="Number of neighbors for kNN")
+    args = parser.parse_args()
+
+    # Check venvs exist
+    for p in (RAPIDS_PYTHON, METRICS_PYTHON):
+        if not Path(p).exists():
+            print(f"ERROR: interpreter not found: {p}", file=sys.stderr)
+            sys.exit(1)
+
+    # Generate data
+    print(f"Generating test data: {args.n_cells} cells, {args.n_features} features, {args.n_neighbors} neighbors ...")
+    rng = np.random.default_rng(42)
+    n_labels, n_batches = 5, 3
+
+    X = rng.normal(size=(args.n_cells, args.n_features)).astype(np.float64)
+    labels = rng.integers(0, n_labels, size=(args.n_cells,))
+    batch = rng.integers(0, n_batches, size=(args.n_cells,))
+
+    nbrs = NearestNeighbors(n_neighbors=args.n_neighbors, algorithm="kd_tree").fit(X)
+    dists, inds = nbrs.kneighbors(X)
+    X_post = rng.normal(size=X.shape).astype(np.float64)
+
+    tmpdir = tempfile.mkdtemp(prefix="scib_bench_")
+    for name, arr in [("X", X), ("labels", labels), ("batch", batch),
+                      ("dists", dists), ("inds", inds), ("X_post", X_post)]:
+        np.save(f"{tmpdir}/{name}.npy", arr)
+
+    print(f"Data saved to {tmpdir}\n")
+
+    # Run benchmarks
+    results: list[MetricResult] = []
+    n = len(METRICS)
+    for i, (name, code) in enumerate(METRICS, 1):
+        print(f"[{i}/{n}] {name} ...", end=" ", flush=True)
+        try:
+            t_rapids, r_rapids = _run_metric_timed(RAPIDS_PYTHON, "scib_rapids", tmpdir, code)
+            t_metrics, r_metrics = _run_metric_timed(METRICS_PYTHON, "scib_metrics", tmpdir, code)
+            max_abs, max_rel = _compute_diffs(r_rapids, r_metrics)
+            mr = MetricResult(
+                name=name,
+                rapids_time=t_rapids,
+                metrics_time=t_metrics,
+                rapids_value=r_rapids,
+                metrics_value=r_metrics,
+                max_abs_diff=max_abs,
+                max_rel_diff=max_rel,
+            )
+            results.append(mr)
+            print(f"done ({mr.speedup:.1f}x)")
+        except Exception as e:
+            print(f"FAILED: {e}")
+
+    # ---------------------------------------------------------------------------
+    # Print report
+    # ---------------------------------------------------------------------------
+    print("\n")
+    print("=" * 100)
+    print(f"  BENCHMARK REPORT: scib-rapids vs scib-metrics")
+    print(f"  Dataset: {args.n_cells} cells, {args.n_features} features, {args.n_neighbors} neighbors")
+    print(f"  Timing: best of {REPEATS} runs after {WARMUP} warmup")
+    print("=" * 100)
+
+    # Table header
+    hdr = f"{'Metric':<22} {'rapids (s)':>10} {'metrics (s)':>11} {'Speedup':>8} {'Max |Δ|':>10} {'Max |Δ|/|x|':>12} {'rapids val':>16} {'metrics val':>16}"
+    print()
+    print(hdr)
+    print("-" * len(hdr))
+
+    for r in results:
+        rv_str = _representative_value(r.rapids_value)
+        mv_str = _representative_value(r.metrics_value)
+        speedup_str = f"{r.speedup:.1f}x"
+        print(
+            f"{r.name:<22} {r.rapids_time:>10.4f} {r.metrics_time:>11.4f} {speedup_str:>8} "
+            f"{r.max_abs_diff:>10.2e} {r.max_rel_diff:>12.2e} {rv_str:>16} {mv_str:>16}"
+        )
+
+    print("-" * len(hdr))
+
+    # Summary
+    speedups = [r.speedup for r in results if np.isfinite(r.speedup)]
+    if speedups:
+        geo_mean = np.exp(np.mean(np.log(speedups)))
+        print(f"\nGeometric mean speedup: {geo_mean:.2f}x")
+
+    max_abs_all = max(r.max_abs_diff for r in results)
+    max_rel_all = max(r.max_rel_diff for r in results)
+    print(f"Worst absolute diff:   {max_abs_all:.2e}")
+    print(f"Worst relative diff:   {max_rel_all:.2e}")
+
+    all_close = all(r.max_abs_diff < 1e-3 for r in results)
+    if all_close:
+        print("\nAll metrics agree within 1e-3 absolute tolerance.")
+    else:
+        print("\nWARNING: Some metrics differ by more than 1e-3!")
+        for r in results:
+            if r.max_abs_diff >= 1e-3:
+                print(f"  - {r.name}: max |Δ| = {r.max_abs_diff:.2e}")
+
+    print()
+
+
+if __name__ == "__main__":
+    main()