diff --git a/.gitattributes b/.gitattributes
index a6344aac8c09253b3b630fb776ae94478aa0275b..bf39d270e09af3183514bed9d9ef53fe893250d0 100644
--- a/.gitattributes
+++ b/.gitattributes
@@ -33,3 +33,10 @@ saved_model/**/* filter=lfs diff=lfs merge=lfs -text
 *.zip filter=lfs diff=lfs merge=lfs -text
 *.zst filter=lfs diff=lfs merge=lfs -text
 *tfevents* filter=lfs diff=lfs merge=lfs -text
+assets/fonts/Ubuntu-Bold.ttf filter=lfs diff=lfs merge=lfs -text
+assets/fonts/Ubuntu-Regular.ttf filter=lfs diff=lfs merge=lfs -text
+assets/remedi_overview.png filter=lfs diff=lfs merge=lfs -text
+examples/sciplex/sciplex_prototype.h5ad filter=lfs diff=lfs merge=lfs -text
+examples/tahoe/tahoe.h5ad filter=lfs diff=lfs merge=lfs -text
+src/remedi/fonts/Ubuntu-Bold.ttf filter=lfs diff=lfs merge=lfs -text
+src/remedi/fonts/Ubuntu-Regular.ttf filter=lfs diff=lfs merge=lfs -text
diff --git a/LICENSE b/LICENSE
new file mode 100644
index 0000000000000000000000000000000000000000..29fdd8631ecc4e3b69ddfb5dafdf7b3bca967eb7
--- /dev/null
+++ b/LICENSE
@@ -0,0 +1,24 @@
+MIT License
+
+Copyright (c) 2026 ChatterjeeLab
+
+Permission is hereby granted, free of charge, to any person obtaining a copy
+of this software and associated documentation files (the "Software"), to deal
+in the Software without restriction, including without limitation the rights
+to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
+copies of the Software, and to permit persons to whom the Software is
+furnished to do so, subject to the following conditions:
+
+The above copyright notice and this permission notice shall be included in all
+copies or substantial portions of the Software.
+
+THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
+IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
+AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
+OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+SOFTWARE.
+
+The included dataset subsets and Ubuntu fonts retain their respective
+licenses, as described in examples/ATTRIBUTION.md and assets/fonts/UFL.txt.
diff --git a/Makefile b/Makefile
new file mode 100644
index 0000000000000000000000000000000000000000..a8557cc62b42b7e1fcb58234700a70642216a197
--- /dev/null
+++ b/Makefile
@@ -0,0 +1,9 @@
+.PHONY: install test smoke validate
+install:
+	python -m pip install -e '.[data,dev]'
+test:
+	python -m pytest
+smoke:
+	python -m remedi.cli demo --output runs/demo --max-queries 4
+validate:
+	python scripts/validate_release.py --output runs/validation
diff --git a/README.md b/README.md
index 154df8298fab5ecf322016157858e08cd1bccbe1..5d04d68ce54bd40534d6c8ab8a5fe80680b9a41c 100644
--- a/README.md
+++ b/README.md
@@ -1,3 +1,179 @@
 ---
-license: apache-2.0
+license: mit
+tags:
+- biology
+- single-cell
+- molecular-design
+- perturbation
+- dose-response
+language:
+- en
 ---
+
+# ReMEDi
+
+![ReMEDi](assets/remedi_overview.png)
+
+**Replicate-Enabled Minimax Endpoint-Directed Inference**
+
+Generate catalog molecules and concentrations for a desired cellular endpoint. Fit a response ensemble to treated-minus-control features, estimate joint candidate-loss uncertainty from shared draws and replicates, and sample molecule–dose pairs by entropy-regularized minimax regret.
+
+## Installation
+
+Python 3.11 or 3.12 on Linux is recommended. The release was tested with Python 3.12. From the repository root:
+
+```bash
+python -m venv .venv
+source .venv/bin/activate
+python -m pip install --upgrade pip
+pip install -e '.[data,dev]'
+pytest
+remedi demo --output runs/demo --max-queries 4
+```
+
+To reproduce the bundled real-count integration checks for both studies and frozen external transfer:
+
+```bash
+python scripts/validate_release.py --output runs/validation
+```
+
+The measured fixtures are included in `examples/`. This command works without downloading the full studies.
+
+The demo uses synthetic responses and checks the complete computation. CPU execution is sufficient for the included ridge head, calibration, and convex solver. Frozen transformer embedding extraction requires the optional `encoders` dependencies.
+
+## Data
+
+Preparation uses RNA counts, molecular structures, concentrations, cell-line identifiers, and matched control pools from each study.
+
+### Tahoe-100M
+
+Download selected official expression shards and metadata:
+
+```bash
+python scripts/download_data.py tahoe --output data/tahoe --indices 600 601 2600
+python scripts/prepare_tahoe.py \
+  --shards data/tahoe/train-*.parquet --metadata data/tahoe \
+  --output data/tahoe/cohort --cap 128
+```
+
+The selected shards cover plates 6 and 14. They support a small replicate-based check at 5 µM. Download additional plates for a broader dose and compound study. `prepare_tahoe.py` accepts any collection of official local expression shards. It joins concentration by sample, removes the CLS pair, converts gene tokens to Ensembl identifiers, and retains a reproducible reservoir per compound, dose, cell line, and plate.
+
+### sci-Plex3
+
+The corrected scPerturb v1.4 release is public at [Zenodo](https://doi.org/10.5281/zenodo.13350497). The complete file is approximately 2.53 GB.
+
+```bash
+python scripts/download_data.py sciplex --output data/sciplex
+python scripts/fetch_sciplex_subset.py \
+  --source data/sciplex/sciplex3.h5ad --metadata-cache data/sciplex \
+  --output data/sciplex/cohort --cells-per-condition 32 --max-molecules 40
+```
+
+This bounded cohort uses A549 at 24 hours, 0.1 and 1 µM, with both replicates and plate-matched controls. Compound structures are joined by exact case-folded names to Tahoe metadata. The nine exact name matches define a small integration test. A larger experiment requires an independently checked `drug,smiles` mapping for the remaining compounds. The `prepare` command accepts the complete corrected file and a corresponding mapping.
+
+The subset script also accepts the source HTTPS URL and caches byte ranges. The full-file route is preferable when the server limits range requests.
+
+## Train and evaluate
+
+For an end-to-end run on either prepared cohort:
+
+```bash
+python scripts/run_pipeline.py \
+  --h5ad data/sciplex/cohort/sciplex_prototype.h5ad \
+  --mapping configs/sciplex.yaml \
+  --structures data/sciplex/cohort/structures.csv \
+  --output runs/sciplex --min-cells 8 --max-cells 32 \
+  --tolerance 1 --max-queries 6
+```
+
+For Tahoe, replace the three input paths with `data/tahoe/cohort/tahoe.h5ad`, `configs/tahoe.yaml`, and its `structures.csv`. The tolerance is a squared endpoint-distance threshold. Select its value on validation data for a scientific comparison. The command above uses 1 for an integration check.
+
+The output includes molecule assignments, condition features, response checkpoints, calibration, all eight generation rules, catalog-pair samples, per-query measurements, molecule-bootstrap summaries, and Ubuntu PDF/PNG plots. The example generation catalog contains measured cohort compounds. Supply your own REAL catalog for REAL generation.
+
+For larger studies, use the individual commands to retain fixed splits and hyperparameters:
+
+```bash
+remedi split --molecules structures.csv --output splits.csv --seed 0
+remedi prepare --h5ad cohort.h5ad --mapping configs/tahoe.yaml \
+  --structures structures.csv --splits splits.csv --output runs/tahoe/data
+remedi train --data runs/tahoe/data --output runs/tahoe/model --ensemble 5
+remedi calibrate --data runs/tahoe/data --model runs/tahoe/model \
+  --output runs/tahoe/calibration
+remedi evaluate --data runs/tahoe/data --model runs/tahoe/model \
+  --calibration runs/tahoe/calibration --output runs/tahoe/evaluation \
+  --tolerance 1 --max-queries 100
+remedi plot --summary runs/tahoe/evaluation/summary.csv --output runs/tahoe/plots
+```
+
+Splits are assigned by canonical molecular structure, with every dose and replicate of that structure in one split. `--scaffold` creates grouped Murcko-scaffold splits. Calibration needs at least three held-out compounds with repeated measurements. Retrospective queries use separate replicate/control pools for target construction and scoring.
+
+## Representations
+
+The executable CPU configuration uses Morgan fingerprints and training-only log-count PCA. Both are explicit baselines for the frozen-encoder configuration in the manuscript.
+
+For frozen UCE features, store per-cell embeddings in `AnnData.obsm['X_uce']` and record the exact checkpoint identifier as `cell_checkpoint` in the mapping YAML. Use `--representation obsm:X_uce` during preparation. UCE extraction is performed with the [released UCE implementation](https://github.com/snap-stanford/UCE). Per-cell embeddings are aggregated after extraction.
+
+Cache molecular embeddings with a pinned MolFormer checkpoint:
+
+```bash
+pip install -e '.[encoders]'
+remedi encode-molecules --molecules structures.csv --output molecular_features.npz \
+  --revision CHECKPOINT_COMMIT --device cpu
+remedi train --data runs/tahoe/data --output runs/tahoe/molformer \
+  --encoder cached --embeddings molecular_features.npz
+```
+
+The cache must cover training, calibration, test, and generation structures. Pass `--embeddings` after moving a cache to another machine. A small MLP response head is available through `train --head mlp`. Encoder extraction and MLP training have separate computational requirements from the CPU reference configuration.
+
+## Generate molecule–dose pairs
+
+Provide a CSV with `catalog_id,smiles`, using actual catalog records. ReMEDi retains those IDs in every output.
+
+```bash
+remedi generate --data runs/tahoe/data --model runs/tahoe/model \
+  --calibration runs/tahoe/calibration --catalog real_subset.csv \
+  --target-rows 10 11 --doses 0.05 0.5 5 --shortlist 512 \
+  --draws 20 --output runs/design
+```
+
+Target row numbers refer to `conditions.csv` in the prepared data directory. All selected rows must describe the same compound, dose, and cellular context. Their observed endpoint defines the desired state. `distribution.csv` contains probabilities for the retained pairs. `generated_pairs.csv` contains samples with their catalog IDs and concentrations in µM.
+
+Catalog membership is determined from the supplied records. Current stock, synthesis availability, and biological activity require separate confirmation. Catalog preparation is bounded by available memory. Large REAL collections should be partitioned into explicit candidate pools before optimization.
+
+## Formulation
+
+For action $a=(m,d)$ and context $c$, a response head predicts $f_\theta(m,\log d,c)$. Shared head and target/control draws give a loss mean $\bar L$ and covariance $C=BB^\top+\operatorname{diag}(v)$. Generation solves
+
+$$
+\min_{p\in\Delta_N}\max_b\left[(p-e_b)^\top\bar L+
+\beta\sqrt{\|B^\top(p-e_b)\|_2^2+\sum_a v_a(p_a-\mathbf 1[a=b])^2}\right]
++\tau\,\mathrm{KL}(p\|\pi).
+$$
+
+Optimization uses CLARABEL with active comparator constraints. Output records include solver status and constraint residuals. The radius is estimated from held-out replicate panels. The selection bound in the manuscript assumes joint coverage of the candidate losses.
+
+## Experiments and plots
+
+See [experiment commands](docs/EXPERIMENTS.md) for generation baselines, ablations, frozen external transfer, and figure exports. See [executed checks](docs/VALIDATION.md) for the actual validation scope and dataset sizes.
+
+Seven comparison rules are included alongside ReMEDi. These are uniform sampling, minimum predicted loss, Gibbs sampling, marginal uncertainty, absolute robust loss, Thompson sampling, and finite-scenario minimax. All use the same fitted response ensemble and query information.
+
+## Repository layout
+
+- `src/remedi/` contains ingestion, response heads, uncertainty estimation, optimization, metrics, and plotting.
+- `scripts/` contains source download, Tahoe conversion, sci-Plex subsetting, gene alignment, and complete-run commands.
+- `configs/` contains explicit study field mappings.
+- `tests/` contains mathematical, data-integrity, and pipeline checks.
+- `assets/` contains the overview figure and Ubuntu font license.
+
+Upload these files directly to the root of a Hugging Face model repository. Training runs locally or on your compute environment after installation. No server process or Space configuration is required. Downloaded datasets and trained checkpoints are excluded from source control by default.
+
+## Sources
+
+- [Tahoe-100M data and metadata](https://huggingface.co/datasets/tahoebio/Tahoe-100M)
+- [Corrected sci-Plex3 release](https://doi.org/10.5281/zenodo.13350497)
+- [scPerturb publication](https://doi.org/10.1038/s41592-023-02144-y)
+- [MolFormer checkpoint](https://huggingface.co/ibm/MoLFormer-XL-both-10pct)
+- [UCE implementation](https://github.com/snap-stanford/UCE)
+
+Dataset and encoder licenses remain those of their respective releases. Ubuntu is distributed under its included font license.
diff --git a/assets/fonts/UFL.txt b/assets/fonts/UFL.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ae78a8f94eae372cb1ca4e14b67244342fce604e
--- /dev/null
+++ b/assets/fonts/UFL.txt
@@ -0,0 +1,96 @@
+-------------------------------
+UBUNTU FONT LICENCE Version 1.0
+-------------------------------
+
+PREAMBLE
+This licence allows the licensed fonts to be used, studied, modified and
+redistributed freely. The fonts, including any derivative works, can be
+bundled, embedded, and redistributed provided the terms of this licence
+are met. The fonts and derivatives, however, cannot be released under
+any other licence. The requirement for fonts to remain under this
+licence does not require any document created using the fonts or their
+derivatives to be published under this licence, as long as the primary
+purpose of the document is not to be a vehicle for the distribution of
+the fonts.
+
+DEFINITIONS
+"Font Software" refers to the set of files released by the Copyright
+Holder(s) under this licence and clearly marked as such. This may
+include source files, build scripts and documentation.
+
+"Original Version" refers to the collection of Font Software components
+as received under this licence.
+
+"Modified Version" refers to any derivative made by adding to, deleting,
+or substituting -- in part or in whole -- any of the components of the
+Original Version, by changing formats or by porting the Font Software to
+a new environment.
+
+"Copyright Holder(s)" refers to all individuals and companies who have a
+copyright ownership of the Font Software.
+
+"Substantially Changed" refers to Modified Versions which can be easily
+identified as dissimilar to the Font Software by users of the Font
+Software comparing the Original Version with the Modified Version.
+
+To "Propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy. Propagation includes copying,
+distribution (with or without modification and with or without charging
+a redistribution fee), making available to the public, and in some
+countries other activities as well.
+
+PERMISSION & CONDITIONS
+This licence does not grant any rights under trademark law and all such
+rights are reserved.
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of the Font Software, to propagate the Font Software, subject to
+the below conditions:
+
+1) Each copy of the Font Software must contain the above copyright
+notice and this licence. These can be included either as stand-alone
+text files, human-readable headers or in the appropriate machine-
+readable metadata fields within text or binary files as long as those
+fields can be easily viewed by the user.
+
+2) The font name complies with the following:
+(a) The Original Version must retain its name, unmodified.
+(b) Modified Versions which are Substantially Changed must be renamed to
+avoid use of the name of the Original Version or similar names entirely.
+(c) Modified Versions which are not Substantially Changed must be
+renamed to both (i) retain the name of the Original Version and (ii) add
+additional naming elements to distinguish the Modified Version from the
+Original Version. The name of such Modified Versions must be the name of
+the Original Version, with "derivative X" where X represents the name of
+the new work, appended to that name.
+
+3) The name(s) of the Copyright Holder(s) and any contributor to the
+Font Software shall not be used to promote, endorse or advertise any
+Modified Version, except (i) as required by this licence, (ii) to
+acknowledge the contribution(s) of the Copyright Holder(s) or (iii) with
+their explicit written permission.
+
+4) The Font Software, modified or unmodified, in part or in whole, must
+be distributed entirely under this licence, and must not be distributed
+under any other licence. The requirement for fonts to remain under this
+licence does not affect any document created using the Font Software,
+except any version of the Font Software extracted from a document
+created using the Font Software may only be distributed under this
+licence.
+
+TERMINATION
+This licence becomes null and void if any of the above conditions are
+not met.
+
+DISCLAIMER
+THE FONT SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO ANY WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF
+COPYRIGHT, PATENT, TRADEMARK, OR OTHER RIGHT. IN NO EVENT SHALL THE
+COPYRIGHT HOLDER BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+INCLUDING ANY GENERAL, SPECIAL, INDIRECT, INCIDENTAL, OR CONSEQUENTIAL
+DAMAGES, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+FROM, OUT OF THE USE OR INABILITY TO USE THE FONT SOFTWARE OR FROM OTHER
+DEALINGS IN THE FONT SOFTWARE.
diff --git a/assets/fonts/Ubuntu-Bold.ttf b/assets/fonts/Ubuntu-Bold.ttf
new file mode 100644
index 0000000000000000000000000000000000000000..747fecc4350959d29aa23aef31f5c2efa644da81
--- /dev/null
+++ b/assets/fonts/Ubuntu-Bold.ttf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:679b5c1e09cab3156bb8ef529735f9382bf31ca7ac737382ab959297f8d82ad4
+size 331612
diff --git a/assets/fonts/Ubuntu-Regular.ttf b/assets/fonts/Ubuntu-Regular.ttf
new file mode 100644
index 0000000000000000000000000000000000000000..361486dccd2671c377d563ba6dfa179822228b52
--- /dev/null
+++ b/assets/fonts/Ubuntu-Regular.ttf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3128df86a31805618436d0ae5651ba4285d0c9de0a39057d025f64ee33bceb64
+size 351884
diff --git a/assets/remedi_overview.pdf b/assets/remedi_overview.pdf
new file mode 100644
index 0000000000000000000000000000000000000000..a1eefeb4ab8c9917a344511b6e13faad96e3080d
Binary files /dev/null and b/assets/remedi_overview.pdf differ
diff --git a/assets/remedi_overview.png b/assets/remedi_overview.png
new file mode 100644
index 0000000000000000000000000000000000000000..9c553326979294ec3f8221b792ec27108ee119ad
--- /dev/null
+++ b/assets/remedi_overview.png
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:001f90938f5cb0c848fe2622019f5fed68dfc15d23a181ce8598f7c5a3dfad3f
+size 156671
diff --git a/configs/sciplex.yaml b/configs/sciplex.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..894405d5ff6a925959c4fcc95885f99f3b84fb9d
--- /dev/null
+++ b/configs/sciplex.yaml
@@ -0,0 +1,10 @@
+study: sciplex3-corrected-v1.4
+drug: perturbation
+dose: dose_value
+dose_unit: nM
+context: cell_line
+block: replicate
+control_group: plate
+control_names: [control]
+# Supply an audited drug,smiles mapping via --structures.
+# For frozen UCE input also record cell_checkpoint: <model-id@commit>.
diff --git a/configs/tahoe.yaml b/configs/tahoe.yaml
new file mode 100644
index 0000000000000000000000000000000000000000..1802022f3c59b75212cf9b8489c9579dc85f1713
--- /dev/null
+++ b/configs/tahoe.yaml
@@ -0,0 +1,11 @@
+study: Tahoe-100M
+drug: drug
+dose: dose_um
+dose_unit: uM
+context: cell_line_id
+block: plate
+control_group: plate
+control_names: [DMSO_TF]
+smiles: canonical_smiles
+# The ingestion script must join sample_metadata to obtain the numeric dose_um.
+# Plate is the declared independent block here; audit the experimental design.
diff --git a/docs/EXPERIMENTS.md b/docs/EXPERIMENTS.md
new file mode 100644
index 0000000000000000000000000000000000000000..1d376e27d7072a7ccc62c9d24bcf1ce0f2f0739c
--- /dev/null
+++ b/docs/EXPERIMENTS.md
@@ -0,0 +1,71 @@
+# Experiments
+
+Use fixed molecular splits and record all seeds. Tune response-head parameters, temperature, and success thresholds on the `tune` partition. Fit the uncertainty radius on `calibration` after choosing those parameters.
+
+## Matched generation comparisons
+
+`remedi evaluate` runs uniform, mean, Gibbs, marginal, absolute, Thompson, finite-scenario minimax, and ReMEDi on identical candidate pools. `per_query.csv` contains endpoint loss, regret, success, best-of-10, dose regret, and timing. `summary.csv` contains molecule-bootstrap intervals. The query compound is excluded from the candidate pool at every dose.
+
+## Uncertainty ablations
+
+Calibrate each uncertainty variant before evaluation:
+
+```bash
+for variant in none permuted diagonal no-floor; do
+  remedi calibrate --data runs/tahoe/data --model runs/tahoe/model \
+    --output runs/tahoe/cal-$variant --ablation $variant
+  remedi evaluate --data runs/tahoe/data --model runs/tahoe/model \
+    --calibration runs/tahoe/cal-$variant --output runs/tahoe/eval-$variant \
+    --ablation $variant --tolerance 1
+done
+```
+
+`permuted` breaks scenario alignment across candidates. `diagonal` retains marginal variances. `no-floor` removes the residual variance term. `no-target-sampling` is an evaluation ablation with fixed target/control means. Compare `tau` values and ensemble sizes through separate calibration/evaluation runs.
+
+## Representations and response heads
+
+Prepare PCA and frozen per-cell UCE features under the same molecular split. Train fingerprint and cached MolFormer inputs with ridge and MLP heads. Keep an independent expression-space metric for comparisons across cell encoders. `evaluate --metric-data PATH` accepts aligned condition records with that representation and training-derived metric scaling.
+
+## Frozen Tahoe-to-sci-Plex transfer
+
+For PCA, align the raw cohorts on common gene identifiers before fitting:
+
+```bash
+python scripts/align_genes.py --inputs tahoe.h5ad sciplex.h5ad \
+  --gene-columns gene_symbol index --drop-ambiguous --output data/aligned
+```
+
+Fit the Tahoe projector on the aligned Tahoe training cells. Prepare sci-Plex using `--feature-model runs/tahoe/data/cell_feature_model.joblib`. Supply a sci-Plex split CSV containing `smiles,split` with `split=test`. Use an explicit structure mapping and harmonized gene IDs. Molecular overlaps are determined from canonical structures.
+
+```bash
+remedi prepare --h5ad data/aligned/1_sciplex.h5ad \
+  --mapping configs/sciplex.yaml --structures sciplex_structures.csv \
+  --splits sciplex_test.csv --feature-model runs/tahoe/data/cell_feature_model.joblib \
+  --output runs/external/data
+remedi evaluate --data runs/external/data --model runs/tahoe/model \
+  --calibration runs/tahoe/calibration --external-track unseen \
+  --output runs/external/evaluation --tolerance 1
+```
+
+For `unseen`, every test structure must be absent from Tahoe response-head training. For `shared`, every test structure must be present in that training set. Construct separate cohort files for the two tracks. Record exposure time and cell-line overlap. The scalar-dose head is fitted at the Tahoe exposure time. Transfer across a different exposure time measures an additional distribution shift.
+
+Frozen UCE preparation requires the same checkpoint for both studies. External treated cells are used for queries and final scoring. Fit response parameters, feature transforms, and calibration on the source study.
+
+## Dose, pool-size, and learning curves
+
+The current API predicts at any positive concentration. Scientific validation should begin with the observed dose grids. Middle-dose holdout requires a dedicated condition split in the training data, while keeping outer-dose measurements of those training compounds. The default molecular split is intentionally stricter and does not implement this condition-level experiment.
+
+For pool-size experiments, supply nested catalog CSVs and vary `--shortlist`. Record the candidate IDs, runtime, and solver residual for each run. For learning curves, reduce training compounds while preserving calibration and test compounds, then refit the response head and projector.
+
+```bash
+remedi plot-curve --csv scaling.csv --output figures/scaling.pdf \
+  --x candidates --y seconds --hue method
+```
+
+## Additional manuscript comparisons
+
+The manuscript also proposes CPA/chemCPA, scGen, scVI, transport models, COMs, and molecular-generator comparisons. Those models require their released implementations and compatible study splits. This package implements the ReMEDi computation and its matched generation rules. It does not include substitute implementations under those method names. Keep their measured predictions and resource usage in separate result files, with their checkpoint, version, and training split documented.
+
+## Interpretation
+
+A source-data smoke test establishes execution on measured RNA and dataset metadata. Performance conclusions require the full held-out study comparison. REAL candidates have predicted endpoints until those compounds have independent measured responses.
diff --git a/docs/VALIDATION.md b/docs/VALIDATION.md
new file mode 100644
index 0000000000000000000000000000000000000000..764765b86801a148f76f572a94abfc273b017450
--- /dev/null
+++ b/docs/VALIDATION.md
@@ -0,0 +1,49 @@
+# Executed validation
+
+The release was exercised on 12 September 2026 with Python 3.12 and CPU execution. Installed dependency versions are recorded in `requirements-tested.txt`.
+
+## Tests
+
+All 15 tests passed. These cover convex optimization, common-shift invariance, covariance calculations, molecular split integrity, independent control pools, concentration conversion, sparse counts, gene reordering, missing genes, and serialization. The synthetic end-to-end demo also completed.
+
+## Measured-data runs
+
+The following runs completed through the repository's `validate_release.py` command.
+
+| Run | Cells used | Compounds | Treated condition summaries | Queries | Candidate pairs per query |
+|---|---:|---:|---:|---:|---:|
+| Tahoe A549, plates 6 and 14 | 387 | 48 | 96 | 4 | 6 |
+| sci-Plex3 A549, both replicates | 120 | 9 | 36 | 4 | 2 |
+| Frozen Tahoe → sci-Plex3 | Source 387, external 120 | Source 48, external 9 | Source 96, external 36 | 4 | 16 |
+
+The Tahoe response ensemble was fitted on 26 compounds. The small sci-Plex ensemble was fitted on three compounds, with three others reserved for calibration and two for testing. Both study-specific runs completed all eight generation rules, CSV exports, catalog sampling, and Ubuntu plots. The catalog sampler was also exercised with a five-pair shortlist.
+
+For external transfer, the studies were aligned on 35,460 unique gene symbols. The PCA transform, response ensemble, and calibration were fitted using Tahoe. All nine external structures were absent from response-head training. Every ReMEDi query solve returned `optimal`. Maximum reported comparator-constraint violations were below 1e-6.
+
+Five uncertainty configurations completed on sci-Plex: the main configuration, permuted scenarios, diagonal covariance, omitted residual floor, and fixed target/control means. The forward-evaluation command completed the zero, training-mean, chemical-neighbor, and fitted-head comparisons on Tahoe.
+
+Measured result CSVs are included under `examples/results/`. These small cohorts establish execution on the released schemas and counts. Full-study accuracy, encoder comparisons, dose interpolation, and REAL candidate activity remain experimental work. The tests used fingerprints and eight-dimensional training-only PCA features. UCE and MolFormer checkpoint extraction were not executed in this environment.
+
+## Reproduction
+
+```bash
+pip install -e '.[data,dev]'
+python -m pytest
+python scripts/validate_release.py --output runs/validation
+python scripts/run_ablations.py \
+  --data runs/validation/sciplex/data --model runs/validation/sciplex/model \
+  --output runs/ablations --tolerance 1 --max-queries 2 --scenarios 16
+python scripts/evaluate_forward.py \
+  --data runs/validation/tahoe/data --model runs/validation/tahoe/model \
+  --output runs/forward
+```
+
+The fixtures include measured cells and source identifiers. The Tahoe fixture contains 568 cells before the two-block eligibility filter. The sci-Plex fixture uses two treated cells per condition and four controls per plate. Increase these counts for scientific evaluation.
+
+## Source acquisition
+
+Tahoe shards `train-00600-of-03388.parquet`, `train-00601-of-03388.parquet`, and `train-02600-of-03388.parquet` were downloaded from the official release and decoded with its sample and gene metadata. Reservoir sampling was restricted to A549 (`CVCL_0023`) with a cap of four cells per drug, concentration, and plate.
+
+sci-Plex was read from the corrected scPerturb v1.4 H5AD through cached HTTP byte ranges. The published full-file MD5 is recorded with the subset. Only subset bytes were retrieved for the successful range-based run, so a full-file checksum was not recomputed. The full-download command verifies the published MD5 after download.
+
+The standalone wheel built successfully and contains both Ubuntu fonts. The source archive excludes temporary downloads, caches, environments, and build products.
diff --git a/examples/ATTRIBUTION.md b/examples/ATTRIBUTION.md
new file mode 100644
index 0000000000000000000000000000000000000000..3fb6a6a3e78b67c7caa84106b59fc58520340c6d
--- /dev/null
+++ b/examples/ATTRIBUTION.md
@@ -0,0 +1,16 @@
+# Measured-data fixtures
+
+The included H5AD files contain subsets of released measured RNA counts.
+
+`examples/tahoe/tahoe.h5ad` was prepared from the Tahoe-100M expression shards 600, 601, and 2600, with the official sample and gene metadata. The subset contains A549 cells (CVCL_0023) from plates 6 and 14. Source records are retained in the observation fields.
+
+Source: https://huggingface.co/datasets/tahoebio/Tahoe-100M
+License listed by the source dataset: CC0-1.0.
+
+`examples/sciplex/sciplex_prototype.h5ad` contains 120 A549 cells from the corrected sci-Plex3 object in scPerturb v1.4, with 24-hour exposures at 0.1 and 1 µM and matched controls. `source_row` identifies the original row. The subset includes nine compound names with exact matches to Tahoe structure metadata.
+
+Source: Peidli, Green, Shen, Gross, Min, Garda, Yuan, Schumacher, Taylor-King, Marks, Luna, Blüthgen, and Sander. scPerturb Single-Cell Perturbation Data, version 1.4. https://doi.org/10.5281/zenodo.13350497
+License: Creative Commons Attribution 4.0 International. https://creativecommons.org/licenses/by/4.0/
+Original sci-Plex publication: Srivatsan et al. Massively multiplex chemical transcriptomics at single-cell resolution. Science (2020).
+
+The fixture sizes are intended for execution checks. They do not support full-study performance conclusions.
diff --git a/examples/results/ablations/diagonal/per_query.csv b/examples/results/ablations/diagonal/per_query.csv
new file mode 100644
index 0000000000000000000000000000000000000000..b819a3db2604ca9d16cf387bfbdd4b1db862fb4e
--- /dev/null
+++ b/examples/results/ablations/diagonal/per_query.csv
@@ -0,0 +1,17 @@
+query_id,target_molecule,target_dose_um,context,method,ablation,candidates,dose_regret,seconds,endpoint_loss,regret,success,best_of_10,entropy,effective_actions,unique_molecules_10,solver_status,solver_violation,estimated_worst_regret
+0,14c7df7d1e83fa4e24ed,0.1,A549,uniform,diagonal,2,0.6210683233378844,0.0010002739982155617,2.3458614971911134,0.6210683233378842,0.0,1.724793173853229,0.6931471805599453,2.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,mean,diagonal,2,1.2421366466757688,0.0009664329991210252,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,-0.0,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,gibbs,diagonal,2,1.2421366466757688,0.0010061940010928083,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,2.939500731829791e-22,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,marginal,diagonal,2,1.2421366466757688,0.0009615259987185709,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,2.0115195926483795e-17,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,absolute,diagonal,2,0.7948606990822286,0.007400439000775805,2.5196538729354576,0.7948606990822285,0.0,1.7434252235533654,0.6534676281378349,1.9221947422727113,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,thompson,diagonal,2,1.1645031062585334,0.0009758669984876178,2.8892962801117625,1.1645031062585334,0.0,2.3893362798247653,0.2337916587064593,1.2633812503027273,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,finite-minimax,diagonal,2,1.1837323893846905,0.004732806999527384,2.9085255632379194,1.1837323893846903,0.0,2.4887072115588267,0.18964308195639096,1.208818071670907,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,remedi,diagonal,2,0.8507494300861778,0.01236355199944228,2.575542603939407,0.8507494300861778,0.0,1.762057273253502,0.6231143789707472,1.8647264728169863,1.0,optimal,2.0418136514877006e-08,3.125500450998901
+1,08427e2fb6b299310c13,0.1,A549,uniform,diagonal,2,0.41460707956381565,0.0013860850012861192,1.3731752055534014,0.4146070795638156,0.5,0.9627141967852242,0.6931471805599453,2.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,mean,diagonal,2,0.8292141591276313,0.0009198529987770598,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,gibbs,diagonal,2,0.8292141591276313,0.0010094770004798193,1.7877822851172172,0.8292141591276313,4.257371470231585e-22,1.7877822851172172,2.0949767236225285e-20,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,marginal,diagonal,2,0.8292141591276313,0.0010812059990712442,1.7877822851172172,0.8292141591276313,7.472519673743183e-17,1.7877822851172172,2.7747493839213196e-15,1.0000000000000027,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,absolute,diagonal,2,0.5184445801103514,0.006109147001552628,1.4770127060999374,0.5184445801103515,0.37477601605865324,0.9627141967852242,0.6614487143861156,1.9375973272171612,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,thompson,diagonal,2,0.8292141591276313,0.0009567379966028966,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,finite-minimax,diagonal,2,0.8292141591276313,0.004650492999644484,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,remedi,diagonal,2,0.562287207703702,0.012452276001567952,1.5208553336932877,0.5622872077037019,0.32190351368909065,0.9627141967852242,0.6282959554744223,1.8744137717603098,1.0,optimal,4.236280215508259e-09,3.0155707541813337
diff --git a/examples/results/ablations/no-floor/per_query.csv b/examples/results/ablations/no-floor/per_query.csv
new file mode 100644
index 0000000000000000000000000000000000000000..58b3eafea6389593723564a14a5deb6d468cf2b0
--- /dev/null
+++ b/examples/results/ablations/no-floor/per_query.csv
@@ -0,0 +1,17 @@
+query_id,target_molecule,target_dose_um,context,method,ablation,candidates,dose_regret,seconds,endpoint_loss,regret,success,best_of_10,entropy,effective_actions,unique_molecules_10,solver_status,solver_violation,estimated_worst_regret
+0,14c7df7d1e83fa4e24ed,0.1,A549,uniform,no-floor,2,0.6210683233378844,0.0009935039997799322,2.3458614971911134,0.6210683233378842,0.0,1.724793173853229,0.6931471805599453,2.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,mean,no-floor,2,1.2421366466757688,0.0008672539988765493,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,-0.0,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,gibbs,no-floor,2,1.2421366466757688,0.001022868000291055,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,2.939500731829791e-22,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,marginal,no-floor,2,1.2421366466757688,0.0011904999992111698,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,3.464124802637856e-78,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,absolute,no-floor,2,1.2421366466757688,0.007212044998595957,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,-0.0,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,thompson,no-floor,2,1.1645031062585334,0.0008871739992173389,2.8892962801117625,1.1645031062585334,0.0,2.3893362798247653,0.2337916587064593,1.2633812503027273,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,finite-minimax,no-floor,2,1.1837323893846905,0.004739247000543401,2.9085255632379194,1.1837323893846903,0.0,2.4887072115588267,0.18964308195639096,1.208818071670907,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,remedi,no-floor,2,0.7227540312199966,0.012345715000265045,2.4475472050732257,0.7227540312199965,0.0,1.7310038570866078,0.6796833668146556,1.9732528359732782,1.0,optimal,6.688992737480248e-08,7.958927426358219
+1,08427e2fb6b299310c13,0.1,A549,uniform,no-floor,2,0.41460707956381565,0.0009695680018921848,1.3731752055534014,0.4146070795638156,0.5,0.9627141967852242,0.6931471805599453,2.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,mean,no-floor,2,0.8292141591276313,0.0008755570015637204,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,gibbs,no-floor,2,0.8292141591276313,0.0010444709987496026,1.7877822851172172,0.8292141591276313,4.257371470231585e-22,1.7877822851172172,2.0949767236225285e-20,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,marginal,no-floor,2,0.8292141591276313,0.0009344549980596639,1.7877822851172172,0.8292141591276313,8.144671043221572e-42,1.7877822851172172,7.705771828753628e-40,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,absolute,no-floor,2,0.8292141370902666,0.005874844999198103,1.7877822630798523,0.8292141370902665,2.65762039445289e-08,1.7877822851172172,4.901515627552698e-07,1.0000004901516828,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,thompson,no-floor,2,0.8292141591276313,0.0009403839976585004,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,finite-minimax,no-floor,2,0.8292141591276313,0.004586876999383094,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,remedi,no-floor,2,0.49497170229228493,0.011339101001794916,1.4535398282818708,0.49497170229228493,0.40308339306094787,0.9627141967852242,0.6742420852208934,1.9625449703126765,1.0,optimal,9.971349523851814e-08,6.108267267141067
diff --git a/examples/results/ablations/no-target-sampling/per_query.csv b/examples/results/ablations/no-target-sampling/per_query.csv
new file mode 100644
index 0000000000000000000000000000000000000000..5519629cebe18335a76c117fc1d8fbd5839263ce
--- /dev/null
+++ b/examples/results/ablations/no-target-sampling/per_query.csv
@@ -0,0 +1,17 @@
+query_id,target_molecule,target_dose_um,context,method,ablation,candidates,dose_regret,seconds,endpoint_loss,regret,success,best_of_10,entropy,effective_actions,unique_molecules_10,solver_status,solver_violation,estimated_worst_regret
+0,14c7df7d1e83fa4e24ed,0.1,A549,uniform,no-target-sampling,2,0.6210683233378844,0.0010524429999350104,2.3458614971911134,0.6210683233378842,0.0,1.724793173853229,0.6931471805599453,2.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,mean,no-target-sampling,2,1.2421366466757688,0.000918961002753349,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,-0.0,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,gibbs,no-target-sampling,2,1.2421366466757688,0.0010767389976535924,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,1.349209055803046e-25,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,marginal,no-target-sampling,2,1.2421366466757688,0.0009822470019571483,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,3.393817457709605e-15,1.0000000000000033,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,absolute,no-target-sampling,2,0.841146384709467,0.0062025279985391535,2.5659395585626963,0.8411463847094671,0.0,1.7558465900201232,0.6289790618495088,1.8756946331898454,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,thompson,no-target-sampling,2,1.2421366466757688,0.0008956370002124459,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,-0.0,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,finite-minimax,no-target-sampling,2,1.2421366466757688,0.004830543999560177,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,-0.0,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,remedi,no-target-sampling,2,0.8633844258074854,0.011631882000074256,2.5881775996607144,0.8633844258074852,0.0,1.762057273253502,0.6149754946725026,1.8496112736729622,1.0,optimal,1.2256626824580508e-09,3.334420846171371
+1,08427e2fb6b299310c13,0.1,A549,uniform,no-target-sampling,2,0.41460707956381565,0.0009755670034792274,1.3731752055534014,0.4146070795638156,0.5,0.9627141967852242,0.6931471805599453,2.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,mean,no-target-sampling,2,0.8292141591276313,0.0008587510019424371,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,gibbs,no-target-sampling,2,0.8292141591276313,0.0010410750001028646,1.7877822851172172,0.8292141591276313,1.931629832771395e-24,1.7877822851172172,1.0547409362089728e-22,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,marginal,no-target-sampling,2,0.8292141590832298,0.0009638890005589928,1.787782285072823,0.8292141590832373,5.3554332307996363e-11,1.7877822851172172,1.3201238337566032e-09,1.0000000013201238,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,absolute,no-target-sampling,2,0.488046143741546,0.005346689998987131,1.4466142697311317,0.48804614374154587,0.41143534710623925,0.9627141967852242,0.6773767064798296,1.968706457417316,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,thompson,no-target-sampling,2,0.8292141591276313,0.0009280560007027816,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,finite-minimax,no-target-sampling,2,0.8292141591206273,0.005026388997066533,1.7877822851102132,0.8292141591206273,8.44645767568175e-12,1.7877822851172172,2.238081112024106e-10,1.000000000223808,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,remedi,no-target-sampling,2,0.5717797137135874,0.011306774002150632,1.5303478397031733,0.5717797137135875,0.310455920922619,0.9710063383765005,0.6194649697958107,1.8579337252149741,1.0,optimal,0.0,3.0834980282590756
diff --git a/examples/results/ablations/none/per_query.csv b/examples/results/ablations/none/per_query.csv
new file mode 100644
index 0000000000000000000000000000000000000000..5216ab4edd95d5fe21be161009d605f61f0bc99b
--- /dev/null
+++ b/examples/results/ablations/none/per_query.csv
@@ -0,0 +1,17 @@
+query_id,target_molecule,target_dose_um,context,method,ablation,candidates,dose_regret,seconds,endpoint_loss,regret,success,best_of_10,entropy,effective_actions,unique_molecules_10,solver_status,solver_violation,estimated_worst_regret
+0,14c7df7d1e83fa4e24ed,0.1,A549,uniform,none,2,0.6210683233378844,0.001213707000715658,2.3458614971911134,0.6210683233378842,0.0,1.724793173853229,0.6931471805599453,2.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,mean,none,2,1.2421366466757688,0.0009324540005763993,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,-0.0,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,gibbs,none,2,1.2421366466757688,0.001097262000257615,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,2.939500731829791e-22,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,marginal,none,2,1.2421366466757688,0.0009660049981903285,2.9669298205289985,1.2421366466757693,0.0,2.9669298205289976,4.879519202998633e-15,1.0000000000000049,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,absolute,none,2,0.9272808925276437,0.15478037399952882,2.652074066380873,0.9272808925276439,0.0,1.805532055887154,0.5661252130573241,1.7614286505808048,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,thompson,none,2,1.1645031062585334,0.0012549089988169726,2.8892962801117625,1.1645031062585334,0.0,2.3893362798247653,0.2337916587064593,1.2633812503027273,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,finite-minimax,none,2,1.1837323893846905,0.005429281998658553,2.9085255632379194,1.1837323893846903,0.0,2.4887072115588267,0.18964308195639096,1.208818071670907,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,remedi,none,2,0.8549315315591413,0.012993177999305772,2.5797247054123704,0.8549315315591413,0.0,1.762057273253502,0.6204739656066904,1.8598093186344355,1.0,optimal,3.566613671068808e-11,3.051736597002783
+1,08427e2fb6b299310c13,0.1,A549,uniform,none,2,0.41460707956381565,0.0012252750020707026,1.3731752055534014,0.4146070795638156,0.5,0.9627141967852242,0.6931471805599453,2.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,mean,none,2,0.8292141591276313,0.0010998550023941789,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,gibbs,none,2,0.8292141591276313,0.0010322340021957643,1.7877822851172172,0.8292141591276313,4.257371470231585e-22,1.7877822851172172,2.0949767236225285e-20,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,marginal,none,2,0.8292141591276313,0.0009672560008766595,1.7877822851172456,0.8292141591276597,2.9642293057265165e-14,1.7877822851172172,9.233448080494037e-13,1.0000000000009233,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,absolute,none,2,0.5946627052722395,0.006787717000406701,1.5532308312618253,0.5946627052722394,0.2828599237887471,0.9792984799677769,0.5956341619845089,1.8141810649048173,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,thompson,none,2,0.8292141591276313,0.0009179420012515038,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,finite-minimax,none,2,0.8292141591276313,0.004703851001977455,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,remedi,none,2,0.5653046608129032,0.012059862998285098,1.5238727868024893,0.5653046608129034,0.3182645826891961,0.9627141967852242,0.6255544220847281,1.8692820414381144,1.0,optimal,9.820026125595405e-09,2.937462120748836
diff --git a/examples/results/ablations/permuted/per_query.csv b/examples/results/ablations/permuted/per_query.csv
new file mode 100644
index 0000000000000000000000000000000000000000..0a30243899c27a5f75dec9a5ae76c22084d08220
--- /dev/null
+++ b/examples/results/ablations/permuted/per_query.csv
@@ -0,0 +1,17 @@
+query_id,target_molecule,target_dose_um,context,method,ablation,candidates,dose_regret,seconds,endpoint_loss,regret,success,best_of_10,entropy,effective_actions,unique_molecules_10,solver_status,solver_violation,estimated_worst_regret
+0,14c7df7d1e83fa4e24ed,0.1,A549,uniform,permuted,2,0.6210683233378844,0.0016888719983398914,2.3458614971911134,0.6210683233378842,0.0,1.724793173853229,0.6931471805599453,2.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,mean,permuted,2,1.2421366466757688,0.0014768729997740593,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,-0.0,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,gibbs,permuted,2,1.2421366466757688,0.0016587919999437872,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,2.939500731829627e-22,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,marginal,permuted,2,1.2421366466757688,0.0015245649992721155,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,7.219710493680947e-28,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,absolute,permuted,2,0.8227257185467983,0.011945320999075193,2.5475188924000274,0.8227257185467982,0.0,1.7496359067867444,0.6394663103755487,1.8954690171708741,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,thompson,permuted,2,0.69870186375512,0.0014056660002097487,2.423495037608349,0.6987018637551197,0.0,1.724793173853229,0.6853142072764582,1.9843952490181156,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,finite-minimax,permuted,2,0.6467203757892618,0.0073726769987843,2.371513549642491,0.6467203757892617,0.0,1.724793173853229,0.6922939644801263,1.998294295611044,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,remedi,permuted,2,0.7797306792097395,0.027244690001680283,2.5045238530629685,0.7797306792097394,0.0,1.7372145403199868,0.6601509847739722,1.9350844806395862,1.0,optimal,0.0,4.899262250810139
+1,08427e2fb6b299310c13,0.1,A549,uniform,permuted,2,0.41460707956381565,0.001529922999907285,1.3731752055534014,0.4146070795638156,0.5,0.9627141967852242,0.6931471805599453,2.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,mean,permuted,2,0.8292141591276313,0.0010190739994868636,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,gibbs,permuted,2,0.8292141591276313,0.0013229409996711183,1.7877822851172172,0.8292141591276313,4.257371470232068e-22,1.7877822851172172,2.0949767236227617e-20,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,marginal,permuted,2,0.8292141591276313,0.001460939001844963,1.7877822851172172,0.8292141591276313,5.593146755023605e-29,1.7877822851172172,3.638533565342826e-27,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,absolute,permuted,2,0.5552809701454036,0.007075068999256473,1.5138490961349893,0.5552809701454035,0.33035276347719045,0.9627141967852242,0.6344281761707704,1.8859434056847724,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,thompson,permuted,2,0.3627811946183387,0.0009871460024442058,1.3213493206079245,0.3627811946183387,0.5625,0.958568125989586,0.6853142072764582,1.9843952490181156,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,finite-minimax,permuted,2,0.5023501576546154,0.005492406999110244,1.4609182836442014,0.5023501576546155,0.39418526308920104,0.9627141967852242,0.6705834384362565,1.9553778304786524,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,remedi,permuted,2,0.5243882113935283,0.01714433699817164,1.4829563373831143,0.5243882113935284,0.3676082280780069,0.9627141967852242,0.6576704578335897,1.930290399834909,1.0,optimal,0.0,4.3203418079097275
diff --git a/examples/results/forward/forward_summary.csv b/examples/results/forward/forward_summary.csv
new file mode 100644
index 0000000000000000000000000000000000000000..91517a69198faeae421812780a3fe55f902def02
--- /dev/null
+++ b/examples/results/forward/forward_summary.csv
@@ -0,0 +1,5 @@
+method,mse,pearson
+chemical-neighbor,3.6675466897442717,0.40206437900880526
+response-head,1.1973462467667448,0.6604562483673746
+training-mean,1.3487781295516517,0.5797383696951393
+zero,2.1842505254797233,
diff --git a/examples/results/sciplex/per_query.csv b/examples/results/sciplex/per_query.csv
new file mode 100644
index 0000000000000000000000000000000000000000..03b163ca185f4a68831e36482c6cd99a922f4915
--- /dev/null
+++ b/examples/results/sciplex/per_query.csv
@@ -0,0 +1,33 @@
+query_id,target_molecule,target_dose_um,context,method,ablation,candidates,dose_regret,seconds,endpoint_loss,regret,success,best_of_10,entropy,effective_actions,unique_molecules_10,solver_status,solver_violation,estimated_worst_regret
+0,14c7df7d1e83fa4e24ed,0.1,A549,uniform,none,2,0.6210683233378844,0.0010568409998086281,2.3458614971911134,0.6210683233378842,0.0,1.724793173853229,0.6931471805599453,2.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,mean,none,2,1.2421366466757688,0.0009057839997694828,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,-0.0,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,gibbs,none,2,1.2421366466757688,0.0010445929983688984,2.966929820528998,1.2421366466757688,0.0,2.9669298205289976,2.939500731829791e-22,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,marginal,none,2,1.2421366466757688,0.0012421399987943005,2.9669298205289985,1.2421366466757693,0.0,2.9669298205289976,4.879519202998633e-15,1.0000000000000049,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,absolute,none,2,0.9272808925276437,0.1787452469980053,2.652074066380873,0.9272808925276439,0.0,1.805532055887154,0.5661252130573241,1.7614286505808048,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,thompson,none,2,1.1645031062585334,0.0012781340010405984,2.8892962801117625,1.1645031062585334,0.0,2.3893362798247653,0.2337916587064593,1.2633812503027273,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,finite-minimax,none,2,1.1837323893846905,0.005864108003152069,2.9085255632379194,1.1837323893846903,0.0,2.4887072115588267,0.18964308195639096,1.208818071670907,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,remedi,none,2,0.8549315315591413,0.01353932299753069,2.5797247054123704,0.8549315315591413,0.0,1.762057273253502,0.6204739656066904,1.8598093186344355,1.0,optimal,3.566613671068808e-11,3.051736597002783
+1,08427e2fb6b299310c13,0.1,A549,uniform,none,2,0.41460707956381565,0.0010660150001058355,1.3731752055534014,0.4146070795638156,0.5,0.9627141967852242,0.6931471805599453,2.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,mean,none,2,0.8292141591276313,0.0011393759996281005,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,gibbs,none,2,0.8292141591276313,0.002222175000497373,1.7877822851172172,0.8292141591276313,4.257371470231585e-22,1.7877822851172172,2.0949767236225285e-20,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,marginal,none,2,0.8292141591276313,0.0012320749992795754,1.7877822851172456,0.8292141591276597,2.9642293057265165e-14,1.7877822851172172,9.233448080494037e-13,1.0000000000009233,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,absolute,none,2,0.5946627052722395,0.007546040000306675,1.5532308312618253,0.5946627052722394,0.2828599237887471,0.9792984799677769,0.5956341619845089,1.8141810649048173,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,thompson,none,2,0.8292141591276313,0.0009839109989115968,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,finite-minimax,none,2,0.8292141591276313,0.005074109001725446,1.7877822851172172,0.8292141591276313,0.0,1.7877822851172172,-0.0,1.0,1.0,,,
+1,08427e2fb6b299310c13,0.1,A549,remedi,none,2,0.5653046608129032,0.013257608999992954,1.5238727868024893,0.5653046608129034,0.3182645826891961,0.9627141967852242,0.6255544220847281,1.8692820414381144,1.0,optimal,9.820026125595405e-09,2.937462120748836
+2,08427e2fb6b299310c13,1.0,A549,uniform,none,2,0.037573553499094414,0.001326116998825455,1.0837318015315083,0.037573553499094414,0.0,1.0461582480324143,0.6931471805599453,2.0,1.0,,,
+2,08427e2fb6b299310c13,1.0,A549,mean,none,2,0.0,0.0009192240031552501,1.0461582480324139,0.0,0.0,1.0461582480324143,-0.0,1.0,1.0,,,
+2,08427e2fb6b299310c13,1.0,A549,gibbs,none,2,6.325270427900244e-36,0.0011685799981933087,1.0461582480324139,0.0,0.0,1.0461582480324143,6.604139813057347e-33,1.0,1.0,,,
+2,08427e2fb6b299310c13,1.0,A549,marginal,none,2,6.721598183927344e-35,0.001084221999917645,1.0461582480324139,0.0,0.0,1.0461582480324143,6.806548291629058e-32,1.0,1.0,,,
+2,08427e2fb6b299310c13,1.0,A549,absolute,none,2,3.6851652224222005e-09,0.005545459000131814,1.0461582517175791,3.6851652840397264e-09,0.0,1.0461582480324143,8.744031185739841e-07,1.000000874403501,1.0,,,
+2,08427e2fb6b299310c13,1.0,A549,thompson,none,2,0.004696694187386802,0.001042159001372056,1.0508549422198006,0.0046966941873867185,0.0,1.0461582480324143,0.2337916587064593,1.2633812503027273,1.0,,,
+2,08427e2fb6b299310c13,1.0,A549,finite-minimax,none,2,0.002287350492080001,0.005708774999220623,1.048445598524494,0.0022873504920801757,0.0,1.0461582480324143,0.1362624709487634,1.1459826412754928,1.0,,,
+2,08427e2fb6b299310c13,1.0,A549,remedi,none,2,0.014531819461003686,0.01251371000034851,1.0606900674934177,0.014531819461003792,0.0,1.0461582480324143,0.4910846152148521,1.63408761543058,1.0,optimal,6.64557964213941e-10,1.9956950059953673
+3,14c7df7d1e83fa4e24ed,1.0,A549,uniform,none,2,0.5223619289032562,0.0010797450013342313,2.428722642442487,0.5223619289032562,0.0,1.90636071353923,0.6931471805599453,2.0,1.0,,,
+3,14c7df7d1e83fa4e24ed,1.0,A549,mean,none,2,1.0447238578065123,0.0009034300019266084,2.951084571345743,1.0447238578065123,0.0,2.951084571345743,-0.0,1.0,1.0,,,
+3,14c7df7d1e83fa4e24ed,1.0,A549,gibbs,none,2,1.0447238578065123,0.001108239001041511,2.951084571345743,1.0447238578065123,0.0,2.951084571345743,3.7090613682420546e-21,1.0,1.0,,,
+3,14c7df7d1e83fa4e24ed,1.0,A549,marginal,none,2,1.0447238578065123,0.000993975998426322,2.951084571345743,1.0447238578065123,0.0,2.951084571345743,1.6365442547476746e-29,1.0,1.0,,,
+3,14c7df7d1e83fa4e24ed,1.0,A549,absolute,none,2,1.0447237959520446,0.005720912999095162,2.951084509491275,1.0447237959520446,0.0,2.951084571345743,1.044535347615547e-06,1.0000010445358931,1.0,,,
+3,14c7df7d1e83fa4e24ed,1.0,A549,thompson,none,2,0.9794286166936053,0.0008423079998465255,2.8857893302328357,0.979428616693605,0.0,2.4339462617315193,0.2337916587064593,1.2633812503027273,1.0,,,
+3,14c7df7d1e83fa4e24ed,1.0,A549,finite-minimax,none,2,0.9723551414505388,0.004529568999714684,2.87871585498977,0.9723551414505391,0.0,2.386933688130226,0.25174828513945446,1.2862722226989185,1.0,,,
+3,14c7df7d1e83fa4e24ed,1.0,A549,remedi,none,2,0.8333596309969656,0.01226573599706171,2.739720344536196,0.8333596309969655,0.0,2.0473984343431093,0.5035962030371717,1.65466108112194,1.0,optimal,0.0,1.3817553061677006
diff --git a/examples/results/sciplex/summary.csv b/examples/results/sciplex/summary.csv
new file mode 100644
index 0000000000000000000000000000000000000000..92fb679e4da47aa8163ea9652a9fbb48648f1444
--- /dev/null
+++ b/examples/results/sciplex/summary.csv
@@ -0,0 +1,57 @@
+method,metric,mean,ci_low,ci_high,molecules
+absolute,endpoint_loss,2.0506369147128884,1.299694541489702,2.8015792879360744,2
+finite-minimax,endpoint_loss,2.1558673254673497,1.4181139418208555,2.8936207091138444,2
+gibbs,endpoint_loss,2.187988731256093,1.4169702665748156,2.9590071959373705,2
+marginal,endpoint_loss,2.1879887312561004,1.4169702665748298,2.959007195937371,2
+mean,endpoint_loss,2.187988731256093,1.4169702665748156,2.9590071959373705,2
+remedi,endpoint_loss,1.9760019760611183,1.2922814271479535,2.659722524974283,2
+thompson,endpoint_loss,2.1534307094204044,1.419318613668509,2.8875428051722993,2
+uniform,endpoint_loss,1.8078727866796274,1.2284535035424549,2.3872920698168,2
+absolute,regret,0.6416668493592733,0.29733135447870235,0.9860023442398442,2
+finite-minimax,regret,0.7468972601137353,0.41575075480985574,1.0780437654176147,2
+gibbs,regret,0.7790186659024781,0.41460707956381565,1.1434302522411406,2
+marginal,regret,0.7790186659024854,0.41460707956382986,1.1434302522411408,2
+mean,regret,0.7790186659024781,0.41460707956381565,1.1434302522411406,2
+remedi,regret,0.5670319107075035,0.2899182401369536,0.8441455812780534,2
+thompson,regret,0.7444606440667891,0.416955426657509,1.0719658614760692,2
+uniform,regret,0.39890272132601257,0.226090316531455,0.5717151261205702,2
+absolute,success,0.07071498094718677,0.0,0.14142996189437354,2
+finite-minimax,success,0.0,0.0,0.0,2
+gibbs,success,1.0643428675578962e-22,0.0,2.1286857351157924e-22,2
+marginal,success,7.410573264316291e-15,0.0,1.4821146528632583e-14,2
+mean,success,0.0,0.0,0.0,2
+remedi,success,0.07956614567229903,0.0,0.15913229134459805,2
+thompson,success,0.0,0.0,0.0,2
+uniform,success,0.125,0.0,0.25,2
+absolute,best_of_10,1.695518338808272,1.0127283640000955,2.3783083136164485,2
+finite-minimax,best_of_10,1.927395358209671,1.4169702665748156,2.4378204498445264,2
+gibbs,best_of_10,2.187988731256093,1.4169702665748156,2.95900719593737,2
+marginal,best_of_10,2.187988731256093,1.4169702665748156,2.95900719593737,2
+mean,best_of_10,2.187988731256093,1.4169702665748156,2.95900719593737,2
+remedi,best_of_10,1.4545820381035623,1.0044362224088192,1.9047278537983057,2
+thompson,best_of_10,1.914305768676479,1.4169702665748156,2.4116412707781425,2
+uniform,best_of_10,1.4100065830525244,1.0044362224088192,1.8155769436962295,2
+absolute,dose_regret,0.6416668493592732,0.29733135447870235,0.9860023442398441,2
+finite-minimax,dose_regret,0.7468972601137351,0.4157507548098556,1.0780437654176147,2
+gibbs,dose_regret,0.7790186659024781,0.41460707956381565,1.1434302522411406,2
+marginal,dose_regret,0.7790186659024781,0.41460707956381565,1.1434302522411406,2
+mean,dose_regret,0.7790186659024781,0.41460707956381565,1.1434302522411406,2
+remedi,dose_regret,0.5670319107075035,0.28991824013695344,0.8441455812780534,2
+thompson,dose_regret,0.7444606440667891,0.41695542665750907,1.0719658614760692,2
+uniform,dose_regret,0.3989027213260127,0.22609031653145503,0.5717151261205703,2
+absolute,entropy,0.2904403234950748,0.2830631287963359,0.29781751819381375,2
+finite-minimax,entropy,0.1444134595111522,0.0681312354743817,0.2206956835479227,2
+gibbs,entropy,6.238194669414231e-21,2.001505720712517e-21,1.0474883618115945e-20,2
+marginal,entropy,2.320560818131006e-13,2.439759601499325e-15,4.616724040247019e-13,2
+mean,entropy,0.0,0.0,0.0,2
+remedi,entropy,0.5601773014858606,0.5583195186497901,0.562035084321931,2
+thompson,entropy,0.17534374402984448,0.11689582935322965,0.2337916587064593,2
+uniform,entropy,0.6931471805599453,0.6931471805599453,0.6931471805599453,2
+absolute,unique_molecules_10,1.0,1.0,1.0,2
+finite-minimax,unique_molecules_10,1.0,1.0,1.0,2
+gibbs,unique_molecules_10,1.0,1.0,1.0,2
+marginal,unique_molecules_10,1.0,1.0,1.0,2
+mean,unique_molecules_10,1.0,1.0,1.0,2
+remedi,unique_molecules_10,1.0,1.0,1.0,2
+thompson,unique_molecules_10,1.0,1.0,1.0,2
+uniform,unique_molecules_10,1.0,1.0,1.0,2
diff --git a/examples/results/tahoe/per_query.csv b/examples/results/tahoe/per_query.csv
new file mode 100644
index 0000000000000000000000000000000000000000..d5e790c9371473ec8c574e333fea30f584c9b23c
--- /dev/null
+++ b/examples/results/tahoe/per_query.csv
@@ -0,0 +1,33 @@
+query_id,target_molecule,target_dose_um,context,method,ablation,candidates,dose_regret,seconds,endpoint_loss,regret,success,best_of_10,entropy,effective_actions,unique_molecules_10,solver_status,solver_violation,estimated_worst_regret
+0,4d2dddf7eaaab030af74,5.0,CVCL_0023,uniform,none,6,0.0,0.0013273689983179793,1.4797382676103594,0.8229970926981272,0.3333333333333333,0.6818549554841797,1.7917594692280547,5.999999999999998,5.030966502660926,,,
+0,4d2dddf7eaaab030af74,5.0,CVCL_0023,mean,none,6,0.0,0.0010496300019440241,0.8361253218547141,0.17938414694248184,1.0,0.8361253218547142,-0.0,1.0,1.0,,,
+0,4d2dddf7eaaab030af74,5.0,CVCL_0023,gibbs,none,6,0.0,0.0012023710005450994,0.8437882195275845,0.18704704461535226,0.904385896615232,0.660328857851082,0.964014801904684,2.622202994086098,2.7102352149497677,,,
+0,4d2dddf7eaaab030af74,5.0,CVCL_0023,marginal,none,6,0.0,0.0011351690009178128,0.7723094618408193,0.115568286928587,0.9877254010412058,0.6576380956469449,0.7526491311714294,2.122615662353522,2.115452159537401,,,
+0,4d2dddf7eaaab030af74,5.0,CVCL_0023,absolute,none,6,0.0,0.1759201889981341,0.8938512037512288,0.23711002883899657,0.8510410125018211,0.6576380956469446,1.151514064974648,3.162978240694389,3.1071433102924257,,,
+0,4d2dddf7eaaab030af74,5.0,CVCL_0023,thompson,none,6,0.0,0.001231784997798968,0.9635401592192263,0.30679898430699404,0.6875,0.6576380956469446,1.2511559833636317,3.4943800688212554,3.386541112971827,,,
+0,4d2dddf7eaaab030af74,5.0,CVCL_0023,finite-minimax,none,6,0.0,0.005598998999630567,1.0640742220456179,0.4073330471333856,0.4650288000851047,0.6773703518106177,1.1152006558760126,3.0501801543532907,3.0690752592099746,,,
+0,4d2dddf7eaaab030af74,5.0,CVCL_0023,remedi,none,6,0.0,0.021281839999574004,0.9479419169923252,0.29120074208009294,0.6918164639744862,0.6719888274023432,1.1353011500543702,3.1121106140791635,3.0516120561961366,optimal,0.0,0.1805705751537281
+1,952d4eb68fe7ff02293e,5.0,CVCL_0023,uniform,none,6,0.0,0.001499417001468828,2.3186385688646336,1.863473864868645,0.16666666666666666,0.7939951969861689,1.7917594692280547,5.999999999999998,5.030966502660926,,,
+1,952d4eb68fe7ff02293e,5.0,CVCL_0023,mean,none,6,0.0,0.001110150999011239,0.45516470399598863,0.0,1.0,0.4551647039959886,-0.0,1.0,1.0,,,
+1,952d4eb68fe7ff02293e,5.0,CVCL_0023,gibbs,none,6,0.0,0.0012694809993263334,1.1493566627048921,0.6941919587089035,0.7150570849791218,0.4551647039959886,0.6130092370519638,1.8459780345237102,1.9888700487191797,,,
+1,952d4eb68fe7ff02293e,5.0,CVCL_0023,marginal,none,6,0.0,0.0012514839982031845,0.814321378458233,0.3591566744622444,0.8527255077375564,0.4551647039959886,0.4225841351547329,1.5258995960015183,1.801780461663748,,,
+1,952d4eb68fe7ff02293e,5.0,CVCL_0023,absolute,none,6,0.0,0.006859246997919399,1.3636923643606615,0.9085276603646729,0.6227596603882732,0.4551647039959886,0.7602652466647268,2.138843466359281,2.204928400110453,,,
+1,952d4eb68fe7ff02293e,5.0,CVCL_0023,thompson,none,6,0.0,0.0010454939983901568,2.23359208609533,1.7784273820993417,0.25,0.5751017060380507,1.2130075659799042,3.363585661014858,3.4165587704628706,,,
+1,952d4eb68fe7ff02293e,5.0,CVCL_0023,finite-minimax,none,6,0.0,0.005003091999242315,1.4319966741921175,0.9768319701961289,0.5958061003331557,0.4551647039959886,0.7667494824698613,2.1527572932618186,2.1856554603285705,,,
+1,952d4eb68fe7ff02293e,5.0,CVCL_0023,remedi,none,6,0.0,0.012374489997455385,1.5029403424617436,1.047775638465755,0.5702032760750585,0.4551647039959886,0.6995428689912393,2.012832369037561,2.020686665733723,optimal,0.0,0.15921518655635364
+2,7e32a4923f53cf86b1ac,5.0,CVCL_0023,uniform,none,6,0.0,0.0012619500012078788,1.843886655203399,1.1871454802911667,0.16666666666666666,0.750352056712608,1.7917594692280547,5.999999999999998,5.030966502660926,,,
+2,7e32a4923f53cf86b1ac,5.0,CVCL_0023,mean,none,6,0.0,0.0010836020010174252,1.1844419840640414,0.5277008091518092,0.0,1.1844419840640412,-0.0,1.0,1.0,,,
+2,7e32a4923f53cf86b1ac,5.0,CVCL_0023,gibbs,none,6,0.0,0.0011860160011565313,1.811495327413881,1.1547541525016487,0.00018426967577716134,1.1844419840640412,0.9087260139024683,2.4811595562377593,2.579600967875729,,,
+2,7e32a4923f53cf86b1ac,5.0,CVCL_0023,marginal,none,6,0.0,0.0012810789994546212,2.6627841000681203,2.006042925155888,6.556237635930052e-05,1.7707369188718212,0.48754927434281414,1.6283207586040636,1.947327568468316,,,
+2,7e32a4923f53cf86b1ac,5.0,CVCL_0023,absolute,none,6,0.0,0.006626345999393379,2.2712622465978116,1.6145210716855793,0.0020254421148502037,1.2018175342903652,1.3729351889797083,3.946918662242241,3.7819180616639283,,,
+2,7e32a4923f53cf86b1ac,5.0,CVCL_0023,thompson,none,6,0.0,0.0011304519975965377,1.738914945707029,1.0821737707947967,0.125,0.7889261754258163,1.461472129304024,4.31230312570097,4.021673919129171,,,
+2,7e32a4923f53cf86b1ac,5.0,CVCL_0023,finite-minimax,none,6,0.0,0.005141850997461006,2.1717990799792997,1.5150579050670674,0.07573218629383602,0.9316162860861898,1.5381906792239834,4.6561581613950676,4.3312157166056995,,,
+2,7e32a4923f53cf86b1ac,5.0,CVCL_0023,remedi,none,6,0.0,0.019209931000659708,2.0099986658548454,1.3532574909426132,0.00022779351908091826,1.1791649759725231,1.1190650129027855,3.0619899433897673,3.031866728797554,optimal,0.0,0.2642535768314859
+3,c2b4eac6bdcc4471946b,5.0,CVCL_0023,uniform,none,6,0.0,0.0013346289997571148,1.7108286910786454,1.2556639870826567,0.16666666666666666,0.6387157319132072,1.7917594692280547,5.999999999999998,5.030966502660926,,,
+3,c2b4eac6bdcc4471946b,5.0,CVCL_0023,mean,none,6,0.0,0.0010698209989641327,1.6686370480157346,1.213472344019746,0.0,1.6686370480157349,-0.0,1.0,1.0,,,
+3,c2b4eac6bdcc4471946b,5.0,CVCL_0023,gibbs,none,6,0.0,0.001193427000544034,1.5576594632468097,1.102494759250821,0.1727505918538528,0.6553876407592466,1.2242500690614362,3.4016141503314765,3.493234091303685,,,
+3,c2b4eac6bdcc4471946b,5.0,CVCL_0023,marginal,none,6,0.0,0.0011461649992270395,1.6951754027700967,1.2400106987741082,0.022429546181702533,1.4684141112524762,0.613016062151079,1.845990633549755,2.2323819363202313,,,
+3,c2b4eac6bdcc4471946b,5.0,CVCL_0023,absolute,none,6,0.0,0.0064251929979946,1.6275169218522685,1.1723522178562797,0.13074741532920992,0.78280223688132,1.3085825014018688,3.7009239378216603,3.731772125370696,,,
+3,c2b4eac6bdcc4471946b,5.0,CVCL_0023,thompson,none,6,0.0,0.0011979139999311883,1.8060006554754133,1.3508359514794246,0.0625,1.0877004992005133,1.6843729712135975,5.389070770225885,4.675035286509228,,,
+3,c2b4eac6bdcc4471946b,5.0,CVCL_0023,finite-minimax,none,6,0.0,0.004858926000451902,1.7321354866050187,1.27697078260903,0.15914354634750294,0.7039265345200366,1.5679213874531976,4.796667411215066,4.249188314216482,,,
+3,c2b4eac6bdcc4471946b,5.0,CVCL_0023,remedi,none,6,0.0,0.02921563400013838,1.3783086766733001,0.9231439726773115,0.33147565921296035,0.4733667891562848,1.4538691089986715,4.2796409205169565,3.974118065400982,optimal,0.0,0.20126312730353998
diff --git a/examples/results/tahoe/summary.csv b/examples/results/tahoe/summary.csv
new file mode 100644
index 0000000000000000000000000000000000000000..a845b99483058af314689bbaa5c8a5e9efb5b489
--- /dev/null
+++ b/examples/results/tahoe/summary.csv
@@ -0,0 +1,57 @@
+method,metric,mean,ci_low,ci_high,molecules
+absolute,endpoint_loss,1.5390806841404927,1.0772676332764888,2.0443697760385238,4
+finite-minimax,endpoint_loss,1.6000013657055137,1.231089538185468,1.986848478532504,4
+gibbs,endpoint_loss,1.3405749182232918,0.9965724411162383,1.6845773953303453,4
+marginal,endpoint_loss,1.4861475857843174,0.7933154201495263,2.2006684196656483,4
+mean,endpoint_loss,1.0360922644826196,0.6374840240130017,1.4605091164754795,4
+remedi,endpoint_loss,1.4597974004955536,1.08669152335968,1.8520761685594591,4
+thompson,endpoint_loss,1.6855119616242498,1.1741552832832731,2.109922800998255,4
+uniform,endpoint_loss,1.8382730456892593,1.5707753645086195,2.1666860994181363,4
+absolute,regret,0.9831277446863822,0.4709205760933174,1.4380227188553527,4
+finite-minimax,regret,1.044048426251403,0.6247424810022967,1.3960143438380488,4
+gibbs,regret,0.7846219787691814,0.41590897327421944,1.1286244558762348,4
+marginal,regret,0.9301946463302069,0.2373624806954157,1.6230268119649982,4
+mean,regret,0.4801393250285092,0.08969207347124092,0.9549502947504299,4
+remedi,regret,0.9038444610414432,0.4803444661765085,1.2457291113762876,4
+thompson,regret,1.1295590221701393,0.5678082261001016,1.6043639792732056,4
+uniform,regret,1.282320106235149,0.9311638162942597,1.6943917687242755,4
+absolute,success,0.4016433825835386,0.06638642872203006,0.7369003364450472,4
+finite-minimax,success,0.3239276582648999,0.11743786632066948,0.5304174502091302,4
+gibbs,success,0.4480944607809959,0.08646743076481499,0.8097214907971769,4
+marginal,success,0.46573650433420605,0.011247554279030916,0.9202254543893811,4
+mean,success,0.5,0.0,1.0,4
+remedi,success,0.3984307981953965,0.1427216641580753,0.6310098700247724,4
+thompson,success,0.28125,0.09375,0.546875,4
+uniform,success,0.20833333333333331,0.16666666666666666,0.29166666666666663,4
+absolute,best_of_10,0.7743556427036546,0.5370740872173214,1.06577267462951,4
+finite-minimax,best_of_10,0.6920194691032082,0.5173551616270006,0.8680548025172968,4
+gibbs,best_of_10,0.7388307966675897,0.5064557424597619,1.0521783982378425,4
+marginal,best_of_10,1.0879884574418077,0.5564013998214667,1.6195755150621487,4
+mean,best_of_10,1.0360922644826198,0.6374840240130017,1.4605091164754798,4
+remedi,best_of_10,0.6949213241317849,0.4642657465761367,1.0027154292684635,4
+thompson,best_of_10,0.7773416190778313,0.6163699008424977,0.9801848983121211,4
+uniform,best_of_10,0.716229485274041,0.6602853436986935,0.7721736268493884,4
+absolute,dose_regret,0.0,0.0,0.0,4
+finite-minimax,dose_regret,0.0,0.0,0.0,4
+gibbs,dose_regret,0.0,0.0,0.0,4
+marginal,dose_regret,0.0,0.0,0.0,4
+mean,dose_regret,0.0,0.0,0.0,4
+remedi,dose_regret,0.0,0.0,0.0,4
+thompson,dose_regret,0.0,0.0,0.0,4
+uniform,dose_regret,0.0,0.0,0.0,4
+absolute,entropy,1.148324250505238,0.8973445603490123,1.3407588451907886,4
+finite-minimax,entropy,1.2470155512557637,0.9409750691729369,1.5530560333385908,4
+gibbs,entropy,0.9275000304801381,0.7007606282651438,1.1453690552716942,4
+marginal,entropy,0.5689496507050138,0.45506670474877353,0.6863741669642756,4
+mean,entropy,0.0,0.0,0.0,4
+remedi,entropy,1.1019445352367665,0.808482439257022,1.3701680849747,4
+thompson,entropy,1.4025021624652894,1.232081774671768,1.576068724251106,4
+uniform,entropy,1.7917594692280547,1.7917594692280547,1.7917594692280547,4
+absolute,unique_molecules_10,3.2064404743593755,2.5866393314255136,3.7568450935173123,4
+finite-minimax,unique_molecules_10,3.4587836875901816,2.6273653597692723,4.290202015411091,4
+gibbs,unique_molecules_10,2.6929850807120905,2.169211340276827,3.2648258104466956,4
+marginal,unique_molecules_10,2.024235531497424,1.8745540150660318,2.173917047928816,4
+mean,unique_molecules_10,1.0,1.0,1.0,4
+remedi,unique_molecules_10,3.019570879032099,2.2784180133493264,3.7385552312501247,4
+thompson,unique_molecules_10,3.874952272268274,3.4015499417173487,4.360416157497639,4
+uniform,unique_molecules_10,5.030966502660926,5.030966502660926,5.030966502660926,4
diff --git a/examples/results/transfer/per_query.csv b/examples/results/transfer/per_query.csv
new file mode 100644
index 0000000000000000000000000000000000000000..27488142d33a80a386c0ff151dda6cebc7d460b7
--- /dev/null
+++ b/examples/results/transfer/per_query.csv
@@ -0,0 +1,33 @@
+query_id,target_molecule,target_dose_um,context,method,ablation,candidates,dose_regret,seconds,endpoint_loss,regret,success,best_of_10,entropy,effective_actions,unique_molecules_10,solver_status,solver_violation,estimated_worst_regret
+0,14c7df7d1e83fa4e24ed,0.1,A549,uniform,none,16,0.6443150725027152,0.0016497040014655795,2.876180395105572,2.0354774023754407,0.0625,1.4258732694102405,2.772588722239781,15.999999999999998,5.895395390689373,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,mean,none,16,0.10511493374206449,0.001140697997470852,3.1576328120475896,2.3169298193174583,0.0,3.1576328120475887,-0.0,1.0,1.0,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,gibbs,none,16,0.2504362451458414,0.0013304339991009329,2.2138232327131813,1.3731202399830498,0.4322252267000464,0.8407029927301315,1.1587166257229957,3.1858420236698737,2.8662580547796024,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,marginal,none,16,0.16369893363007443,0.0016338409986929037,1.8278931643953367,0.9871901716652052,0.5880017928332366,0.8407029927301315,1.053234435429744,2.866908969976418,2.6380191724114903,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,absolute,none,16,0.49818470690335037,0.18364997699973173,2.421919009464722,1.5812160167345903,0.36966273836233143,0.8633467162545958,1.5934680262659167,4.92078478173058,3.636305280323022,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,thompson,none,16,0.29420297344580293,0.0015975660026015248,1.9485933791146803,1.1078903863845488,0.5625,0.8407029927301315,0.9470571522825082,2.578111495150432,2.71307756384158,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,finite-minimax,none,16,0.5667743439515941,0.007063968001602916,2.7963018856341475,1.955598892904016,0.23365399350960653,1.0309571887233855,1.1963947255260805,3.3081685414836346,3.0957105365538906,,,
+0,14c7df7d1e83fa4e24ed,0.1,A549,remedi,none,16,0.6045098119714001,0.03465827499894658,2.5979788808305795,1.757275888100448,0.3225126609192627,0.8975750888756471,1.3371756917275883,3.8082725665869046,3.2998691142107868,optimal,0.0,0.35300882077658063
+1,8dbce83dd566935f9221,0.1,A549,uniform,none,16,0.2454635219465067,0.0016169150003406685,0.8038805115379104,0.5836606974391038,0.625,0.2916292068995025,2.772588722239781,15.999999999999998,5.895395390689373,,,
+1,8dbce83dd566935f9221,0.1,A549,mean,none,16,0.45695280459714516,0.0011526560010679532,1.226369022736875,1.0061492086380683,0.0,1.2263690227368753,-0.0,1.0,1.0,,,
+1,8dbce83dd566935f9221,0.1,A549,gibbs,none,16,0.4527858658385084,0.001237653999851318,1.2236000063940236,1.003380192295217,0.005369706395543143,1.2032869263272288,0.06885305770534886,1.0712787814198808,1.055341303462745,,,
+1,8dbce83dd566935f9221,0.1,A549,marginal,none,16,0.4504446850282667,0.0012615400009963196,1.2215052827582675,1.0012854686594608,0.009707488896083524,1.1845774412880419,0.09694461284039614,1.1017993463947136,1.0698134153718617,,,
+1,8dbce83dd566935f9221,0.1,A549,absolute,none,16,0.3746279255037778,0.007307625000976259,1.1750526655715465,0.9548328514727399,0.0967738181704323,0.8941772572348975,0.7609119955383721,2.140227208380616,1.9309946226348478,,,
+1,8dbce83dd566935f9221,0.1,A549,thompson,none,16,0.37127415373518047,0.0010932459990726784,1.1609862730748464,0.9407664589760397,0.0625,0.7010067188164794,0.6019239717322894,1.8256278800579229,2.212463895082692,,,
+1,8dbce83dd566935f9221,0.1,A549,finite-minimax,none,16,0.31048519390331225,0.005529177000425989,1.0349158255611004,0.8146960114622938,0.2088692283666835,0.34848480867699905,1.1702291210647875,3.222730949490644,2.9839951599355548,,,
+1,8dbce83dd566935f9221,0.1,A549,remedi,none,16,0.3436655481273693,0.024280376001115656,1.2013626283104024,0.9811428142115958,0.04227613675176224,1.0222313692600518,0.7258047479422886,2.06639335682975,1.9927456896160085,optimal,1.1605106042811286e-09,0.18988115965667157
+2,14c7df7d1e83fa4e24ed,1.0,A549,uniform,none,16,0.2648465630398466,0.0015727989994047675,1.149269713810575,0.6108295143745597,0.5625,0.6173246359355173,2.772588722239781,15.999999999999998,5.895395390689373,,,
+2,14c7df7d1e83fa4e24ed,1.0,A549,mean,none,16,0.0,0.001173176999145653,1.2510109112224548,0.7125707117864395,0.0,1.2510109112224548,-0.0,1.0,1.0,,,
+2,14c7df7d1e83fa4e24ed,1.0,A549,gibbs,none,16,0.17506984351246402,0.0014063979979255237,1.1438039996426688,0.6053638002066535,0.4370569180616298,0.7534218278405237,1.7266779077545475,5.621946227346267,3.6376181528455467,,,
+2,14c7df7d1e83fa4e24ed,1.0,A549,marginal,none,16,0.22342806963786846,0.0014775050003663637,1.137020099408348,0.5985798999723326,0.5336989151840801,0.776155620908845,1.8920300921714974,6.632820258934401,3.9107984009498065,,,
+2,14c7df7d1e83fa4e24ed,1.0,A549,absolute,none,16,0.2411704018692673,0.00853000699862605,1.1456644231716515,0.6072242237356362,0.5318091290640261,0.6954067674309812,2.1368988018754083,8.4731200180064,4.3173870020044385,,,
+2,14c7df7d1e83fa4e24ed,1.0,A549,thompson,none,16,0.2131776081689259,0.001410534001479391,1.052387535391039,0.5139473359550237,0.5625,0.691584846136955,1.511086176073611,4.531650292231521,4.130858101420017,,,
+2,14c7df7d1e83fa4e24ed,1.0,A549,finite-minimax,none,16,0.18761398298570467,0.0065819929986901116,1.1572497852094266,0.6188095857734113,0.3998546403727868,0.7666257754979958,1.7230832374250509,5.601773462989668,3.3098655100339283,,,
+2,14c7df7d1e83fa4e24ed,1.0,A549,remedi,none,16,0.16269190929589744,0.03506866299721878,1.0718627449366478,0.5334225455006325,0.5292052355513774,0.6247198431606074,1.808009872567144,6.098298955686534,3.5796700161733987,optimal,0.0,0.3015974926731174
+3,97df6318d37dbefefcde,0.1,A549,uniform,none,16,0.3166748893064085,0.0017466999997850507,1.6788635528842226,0.9939433520463352,0.1875,0.7684386466775295,2.772588722239781,15.999999999999998,5.895395390689373,,,
+3,97df6318d37dbefefcde,0.1,A549,mean,none,16,2.0323988070382577,0.0011869570007547736,3.2819705191387216,2.597050318300834,0.0,3.281970519138722,-0.0,1.0,1.0,,,
+3,97df6318d37dbefefcde,0.1,A549,gibbs,none,16,1.9049490287839812,0.0012808690007659607,3.171930412774961,2.487010211937074,0.0024842742732633614,2.338255003957254,0.3200948817026613,1.3772584347620376,1.2132076402961824,,,
+3,97df6318d37dbefefcde,0.1,A549,marginal,none,16,1.8645034168167547,0.0012675190009758808,3.157331453186382,2.4724112523484942,0.004972355282530586,2.298914304202624,0.37412574673536625,1.453719939500578,1.3278928649068344,,,
+3,97df6318d37dbefefcde,0.1,A549,absolute,none,16,1.2308205638971343,0.007664844000828452,2.6483804809172216,1.9634602800793342,0.02812314176056409,1.21416553320923,1.3608890443887045,3.8996587308193416,2.596523998664006,,,
+3,97df6318d37dbefefcde,0.1,A549,thompson,none,16,1.8241956042266425,0.0011864170010085218,3.145781156588227,2.4608609557503396,0.0,2.313898842079318,0.46341355882643,1.589490554357095,1.9510790489712235,,,
+3,97df6318d37dbefefcde,0.1,A549,finite-minimax,none,16,1.7259342565487175,0.006497856000351021,3.1925284697490137,2.5076082689111265,0.00018405949512587726,2.781382202303947,0.488699396972656,1.6301946045269182,1.8368115771998015,,,
+3,97df6318d37dbefefcde,0.1,A549,remedi,none,16,1.2941521372160594,0.0341392629998154,2.724803537761057,2.03988333692317,0.00626467322507552,1.459156075361722,1.182897060305342,3.2638159913231215,2.664516951814567,optimal,0.0,0.25189030894311104
diff --git a/examples/results/transfer/summary.csv b/examples/results/transfer/summary.csv
new file mode 100644
index 0000000000000000000000000000000000000000..acb10ffdc0dc8e228bc2bc006dc9846a1e317827
--- /dev/null
+++ b/examples/results/transfer/summary.csv
@@ -0,0 +1,57 @@
+method,metric,mean,ci_low,ci_high,molecules
+absolute,endpoint_loss,1.869074954268985,1.1750526655715465,2.6483804809172216,3
+finite-minimax,endpoint_loss,2.068073376910634,1.0349158255611004,3.1925284697490137,3
+gibbs,endpoint_loss,2.0247813451156365,1.2236000063940236,3.171930412774961,3
+marginal,endpoint_loss,1.9537644559488305,1.2215052827582675,3.157331453186382,3
+mean,endpoint_loss,2.2375538011702063,1.226369022736875,3.281970519138722,3
+remedi,endpoint_loss,1.9203623263183578,1.2013626283104024,2.724803537761057,3
+thompson,endpoint_loss,1.9357526289719778,1.1609862730748464,3.1457811565882268,3
+uniform,endpoint_loss,1.498489706293402,0.8038805115379105,2.0127250544580733,3
+absolute,regret,1.337504417262396,0.9548328514727399,1.9634602800793342,3
+finite-minimax,regret,1.5365028399040448,0.8146960114622938,2.5076082689111265,3
+gibbs,regret,1.4932108081090476,0.9892420200948516,2.487010211937074,3
+marginal,regret,1.4221939189422415,0.7928850358187689,2.4724112523484942,3
+mean,regret,1.705983264163617,1.0061492086380683,2.597050318300834,3
+remedi,regret,1.3887917893117685,0.9811428142115958,2.03988333692317,3
+thompson,regret,1.4041820919653887,0.8109188611697863,2.4608609557503396,3
+uniform,regret,0.9669191692868132,0.5836606974391038,1.3231534583750002,3
+absolute,success,0.19187763121472504,0.02812314176056409,0.4507359337131788,3
+finite-minimax,success,0.175269201601002,0.00018405949512587726,0.31675431694119666,3
+gibbs,success,0.1474983510165482,0.0024842742732633614,0.4346410723808381,3
+marginal,success,0.19184339939575748,0.004972355282530586,0.5608503540086583,3
+mean,success,0.0,0.0,0.0,3
+remedi,success,0.15813325273738593,0.00626467322507552,0.42585894823532006,3
+thompson,success,0.20833333333333334,0.0,0.5625,3
+uniform,success,0.375,0.1875,0.625,3
+absolute,best_of_10,0.962573177428972,0.7793767418427885,1.21416553320923,3
+finite-minimax,best_of_10,1.342886164363879,0.348484808676999,2.781382202303947,3
+gibbs,best_of_10,1.4462014468566036,0.7970624102853275,2.338255003957254,3
+marginal,best_of_10,1.4306403507700514,0.8084293068194883,2.298914304202624,3
+mean,best_of_10,2.2375538011702063,1.2263690227368753,3.281970519138722,3
+remedi,best_of_10,1.0808449702133005,0.7611474660181271,1.459156075361722,3
+thompson,best_of_10,1.2603498267764468,0.7010067188164794,2.313898842079318,3
+uniform,best_of_10,0.6938889354166369,0.2916292068995025,1.0215989526728788,3
+absolute,dose_regret,0.6583753479290736,0.36967755438630884,1.2308205638971343,3
+finite-minimax,dose_regret,0.804537871306893,0.31048519390331225,1.7259342565487177,3
+gibbs,dose_regret,0.8568293129838808,0.21275304432915268,1.904949028783981,3
+marginal,dose_regret,0.8361705344929976,0.19356350163397149,1.8645034168167547,3
+mean,dose_regret,0.8473030261688117,0.05255746687103224,2.0323988070382577,3
+remedi,dose_regret,0.6738061819923592,0.3436655481273693,1.2941521372160594,3
+thompson,dose_regret,0.8163866829230625,0.2536902908073644,1.8241956042266425,3
+uniform,dose_regret,0.33890640967473207,0.2454635219465067,0.4545808177712809,3
+absolute,entropy,1.3289948179992463,0.7609119955383722,1.8651834140706625,3
+finite-minimax,entropy,1.039555833171003,0.488699396972656,1.4597389814755657,3
+gibbs,entropy,0.6105484020489272,0.06885305770534886,1.4426972667387716,3
+marginal,entropy,0.6479008744587944,0.09694461284039614,1.4726322638006206,3
+mean,entropy,0.0,0.0,0.0,3
+remedi,entropy,1.1604315301316657,0.7258047479422886,1.5725927821473664,3
+thompson,entropy,0.7648030649122597,0.4634135588264299,1.2290716641780595,3
+uniform,entropy,2.772588722239781,2.772588722239781,2.772588722239781,3
+absolute,unique_molecules_10,2.8347882541541947,1.9309946226348478,3.97684614116373,3
+finite-minimax,unique_molecules_10,2.674531586809755,1.8368115771998015,3.2027880232939094,3
+gibbs,unique_molecules_10,1.8401623491905006,1.055341303462745,3.251938103812575,3
+marginal,unique_molecules_10,1.8907050223197814,1.0698134153718617,3.2744087866806484,3
+mean,unique_molecules_10,1.0,1.0,1.0,3
+remedi,unique_molecules_10,2.6990107355408894,1.9927456896160087,3.4397695651920928,3
+thompson,unique_molecules_10,2.528503592228238,1.9510790489712235,3.4219678326307985,3
+uniform,unique_molecules_10,5.895395390689373,5.895395390689373,5.895395390689373,3
diff --git a/examples/sciplex/sciplex_prototype.h5ad b/examples/sciplex/sciplex_prototype.h5ad
new file mode 100644
index 0000000000000000000000000000000000000000..5f126dceb456784cb45c7f76908ee4ff1e8db51a
--- /dev/null
+++ b/examples/sciplex/sciplex_prototype.h5ad
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:0c89358e2ae92614a170c30adf38758f39f58c82b2b1160023a14d500bc3c131
+size 3593503
diff --git a/examples/sciplex/source_manifest.json b/examples/sciplex/source_manifest.json
new file mode 100644
index 0000000000000000000000000000000000000000..9390bd745f24d01ebe05347805224d6cd7ac4446
--- /dev/null
+++ b/examples/sciplex/source_manifest.json
@@ -0,0 +1,21 @@
+{
+  "source_url": "https://zenodo.org/records/13350497/files/SrivatsanTrapnell2020_sciplex3.h5ad?download=1",
+  "source_record": "scPerturb corrected v1.4 / Zenodo 13350497",
+  "source_md5": "c9e70629505d98c7ca1a837f62b14e89",
+  "full_source_checksum_recomputed": false,
+  "smiles_source": "https://huggingface.co/datasets/tahoebio/Tahoe-100M/resolve/main/metadata/drug_metadata.parquet",
+  "molecules": 9,
+  "cells": 120,
+  "gene_columns": 110983,
+  "nonzero_counts": 118728,
+  "cell_line": "A549",
+  "exposure_hours": 24,
+  "dose_um": [
+    0.1,
+    1.0
+  ],
+  "seed": 0,
+  "prototype_only": true,
+  "elapsed_seconds": 192.63887066200186,
+  "selected_rows_sha256": "1f01fb9dbb95b637bf2e55fefebd48d4f8f0fadce68c3d4a9f16b0c674313b94"
+}
diff --git a/examples/sciplex/structures.csv b/examples/sciplex/structures.csv
new file mode 100644
index 0000000000000000000000000000000000000000..e4a19e3e03b8e15e261081a8d12d0e54d43d740c
--- /dev/null
+++ b/examples/sciplex/structures.csv
@@ -0,0 +1,10 @@
+drug,smiles,tahoe_name,matching_rule
+Altretamine,CN(C)C1=NC(=NC(=N1)N(C)C)N(C)C,Altretamine,exact case-folded name after outer whitespace removal
+Busulfan ,CS(=O)(=O)OCCCCOS(=O)(=O)C,Busulfan,exact case-folded name after outer whitespace removal
+Celecoxib,CC1=CC=C(C=C1)C2=CC(=NN2C3=CC=C(C=C3)S(=O)(=O)N)C(F)(F)F,Celecoxib,exact case-folded name after outer whitespace removal
+Decitabine,C1C(C(OC1N2C=NC(=NC2=O)N)CO)O,Decitabine,exact case-folded name after outer whitespace removal
+Fulvestrant,CC12CCC3C(C1CCC2O)C(CC4=C3C=CC(=C4)O)CCCCCCCCCS(=O)CCCC(C(F)(F)F)(F)F,Fulvestrant,exact case-folded name after outer whitespace removal
+Raltitrexed,CC1=NC2=C(C=C(C=C2)CN(C)C3=CC=C(S3)C(=O)NC(CCC(=O)O)C(=O)O)C(=O)N1,Raltitrexed,exact case-folded name after outer whitespace removal
+Resveratrol,C1=CC(=CC=C1C=CC2=CC(=CC(=C2)O)O)O,Resveratrol,exact case-folded name after outer whitespace removal
+TAK-901,CCS(=O)(=O)C1=CC=CC(=C1)C2=CC(=C(C3=C2C4=C(N3)N=CC(=C4)C)C)C(=O)NC5CCN(CC5)C,TAK-901,exact case-folded name after outer whitespace removal
+Zileuton,CC(C1=CC2=CC=CC=C2S1)N(C(=O)N)O,Zileuton,exact case-folded name after outer whitespace removal
diff --git a/examples/tahoe/source.json b/examples/tahoe/source.json
new file mode 100644
index 0000000000000000000000000000000000000000..f80e08a0509bdaaf4a5868bea5f89fe16d0b4a62
--- /dev/null
+++ b/examples/tahoe/source.json
@@ -0,0 +1,12 @@
+{
+  "dataset": "tahoebio/Tahoe-100M",
+  "shards": [
+    "data/train-00600-of-03388.parquet",
+    "data/train-00601-of-03388.parquet",
+    "data/train-02600-of-03388.parquet"
+  ],
+  "cells": 568,
+  "conditions": 144,
+  "reservoir_cap": 4,
+  "seed": 0
+}
diff --git a/examples/tahoe/structures.csv b/examples/tahoe/structures.csv
new file mode 100644
index 0000000000000000000000000000000000000000..6e83e823ccfc532c699ca4b1027f50aca7176299
--- /dev/null
+++ b/examples/tahoe/structures.csv
@@ -0,0 +1,95 @@
+drug,smiles
+(R)-Verapamil (hydrochloride),CC(C)C(CCCN(C)CCC1=CC(=C(C=C1)OC)OC)(C#N)C2=CC(=C(C=C2)OC)OC.Cl
+18β-Glycyrrhetinic acid,C[C@]12CC[C@](C[C@H]1C3=CC(=O)[C@@H]4[C@]5(CC[C@@H](C([C@@H]5CC[C@]4([C@@]3(CC2)C)C)(C)C)O)C)(C)C(=O)O
+8-Hydroxyquinoline,C1=CC2=C(C(=C1)O)N=CC=C2
+Acetohexamide,CC(=O)C1=CC=C(C=C1)S(=O)(=O)NC(=O)NC2CCCCC2
+Adagrasib,CN1CCCC1COC2=NC3=C(CCN(C3)C4=CC=CC5=C4C(=CC=C5)Cl)C(=N2)N6CCN(C(C6)CC#N)C(=O)C(=C)F
+Adenine,C1=NC2=NC=NC(=C2N1)N
+Adenosine,C1=NC(=C2C(=N1)N(C=N2)C3C(C(C(O3)CO)O)O)N
+Aliskiren,CC(C)C(CC1=CC(=C(C=C1)OC)OCCCOC)CC(C(CC(C(C)C)C(=O)NCC(C)(C)C(=O)N)O)N
+Allantoin,C1(C(=O)NC(=O)N1)NC(=O)N
+Almonertinib (hydrochloride),CN(C)CCN(C)C1=CC(=C(C=C1NC(=O)C=C)NC2=NC=CC(=N2)C3=CN(C4=CC=CC=C43)C5CC5)OC.Cl
+Amsacrine,COC1=C(C=CC(=C1)NS(=O)(=O)C)NC2=C3C=CC=CC3=NC4=CC=CC=C42
+Anethole trithione,COC1=CC=C(C=C1)C2=CC(=S)SS2
+Artesunate,CC1CCC2C(C(OC3C24C1CCC(O3)(OO4)C)OC(=O)CCC(=O)O)C
+Azithromycin (hydrate),CCC1C(C(C(N(CC(CC(C(C(C(C(C(=O)O1)C)OC2CC(C(C(O2)C)O)(C)OC)C)OC3C(C(CC(O3)C)N(C)C)O)(C)O)C)C)C)O)(C)O.O.O
+Baicalin,C1=CC=C(C=C1)C2=CC(=O)C3=C(C(=C(C=C3O2)OC4C(C(C(C(O4)C(=O)O)O)O)O)O)O
+Balsalazide (sodium hydrate),C1=CC(=CC=C1C(=O)NCCC(=O)[O-])N=NC2=CC(=C(C=C2)[O-])C(=O)O.O.O.[Na+].[Na+]
+Berbamine,CN1CCC2=CC(=C3C=C2C1CC4=CC=C(C=C4)OC5=C(C=CC(=C5)CC6C7=C(O3)C(=C(C=C7CCN6C)OC)OC)O)OC
+Berberine (chloride hydrate),COC1=C(C2=C[N+]3=C(C=C2C=C1)C4=CC5=C(C=C4CC3)OCO5)OC.O.[Cl-]
+Bergenin,COC1=C(C=C2C(=C1O)C3C(C(C(C(O3)CO)O)O)OC2=O)O
+Bimatoprost,CCNC(=O)CCCC=CCC1C(CC(C1C=CC(CCC2=CC=CC=C2)O)O)O
+Busulfan,CS(=O)(=O)OCCCCOS(=O)(=O)C
+Cepharanthine,CN1CCC2=CC3=C(C4=C2C1CC5=CC=C(C=C5)OC6=C(C=CC(=C6)CC7C8=CC(=C(C=C8CCN7C)OC)O4)OC)OCO3
+Chlorhexidine (diacetate),CC(=O)O.CC(=O)O.C1=CC(=CC=C1NC(=NC(=NCCCCCCN=C(N)N=C(N)NC2=CC=C(C=C2)Cl)N)N)Cl
+Ciclopirox,CC1=CC(=O)N(C(=C1)C2CCCCC2)O
+Clonidine (hydrochloride),C1CN=C(N1)NC2=C(C=CC=C2Cl)Cl.Cl
+Clotrimazole,C1=CC=C(C=C1)C(C2=CC=CC=C2)(C3=CC=CC=C3Cl)N4C=CN=C4
+Cloxacillin (sodium),CC1=C(C(=NO1)C2=CC=CC=C2Cl)C(=O)NC3C4N(C3=O)C(C(S4)(C)C)C(=O)[O-].[Na+]
+Cysteamine (hydrochloride),C(CS)N.Cl
+Darinaparsin,C[As](C)SCC(C(=O)NCC(=O)O)NC(=O)CCC(C(=O)O)N
+Demeclocycline,CN(C)C1C2CC3C(C4=C(C=CC(=C4C(=C3C(=O)C2(C(=C(C1=O)C(=O)N)O)O)O)O)Cl)O
+Diammonium Glycyrrhizinate,CC1(C2CCC3(C(C2(CCC1OC4C(C(C(C(O4)C(=O)O)O)O)OC5C(C(C(C(O5)C(=O)O)O)O)O)C)C(=O)C=C6C3(CCC7(C6CC(CC7)(C)C(=O)O)C)C)C)C.N.N
+Docetaxel,CC1=C2C(C(=O)C3(C(CC4C(C3C(C(C2(C)C)(CC1OC(=O)C(C(C5=CC=CC=C5)NC(=O)OC(C)(C)C)O)O)OC(=O)C6=CC=CC=C6)(CO4)OC(=O)C)O)C)O
+Elagolix sodium,CC1=C(C(=O)N(C(=O)N1CC2=C(C=CC=C2F)C(F)(F)F)CC(C3=CC=CC=C3)NCCCC(=O)[O-])C4=C(C(=CC=C4)OC)F.[Na+]
+Fluvoxamine (maleate),COCCCCC(=NOCCN)C1=CC=C(C=C1)C(F)(F)F.C(=CC(=O)O)C(=O)O
+Folic acid,C1=CC(=CC=C1C(=O)NC(CCC(=O)O)C(=O)O)NCC2=CN=C3C(=N2)C(=O)NC(=N3)N
+Furosemide,C1=COC(=C1)CNC2=CC(=C(C=C2C(=O)O)S(=O)(=O)N)Cl
+Fusidic acid,CC1C2CCC3(C(C2(CCC1O)C)C(CC4C3(CC(C4=C(CCC=C(C)C)C(=O)O)OC(=O)C)C)O)C
+Gallic acid,C1=C(C=C(C(=C1O)O)O)C(=O)O
+Gallic acid (hydrate),C1=C(C=C(C(=C1O)O)O)C(=O)O.O
+Gemfibrozil,CC1=CC(=C(C=C1)C)OCCCC(C)(C)C(=O)O
+Goserelin (acetate),CC(C)CC(C(=O)NC(CCCN=C(N)N)C(=O)N1CCCC1C(=O)NNC(=O)N)NC(=O)C(COC(C)(C)C)NC(=O)C(CC2=CC=C(C=C2)O)NC(=O)C(CO)NC(=O)C(CC3=CNC4=CC=CC=C43)NC(=O)C(CC5=CN=CN5)NC(=O)C6CCC(=O)N6.CC(=O)O
+Hexylresorcinol,CCCCCCC1=C(C=C(C=C1)O)O
+Homoharringtonine,CC(C)(CCCC(CC(=O)OC)(C(=O)OC1C2C3=CC4=C(C=C3CCN5C2(CCC5)C=C1OC)OCO4)O)O
+Idarubicin (hydrochloride),CC1C(C(CC(O1)OC2CC(CC3=C2C(=C4C(=C3O)C(=O)C5=CC=CC=C5C4=O)O)(C(=O)C)O)N)O.Cl
+Isocorydine,CN1CCC2=CC(=C(C3=C2C1CC4=C3C(=C(C=C4)OC)O)OC)OC
+L-Thyroxine (sodium salt pentahydrate),C1=C(C=C(C(=C1I)OC2=CC(=C(C(=C2)I)O)I)I)CC(C(=O)[O-])N.O.O.O.O.O.[Na+]
+Larotrectinib,C1CC(N(C1)C2=NC3=C(C=NN3C=C2)NC(=O)N4CCC(C4)O)C5=C(C=CC(=C5)F)F
+Lenalidomide (hemihydrate),C1CC(=O)NC(=O)C1N2CC3=C(C2=O)C=CC=C3N.C1CC(=O)NC(=O)C1N2CC3=C(C2=O)C=CC=C3N.O
+Lidocaine (hydrochloride),CCN(CC)CC(=O)NC1=C(C=CC=C1C)C.Cl
+Lipoic acid,C1CSSC1CCCCC(=O)O
+Lumateperone (tosylate),CC1=CC=C(C=C1)S(=O)(=O)O.CN1CCN2C3CCN(CC3C4=C2C1=CC=C4)CCCC(=O)C5=CC=C(C=C5)F
+Medroxyprogesterone acetate,CC1CC2C(CCC3(C2CCC3(C(=O)C)OC(=O)C)C)C4(C1=CC(=O)CC4)C
+Megestrol,CC1=CC2C(CCC3(C2CCC3(C(=O)C)O)C)C4(C1=CC(=O)CC4)C
+Mitoxantrone (dihydrochloride),C1=CC(=C2C(=C1NCCNCCO)C(=O)C3=C(C=CC(=C3C2=O)O)O)NCCNCCO.Cl.Cl
+Neratinib (maleate),CCOC1=C(C=C2C(=C1)N=CC(=C2NC3=CC(=C(C=C3)OCC4=CC=CC=N4)Cl)C#N)NC(=O)C=CCN(C)C.C(=CC(=O)O)C(=O)O
+Norepinephrine (hydrochloride),C1=CC(=C(C=C1C(CN)O)O)O.Cl
+Olanzapine,CC1=CC2=C(S1)NC3=CC=CC=C3N=C2N4CCN(CC4)C
+Palmatine (chloride),COC1=C(C2=C[N+]3=C(C=C2C=C1)C4=CC(=C(C=C4CC3)OC)OC)OC.[Cl-]
+Pasireotide (acetate),CC(=O)O.C1C(CN2C1C(=O)NC(C(=O)NC(C(=O)NC(C(=O)NC(C(=O)NC(C2=O)CC3=CC=CC=C3)CC4=CC=C(C=C4)OCC5=CC=CC=C5)CCCCN)CC6=CNC7=CC=CC=C76)C8=CC=CC=C8)OC(=O)NCCN
+Pentagastrin,CC(C)(C)OC(=O)NCCC(=O)NC(CC1=CNC2=CC=CC=C21)C(=O)NC(CCSC)C(=O)NC(CC(=O)O)C(=O)NC(CC3=CC=CC=C3)C(=O)N
+Pentamidine (isethionate),C1=CC(=CC=C1C(=N)N)OCCCCCOC2=CC=C(C=C2)C(=N)N.C(CS(=O)(=O)O)O.C(CS(=O)(=O)O)O
+Pyridoxine,CC1=NC=C(C(=C1O)CO)CO
+Pyridoxine (hydrochloride),CC1=NC=C(C(=C1O)CO)CO.Cl
+Quinestrol,CC12CCC3C(C1CCC2(C#C)O)CCC4=C3C=CC(=C4)OC5CCCC5
+Radotinib,CC1=C(C=C(C=C1)C(=O)NC2=CC(=CC(=C2)C(F)(F)F)N3C=C(N=C3)C)NC4=NC=CC(=N4)C5=NC=CN=C5
+Raltitrexed,CC1=NC2=C(C=C(C=C2)CN(C)C3=CC=C(S3)C(=O)NC(CCC(=O)O)C(=O)O)C(=O)N1
+Rapamycin,CC1CCC2CC(C(=CC=CC=CC(CC(C(=O)C(C(C(=CC(C(=O)CC(OC(=O)C3CCCCN3C(=O)C(=O)C1(O2)O)C(C)CC4CCC(C(C4)OC)O)C)C)O)OC)C)C)C)OC
+Relugolix,CN(C)CC1=C(SC2=C1C(=O)N(C(=O)N2CC3=C(C=CC=C3F)F)C4=NN=C(C=C4)OC)C5=CC=C(C=C5)NC(=O)NOC
+Resveratrol,C1=CC(=CC=C1C=CC2=CC(=CC(=C2)O)O)O
+Retinoic acid,CC1=C(C(CCC1)(C)C)C=CC(=CC=CC(=CC(=O)O)C)C
+Riluzole hydrochloride,C1=CC2=C(C=C1OC(F)(F)F)SC(=N2)N.Cl
+Rimonabant,CC1=C(N(N=C1C(=O)NN2CCCCC2)C3=C(C=C(C=C3)Cl)Cl)C4=CC=C(C=C4)Cl
+Ritonavir,CC(C)C1=NC(=CS1)CN(C)C(=O)NC(C(C)C)C(=O)NC(CC2=CC=CC=C2)CC(C(CC3=CC=CC=C3)NC(=O)OCC4=CN=CS4)O
+Ropivacaine (hydrochloride monohydrate),CCCN1CCCCC1C(=O)NC2=C(C=CC=C2C)C.O.Cl
+S-Adenosyl-L-methionine (disulfate tosylate),CC1=CC=C(C=C1)S(=O)(=O)O.C[S+](CCC(C(=O)O)N)CC1C(C(C(O1)N2C=NC3=C(N=CN=C32)N)O)O.OS(=O)(=O)O.OS(=O)(=O)[O-]
+Salicylic acid,C1=CC=C(C(=C1)C(=O)O)O
+Saquinavir,CC(C)(C)NC(=O)C1CC2CCCCC2CN1CC(C(CC3=CC=CC=C3)NC(=O)C(CC(=O)N)NC(=O)C4=NC5=CC=CC=C5C=C4)O
+Selinexor,C1=CN=C(C=N1)NNC(=O)C=CN2C=NC(=N2)C3=CC(=CC(=C3)C(F)(F)F)C(F)(F)F
+Sinomenine,CN1CCC23CC(=O)C(=CC2C1CC4=C3C(=C(C=C4)OC)O)OC
+Sodium Salicylate,C1=CC=C(C(=C1)C(=O)[O-])O.[Na+]
+Sulfisoxazole,CC1=C(ON=C1C)NS(=O)(=O)C2=CC=C(C=C2)N
+Temsirolimus,CC1CCC2CC(C(=CC=CC=CC(CC(C(=O)C(C(C(=CC(C(=O)CC(OC(=O)C3CCCCN3C(=O)C(=O)C1(O2)O)C(C)CC4CCC(C(C4)OC)OC(=O)C(C)(CO)CO)C)C)O)OC)C)C)C)OC
+Terfenadine,CC(C)(C)C1=CC=C(C=C1)C(CCCN2CCC(CC2)C(C3=CC=CC=C3)(C4=CC=CC=C4)O)O
+Thymol,CC1=CC(=C(C=C1)C(C)C)O
+Tofacitinib,CC1CCN(CC1N(C)C2=NC=NC3=C2C=CN3)C(=O)CC#N
+Tofacitinib (citrate),CC1CCN(CC1N(C)C2=NC=NC3=C2C=CN3)C(=O)CC#N.C(C(=O)O)C(CC(=O)O)(C(=O)O)O
+Tolcapone,CC1=CC=C(C=C1)C(=O)C2=CC(=C(C(=C2)O)O)[N+](=O)[O-]
+Triamcinolone,CC12CC(C3(C(C1CC(C2(C(=O)CO)O)O)CCC4=CC(=O)C=CC43C)F)O
+Triclosan,C1=CC(=C(C=C1Cl)O)OC2=C(C=C(C=C2)Cl)Cl
+Trifluridine,C1C(C(OC1N2C=C(C(=O)NC2=O)C(F)(F)F)CO)O
+Trimetrexate,CC1=C(C=CC2=C1C(=NC(=N2)N)N)CNC3=CC(=C(C(=C3)OC)OC)OC
+Tucidinostat,C1=CC(=CN=C1)C=CC(=O)NCC2=CC=C(C=C2)C(=O)NC3=C(C=C(C=C3)F)N
+Vilanterol,C1=CC(=C(C(=C1)Cl)COCCOCCCCCCNCC(C2=CC(=C(C=C2)O)CO)O)Cl
+Vinblastine (sulfate),CCC1(CC2CC(C3=C(CCN(C2)C1)C4=CC=CC=C4N3)(C5=C(C=C6C(=C5)C78CCN9C7C(C=CC9)(C(C(C8N6C)(C(=O)OC)O)OC(=O)C)CC)OC)C(=O)OC)O.OS(=O)(=O)O
diff --git a/examples/tahoe/tahoe.h5ad b/examples/tahoe/tahoe.h5ad
new file mode 100644
index 0000000000000000000000000000000000000000..dc927c4a9984ef4e9700a22d6a2251f8a5c188ac
--- /dev/null
+++ b/examples/tahoe/tahoe.h5ad
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:f1482038edcf19663c8e48d84e743328999444c02b60af264abf39407f07ae63
+size 6067969
diff --git a/pyproject.toml b/pyproject.toml
new file mode 100644
index 0000000000000000000000000000000000000000..7e8fdade837b6c8a5959842788fe708f43982573
--- /dev/null
+++ b/pyproject.toml
@@ -0,0 +1,33 @@
+[build-system]
+requires = ["setuptools>=68", "wheel"]
+build-backend = "setuptools.build_meta"
+
+[project]
+name = "remedi-cell"
+version = "0.1.0"
+description = "Replicate-enabled minimax endpoint-directed inference for molecule-dose selection"
+readme = "README.md"
+requires-python = ">=3.11"
+dependencies = ["numpy>=1.24", "scipy>=1.10", "pandas>=2.0", "scikit-learn>=1.7", "cvxpy>=1.5", "clarabel>=0.9", "joblib>=1.3", "matplotlib>=3.7", "requests>=2.31", "PyYAML>=6", "rdkit>=2023.9"]
+
+[project.optional-dependencies]
+data = ["anndata>=0.10", "h5py>=3.9", "rdkit>=2023.9", "huggingface_hub>=0.25", "pyarrow>=14", "fsspec>=2024.2", "aiohttp>=3.9"]
+encoders = ["torch>=2.2", "transformers>=4.40", "rdkit>=2023.9"]
+dev = ["pytest>=8", "ruff>=0.5"]
+
+[project.scripts]
+remedi = "remedi.cli:main"
+
+[tool.setuptools.packages.find]
+where = ["src"]
+
+[tool.pytest.ini_options]
+testpaths = ["tests"]
+addopts = "-q"
+
+[tool.ruff]
+line-length = 100
+
+
+[tool.setuptools.package-data]
+remedi = ["fonts/*.ttf", "fonts/UFL.txt"]
diff --git a/requirements-tested.txt b/requirements-tested.txt
new file mode 100644
index 0000000000000000000000000000000000000000..9fe2529bfceaba8088f092ef3ebc92472753cac2
--- /dev/null
+++ b/requirements-tested.txt
@@ -0,0 +1,19 @@
+# Versions used for the Python 3.12 CPU integration checks.
+numpy==2.3.5
+scipy==1.18.1
+pandas==3.0.5
+scikit-learn==1.8.0
+cvxpy==1.9.2
+clarabel==0.11.1
+joblib==1.5.3
+matplotlib==3.10.8
+requests==2.34.2
+PyYAML==6.0.3
+rdkit==2026.3.6
+anndata==0.13.3.post0
+h5py==3.16.0
+huggingface-hub==1.31.0
+pyarrow==25.0.1
+fsspec==2026.7.0
+aiohttp==3.13.5
+pytest==9.1.1
diff --git a/requirements.txt b/requirements.txt
new file mode 100644
index 0000000000000000000000000000000000000000..a73fdfd8539275564040864af85df8c0c5333e7b
--- /dev/null
+++ b/requirements.txt
@@ -0,0 +1 @@
+-e .[data]
diff --git a/scripts/align_genes.py b/scripts/align_genes.py
new file mode 100644
index 0000000000000000000000000000000000000000..a882f6f7456b2caa64292f54c5e4def1c9373c35
--- /dev/null
+++ b/scripts/align_genes.py
@@ -0,0 +1,31 @@
+"""Align raw study cohorts to shared gene identifiers before fitting PCA."""
+import argparse
+from pathlib import Path
+import anndata as ad
+
+p = argparse.ArgumentParser()
+p.add_argument('--inputs', nargs='+', required=True)
+p.add_argument('--output', required=True)
+p.add_argument('--gene-columns', nargs='+', help='One var column per input. Use index for var_names.')
+p.add_argument('--drop-ambiguous', action='store_true', help='Drop all duplicated gene identifiers')
+a = p.parse_args()
+cohorts = [ad.read_h5ad(f) for f in a.inputs]
+if a.gene_columns and len(a.gene_columns) != len(cohorts): raise ValueError('Supply one gene-column choice per input')
+for i, cohort in enumerate(cohorts):
+    field = a.gene_columns[i] if a.gene_columns else 'index'
+    if field != 'index':
+        if cohort.var[field].isna().any():
+            cohort = cohort[:, cohort.var[field].notna()].copy()
+        cohort.var_names = cohort.var[field].astype(str)
+    if not cohort.var_names.is_unique:
+        if not a.drop_ambiguous: raise ValueError('Duplicate gene identifiers. Use --drop-ambiguous to discard these genes.')
+        cohort = cohort[:, ~cohort.var_names.duplicated(keep=False)].copy()
+    cohorts[i] = cohort
+shared = sorted(set.intersection(*(set(c.var_names) for c in cohorts)))
+if len(shared) < 100:
+    raise ValueError('Fewer than 100 shared genes. Check identifier conventions.')
+out = Path(a.output); out.mkdir(parents=True, exist_ok=True)
+for i, (path, cohort) in enumerate(zip(a.inputs, cohorts)):
+    cohort[:, shared].copy().write_h5ad(out / f'{i}_{Path(path).stem}.h5ad', compression='gzip')
+(out / 'genes.txt').write_text('\n'.join(shared)+'\n')
+print(f'Aligned {len(cohorts)} cohorts on {len(shared)} genes')
diff --git a/scripts/download_data.py b/scripts/download_data.py
new file mode 100644
index 0000000000000000000000000000000000000000..a98d348115bf33bc50fc2fe781d2e1f2a601029d
--- /dev/null
+++ b/scripts/download_data.py
@@ -0,0 +1,39 @@
+"""Download released source files with retries and streamed checksum verification."""
+import argparse
+import hashlib
+from pathlib import Path
+import requests
+from requests.adapters import HTTPAdapter
+from urllib3.util.retry import Retry
+
+p=argparse.ArgumentParser()
+p.add_argument('dataset',choices=['tahoe','sciplex'])
+p.add_argument('--output',required=True)
+p.add_argument('--shards',type=int,default=2)
+p.add_argument('--indices',type=int,nargs='+',help='Explicit official Tahoe shard indices')
+a=p.parse_args()
+out=Path(a.output);out.mkdir(parents=True,exist_ok=True)
+s=requests.Session()
+s.mount('https://',HTTPAdapter(max_retries=Retry(total=5,backoff_factor=2,status_forcelist=[429,500,502,503,504],respect_retry_after_header=True)))
+base='https://huggingface.co/datasets/tahoebio/Tahoe-100M/resolve/main/'
+files=[(base+'metadata/'+n+'.parquet',n+'.parquet',None) for n in ['sample_metadata','gene_metadata','drug_metadata']]
+if a.dataset=='tahoe':
+    files += [(base+f'data/train-{i:05d}-of-03388.parquet',f'train-{i:05d}-of-03388.parquet',None) for i in (a.indices if a.indices is not None else range(a.shards))]
+else:
+    files += [('https://zenodo.org/records/13350497/files/SrivatsanTrapnell2020_sciplex3.h5ad?download=1','sciplex3.h5ad','c9e70629505d98c7ca1a837f62b14e89')]
+for url,name,expected in files:
+    path=out/name
+    if path.exists():
+        if expected and hashlib.file_digest(path.open('rb'),'md5').hexdigest()!=expected:
+            raise ValueError(f'Checksum mismatch for {name}')
+        print(f'Using {path}',flush=True);continue
+    temp=path.with_suffix(path.suffix+'.part')
+    with s.get(url,stream=True,timeout=(30,180)) as r:
+        r.raise_for_status()
+        with temp.open('wb') as f:
+            for block in r.iter_content(8*1024*1024):f.write(block)
+    if expected:
+        with temp.open('rb') as f:digest=hashlib.file_digest(f,'md5').hexdigest()
+        if digest!=expected:raise ValueError(f'Checksum mismatch for {name}')
+    temp.replace(path)
+    print(f'Downloaded {path} ({path.stat().st_size:,} bytes)',flush=True)
diff --git a/scripts/evaluate_forward.py b/scripts/evaluate_forward.py
new file mode 100644
index 0000000000000000000000000000000000000000..4c89a68aae0e696b00074d12022d474dfff81604
--- /dev/null
+++ b/scripts/evaluate_forward.py
@@ -0,0 +1,40 @@
+"""Compare response-head predictions with zero, training mean, and chemical neighbors."""
+import argparse
+from pathlib import Path
+import joblib
+import numpy as np
+import pandas as pd
+from remedi.io import ResponseData
+from remedi.models import predict,check_feature_space
+from remedi.chemistry import fingerprints
+
+p=argparse.ArgumentParser()
+for name in ['data','model','output']:p.add_argument('--'+name,required=True)
+p.add_argument('--split',default='test');p.add_argument('--embeddings')
+a=p.parse_args();d=ResponseData.load(a.data);m=joblib.load(Path(a.model)/'model.joblib');check_feature_space(m,d)
+tr=np.flatnonzero(d.obs.split.to_numpy()=='train');te=np.flatnonzero(d.obs.split.to_numpy()==a.split)
+if not len(tr) or not len(te):raise ValueError('Training and evaluation conditions are required')
+o=d.obs.iloc[te];y=d.response[te]
+preds={'zero':np.zeros_like(y),'training-mean':np.repeat(d.response[tr].mean(0)[None,:],len(te),axis=0),
+       'response-head':predict(m,o.smiles,o.dose_um,d.control[te],a.embeddings).mean(0)}
+fp=fingerprints(d.obs.smiles);nearest=[]
+for i in te:
+    candidates=tr[d.obs.iloc[tr].context.to_numpy()==d.obs.iloc[i].context]
+    if not len(candidates):candidates=tr
+    overlap=fp[candidates]@fp[i];union=fp[candidates].sum(1)+fp[i].sum()-overlap
+    similarity=np.divide(overlap,union,out=np.zeros_like(overlap),where=union>0)
+    # Select chemical neighbors first and the closest observed log dose within ties.
+    top=candidates[np.isclose(similarity,similarity.max())]
+    j=top[np.argmin(abs(np.log(d.obs.iloc[top].dose_um.to_numpy()/d.obs.iloc[i].dose_um)))]
+    nearest.append(d.response[j])
+preds['chemical-neighbor']=np.asarray(nearest)
+rows=[]
+for name,values in preds.items():
+    for k,i in enumerate(te):
+        truth=y[k];estimate=values[k]
+        corr=float(np.corrcoef(truth,estimate)[0,1]) if truth.std()>0 and estimate.std()>0 else np.nan
+        rows.append({'condition':int(i),'molecule_id':d.obs.iloc[i].molecule_id,'dose_um':d.obs.iloc[i].dose_um,
+          'method':name,'mse':float(np.mean((estimate-truth)**2)),'pearson':corr})
+out=Path(a.output);out.mkdir(parents=True,exist_ok=True)
+r=pd.DataFrame(rows);r.to_csv(out/'forward_conditions.csv',index=False)
+r.groupby(['method','molecule_id'])[['mse','pearson']].mean().groupby('method').mean().to_csv(out/'forward_summary.csv')
diff --git a/scripts/fetch_sciplex_subset.py b/scripts/fetch_sciplex_subset.py
new file mode 100644
index 0000000000000000000000000000000000000000..2f6250bef11e119eb7d50c72de1fe78f0be0061a
--- /dev/null
+++ b/scripts/fetch_sciplex_subset.py
@@ -0,0 +1,106 @@
+"""Fetch a bounded sci-Plex3 prototype by HTTP ranges, with exact-name Tahoe SMILES.
+
+The source is the corrected scPerturb object. The selected cohort is intentionally
+small and is an engineering prototype, not the external ICML benchmark.
+"""
+import argparse
+from remedi.download import RangeReader
+from pathlib import Path
+import json
+import hashlib
+import time
+import numpy as np
+import pandas as pd
+import requests
+import fsspec
+import h5py
+import anndata as ad
+import pyarrow.parquet as pq
+from scipy import sparse
+
+URL = "https://zenodo.org/records/13350497/files/SrivatsanTrapnell2020_sciplex3.h5ad?download=1"
+DRUG_URL = "https://huggingface.co/datasets/tahoebio/Tahoe-100M/resolve/main/metadata/drug_metadata.parquet"
+
+def read_column(node):
+    if isinstance(node, h5py.Group):
+        categories = read_column(node["categories"])
+        codes = node["codes"][:]
+        out = categories[np.maximum(codes, 0)].astype(object)
+        out[codes < 0] = None
+        return out
+    values = node[:]
+    return np.array([v.decode() if isinstance(v, bytes) else v for v in values]) if values.dtype.kind in "SO" else values
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument("--output", required=True)
+    parser.add_argument("--source", default=URL, help="Local corrected h5ad or source URL")
+    parser.add_argument("--cells-per-condition", type=int, default=8)
+    parser.add_argument("--cell-line", default="A549")
+    parser.add_argument("--max-molecules", type=int, default=40)
+    parser.add_argument("--seed", type=int, default=0)
+    parser.add_argument("--metadata-cache", help="Optional locally generated metadata cache directory")
+    args = parser.parse_args()
+    output = Path(args.output); output.mkdir(parents=True, exist_ok=True)
+    cache = Path(args.metadata_cache) if args.metadata_cache else output
+    drug_file = cache/"drug_metadata.parquet"
+    if not drug_file.exists():
+        response = requests.get(DRUG_URL, timeout=60); response.raise_for_status()
+        drug_file.parent.mkdir(parents=True, exist_ok=True); drug_file.write_bytes(response.content)
+    # Arrow conversion avoids dependence on a pandas-specific parquet extension registry.
+    drugs = pd.DataFrame(pq.read_table(drug_file).to_pylist())
+    lookup = {s.strip().casefold(): i for i, s in enumerate(drugs.drug) if isinstance(s, str)}
+    rng = np.random.default_rng(args.seed)
+    start = time.perf_counter()
+    handle = RangeReader(args.source, output/'range_cache', 2526631614) if args.source.startswith('https://') else open(args.source, 'rb')
+    with handle:
+        with h5py.File(handle, "r") as h:
+            names = ["cell_line","dose_value","perturbation","replicate","plate","well","time"]
+            obs = pd.DataFrame({name: read_column(h["obs"][name]) for name in names})
+            print(f"Loaded metadata for {len(obs)} cells", flush=True)
+            matches = sorted(s for s in obs.perturbation.dropna().unique() if s.strip().casefold() in lookup)[:args.max_molecules]
+            mappings = [{"drug": s, "smiles": drugs.iloc[lookup[s.strip().casefold()]].canonical_smiles,
+                         "tahoe_name": drugs.iloc[lookup[s.strip().casefold()]].drug,
+                         "matching_rule": "exact case-folded name after outer whitespace removal"} for s in matches]
+            pd.DataFrame(mappings).to_csv(output/"structures.csv", index=False)
+            eligible = obs[(obs.cell_line == args.cell_line) & (obs.time == 24.)]
+            treated = eligible[eligible.perturbation.isin(matches) & eligible.dose_value.isin([100.,1000.])]
+            selected = []
+            for _, group in treated.groupby(["perturbation","dose_value","replicate","plate"], sort=True):
+                if len(group) >= args.cells_per_condition:
+                    selected.extend(rng.choice(group.index, args.cells_per_condition, replace=False))
+            plates = set(obs.loc[selected, "plate"])
+            controls = eligible[(eligible.perturbation == "control") & eligible.plate.isin(plates)]
+            for _, group in controls.groupby(["replicate","plate"], sort=True):
+                selected.extend(rng.choice(group.index, min(args.cells_per_condition*2,len(group)), replace=False))
+            selected = np.sort(np.unique(selected))
+            x = h["X"]
+            if x.attrs.get("encoding-type") != "csr_matrix": raise ValueError("Expected source CSR matrix")
+            pointers = x["indptr"][:]
+            values, columns, pointer = [], [], [0]
+            for j, index in enumerate(selected):
+                begin, end = int(pointers[index]), int(pointers[index+1])
+                values.append(x["data"][begin:end]); columns.append(x["indices"][begin:end])
+                pointer.append(pointer[-1]+end-begin)
+                if j % 100 == 0: print(f"Read {j}/{len(selected)} selected cell rows", flush=True)
+            matrix = sparse.csr_matrix((np.concatenate(values),np.concatenate(columns),np.asarray(pointer)),
+                                       shape=(len(selected),int(x.attrs["shape"][1])))
+            var_index = h["var"].attrs.get("_index", "_index")
+            if isinstance(var_index, bytes): var_index=var_index.decode()
+            gene_names = read_column(h["var"][var_index]).astype(str)
+    subset_obs = obs.iloc[selected].copy()
+    subset_obs["source_row"] = selected
+    subset_obs.index = [f"sciplex3_row_{i}" for i in selected]
+    result = ad.AnnData(matrix, obs=subset_obs, var=pd.DataFrame(index=gene_names))
+    result.write_h5ad(output/"sciplex_prototype.h5ad", compression="gzip")
+    manifest = {"source_url":URL,"source_record":"scPerturb corrected v1.4 / Zenodo 13350497",
+                "source_md5":"c9e70629505d98c7ca1a837f62b14e89", "full_source_checksum_recomputed":False,
+                "smiles_source":DRUG_URL,"molecules":len(matches),"cells":len(selected),
+                "gene_columns":matrix.shape[1],"nonzero_counts":matrix.nnz,"cell_line":args.cell_line,
+                "exposure_hours":24,"dose_um":[.1,1.],"seed":args.seed,
+                "prototype_only":True,"elapsed_seconds":time.perf_counter()-start,
+                "selected_rows_sha256":hashlib.sha256(selected.tobytes()).hexdigest()}
+    (output/"source_manifest.json").write_text(json.dumps(manifest, indent=2)+"\n")
+    print(json.dumps(manifest,indent=2))
+
+if __name__ == "__main__": main()
diff --git a/scripts/prepare_tahoe.py b/scripts/prepare_tahoe.py
new file mode 100644
index 0000000000000000000000000000000000000000..35db6b941033b293be8dda7062571e551db12c52
--- /dev/null
+++ b/scripts/prepare_tahoe.py
@@ -0,0 +1,78 @@
+"""Convert official Tahoe expression shards to a bounded raw-count AnnData cohort."""
+import argparse
+import ast
+from collections import Counter
+from pathlib import Path
+import json
+import anndata as ad
+import numpy as np
+import pandas as pd
+import pyarrow.parquet as pq
+from scipy import sparse
+from remedi.data import UNIT_TO_UM
+
+
+def concentration(value, drug):
+    terms = ast.literal_eval(value)
+    if len(terms) != 1:
+        raise ValueError('Only single-molecule perturbations are supported')
+    name, dose, unit = terms[0]
+    if name.strip() != drug.strip():
+        raise ValueError('Sample compound differs from expression compound')
+    return float(dose) * UNIT_TO_UM[unit]
+
+
+def convert(shards, metadata, output, cell_line=None, cap=128, seed=0):
+    meta = Path(metadata)
+    samples = pd.DataFrame(pq.read_table(meta/'sample_metadata.parquet').to_pylist()).set_index('sample')
+    genes = pd.DataFrame(pq.read_table(meta/'gene_metadata.parquet').to_pylist())
+    gene_ids = genes.ensembl_id.astype(str).tolist()
+    if len(set(gene_ids)) != len(gene_ids): raise ValueError('Duplicate gene identifiers')
+    token = dict(zip(genes.token_id.astype(int), range(len(genes))))
+    rng = np.random.default_rng(seed)
+    pools, seen = {}, Counter()
+    for path in shards:
+        for batch in pq.ParquetFile(path).iter_batches(batch_size=1024):
+            for row in batch.to_pylist():
+                if cell_line and row['cell_line_id'] != cell_line: continue
+                sample = samples.loc[row['sample']]
+                if isinstance(sample, pd.DataFrame): raise ValueError('Duplicate sample identifiers')
+                if sample.plate != row['plate']: raise ValueError('Plate metadata mismatch')
+                dose = 0. if row['drug']=='DMSO_TF' else concentration(sample.drugname_drugconc, row['drug'])
+                key = (row['drug'], dose, row['cell_line_id'], row['plate'])
+                seen[key] += 1
+                pool = pools.setdefault(key, [])
+                j = len(pool) if len(pool) < cap else int(rng.integers(seen[key]))
+                if j >= cap: continue
+                # The first gene/expression pair is the released CLS token.
+                ids, values = row['genes'][1:], row['expressions'][1:]
+                if len(ids) != len(values): raise ValueError('Misaligned genes and counts')
+                cols = [token[int(i)] for i in ids]
+                record = ({k: row[k] for k in ['drug','sample','cell_line_id','plate','canonical_smiles','BARCODE_SUB_LIB_ID']}, cols, values, dose)
+                if j == len(pool): pool.append(record)
+                else: pool[j] = record
+        print(f'Read {path}. Retained {sum(map(len,pools.values()))} cells', flush=True)
+    observations, columns, values, pointers = [], [], [], [0]
+    for key in sorted(pools):
+        for row, cols, vals, dose in pools[key]:
+            row['dose_um'] = dose
+            observations.append(row); columns.extend(cols); values.extend(vals); pointers.append(len(values))
+    if not observations: raise ValueError('No matching cells')
+    obs = pd.DataFrame(observations)
+    obs.index = obs.BARCODE_SUB_LIB_ID.astype(str)
+    if not obs.index.is_unique: raise ValueError('Duplicate source cells across shards')
+    x = sparse.csr_matrix((np.asarray(values, dtype=np.float32), columns, pointers), shape=(len(obs),len(genes)))
+    out = Path(output); out.mkdir(parents=True,exist_ok=True)
+    ad.AnnData(x, obs=obs, var=genes.set_index('ensembl_id')).write_h5ad(out/'tahoe.h5ad',compression='gzip')
+    obs.loc[obs.drug!='DMSO_TF',['drug','canonical_smiles']].drop_duplicates().rename(columns={'canonical_smiles':'smiles'}).to_csv(out/'structures.csv',index=False)
+    (out/'source.json').write_text(json.dumps({'dataset':'tahoebio/Tahoe-100M','shards':[str(p) for p in shards], 'cells':len(obs),'conditions':len(pools),'reservoir_cap':cap,'seed':seed},indent=2)+'\n')
+
+if __name__ == '__main__':
+    p=argparse.ArgumentParser()
+    p.add_argument('--shards',nargs='+',required=True)
+    p.add_argument('--metadata',required=True)
+    p.add_argument('--output',required=True)
+    p.add_argument('--cell-line')
+    p.add_argument('--cap',type=int,default=128)
+    p.add_argument('--seed',type=int,default=0)
+    convert(**vars(p.parse_args()))
diff --git a/scripts/run_ablations.py b/scripts/run_ablations.py
new file mode 100644
index 0000000000000000000000000000000000000000..e32cecdb73c9aebe1e680c6c27913460af394ee8
--- /dev/null
+++ b/scripts/run_ablations.py
@@ -0,0 +1,25 @@
+"""Run uncertainty ablations with variant-specific calibration."""
+import argparse
+from pathlib import Path
+import joblib
+from remedi.io import ResponseData
+from remedi.pipeline import calibrate,evaluate
+from remedi.plotting import plot_summary
+
+p=argparse.ArgumentParser()
+for name in ['data','model','output']:p.add_argument('--'+name,required=True)
+p.add_argument('--tolerance',type=float,required=True)
+p.add_argument('--max-queries',type=int,default=20)
+p.add_argument('--scenarios',type=int,default=32)
+p.add_argument('--seed',type=int,default=0)
+p.add_argument('--embeddings')
+a=p.parse_args();d=ResponseData.load(a.data);m=joblib.load(Path(a.model)/'model.joblib')
+for variant in ['none','permuted','diagonal','no-floor','no-target-sampling']:
+    out=Path(a.output)/variant
+    cal_variant='none' if variant=='no-target-sampling' else variant
+    c=calibrate(d,m,out/'calibration',ablation=cal_variant,scenarios=a.scenarios,
+                seed=a.seed,embeddings=a.embeddings,panel_molecules=3)
+    evaluate(d,m,c,out/'evaluation',ablation=variant,scenarios=a.scenarios,
+             seed=a.seed,embeddings=a.embeddings,tolerance=a.tolerance,max_queries=a.max_queries)
+    plot_summary(out/'evaluation/summary.csv',out/'plots',title=variant)
+    print(f'Completed {variant}',flush=True)
diff --git a/scripts/run_pipeline.py b/scripts/run_pipeline.py
new file mode 100644
index 0000000000000000000000000000000000000000..f26bbc9a8954d143b93d7e391fe1514f9e486fc9
--- /dev/null
+++ b/scripts/run_pipeline.py
@@ -0,0 +1,68 @@
+"""Run preparation, fitting, calibration, evaluation, generation, and figures."""
+import argparse
+from pathlib import Path
+import json
+import anndata as ad
+import pandas as pd
+import numpy as np
+import yaml
+from remedi.data import standardize_obs, prepare_h5ad
+from remedi.splits import make_splits
+from remedi.models import train
+from remedi.pipeline import calibrate, evaluate, generate_catalog
+from remedi.plotting import plot_summary
+
+p=argparse.ArgumentParser()
+p.add_argument('--h5ad',required=True)
+p.add_argument('--mapping',required=True)
+p.add_argument('--structures',required=True)
+p.add_argument('--output',required=True)
+p.add_argument('--seed',type=int,default=0)
+p.add_argument('--max-queries',type=int,default=6)
+p.add_argument('--min-cells',type=int,default=2)
+p.add_argument('--max-cells',type=int,default=32)
+p.add_argument('--dimensions',type=int,default=16)
+p.add_argument('--genes',type=int,default=500)
+p.add_argument('--ensemble',type=int,default=3)
+p.add_argument('--scenarios',type=int,default=16)
+p.add_argument('--tolerance',type=float,required=True,help='Preselected endpoint-success threshold')
+p.add_argument('--representation',default='pca')
+p.add_argument('--embeddings')
+a=p.parse_args()
+out=Path(a.output);out.mkdir(parents=True,exist_ok=True)
+mapping=yaml.safe_load(Path(a.mapping).read_text())
+structures=pd.read_csv(a.structures)
+x=ad.read_h5ad(a.h5ad)
+obs=standardize_obs(x.obs,mapping,structures)
+# Retain compounds with an eligible two-block query at one or more doses.
+counts=obs[~obs.is_control].groupby(['smiles','dose_um','context','block','control_id']).size()
+valid=counts[counts>=a.min_cells].reset_index()
+valid=valid.groupby(['smiles','dose_um','context']).block.nunique()
+eligible=set(valid[valid>=2].reset_index().smiles)
+keep=obs.is_control|obs.smiles.isin(eligible)
+x=x[keep.to_numpy()].copy();x.write_h5ad(out/'cohort.h5ad',compression='gzip')
+smiles=sorted(eligible)
+if len(smiles)<9:raise ValueError(f'Need at least 9 compounds with two blocks. Found {len(smiles)}. Increase input cohort.')
+# A small engineering cohort needs three calibration and two test compounds.
+fractions=((len(smiles)-6)/len(smiles),1/len(smiles),3/len(smiles),2/len(smiles)) if len(smiles)<30 else (.55,.1,.2,.15)
+splits=make_splits(pd.DataFrame({'smiles':smiles}),a.seed,fractions)
+splits.to_csv(out/'splits.csv',index=False)
+d=prepare_h5ad(out/'cohort.h5ad',out/'data',mapping,splits,structures,
+               representation=a.representation,max_cells=a.max_cells,min_cells=a.min_cells,
+               genes=a.genes,dimensions=a.dimensions,seed=a.seed)
+m=train(d,out/'model',ensemble=a.ensemble,folds=3,seed=a.seed,
+        encoder='cached' if a.embeddings else 'fingerprint',embeddings=a.embeddings)
+c=calibrate(d,m,out/'calibration',scenarios=a.scenarios,panel_molecules=3,seed=a.seed,embeddings=a.embeddings)
+e=evaluate(d,m,c,out/'evaluation',scenarios=a.scenarios,max_queries=a.max_queries,
+           tolerance=a.tolerance,seed=a.seed,embeddings=a.embeddings)
+plot_summary(out/'evaluation/summary.csv',out/'plots')
+query=d.obs[d.obs.split=='test'].groupby(['molecule_id','dose_um','context'],sort=True)
+indices=next(g.index.tolist() for _,g in query if len(g)>=2)
+catalog=pd.DataFrame({'smiles':smiles,'catalog_id':[f'cohort_{i}' for i in range(len(smiles))]})
+catalog.to_csv(out/'example_catalog.csv',index=False)
+generate_catalog(d,m,c,catalog,indices,out/'generation',scenarios=a.scenarios,
+                 seed=a.seed,embeddings=a.embeddings,doses=sorted(d.obs.dose_um.unique()))
+(out/'run.json').write_text(json.dumps({'source':a.h5ad,'cells':x.n_obs,'genes':x.n_vars,
+ 'compounds':len(smiles),'conditions':len(d.obs),'queries':int(e.query_id.nunique()),
+ 'representation':a.representation,'seed':a.seed,'example_catalog':'measured cohort, not Enamine REAL'},indent=2)+'\n')
+print((out/'run.json').read_text())
diff --git a/scripts/validate_release.py b/scripts/validate_release.py
new file mode 100644
index 0000000000000000000000000000000000000000..61a457c6d03a852a73c10a23b271a5f5fbcfb039
--- /dev/null
+++ b/scripts/validate_release.py
@@ -0,0 +1,42 @@
+"""Reproduce the bundled measured-data checks, including frozen external transfer."""
+import argparse
+import json
+from pathlib import Path
+import subprocess
+import sys
+import joblib
+import pandas as pd
+import yaml
+from remedi.data import prepare_h5ad
+from remedi.pipeline import evaluate
+from remedi.plotting import plot_summary
+
+p=argparse.ArgumentParser()
+p.add_argument('--output',default='runs/validation')
+a=p.parse_args()
+root=Path(__file__).resolve().parents[1]
+out=Path(a.output).resolve();out.mkdir(parents=True,exist_ok=True)
+
+def run(*args):
+    subprocess.run([sys.executable,*map(str,args)],check=True,cwd=root)
+
+for study,filename in [('tahoe','tahoe.h5ad'),('sciplex','sciplex_prototype.h5ad')]:
+    run(root/'scripts/run_pipeline.py','--h5ad',root/'examples'/study/filename,
+        '--mapping',root/'configs'/f'{study}.yaml','--structures',root/'examples'/study/'structures.csv',
+        '--output',out/study,'--min-cells',2,'--dimensions',8,'--genes',250,'--max-queries',4,'--tolerance',1)
+run(root/'scripts/align_genes.py','--inputs',root/'examples/tahoe/tahoe.h5ad',
+    root/'examples/sciplex/sciplex_prototype.h5ad','--gene-columns','gene_symbol','index',
+    '--drop-ambiguous','--output',out/'aligned')
+run(root/'scripts/run_pipeline.py','--h5ad',out/'aligned/0_tahoe.h5ad',
+    '--mapping',root/'configs/tahoe.yaml','--structures',root/'examples/tahoe/structures.csv',
+    '--output',out/'transfer-source','--min-cells',2,'--dimensions',8,'--genes',250,'--max-queries',4,'--tolerance',1)
+s=pd.read_csv(root/'examples/sciplex/structures.csv')
+splits=s[['smiles']].copy();splits['split']='test'
+d=prepare_h5ad(out/'aligned/1_sciplex_prototype.h5ad',out/'transfer/data',
+    yaml.safe_load((root/'configs/sciplex.yaml').read_text()),splits,s,
+    feature_model=out/'transfer-source/data/cell_feature_model.joblib',min_cells=2,max_cells=32)
+m=joblib.load(out/'transfer-source/model/model.joblib')
+c=json.loads((out/'transfer-source/calibration/calibration.json').read_text())
+evaluate(d,m,c,out/'transfer/evaluation',external_track='unseen',scenarios=16,max_queries=4,tolerance=1)
+plot_summary(out/'transfer/evaluation/summary.csv',out/'transfer/plots')
+print(f'Completed measured-data integration checks in {out}')
diff --git a/src/remedi/__init__.py b/src/remedi/__init__.py
new file mode 100644
index 0000000000000000000000000000000000000000..20fe1ebc7005e03143c6e5e146e5a56356f704ad
--- /dev/null
+++ b/src/remedi/__init__.py
@@ -0,0 +1,3 @@
+"""ReMEDi reference implementation. No pretrained response weights are distributed."""
+__version__ = "0.1.0"
+
diff --git a/src/remedi/chemistry.py b/src/remedi/chemistry.py
new file mode 100644
index 0000000000000000000000000000000000000000..03042b18c8398006f898bfa8d91977f7ae472075
--- /dev/null
+++ b/src/remedi/chemistry.py
@@ -0,0 +1,90 @@
+"""Canonical structure identity, fingerprints, and frozen embedding adapters."""
+import hashlib
+import numpy as np
+import pandas as pd
+
+def canonical_smiles(smiles):
+    from rdkit import Chem
+    mol = Chem.MolFromSmiles(str(smiles))
+    if mol is None:
+        raise ValueError(f"Invalid SMILES: {smiles}")
+    # Preserve stereochemistry and disconnected components; do not silently desalt.
+    return Chem.MolToSmiles(mol, canonical=True, isomericSmiles=True)
+
+def molecule_id(smiles):
+    return hashlib.sha256(canonical_smiles(smiles).encode()).hexdigest()[:20]
+
+def fingerprints(smiles, bits=1024, radius=2):
+    from rdkit import Chem
+    from rdkit.Chem import rdFingerprintGenerator
+    gen = rdFingerprintGenerator.GetMorganGenerator(radius=radius, fpSize=bits)
+    out = []
+    for s in smiles:
+        mol = Chem.MolFromSmiles(str(s))
+        if mol is None:
+            raise ValueError(f"Invalid SMILES: {s}")
+        out.append(gen.GetFingerprintAsNumPy(mol))
+    return np.asarray(out, dtype=float)
+
+def max_tanimoto(query, reference, batch_size=1024):
+    query, reference = np.asarray(query, float), np.asarray(reference, float)
+    if reference.shape[0] == 0:
+        raise ValueError("Chemical-distance calculation needs training structures")
+    out = []
+    for start in range(0, len(query), batch_size):
+        q = query[start:start+batch_size]
+        overlap = q @ reference.T
+        union = q.sum(1)[:, None] + reference.sum(1)[None, :] - overlap
+        similarity = np.divide(overlap, union, out=np.zeros_like(overlap), where=union > 0)
+        out.extend(similarity.max(1))
+    return np.asarray(out)
+
+def encode(smiles, kind="fingerprint", embedding_file=None):
+    if kind == "fingerprint":
+        return fingerprints(smiles)
+    if kind != "cached" or embedding_file is None:
+        raise ValueError("Use fingerprint or cached molecular embeddings")
+    with np.load(embedding_file, allow_pickle=False) as a:
+        keys = [canonical_smiles(s) for s in a["smiles"].astype(str)]
+        values = a["embeddings"]
+    if len(keys) != len(set(keys)) or values.shape[0] != len(keys):
+        raise ValueError("Embedding cache has duplicate or misaligned structures")
+    lookup = dict(zip(keys, values))
+    missing = [s for s in smiles if canonical_smiles(s) not in lookup]
+    if missing:
+        raise ValueError(f"Embedding cache misses {len(missing)} input structures")
+    return np.stack([lookup[canonical_smiles(s)] for s in smiles]).astype(float)
+
+def cache_molformer(smiles, output, model_id, revision, batch_size=32, device="cpu"):
+    """Load a revision-pinned released encoder; keep parameters frozen."""
+    import torch
+    from transformers import AutoModel, AutoTokenizer
+    if not revision or revision == "main":
+        raise ValueError("Supply a concrete checkpoint commit hash via --revision")
+    tokenizer = AutoTokenizer.from_pretrained(model_id, revision=revision, trust_remote_code=True)
+    model = AutoModel.from_pretrained(model_id, revision=revision, trust_remote_code=True,
+                                      deterministic_eval=True).eval().to(device)
+    unique = sorted({canonical_smiles(s) for s in smiles})
+    vectors = []
+    with torch.inference_mode():
+        for start in range(0, len(unique), batch_size):
+            tokens = tokenizer(unique[start:start+batch_size], padding=True, return_tensors="pt")
+            result = model(**{k: v.to(device) for k, v in tokens.items()})
+            if getattr(result, "pooler_output", None) is None:
+                raise ValueError("Checkpoint lacks the documented MolFormer pooled output")
+            vectors.append(result.pooler_output.cpu().numpy())
+    np.savez_compressed(output, smiles=np.asarray(unique), embeddings=np.concatenate(vectors),
+                        model_id=np.asarray(model_id), revision=np.asarray(revision))
+
+def verify_catalog(frame):
+    required = {"catalog_id", "smiles"}
+    if not required.issubset(frame):
+        raise ValueError("Catalog needs catalog_id and smiles columns")
+    frame = frame.copy()
+    if frame.catalog_id.isna().any() or frame.catalog_id.duplicated().any():
+        raise ValueError("Catalog IDs must be present and unique")
+    frame["smiles"] = frame.smiles.map(canonical_smiles)
+    frame["molecule_id"] = frame.smiles.map(molecule_id)
+    # One deterministic representative ID per exact canonical structure.
+    return frame.sort_values("catalog_id").drop_duplicates("smiles").reset_index(drop=True)
+
diff --git a/src/remedi/cli.py b/src/remedi/cli.py
new file mode 100644
index 0000000000000000000000000000000000000000..72da23a365cfc3b05e2ff020e1a4d5d562f86696
--- /dev/null
+++ b/src/remedi/cli.py
@@ -0,0 +1,112 @@
+"""Reproducible command-line entry points for the complete reference pipeline."""
+import argparse
+import json
+from pathlib import Path
+import joblib
+import pandas as pd
+import yaml
+from .io import ResponseData
+
+def parser():
+    p = argparse.ArgumentParser(prog="remedi", description="Molecule-dose endpoint generation")
+    sub = p.add_subparsers(dest="command", required=True)
+    d = sub.add_parser("demo", help="Run a synthetic end-to-end engineering check")
+    d.add_argument("--output", required=True); d.add_argument("--seed", type=int, default=0)
+    d.add_argument("--max-queries", type=int, default=8)
+    s = sub.add_parser("split", help="Create grouped structure/scaffold assignments")
+    s.add_argument("--molecules", required=True); s.add_argument("--output", required=True)
+    s.add_argument("--seed", type=int, default=0); s.add_argument("--scaffold", action="store_true")
+    s.add_argument("--fractions", type=float, nargs=4, default=[.6,.15,.1,.15])
+    q = sub.add_parser("prepare", help="Aggregate a local RNA/embedding AnnData cohort")
+    q.add_argument("--h5ad", required=True); q.add_argument("--mapping", required=True)
+    q.add_argument("--splits", required=True); q.add_argument("--structures")
+    q.add_argument("--output", required=True); q.add_argument("--representation", default="pca")
+    q.add_argument("--feature-model"); q.add_argument("--max-cells", type=int, default=128)
+    q.add_argument("--min-cells", type=int, default=8); q.add_argument("--genes", type=int, default=1000)
+    q.add_argument("--dimensions", type=int, default=32); q.add_argument("--seed", type=int, default=0)
+    q.add_argument("--layer")
+    t = sub.add_parser("train", help="Fit grouped response ensemble and cross-fitted residuals")
+    t.add_argument("--data", required=True); t.add_argument("--output", required=True)
+    t.add_argument("--head", choices=["ridge","mlp"], default="ridge")
+    t.add_argument("--alpha", type=float, default=10.); t.add_argument("--hidden", type=int, default=128)
+    t.add_argument("--ensemble", type=int, default=5); t.add_argument("--folds", type=int, default=3)
+    t.add_argument("--encoder", choices=["fingerprint","cached"], default="fingerprint")
+    t.add_argument("--embeddings"); t.add_argument("--seed", type=int, default=0)
+    c = sub.add_parser("calibrate", help="Estimate an empirical held-out panel radius")
+    for name in ["data","model","output"]: c.add_argument("--"+name, required=True)
+    c.add_argument("--coverage", type=float, default=.9); c.add_argument("--panel-molecules", type=int, default=6)
+    c.add_argument("--scenarios", type=int, default=32); c.add_argument("--seed", type=int, default=0)
+    c.add_argument("--ablation", choices=["none","permuted","diagonal","no-floor"], default="none")
+    c.add_argument("--embeddings")
+    e = sub.add_parser("evaluate", help="Compare generators on hidden measured candidate responses")
+    for name in ["data","model","calibration","output"]: e.add_argument("--"+name, required=True)
+    e.add_argument("--split", default="test"); e.add_argument("--scenarios", type=int, default=32)
+    e.add_argument("--tau", type=float, default=.05); e.add_argument("--tolerance", type=float, required=True,
+                       help="Endpoint success threshold fixed on validation data")
+    e.add_argument("--seed", type=int, default=0); e.add_argument("--max-queries", type=int, default=30)
+    e.add_argument("--methods", nargs="+"); e.add_argument("--embeddings"); e.add_argument("--metric-data")
+    e.add_argument("--external-track", choices=["none","shared","unseen"], default="none")
+    e.add_argument("--ablation", choices=["none","permuted","diagonal","no-floor","no-target-sampling"], default="none")
+    g = sub.add_parser("generate", help="Sample exact members of a supplied verified chemical catalog")
+    for name in ["data","model","calibration","catalog","output"]: g.add_argument("--"+name, required=True)
+    g.add_argument("--target-rows", type=int, nargs="+", required=True)
+    g.add_argument("--doses", type=float, nargs="+", default=[.05,.5,5.])
+    g.add_argument("--scenarios", type=int, default=32); g.add_argument("--tau", type=float, default=.05)
+    g.add_argument("--shortlist", type=int, default=512); g.add_argument("--draws", type=int, default=20)
+    g.add_argument("--seed", type=int, default=0); g.add_argument("--embeddings")
+    z = sub.add_parser("plot", help="Draw Ubuntu plots from measured result CSVs")
+    z.add_argument("--summary", required=True); z.add_argument("--output", required=True)
+    z.add_argument("--font-dir"); z.add_argument("--title", default="Retrospective endpoint evaluation")
+    r = sub.add_parser("plot-curve", help="Plot calibration, scaling, or runtime records")
+    for name in ["csv","output","x","y"]: r.add_argument("--"+name, required=True)
+    r.add_argument("--hue", default="method"); r.add_argument("--lower"); r.add_argument("--upper"); r.add_argument("--font-dir")
+    m = sub.add_parser("encode-molecules", help="Cache a revision-pinned frozen MolFormer checkpoint")
+    m.add_argument("--molecules", required=True); m.add_argument("--output", required=True)
+    m.add_argument("--model-id", default="ibm/MoLFormer-XL-both-10pct"); m.add_argument("--revision", required=True)
+    m.add_argument("--batch-size", type=int, default=32); m.add_argument("--device", default="cpu")
+    return p
+
+def main(argv=None):
+    a = vars(parser().parse_args(argv)); command = a.pop("command")
+    if command == "demo":
+        from .demo import run_demo
+        print(run_demo(**a)); return
+    if command == "split":
+        from .splits import make_splits
+        frame = make_splits(pd.read_csv(a.pop("molecules")), seed=a["seed"], fractions=a["fractions"], scaffold=a["scaffold"])
+        Path(a["output"]).parent.mkdir(parents=True, exist_ok=True); frame.to_csv(a["output"], index=False); return
+    if command == "prepare":
+        from .data import prepare_h5ad
+        a["mapping"] = yaml.safe_load(Path(a["mapping"]).read_text())
+        a["splits"] = pd.read_csv(a["splits"])
+        a["structures"] = pd.read_csv(a["structures"]) if a["structures"] else None
+        a["path"] = a.pop("h5ad"); data = prepare_h5ad(**a)
+        print(f"Prepared {len(data.obs)} conditions in {data.response.shape[1]} dimensions"); return
+    if command == "train":
+        from .models import train
+        a["data"] = ResponseData.load(a["data"]); a["kind"] = a.pop("head")
+        train(**a); print("Saved fitted ensemble and cross-fitted residuals"); return
+    if command in ["calibrate","evaluate","generate"]:
+        from .pipeline import calibrate, evaluate, generate_catalog
+        a["data"] = ResponseData.load(a["data"])
+        a["model"] = joblib.load(Path(a["model"])/"model.joblib")
+        if "calibration" in a: a["calibration"] = json.loads((Path(a["calibration"])/"calibration.json").read_text())
+        if command == "calibrate": print(json.dumps(calibrate(**a), indent=2)); return
+        if command == "evaluate":
+            if a["metric_data"]: a["metric_data"] = ResponseData.load(a["metric_data"])
+            frame = evaluate(**a); print(f"Evaluated {frame.query_id.nunique()} queries"); return
+        a["catalog"] = pd.read_csv(a["catalog"]); a["target_indices"] = a.pop("target_rows")
+        frame = generate_catalog(**a); print(f"Generated a distribution on {len(frame)} verified pairs"); return
+    if command == "plot":
+        from .plotting import plot_summary
+        plot_summary(**a); return
+    if command == "plot-curve":
+        from .plotting import plot_curve
+        plot_curve(**a); return
+    if command == "encode-molecules":
+        from .chemistry import cache_molformer
+        a["smiles"] = pd.read_csv(a.pop("molecules")).smiles
+        cache_molformer(**a)
+
+if __name__ == "__main__":
+    main()
diff --git a/src/remedi/data.py b/src/remedi/data.py
new file mode 100644
index 0000000000000000000000000000000000000000..e9b6abf13077edc3f1145673377ace68bc7d7845
--- /dev/null
+++ b/src/remedi/data.py
@@ -0,0 +1,166 @@
+"""AnnData ingestion with explicit units, control matching, and train-only features."""
+from pathlib import Path
+import json
+import joblib
+import numpy as np
+import pandas as pd
+from scipy import sparse
+from sklearn.decomposition import PCA
+from .chemistry import canonical_smiles, molecule_id
+from .io import ResponseData, fingerprint_config
+
+UNIT_TO_UM = {"uM": 1., "µM": 1., "μM": 1., "nM": .001, "mM": 1000., "M": 1e6}
+
+def standardize_obs(obs, mapping, structures=None):
+    """Map original fields without guessing dose units or replicate identities."""
+    required = ["drug", "dose", "context", "block", "control_group"]
+    for name in required:
+        if name not in mapping or mapping[name] not in obs:
+            raise ValueError(f"Mapping requires an existing {name} column")
+    unit = mapping.get("dose_unit")
+    if unit not in UNIT_TO_UM:
+        raise ValueError("Specify dose_unit explicitly as uM, nM, mM, or M")
+    result = pd.DataFrame(index=obs.index)
+    for name in required:
+        result[name] = obs[mapping[name]].to_numpy()
+    result["is_control"] = result.drug.astype(str).isin(mapping.get("control_names", ["DMSO"]))
+    result["dose_um"] = pd.to_numeric(result.dose, errors="raise")*UNIT_TO_UM[unit]
+    result.loc[result.is_control, "dose_um"] = 0.
+    for field in ["drug", "context", "block", "control_group"]:
+        if result[field].isna().any(): raise ValueError(f"Missing {field} identifiers")
+        result[field] = result[field].astype(str)
+    result["control_id"] = result.context+"|"+result.control_group
+    if "smiles" in mapping:
+        result["smiles"] = obs[mapping["smiles"]].to_numpy()
+    elif structures is not None:
+        if structures.drug.duplicated().any(): raise ValueError("Ambiguous drug-to-structure mapping")
+        result["smiles"] = result.drug.map(structures.set_index("drug").smiles)
+    else:
+        raise ValueError("Provide a SMILES column or a drug,smiles mapping file")
+    result.loc[result.is_control, "smiles"] = "CS(C)=O"
+    if result.smiles.isna().any():
+        missing = result.loc[result.smiles.isna(), "drug"].unique().tolist()
+        raise ValueError(f"Unmapped structures: {missing[:10]}")
+    canonical = {s: canonical_smiles(s) for s in result.smiles.unique()}
+    result["smiles"] = result.smiles.map(canonical)
+    result["molecule_id"] = result.smiles.map(molecule_id)
+    if not np.isfinite(result.dose_um).all(): raise ValueError("Concentrations must be finite")
+    if ((result.dose_um <= 0) & ~result.is_control).any():
+        raise ValueError("Non-control records have nonpositive concentration")
+    result["source_cell_id"] = result.index.astype(str)
+    return result.reset_index(drop=True)
+
+def normalize_counts(x):
+    if sparse.issparse(x):
+        x = x.astype(np.float64).tocsr()
+        if not np.isfinite(x.data).all() or (x.data < 0).any():
+            raise ValueError("RNA input must contain finite nonnegative counts")
+        total = np.asarray(x.sum(1)).ravel()
+        x = sparse.diags(np.divide(1e4, total, out=np.zeros_like(total), where=total > 0))@x
+        x.data = np.log1p(x.data)
+        return x.tocsr()
+    x = np.asarray(x)
+    x = np.asarray(x, dtype=np.float64)
+    if not np.isfinite(x).all() or (x < 0).any():
+        raise ValueError("RNA input must contain finite nonnegative counts")
+    total = x.sum(axis=1, keepdims=True)
+    return np.log1p(np.divide(x*1e4, total, out=np.zeros_like(x), where=total > 0))
+
+def prepare_h5ad(path, output, mapping, splits, structures=None, representation="pca",
+                 feature_model=None, max_cells=128, min_cells=8, genes=1000, dimensions=32,
+                 seed=0, layer=None, fit_cell_cap=4096):
+    import anndata as ad
+    a = ad.read_h5ad(path, backed="r")
+    obs = standardize_obs(a.obs, mapping, structures)
+    splits = splits.copy()
+    splits["smiles"] = splits.smiles.map(canonical_smiles)
+    if splits.smiles.duplicated().any(): raise ValueError("Split manifest has duplicate structures")
+    obs["split"] = obs.smiles.map(splits.set_index("smiles").split)
+    obs.loc[obs.is_control, "split"] = "control"
+    if obs.loc[~obs.is_control, "split"].isna().any():
+        raise ValueError("Split manifest is missing treated structures")
+    rng = np.random.default_rng(seed)
+    group_fields = ["molecule_id", "dose_um", "context", "block", "control_id", "is_control"]
+    groups = []
+    for _, group in obs.groupby(group_fields, sort=True):
+        if len(group) < min_cells: continue
+        idx = group.index.to_numpy()
+        groups.append(np.sort(rng.choice(idx, min(max_cells, len(idx)), replace=False)))
+    if not groups: raise ValueError("No groups pass the minimum-cell threshold")
+    selected = np.sort(np.concatenate(groups))
+    position = {int(v): k for k, v in enumerate(selected)}
+    source = a.layers[layer] if layer else a.X
+    if not a.var_names.is_unique: raise ValueError("Gene identifiers must be unique")
+    if min_cells < 2 or max_cells < min_cells: raise ValueError("Require max_cells >= min_cells >= 2")
+    if representation.startswith("obsm:"):
+        name = representation.split(":", 1)[1]
+        if name not in a.obsm: raise ValueError(f"Missing frozen cell embeddings in obsm[{name!r}]")
+        features = np.asarray(a.obsm[name][selected], float)
+        feature_metadata = {"representation": representation, "dimensions": features.shape[1],
+                            "checkpoint": mapping.get("cell_checkpoint", "unrecorded")}
+        if feature_metadata["checkpoint"] == "unrecorded":
+            raise ValueError("Record cell_checkpoint in the mapping for frozen embeddings")
+    elif representation == "pca":
+        # Work in bounded batches. Feature selection sees training treated cells only.
+        if feature_model:
+            projector = joblib.load(feature_model)
+            order = a.var_names.get_indexer(projector["genes"])
+            if (order < 0).any():
+                raise ValueError("External cohort misses fitted genes. Align both raw cohorts to a common gene set before fitting.")
+        else:
+            order = np.arange(a.n_vars)
+            train = selected[obs.iloc[selected].split.to_numpy() == "train"]
+            if len(train) < 3: raise ValueError("PCA requires training cells")
+            train = np.sort(rng.choice(train, min(fit_cell_cap, len(train)), replace=False))
+            x = normalize_counts(source[train])
+            variance = np.asarray(x.power(2).mean(0)).ravel()-np.asarray(x.mean(0)).ravel()**2 if sparse.issparse(x) else x.var(0)
+            keep = np.argsort(variance)[-min(genes, x.shape[1]):]
+            d = min(dimensions, len(train)-1, len(keep))
+            xx = x[:, keep].toarray() if sparse.issparse(x) else x[:, keep]
+            pca = PCA(n_components=d, random_state=seed).fit(xx)
+            projector = {"genes": list(a.var_names), "keep": keep, "pca": pca,
+                         "seed": seed, "fit_cells": len(train)}
+        pieces = []
+        for start in range(0, len(selected), 512):
+            x = normalize_counts(source[selected[start:start+512]][:, order])
+            xx = x[:, projector["keep"]]
+            pieces.append(projector["pca"].transform(xx.toarray() if sparse.issparse(xx) else xx))
+        features = np.concatenate(pieces)
+        Path(output).mkdir(parents=True, exist_ok=True)
+        joblib.dump(projector, Path(output)/"cell_feature_model.joblib")
+        feature_metadata = {"representation": "training-only log-count PCA", "dimensions": features.shape[1],
+                            "feature_genes": [projector["genes"][i] for i in projector["keep"]],
+                            "projection_hash": fingerprint_config({"components": projector["pca"].components_.tolist(), "mean": projector["pca"].mean_.tolist(), "normalization_genes": projector["genes"]})}
+    else:
+        raise ValueError("Cell representation must be pca or obsm:<frozen_embedding_key>")
+    summaries = []
+    controls = {}
+    for idx in groups:
+        row = obs.iloc[idx[0]].to_dict()
+        f = features[[position[int(i)] for i in idx]]
+        summary = (row, f.mean(0), f.var(0, ddof=1)/len(f), len(f))
+        if row["is_control"]:
+            # Multiple control samples in one pool are combined at the cell level.
+            controls.setdefault(row["control_id"], []).append(f)
+        else: summaries.append(summary)
+    controls = {k: np.concatenate(v) for k, v in controls.items()}
+    records, responses, contexts, tv, cv = [], [], [], [], []
+    dropped = 0
+    for row, mean, var, n in summaries:
+        if row["control_id"] not in controls:
+            dropped += 1; continue
+        f = controls[row["control_id"]]
+        c, v = f.mean(0), f.var(0, ddof=1)/len(f)
+        row["n_cells"], row["n_control_cells"] = n, len(f)
+        records.append({k: row[k] for k in ["molecule_id", "smiles", "drug", "dose_um", "context", "block", "control_id", "split", "n_cells", "n_control_cells"]})
+        responses.append(mean-c); contexts.append(c); tv.append(var); cv.append(v)
+    if a.isbacked: a.file.close()
+    if not records: raise ValueError("No treated conditions have matching controls")
+    metadata = {"source": str(path), "study": mapping.get("study"), "mapping": mapping, "seed": seed, "max_cells": max_cells,
+                "min_cells": min_cells, "conditions_without_controls": dropped,
+                "cell_features": feature_metadata, "feature_space_id": fingerprint_config(feature_metadata),
+                "sampling_model": "block bootstrap plus diagonal Gaussian sampling error of population means"}
+    data = ResponseData(pd.DataFrame(records), np.asarray(responses), np.asarray(contexts),
+                        np.asarray(tv), np.asarray(cv), metadata).validate()
+    data.save(output)
+    return data
diff --git a/src/remedi/demo.py b/src/remedi/demo.py
new file mode 100644
index 0000000000000000000000000000000000000000..e30f74dc92cfa1bd470f671de0be3d500577a5c3
--- /dev/null
+++ b/src/remedi/demo.py
@@ -0,0 +1,56 @@
+"""A synthetic engineering fixture. These data are not cellular measurements."""
+from pathlib import Path
+import numpy as np
+import pandas as pd
+from .chemistry import canonical_smiles, molecule_id, fingerprints
+from .splits import make_splits
+from .io import ResponseData, write_json
+from .models import train
+from .pipeline import calibrate, evaluate
+from .plotting import plot_summary
+
+def synthetic_data(seed=0, molecules=48, dimensions=8):
+    rng = np.random.default_rng(seed)
+    smiles = []
+    for chain in range(1, 9):
+        for ring in ["c1ccccc1", "c1ccncc1", "c1ccoc1", "c1ccsc1", "C1CCCCC1", "N1CCCCC1"]:
+            smiles.append(canonical_smiles("C"*chain+ring))
+    smiles = sorted(set(smiles))[:molecules]
+    splits = make_splits(pd.DataFrame({"smiles": smiles}), seed, (.55,.1,.2,.15))
+    fp = fingerprints(smiles)
+    projection = rng.normal(size=(fp.shape[1], dimensions))/np.sqrt(fp.shape[1])
+    effects = fp@projection
+    rows, response, context, tv, cv = [], [], [], [], []
+    for c in ["context_A", "context_B"]:
+        baseline = rng.normal(size=dimensions)*.2
+        for block in ["rep1", "rep2"]:
+            ctrl = baseline+rng.normal(size=dimensions)*.015
+            for m, s in enumerate(smiles):
+                for dose in [.05,.5,5.]:
+                    effect = effects[m]*(1+.25*np.log(dose))+.05*baseline
+                    observed = effect+rng.normal(size=dimensions)*.02-(ctrl-baseline)
+                    rows.append({"molecule_id": molecule_id(s), "smiles": s, "drug": f"synthetic_{m}",
+                                 "dose_um": dose, "context": c, "block": block, "control_id": c+"|"+block,
+                                 "split": splits.set_index("smiles").loc[s,"split"], "n_cells": 64,
+                                 "n_control_cells": 128})
+                    response.append(observed); context.append(ctrl)
+                    tv.append(np.full(dimensions,.02**2)); cv.append(np.full(dimensions,.015**2))
+    return ResponseData(pd.DataFrame(rows), np.array(response), np.array(context), np.array(tv),
+                         np.array(cv), {"source": "synthetic engineering fixture", "feature_space_id": "synthetic-v1",
+                                        "cell_features": {"representation": "synthetic", "dimensions": dimensions},
+                                        "biological_data": False}).validate()
+
+def run_demo(output, seed=0, max_queries=8):
+    output = Path(output); output.mkdir(parents=True, exist_ok=True)
+    data = synthetic_data(seed)
+    data.save(output/"data")
+    model = train(data, output/"model", ensemble=3, folds=3, seed=seed)
+    cal = calibrate(data, model, output/"calibration", panel_molecules=5, scenarios=16, seed=seed)
+    frame = evaluate(data, model, cal, output/"evaluation", scenarios=16, seed=seed,
+                      max_queries=max_queries, methods=["mean","gibbs","marginal","finite-minimax","remedi"])
+    plot_summary(output/"evaluation/summary.csv", output/"plots", title="Synthetic engineering check")
+    write_json(output/"run.json", {"source": "synthetic", "biological_validation": False,
+               "queries": int(frame.query_id.nunique()), "seed": seed,
+               "purpose": "Verify serialization, grouped training, calibration, optimization, evaluation, and plots"})
+    return output
+
diff --git a/src/remedi/download.py b/src/remedi/download.py
new file mode 100644
index 0000000000000000000000000000000000000000..87cd2eaeaae5e79ab4541c40b9ed24cbde0953e1
--- /dev/null
+++ b/src/remedi/download.py
@@ -0,0 +1,44 @@
+"""Seekable HTTP reader with persistent range caching and bounded retries."""
+import io
+import hashlib
+from pathlib import Path
+import requests
+from requests.adapters import HTTPAdapter
+from urllib3.util.retry import Retry
+
+class RangeReader(io.RawIOBase):
+    def __init__(self, url, cache, size, block_size=1024*1024):
+        self.url, self.size, self.block_size = url, size, block_size
+        self.position=0
+        self.cache=Path(cache)/hashlib.sha256(url.encode()).hexdigest()[:16]
+        self.cache.mkdir(parents=True,exist_ok=True)
+        self.session=requests.Session()
+        self.session.mount('https://',HTTPAdapter(max_retries=Retry(total=6,backoff_factor=2,
+          status_forcelist=[429,500,502,503,504],respect_retry_after_header=True)))
+    def readable(self):return True
+    def seekable(self):return True
+    def tell(self):return self.position
+    def seek(self,offset,whence=0):
+        self.position=int(offset if whence==0 else self.position+offset if whence==1 else self.size+offset)
+        if self.position<0:raise ValueError('Negative seek')
+        return self.position
+    def read(self,size=-1):
+        end=self.size if size<0 else min(self.size,self.position+size)
+        result=[]
+        while self.position<end:
+            index=self.position//self.block_size
+            start=index*self.block_size;stop=min(start+self.block_size,self.size)
+            path=self.cache/f'{index}.bin'
+            if not path.exists():
+                r=self.session.get(self.url,headers={'Range':f'bytes={start}-{stop-1}'},timeout=(30,120))
+                r.raise_for_status()
+                if r.status_code!=206 or len(r.content)!=stop-start:
+                    raise IOError('Server did not return the requested byte range. Use the full-file download.')
+                temp=path.with_suffix('.part');temp.write_bytes(r.content);temp.replace(path)
+            block=path.read_bytes();offset=self.position-start;n=min(end-self.position,stop-self.position)
+            result.append(block[offset:offset+n]);self.position+=n
+        return b''.join(result)
+    def readinto(self,b):
+        data=self.read(len(b));b[:len(data)]=data;return len(data)
+    def close(self):
+        self.session.close();super().close()
diff --git a/src/remedi/fonts/UFL.txt b/src/remedi/fonts/UFL.txt
new file mode 100644
index 0000000000000000000000000000000000000000..ae78a8f94eae372cb1ca4e14b67244342fce604e
--- /dev/null
+++ b/src/remedi/fonts/UFL.txt
@@ -0,0 +1,96 @@
+-------------------------------
+UBUNTU FONT LICENCE Version 1.0
+-------------------------------
+
+PREAMBLE
+This licence allows the licensed fonts to be used, studied, modified and
+redistributed freely. The fonts, including any derivative works, can be
+bundled, embedded, and redistributed provided the terms of this licence
+are met. The fonts and derivatives, however, cannot be released under
+any other licence. The requirement for fonts to remain under this
+licence does not require any document created using the fonts or their
+derivatives to be published under this licence, as long as the primary
+purpose of the document is not to be a vehicle for the distribution of
+the fonts.
+
+DEFINITIONS
+"Font Software" refers to the set of files released by the Copyright
+Holder(s) under this licence and clearly marked as such. This may
+include source files, build scripts and documentation.
+
+"Original Version" refers to the collection of Font Software components
+as received under this licence.
+
+"Modified Version" refers to any derivative made by adding to, deleting,
+or substituting -- in part or in whole -- any of the components of the
+Original Version, by changing formats or by porting the Font Software to
+a new environment.
+
+"Copyright Holder(s)" refers to all individuals and companies who have a
+copyright ownership of the Font Software.
+
+"Substantially Changed" refers to Modified Versions which can be easily
+identified as dissimilar to the Font Software by users of the Font
+Software comparing the Original Version with the Modified Version.
+
+To "Propagate" a work means to do anything with it that, without
+permission, would make you directly or secondarily liable for
+infringement under applicable copyright law, except executing it on a
+computer or modifying a private copy. Propagation includes copying,
+distribution (with or without modification and with or without charging
+a redistribution fee), making available to the public, and in some
+countries other activities as well.
+
+PERMISSION & CONDITIONS
+This licence does not grant any rights under trademark law and all such
+rights are reserved.
+
+Permission is hereby granted, free of charge, to any person obtaining a
+copy of the Font Software, to propagate the Font Software, subject to
+the below conditions:
+
+1) Each copy of the Font Software must contain the above copyright
+notice and this licence. These can be included either as stand-alone
+text files, human-readable headers or in the appropriate machine-
+readable metadata fields within text or binary files as long as those
+fields can be easily viewed by the user.
+
+2) The font name complies with the following:
+(a) The Original Version must retain its name, unmodified.
+(b) Modified Versions which are Substantially Changed must be renamed to
+avoid use of the name of the Original Version or similar names entirely.
+(c) Modified Versions which are not Substantially Changed must be
+renamed to both (i) retain the name of the Original Version and (ii) add
+additional naming elements to distinguish the Modified Version from the
+Original Version. The name of such Modified Versions must be the name of
+the Original Version, with "derivative X" where X represents the name of
+the new work, appended to that name.
+
+3) The name(s) of the Copyright Holder(s) and any contributor to the
+Font Software shall not be used to promote, endorse or advertise any
+Modified Version, except (i) as required by this licence, (ii) to
+acknowledge the contribution(s) of the Copyright Holder(s) or (iii) with
+their explicit written permission.
+
+4) The Font Software, modified or unmodified, in part or in whole, must
+be distributed entirely under this licence, and must not be distributed
+under any other licence. The requirement for fonts to remain under this
+licence does not affect any document created using the Font Software,
+except any version of the Font Software extracted from a document
+created using the Font Software may only be distributed under this
+licence.
+
+TERMINATION
+This licence becomes null and void if any of the above conditions are
+not met.
+
+DISCLAIMER
+THE FONT SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND,
+EXPRESS OR IMPLIED, INCLUDING BUT NOT LIMITED TO ANY WARRANTIES OF
+MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT OF
+COPYRIGHT, PATENT, TRADEMARK, OR OTHER RIGHT. IN NO EVENT SHALL THE
+COPYRIGHT HOLDER BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER LIABILITY,
+INCLUDING ANY GENERAL, SPECIAL, INDIRECT, INCIDENTAL, OR CONSEQUENTIAL
+DAMAGES, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+FROM, OUT OF THE USE OR INABILITY TO USE THE FONT SOFTWARE OR FROM OTHER
+DEALINGS IN THE FONT SOFTWARE.
diff --git a/src/remedi/fonts/Ubuntu-Bold.ttf b/src/remedi/fonts/Ubuntu-Bold.ttf
new file mode 100644
index 0000000000000000000000000000000000000000..747fecc4350959d29aa23aef31f5c2efa644da81
--- /dev/null
+++ b/src/remedi/fonts/Ubuntu-Bold.ttf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:679b5c1e09cab3156bb8ef529735f9382bf31ca7ac737382ab959297f8d82ad4
+size 331612
diff --git a/src/remedi/fonts/Ubuntu-Regular.ttf b/src/remedi/fonts/Ubuntu-Regular.ttf
new file mode 100644
index 0000000000000000000000000000000000000000..361486dccd2671c377d563ba6dfa179822228b52
--- /dev/null
+++ b/src/remedi/fonts/Ubuntu-Regular.ttf
@@ -0,0 +1,3 @@
+version https://git-lfs.github.com/spec/v1
+oid sha256:3128df86a31805618436d0ae5651ba4285d0c9de0a39057d025f64ee33bceb64
+size 351884
diff --git a/src/remedi/io.py b/src/remedi/io.py
new file mode 100644
index 0000000000000000000000000000000000000000..65ec5265df5508516166fe3ce1e12800c17069eb
--- /dev/null
+++ b/src/remedi/io.py
@@ -0,0 +1,87 @@
+"""Validated condition-level artifacts and reproducible metadata."""
+from dataclasses import dataclass
+from pathlib import Path
+import hashlib
+import json
+import numpy as np
+import pandas as pd
+
+REQUIRED = ["molecule_id", "smiles", "dose_um", "context", "block", "control_id", "split"]
+
+def write_json(path, value):
+    path = Path(path)
+    path.parent.mkdir(parents=True, exist_ok=True)
+    path.write_text(json.dumps(value, indent=2, allow_nan=False, default=str) + "\n")
+
+def file_hash(path):
+    h = hashlib.sha256()
+    with open(path, "rb") as handle:
+        for chunk in iter(lambda: handle.read(1024 * 1024), b""):
+            h.update(chunk)
+    return h.hexdigest()
+
+def fingerprint_config(value):
+    return hashlib.sha256(json.dumps(value, sort_keys=True).encode()).hexdigest()
+
+@dataclass
+class ResponseData:
+    obs: pd.DataFrame
+    response: np.ndarray
+    control: np.ndarray
+    treated_var: np.ndarray
+    control_var: np.ndarray
+    metadata: dict
+
+    def validate(self):
+        missing = set(REQUIRED) - set(self.obs.columns)
+        if missing:
+            raise ValueError(f"Missing condition fields: {sorted(missing)}")
+        n = len(self.obs)
+        if n == 0:
+            raise ValueError("No eligible treated conditions")
+        if self.response.ndim != 2 or self.response.shape[0] != n:
+            raise ValueError("Response matrix and condition table are misaligned")
+        for name in ["response", "control", "treated_var", "control_var"]:
+            a = getattr(self, name)
+            if a.shape != self.response.shape or not np.isfinite(a).all():
+                raise ValueError(f"Invalid {name} matrix")
+        if (self.treated_var < 0).any() or (self.control_var < 0).any():
+            raise ValueError("Sampling variances must be nonnegative")
+        if not np.isfinite(self.obs.dose_um).all() or (self.obs.dose_um <= 0).any():
+            raise ValueError("Treated concentrations must be positive, finite micromolar values")
+        if self.obs[REQUIRED].isna().any().any():
+            raise ValueError("Condition identifiers and splits cannot be missing")
+        if self.obs.groupby("molecule_id").split.nunique().max() > 1:
+            raise ValueError("A molecule appears in multiple splits")
+        if self.obs.groupby("smiles").split.nunique().max() > 1:
+            raise ValueError("An identical structure appears in multiple splits")
+        for cid, group in self.obs.groupby("control_id", sort=False):
+            idx = group.index.to_numpy()
+            if not np.allclose(self.control[idx], self.control[idx[0]]):
+                raise ValueError(f"Shared control {cid} has inconsistent summaries")
+        return self
+
+    def subset(self, mask):
+        idx = np.flatnonzero(np.asarray(mask))
+        return ResponseData(self.obs.iloc[idx].reset_index(drop=True), self.response[idx],
+                            self.control[idx], self.treated_var[idx], self.control_var[idx],
+                            dict(self.metadata)).validate()
+
+    def save(self, directory):
+        self.validate()
+        directory = Path(directory)
+        directory.mkdir(parents=True, exist_ok=True)
+        self.obs.to_csv(directory / "conditions.csv", index=False)
+        np.savez_compressed(directory / "responses.npz", response=self.response,
+                            control=self.control, treated_var=self.treated_var,
+                            control_var=self.control_var)
+        write_json(directory / "metadata.json", self.metadata)
+
+    @classmethod
+    def load(cls, directory):
+        directory = Path(directory)
+        obs = pd.read_csv(directory / "conditions.csv", dtype={k: str for k in REQUIRED if k != "dose_um"})
+        with np.load(directory / "responses.npz", allow_pickle=False) as a:
+            values = {k: np.array(a[k], dtype=float) for k in ["response", "control", "treated_var", "control_var"]}
+        return cls(obs, **values, metadata=json.loads((directory / "metadata.json").read_text())).validate()
+
diff --git a/src/remedi/metrics.py b/src/remedi/metrics.py
new file mode 100644
index 0000000000000000000000000000000000000000..6165a66f84fb34a6ab1b32a9cd2a07afdb1a4028
--- /dev/null
+++ b/src/remedi/metrics.py
@@ -0,0 +1,46 @@
+"""Measured endpoint metrics and grouped uncertainty summaries."""
+import numpy as np
+import pandas as pd
+
+def intervention_metrics(probabilities, observed_losses, tolerance, seed=0, draws=10, repeats=200):
+    p, loss = np.asarray(probabilities, float), np.asarray(observed_losses, float)
+    if p.shape != loss.shape or (p < 0).any() or not np.isclose(p.sum(), 1):
+        raise ValueError("A normalized action distribution and aligned losses are required")
+    if not np.isfinite(loss).all(): raise ValueError("Observed losses must be finite")
+    selected = np.random.default_rng(seed).choice(len(p), size=(repeats, draws), p=p)
+    expected = float(p@loss)
+    return {"endpoint_loss": expected, "regret": expected-float(loss.min()),
+            "success": float(p@(loss <= tolerance)), "best_of_10": float(loss[selected].min(1).mean()),
+            "entropy": float(-np.sum(p*np.log(np.maximum(p, 1e-300)))),
+            "effective_actions": float(np.exp(-np.sum(p*np.log(np.maximum(p, 1e-300)))))}
+
+def grouped_summary(frame, metric="endpoint_loss", group="target_molecule", bootstrap=2000, seed=0):
+    rng = np.random.default_rng(seed)
+    values = frame.groupby([group, "method"])[metric].mean().unstack()
+    if values.isna().any().any():
+        raise ValueError("Methods must share the same evaluated molecule groups")
+    n = len(values)
+    if n == 0: raise ValueError("No evaluation groups")
+    draws = values.to_numpy()[rng.integers(n, size=(bootstrap, n))].mean(1)
+    rows = []
+    for j, name in enumerate(values.columns):
+        lo, hi = np.quantile(draws[:, j], [.025, .975])
+        rows.append({"method": name, "metric": metric, "mean": float(values.iloc[:, j].mean()),
+                     "ci_low": float(lo), "ci_high": float(hi), "molecules": n})
+    return pd.DataFrame(rows)
+
+def paired_difference(frame, first, second, metric="endpoint_loss", seed=0, bootstrap=2000):
+    pivot = frame.groupby(["target_molecule", "method"])[metric].mean().unstack()
+    diff = (pivot[first]-pivot[second]).dropna().to_numpy()
+    if not len(diff): raise ValueError("No paired molecule groups")
+    rng = np.random.default_rng(seed)
+    values = diff[rng.integers(len(diff), size=(bootstrap, len(diff)))].mean(1)
+    return {"first": first, "second": second, "metric": metric, "difference": float(diff.mean()),
+            "ci_low": float(np.quantile(values, .025)), "ci_high": float(np.quantile(values, .975))}
+
+def rbf_mmd(x, y, bandwidth):
+    from scipy.spatial.distance import cdist
+    if bandwidth <= 0 or len(x) < 2 or len(y) < 2: raise ValueError("MMD needs positive bandwidth and two cells per group")
+    xx, yy, xy = [np.exp(-d/(2*bandwidth**2)) for d in [cdist(x, x, "sqeuclidean"), cdist(y, y, "sqeuclidean"), cdist(x, y, "sqeuclidean")]]
+    return float((xx.sum()-np.trace(xx))/(len(x)*(len(x)-1))+(yy.sum()-np.trace(yy))/(len(y)*(len(y)-1))-2*xy.mean())
+
diff --git a/src/remedi/models.py b/src/remedi/models.py
new file mode 100644
index 0000000000000000000000000000000000000000..8361c2608ec9475ada6ddf17357288473c1e3f9d
--- /dev/null
+++ b/src/remedi/models.py
@@ -0,0 +1,107 @@
+"""Small molecule-dose response heads, grouped ensembles, and cross-fitted residuals."""
+from pathlib import Path
+import joblib
+import numpy as np
+from sklearn.linear_model import Ridge
+from sklearn.neural_network import MLPRegressor
+from sklearn.preprocessing import StandardScaler
+from sklearn.model_selection import GroupKFold
+from .chemistry import encode, fingerprints, max_tanimoto
+from .io import write_json, file_hash
+
+def design_features(chemical, dose, context, mode="ridge"):
+    """Build condition rows from molecular features, dose in uM, and control means.
+
+    chemical is (N, P), dose is (N,), and context is (N, D). Ridge
+    features include molecule-specific linear and quadratic log-dose terms.
+    """
+    dose = np.asarray(dose, float)
+    if (dose <= 0).any() or not np.isfinite(dose).all():
+        raise ValueError("Concentrations must be positive micromolar values")
+    logdose = np.log(dose)[:, None]
+    parts = [chemical, logdose, context]
+    if mode == "ridge":
+        parts += [chemical*logdose, chemical*logdose**2]
+    return np.concatenate(parts, axis=1)
+
+def _fit_head(x, y, weights, kind, alpha, hidden, seed):
+    scale = StandardScaler().fit(x, sample_weight=weights)
+    x = scale.transform(x)
+    if kind == "ridge":
+        head = Ridge(alpha=alpha).fit(x, y, sample_weight=weights)
+    elif kind == "mlp":
+        head = MLPRegressor(hidden_layer_sizes=(hidden, hidden), alpha=alpha,
+                            max_iter=500, early_stopping=False, random_state=seed)
+        head.fit(x, y, sample_weight=weights)
+    else: raise ValueError("head must be ridge or mlp")
+    return scale, head
+
+def _predict_head(head, x):
+    return np.asarray(head[1].predict(head[0].transform(x)))
+
+def train(data, output, kind="ridge", alpha=10., hidden=128, ensemble=5, folds=3,
+           encoder="fingerprint", embeddings=None, seed=0):
+    """Fit response heads and molecule-held-out residuals.
+
+    data.response contains treated-minus-control means, one row per condition.
+    Return the fitted model dictionary and save it with training metadata.
+    """
+    if ensemble < 1: raise ValueError("Ensemble must contain at least one head")
+    mask = data.obs.split.to_numpy() == "train"
+    if not mask.any(): raise ValueError("No training conditions")
+    train_obs = data.obs.loc[mask].reset_index(drop=True)
+    chemical = encode(train_obs.smiles, encoder, embeddings)
+    x = design_features(chemical, train_obs.dose_um, data.control[mask], kind)
+    y = data.response[mask]
+    groups = train_obs.molecule_id.to_numpy()
+    molecules, counts = np.unique(groups, return_counts=True)
+    if len(molecules) < max(folds, 3): raise ValueError("Too few independent training molecules")
+    # Give each training molecule equal total weight across its conditions.
+    base_weights = np.array([1/counts[np.where(molecules == m)[0][0]] for m in groups])
+    base_weights *= len(groups)/base_weights.sum()
+    rng = np.random.default_rng(seed)
+    heads = []
+    for k in range(ensemble):
+        # Resample whole molecules so their doses and replicates stay together.
+        sampled = rng.choice(molecules, len(molecules), replace=True)
+        multiplicity = {m: np.sum(sampled == m) for m in molecules}
+        weights = base_weights*np.array([multiplicity[m] for m in groups])
+        heads.append(_fit_head(x, y, weights, kind, alpha, hidden, seed+k))
+    residuals = np.empty_like(y)
+    residual_distance = np.empty(len(y))
+    row_fp = fingerprints(train_obs.smiles)
+    # Estimate residuals on compounds excluded from each fitted fold.
+    for k, (fit, held) in enumerate(GroupKFold(n_splits=folds).split(x, y, groups)):
+        head = _fit_head(x[fit], y[fit], base_weights[fit], kind, alpha, hidden, seed+100+k)
+        residuals[held] = y[held]-_predict_head(head, x[held])
+        residual_distance[held] = 1-max_tanimoto(row_fp[held], np.unique(row_fp[fit], axis=0))
+    # Set endpoint-coordinate weights using training response variance.
+    weight = 1/np.maximum(np.var(y, axis=0), np.median(np.var(y, axis=0))*.05+1e-6)
+    weight /= weight.sum()
+    unique_smiles = train_obs.drop_duplicates("molecule_id").smiles.tolist()
+    train_fp = fingerprints(unique_smiles)
+    # Cross-fitted residual rows and identities remain available for stratification.
+    model = {"heads": heads, "kind": kind, "alpha": alpha, "hidden": hidden,
+             "encoder": encoder, "embeddings": str(embeddings) if embeddings else None,
+             "residuals": residuals, "residual_distance": residual_distance,
+             "residual_smiles": train_obs.smiles.to_numpy(),
+             "train_molecules": molecules.tolist(), "train_smiles": unique_smiles,
+             "train_fingerprints": train_fp, "weight": weight,
+             "feature_space_id": data.metadata.get("feature_space_id"), "study": data.metadata.get("study"), "seed": seed}
+    output = Path(output); output.mkdir(parents=True, exist_ok=True)
+    joblib.dump(model, output/"model.joblib")
+    write_json(output/"training.json", {"kind": kind, "alpha": alpha, "hidden": hidden,
+               "encoder": encoder, "ensemble": ensemble, "crossfit_folds": folds,
+               "training_molecules": len(molecules), "training_conditions": len(y), "seed": seed,
+               "feature_space_id": model["feature_space_id"], "model_sha256": file_hash(output/"model.joblib")})
+    return model
+
+def predict(model, smiles, dose, context, embeddings=None):
+    """Return response predictions with shape (ensemble heads, actions, features)."""
+    chemical = encode(smiles, model["encoder"], embeddings or model["embeddings"])
+    x = design_features(chemical, dose, np.asarray(context, float), model["kind"])
+    return np.stack([_predict_head(h, x) for h in model["heads"]])
+
+def check_feature_space(model, data):
+    if model.get("feature_space_id") != data.metadata.get("feature_space_id"):
+        raise ValueError("Model and response data use different cell-feature spaces")
diff --git a/src/remedi/optimization.py b/src/remedi/optimization.py
new file mode 100644
index 0000000000000000000000000000000000000000..f1e0d15dbe9d69f252e430307c6bc113a9e5a8db
--- /dev/null
+++ b/src/remedi/optimization.py
@@ -0,0 +1,141 @@
+"""Convex finite-action minimax regret and matched generation rules."""
+from dataclasses import dataclass
+import time
+import numpy as np
+from scipy.special import logsumexp
+
+@dataclass
+class Decision:
+    probabilities: np.ndarray
+    objective: float
+    worst_regret: float
+    constraint_violation: float
+    seconds: float
+    status: str
+
+def validate_problem(mean, factor, diagonal, beta, tau, prior=None):
+    mean = np.asarray(mean, dtype=float)
+    factor = np.asarray(factor, dtype=float)
+    diagonal = np.asarray(diagonal, dtype=float)
+    n = len(mean)
+    if n == 0 or mean.shape != (n,) or factor.ndim != 2 or factor.shape[0] != n or diagonal.shape != (n,):
+        raise ValueError("Incompatible finite-action dimensions")
+    if not all(np.isfinite(x).all() for x in [mean, factor, diagonal]):
+        raise ValueError("Nonfinite loss or uncertainty values")
+    if (diagonal < 0).any() or not np.isfinite(beta) or not np.isfinite(tau) or beta < 0 or tau <= 0:
+        raise ValueError("Need nonnegative uncertainty and positive temperature")
+    prior = np.ones(n)/n if prior is None else np.asarray(prior, dtype=float)
+    if prior.shape != (n,) or not np.isfinite(prior).all() or (prior <= 0).any():
+        raise ValueError("Every eligible action needs a positive finite prior")
+    return mean, factor, diagonal, prior/prior.sum()
+
+def regret_components(p, mean, factor, diagonal, beta):
+    """All comparator constraints in O(N*S), including diagonal cross-terms."""
+    q = factor.T @ p
+    variance = np.sum((factor-q[None, :])**2, axis=1)
+    variance += np.dot(diagonal, p*p)-2*diagonal*p+diagonal
+    return np.dot(p, mean)-mean+beta*np.sqrt(np.maximum(variance, 0))
+
+def softmin(values, tau, prior=None):
+    values = np.asarray(values, float)
+    if tau <= 0 or not np.isfinite(values).all():
+        raise ValueError("Softmin needs finite values and positive temperature")
+    prior = np.ones(len(values))/len(values) if prior is None else np.asarray(prior, float)
+    if (prior <= 0).any():
+        raise ValueError("Prior must have positive support")
+    logits = np.log(prior)-values/tau
+    return np.exp(logits-logsumexp(logits))
+
+def solve_regret(mean, factor, diagonal, beta=1., tau=.05, prior=None,
+                 tolerance=1e-6, constraint_generation=True, max_rounds=100, solver="CLARABEL"):
+    """Optimize a distribution over N candidate molecule-dose pairs.
+
+    mean is (N,), factor is (N, S), and diagonal is (N,). Their covariance
+    is factor @ factor.T + diag(diagonal). beta sets the uncertainty radius
+    and tau sets KL regularization. Return probabilities and solver diagnostics.
+    """
+    import cvxpy as cp
+    start = time.perf_counter()
+    mean, factor, diagonal, prior = validate_problem(mean, factor, diagonal, beta, tau, prior)
+    n = len(mean)
+    # At zero uncertainty, the optimizer is the prior-weighted Gibbs distribution.
+    if beta == 0:
+        p = softmin(mean, tau, prior)
+        regret = float(p@mean-mean.min())
+        objective = regret+tau*float(np.sum(p*np.log(np.maximum(p, 1e-300)/prior)))
+        return Decision(p, objective, regret, 0., time.perf_counter()-start, "analytic")
+    active = set(range(n)) if not constraint_generation else {int(np.argmin(mean))}
+    probability, epigraph = cp.Variable(n), cp.Variable()
+    p = None
+    for _ in range(max_rounds):
+        constraints = [probability >= 0, cp.sum(probability) == 1]
+        for b in sorted(active):
+            e = np.zeros(n); e[b] = 1
+            # Compare the randomized intervention with the single action b.
+            contrast = probability-e
+            norm = cp.norm(cp.hstack([factor.T@contrast, cp.multiply(np.sqrt(diagonal), contrast)]), 2)
+            constraints.append(epigraph >= contrast@mean+beta*norm)
+        problem = cp.Problem(cp.Minimize(epigraph+tau*cp.sum(cp.kl_div(probability, prior))), constraints)
+        options = {"tol_gap_abs": tolerance*.1, "tol_feas": tolerance*.1} if solver == "CLARABEL" else {}
+        problem.solve(solver=solver, warm_start=True, **options)
+        if problem.status not in {"optimal", "optimal_inaccurate"} or probability.value is None:
+            raise RuntimeError(f"Cone solver failed: {problem.status}")
+        p = np.maximum(np.asarray(probability.value).ravel(), 0)
+        p /= p.sum()
+        # Check every comparator, including those outside the active constraint set.
+        values = regret_components(p, mean, factor, diagonal, beta)
+        worst = int(np.argmax(values))
+        violation = max(0., float(values[worst]-epigraph.value))
+        if violation <= tolerance:
+            status = problem.status
+            break
+        # Add the most violated comparator and solve the expanded program.
+        active.add(worst)
+    else:
+        status = "iteration_limit"
+    robust_regret = float(np.max(regret_components(p, mean, factor, diagonal, beta)))
+    kl = float(np.sum(p*np.log(np.maximum(p, 1e-300)/prior)))
+    return Decision(p, robust_regret+tau*kl, robust_regret, violation,
+                    time.perf_counter()-start, status)
+
+def finite_minimax(losses, tau=.05, prior=None, solver="CLARABEL"):
+    import cvxpy as cp
+    losses = np.asarray(losses, float)
+    if losses.ndim != 2 or not np.isfinite(losses).all():
+        raise ValueError("Scenarios must be a finite S-by-N matrix")
+    n = losses.shape[1]
+    prior = np.ones(n)/n if prior is None else np.asarray(prior, float)
+    regrets = losses-losses.min(axis=1, keepdims=True)
+    p, z = cp.Variable(n), cp.Variable()
+    problem = cp.Problem(cp.Minimize(z+tau*cp.sum(cp.kl_div(p, prior))),
+                         [p >= 0, cp.sum(p) == 1, regrets@p <= z])
+    problem.solve(solver=solver)
+    if p.value is None:
+        raise RuntimeError(f"Finite minimax failed: {problem.status}")
+    values = np.maximum(np.asarray(p.value), 0)
+    return values/values.sum()
+
+def generate_baseline(name, mean, factor, diagonal, losses, beta=1., tau=.05, seed=0):
+    n = len(mean)
+    if name == "uniform":
+        return np.ones(n)/n
+    if name == "mean":
+        p = np.zeros(n); p[int(np.argmin(mean))] = 1; return p
+    if name == "gibbs":
+        return softmin(mean, tau)
+    if name == "marginal":
+        return softmin(mean+beta*np.sqrt((factor**2).sum(1)+diagonal), tau)
+    if name == "thompson":
+        # Exact empirical distribution of scenario-wise minimizers.
+        return np.bincount(np.argmin(losses, axis=1), minlength=n)/len(losses)
+    if name == "finite-minimax":
+        return finite_minimax(losses, tau)
+    if name == "absolute":
+        import cvxpy as cp
+        p = cp.Variable(n)
+        norm = cp.norm(cp.hstack([factor.T@p, cp.multiply(np.sqrt(diagonal), p)]), 2)
+        objective = cp.Minimize(mean@p+beta*norm+tau*cp.sum(cp.kl_div(p, np.ones(n)/n)))
+        problem = cp.Problem(objective, [p >= 0, cp.sum(p) == 1]); problem.solve(solver="CLARABEL")
+        if p.value is None: raise RuntimeError("Absolute-loss solver failed")
+        values = np.maximum(p.value, 0); return values/values.sum()
+    raise ValueError(f"Unknown baseline {name}")
diff --git a/src/remedi/pipeline.py b/src/remedi/pipeline.py
new file mode 100644
index 0000000000000000000000000000000000000000..58d030822e383181fa64c52bc386f86dfa10e264
--- /dev/null
+++ b/src/remedi/pipeline.py
@@ -0,0 +1,239 @@
+"""Calibration, retrospective evaluation, and catalog generation with locked inputs."""
+from pathlib import Path
+import json
+import time
+import joblib
+import numpy as np
+import pandas as pd
+from .io import write_json, file_hash
+from .models import check_feature_space, predict
+from .uncertainty import target_draws, joint_losses, cross_block_losses, panel_score, calibration_radius
+from .optimization import solve_regret, generate_baseline
+from .metrics import intervention_metrics, grouped_summary, paired_difference
+from .chemistry import verify_catalog, fingerprints, max_tanimoto
+
+BASELINES = ["uniform", "mean", "gibbs", "marginal", "absolute", "thompson", "finite-minimax"]
+
+def action_groups(data, indices, min_blocks=1):
+    rows, groups = [], []
+    for _, g in data.obs.loc[indices].groupby(["molecule_id", "dose_um", "context"], sort=True):
+        g = g.sort_values("block").drop_duplicates("block")
+        if len(g) < min_blocks: continue
+        rows.append(g.iloc[0][["molecule_id", "smiles", "dose_um", "context"]].to_dict())
+        groups.append(g.index.to_list())
+    return pd.DataFrame(rows), groups
+
+def calibrate(data, model, output, coverage=.9, panel_molecules=6, scenarios=32, seed=0,
+               ablation="none", embeddings=None):
+    """Estimate beta from held-out compound panels with two measurement blocks.
+
+    Pool repeated contexts by their maximum panel score within each compound
+    group. Save the panel scores and the associated model identity.
+    """
+    check_feature_space(model, data)
+    obs = data.obs
+    cal = obs[obs.split == "calibration"]
+    molecules = np.sort(cal.molecule_id.unique())
+    if set(molecules) & set(model["train_molecules"]): raise ValueError("Training/calibration leakage")
+    if len(molecules) < 3: raise ValueError("Calibration needs at least three held-out molecules")
+    rng = np.random.default_rng(seed)
+    rng.shuffle(molecules)
+    records, skips = [], []
+    # A group is disjoint in molecule identity; repeated contexts contribute one maximum score.
+    for group_id, start in enumerate(range(0, len(molecules), panel_molecules)):
+        group = molecules[start:start+panel_molecules]
+        if len(group) < 3: continue
+        for context in sorted(cal.context.unique()):
+            panel = cal[(cal.molecule_id.isin(group)) & (cal.context == context)]
+            if panel.molecule_id.nunique() < 3: continue
+            target_molecule = panel.molecule_id.iloc[0]
+            target_rows = panel[panel.molecule_id == target_molecule]
+            # One dose per group/context prevents dependent target-dose panels being counted independently.
+            available_doses = np.sort(target_rows.dose_um.unique())
+            dose = available_doses[len(available_doses)//2]
+            target_rows = target_rows[target_rows.dose_um == dose].sort_values("block").drop_duplicates("block")
+            actions, groups = action_groups(data, panel[panel.molecule_id != target_molecule].index, 2)
+            if len(actions) < 2 or len(target_rows) < 2: continue
+            try:
+                observed = cross_block_losses(data, groups, target_rows.index.to_list(), model["weight"])
+            except ValueError as e:
+                skips.append({"group": group_id, "context": context, "reason": str(e)}); continue
+            targets, controls = target_draws(data, target_rows.index[:2], scenarios, rng)
+            joint = joint_losses(model, actions, targets, controls, seed+group_id,
+                                  embeddings, floor=ablation != "no-floor", ablation=ablation)
+            records.append({"group": group_id, "context": context, "score": panel_score(joint, observed),
+                             "actions": len(actions), "target_molecule": target_molecule})
+    if not records:
+        raise ValueError("No calibration panel has independent target/candidate blocks and distinct control pools")
+    frame = pd.DataFrame(records)
+    scores = frame.groupby("group").score.max().to_numpy()
+    beta = calibration_radius(scores, coverage)
+    output = Path(output); output.mkdir(parents=True, exist_ok=True)
+    frame.to_csv(output/"calibration_panels.csv", index=False)
+    result = {"beta": beta, "coverage": coverage, "independent_groups": len(scores),
+              "scores": scores.tolist(), "ablation": ablation, "seed": seed,
+              "calibration_molecules": molecules.tolist(), "feature_space_id": model["feature_space_id"],
+              "model_training_molecules": model["train_molecules"], "model_digest": joblib.hash(model),
+              "status": "empirical_panel_quantile",
+              "latent_ellipsoid_coverage_guaranteed": False, "skipped": skips}
+    write_json(output/"calibration.json", result)
+    return result
+
+def evaluate(data, model, calibration, output, split="test", scenarios=32, tau=.05,
+              tolerance=1., seed=0, max_queries=30, methods=None, embeddings=None,
+              metric_data=None, ablation="none", external_track="none"):
+    """Evaluate action distributions against separate measured response blocks.
+
+    Each query excludes its source molecule at every dose from the candidate
+    pool. Return per-query metrics and save probabilities and grouped summaries.
+    """
+    check_feature_space(model, data)
+    held = data.obs[data.obs.split == split]
+    protected = set(model["train_molecules"]) | set(calibration.get("calibration_molecules", []))
+    if external_track == "none" and set(held.molecule_id) & protected:
+        raise ValueError("Evaluation molecules overlap training or calibration")
+    if external_track != "none":
+        source = data.metadata.get("study")
+        if not source or not model.get("study") or source == model["study"]:
+            raise ValueError("External transfer requires distinct recorded study identifiers")
+        overlap = set(held.molecule_id) & set(model["train_molecules"])
+        if external_track == "unseen" and overlap:
+            raise ValueError("Unseen external track includes trained molecular structures")
+        if external_track == "shared" and set(held.molecule_id) != overlap:
+            raise ValueError("Shared external track includes structures absent from response-model training")
+    if calibration.get("model_digest") != joblib.hash(model):
+        raise ValueError("Calibration was produced by a different fitted model")
+    if calibration.get("feature_space_id") != model["feature_space_id"]:
+        raise ValueError("Calibration feature space differs from model")
+    methods = methods or BASELINES+["remedi"]
+    metric_data = data if metric_data is None else metric_data
+    keys = ["molecule_id", "dose_um", "context", "block", "control_id"]
+    if not data.obs[keys].equals(metric_data.obs[keys]): raise ValueError("Evaluation representation must preserve condition alignment")
+    if metric_data is data: weight = model["weight"]
+    else:
+        train = metric_data.response[metric_data.obs.split.to_numpy() == "train"]
+        if not len(train): raise ValueError("Independent evaluation space needs training-derived scaling")
+        weight = 1/np.maximum(train.var(0), 1e-6); weight /= weight.sum()
+    rng = np.random.default_rng(seed)
+    query_groups = list(held.groupby(["molecule_id", "dose_um", "context"], sort=True))
+    rng.shuffle(query_groups)
+    rows, choices, skipped = [], [], []
+    query_count = 0
+    for key, target in query_groups:
+        if query_count >= max_queries: break
+        target = target.sort_values("block").drop_duplicates("block")
+        if len(target) < 2: continue
+        query_idx, eval_idx = int(target.index[0]), int(target.index[1])
+        if target.iloc[0].control_id == target.iloc[1].control_id:
+            skipped.append({"target": key, "reason": "query/evaluation controls shared"}); continue
+        pool = held[(held.context == key[2]) & (held.molecule_id != key[0])]
+        pool = pool[pool.control_id != data.obs.iloc[query_idx].control_id]
+        actions, groups = action_groups(data, pool.index)
+        if len(actions) < 2: continue
+        # All candidate responses below remain evaluation-only; no observed value enters joint_losses.
+        targets, controls = target_draws(data, [query_idx], scenarios, rng,
+                                         sampling=ablation != "no-target-sampling")
+        joint = joint_losses(model, actions, targets, controls, seed+query_count,
+                              embeddings, floor=ablation != "no-floor", ablation=ablation)
+        observed = np.array([np.sum((metric_data.response[g].mean(0)-metric_data.response[eval_idx])**2*weight) for g in groups])
+        for name in methods:
+            begin = time.perf_counter()
+            diagnostics = {}
+            if name == "remedi":
+                decision = solve_regret(joint.mean, joint.factor, joint.diagonal,
+                                         calibration["beta"], tau)
+                p = decision.probabilities
+                diagnostics = {"solver_status": decision.status, "solver_violation": decision.constraint_violation,
+                               "estimated_worst_regret": decision.worst_regret}
+            else:
+                p = generate_baseline(name, joint.mean, joint.factor, joint.diagonal, joint.scenarios,
+                                       calibration["beta"], tau, seed+query_count)
+            scores = intervention_metrics(p, observed, tolerance, seed+query_count)
+            molecule_p = pd.Series(p).groupby(actions.molecule_id.to_numpy()).sum().to_numpy()
+            scores["unique_molecules_10"] = float(np.sum(1-(1-molecule_p)**10))
+            dose_regret = np.zeros(len(actions))
+            for _, g in actions.groupby("molecule_id"):
+                ii = g.index.to_numpy(); dose_regret[ii] = observed[ii]-observed[ii].min()
+            row = {"query_id": query_count, "target_molecule": key[0], "target_dose_um": key[1],
+                   "context": key[2], "method": name, "ablation": ablation, "candidates": len(actions),
+                   "dose_regret": float(p@dose_regret), "seconds": time.perf_counter()-begin,
+                   **scores, **diagnostics}
+            rows.append(row)
+            for i in np.flatnonzero(p > 1e-8):
+                choices.append({"query_id": query_count, "method": name, **actions.iloc[i].to_dict(),
+                                "probability": float(p[i]), "observed_loss": float(observed[i])})
+        query_count += 1
+    if not rows: raise ValueError("No independent query/evaluation pairs with eligible alternatives")
+    output = Path(output); output.mkdir(parents=True, exist_ok=True)
+    frame = pd.DataFrame(rows); frame.to_csv(output/"per_query.csv", index=False)
+    pd.DataFrame(choices).to_csv(output/"action_probabilities.csv", index=False)
+    summary = pd.concat([grouped_summary(frame, m, seed=seed) for m in ["endpoint_loss", "regret", "success", "best_of_10", "dose_regret", "entropy", "unique_molecules_10"]])
+    summary.to_csv(output/"summary.csv", index=False)
+    differences = [paired_difference(frame, "remedi", n, seed=seed) for n in methods if n != "remedi"] if "remedi" in methods else []
+    write_json(output/"evaluation.json", {"queries": query_count, "target_molecules": frame.target_molecule.nunique(),
+               "candidate_split": split, "external_track": external_track, "seed": seed, "tau": tau, "tolerance": tolerance,
+               "beta": calibration["beta"], "scenarios": scenarios, "paired_differences": differences,
+               "metric_space": metric_data.metadata.get("cell_features"), "skipped": skipped,
+               "biological_validation": "retrospective measured endpoints; prototype unless complete registered benchmark"})
+    return frame
+
+def generate_catalog(data, model, calibration, catalog, target_indices, output, doses=(.05,.5,5.),
+                      scenarios=32, tau=.05, shortlist=512, draws=20, seed=0, embeddings=None):
+    """Sample catalog IDs, structures, and concentrations for one target endpoint.
+
+    target_indices identify replicates of one condition in data.obs. Expand
+    catalog records over doses in uM, shortlist if needed, then optimize and
+    save the distribution together with the sampled molecule-dose pairs.
+    """
+    check_feature_space(model, data)
+    if calibration.get("model_digest") != joblib.hash(model):
+        raise ValueError("Calibration was produced by a different fitted model")
+    if not target_indices or min(target_indices) < 0 or max(target_indices) >= len(data.obs):
+        raise ValueError("Target rows must be valid condition indices")
+    query = data.obs.iloc[target_indices]
+    if len(query[["molecule_id", "dose_um", "context"]].drop_duplicates()) != 1:
+        raise ValueError("Target rows must describe one molecule, dose, and context")
+    if shortlist < 2 or draws < 1 or scenarios < 2:
+        raise ValueError("Require shortlist >= 2, draws >= 1, and scenarios >= 2")
+    catalog = verify_catalog(catalog)
+    if not len(catalog): raise ValueError("Verified catalog is empty")
+    doses = np.asarray(doses, float)
+    if not np.isfinite(doses).all() or (doses <= 0).any(): raise ValueError("Doses must be finite positive values")
+    actions = catalog.merge(pd.DataFrame({"dose_um": np.unique(doses)}), how="cross")
+    targets, controls = target_draws(data, target_indices, scenarios, np.random.default_rng(seed))
+    source_count = len(actions)
+    # Shared bounded shortlist. This stage is explicitly outside the finite-pool regret guarantee.
+    if len(actions) > shortlist:
+        mean_target, mean_control = targets.mean(0), controls.mean(0)
+        rough = []
+        for start in range(0, len(actions), 1024):
+            a = actions.iloc[start:start+1024]
+            f = predict(model, a.smiles, a.dose_um, np.repeat(mean_control[None, :], len(a), axis=0), embeddings).mean(0)
+            rough.extend(((f-mean_target)**2*model["weight"]).sum(1))
+        rank = np.argsort(rough, kind="stable")
+        low = rank[:max(1, shortlist*3//4)]
+        remaining = np.setdiff1d(np.arange(len(actions)), low)
+        extra = []
+        candidate_fp = fingerprints(actions.iloc[remaining].smiles)
+        similarity = max_tanimoto(candidate_fp, fingerprints(actions.iloc[low].smiles))
+        for _ in range(min(shortlist-len(low), len(remaining))):
+            j = int(np.argmin(similarity))
+            extra.append(remaining[j])
+            similarity = np.maximum(similarity, max_tanimoto(candidate_fp, candidate_fp[j:j+1]))
+            similarity[j] = np.inf
+        extra = np.asarray(extra, dtype=int)
+        actions = actions.iloc[np.concatenate([low, extra])].reset_index(drop=True)
+    joint = joint_losses(model, actions, targets, controls, seed, embeddings)
+    result = solve_regret(joint.mean, joint.factor, joint.diagonal, calibration["beta"], tau)
+    frame = actions.copy(); frame["probability"] = result.probabilities; frame["predicted_loss"] = joint.mean
+    frame["chemical_distance"] = 1-max_tanimoto(fingerprints(frame.smiles), model["train_fingerprints"])
+    selected = np.random.default_rng(seed).choice(len(frame), draws, p=result.probabilities)
+    output = Path(output); output.mkdir(parents=True, exist_ok=True)
+    frame.to_csv(output/"distribution.csv", index=False)
+    frame.iloc[selected].to_csv(output/"generated_pairs.csv", index=False)
+    write_json(output/"generation.json", {"source_actions": source_count, "shortlist_actions": len(frame),
+               "catalog_membership": "exact canonical structure match to supplied verified records",
+               "seconds": result.seconds, "solver_status": result.status, "constraint_violation": result.constraint_violation,
+               "estimated_worst_regret": result.worst_regret, "biological_effect_measured": False,
+               "dose_um": doses.tolist(), "seed": seed, "tau": tau, "beta": calibration["beta"]})
+    return frame
diff --git a/src/remedi/plotting.py b/src/remedi/plotting.py
new file mode 100644
index 0000000000000000000000000000000000000000..11d9ba84bca6ff83235fc88b2ea7f954715dca4d
--- /dev/null
+++ b/src/remedi/plotting.py
@@ -0,0 +1,58 @@
+"""Ubuntu plots from explicit evaluation records; never invent missing measurements."""
+from pathlib import Path
+import numpy as np
+import pandas as pd
+import matplotlib
+matplotlib.use("Agg")
+import matplotlib.pyplot as plt
+from matplotlib import font_manager
+
+GOLD = "#C68B25"
+INK = "#302B25"
+
+def configure_fonts(font_dir=None):
+    if font_dir is None:
+        font_dir = Path(__file__).resolve().parent/"fonts"
+    for path in Path(font_dir).glob("*.ttf"):
+        font_manager.fontManager.addfont(str(path))
+    available = {f.name for f in font_manager.fontManager.ttflist}
+    if "Ubuntu" not in available:
+        raise ValueError("Ubuntu is unavailable; provide --font-dir pointing to the included fonts")
+    plt.rcParams.update({"font.family": "Ubuntu", "font.size": 10, "mathtext.fontset": "cm",
+                         "pdf.fonttype": 42, "ps.fonttype": 42, "axes.spines.top": False,
+                         "axes.spines.right": False, "axes.labelcolor": INK, "text.color": INK,
+                         "savefig.bbox": "tight", "figure.dpi": 160})
+
+def plot_summary(summary, output, font_dir=None, title="Retrospective endpoint evaluation"):
+    configure_fonts(font_dir)
+    frame = pd.read_csv(summary) if not isinstance(summary, pd.DataFrame) else summary
+    output = Path(output); output.mkdir(parents=True, exist_ok=True)
+    for metric in frame.metric.unique():
+        group = frame[frame.metric == metric].sort_values("mean", ascending=metric != "success")
+        fig, ax = plt.subplots(figsize=(6.6, max(2.8, .38*len(group))))
+        y = np.arange(len(group))
+        colors = [GOLD if m == "remedi" else "#8C8982" for m in group.method]
+        ax.barh(y, group["mean"], color=colors, height=.65)
+        ax.errorbar(group["mean"], y, xerr=np.vstack([np.maximum(0, group["mean"]-group.ci_low), np.maximum(0, group.ci_high-group["mean"])]),
+                    fmt="none", ecolor=INK, capsize=2, linewidth=.8)
+        ax.set_yticks(y, group.method); ax.invert_yaxis()
+        ax.set_xlabel(metric.replace("_", " ")); ax.set_title(title, loc="left", fontsize=11)
+        fig.savefig(output/f"{metric}.pdf"); fig.savefig(output/f"{metric}.png"); plt.close(fig)
+
+def plot_curve(csv, output, x, y, hue="method", lower=None, upper=None, font_dir=None):
+    configure_fonts(font_dir)
+    frame = pd.read_csv(csv)
+    for field in [x, y, hue]:
+        if field not in frame: raise ValueError(f"Missing plot column {field}")
+    fig, ax = plt.subplots(figsize=(5.2, 3.4))
+    palette = [GOLD, "#74716B", "#9E4A35", "#457679", "#756589"]
+    for i, (name, group) in enumerate(frame.groupby(hue, sort=True)):
+        group = group.sort_values(x)
+        ax.plot(group[x], group[y], marker="o", markersize=3.5, label=name,
+                color=GOLD if name == "remedi" else palette[(i+1) % len(palette)])
+        if lower and upper:
+            ax.fill_between(group[x], group[lower], group[upper], alpha=.13)
+    ax.set_xlabel(x.replace("_", " ")); ax.set_ylabel(y.replace("_", " ")); ax.legend(frameon=False)
+    output = Path(output); output.parent.mkdir(parents=True, exist_ok=True)
+    fig.savefig(output); plt.close(fig)
+
diff --git a/src/remedi/splits.py b/src/remedi/splits.py
new file mode 100644
index 0000000000000000000000000000000000000000..e963108f191418a94dbd12ee223141f80df93e34
--- /dev/null
+++ b/src/remedi/splits.py
@@ -0,0 +1,41 @@
+"""Molecule-disjoint or scaffold-disjoint splits before feature fitting."""
+import numpy as np
+import pandas as pd
+from .chemistry import canonical_smiles, molecule_id
+
+def make_splits(frame, seed=0, fractions=(.6, .15, .1, .15), scaffold=False):
+    if len(fractions) != 4 or not np.isclose(sum(fractions), 1) or min(fractions) <= 0:
+        raise ValueError("Four positive split fractions must sum to one")
+    frame = frame[["smiles"]].drop_duplicates().copy()
+    frame["smiles"] = frame.smiles.map(canonical_smiles)
+    frame = frame.drop_duplicates("smiles")
+    frame["molecule_id"] = frame.smiles.map(molecule_id)
+    frame["group"] = frame.molecule_id
+    if scaffold:
+        from rdkit import Chem
+        from rdkit.Chem.Scaffolds import MurckoScaffold
+        def key(s):
+            scaffold = MurckoScaffold.MurckoScaffoldSmiles(mol=Chem.MolFromSmiles(s))
+            # All acyclic structures remain one group in this conservative implementation.
+            return scaffold or "ACYCLIC"
+        frame["group"] = frame.smiles.map(key)
+    groups = frame.groupby("group").size()
+    if len(groups) < 8:
+        raise ValueError("At least eight independent molecule/scaffold groups are required")
+    rng = np.random.default_rng(seed)
+    order = groups.index.to_numpy()[rng.permutation(len(groups))]
+    # Allocate whole groups. Each split receives at least one group.
+    counts = np.maximum(1, np.floor(np.asarray(fractions) * len(order)).astype(int))
+    while counts.sum() < len(order):
+        counts[np.argmax(np.asarray(fractions) * len(order) - counts)] += 1
+    while counts.sum() > len(order):
+        eligible = np.where(counts > 1, counts - np.asarray(fractions)*len(order), -np.inf)
+        counts[np.argmax(eligible)] -= 1
+    assignments = {}
+    pos = 0
+    for name, n in zip(["train", "tune", "calibration", "test"], counts):
+        assignments.update({g: name for g in order[pos:pos+n]})
+        pos += n
+    frame["split"] = frame.group.map(assignments)
+    return frame.reset_index(drop=True)
+
diff --git a/src/remedi/uncertainty.py b/src/remedi/uncertainty.py
new file mode 100644
index 0000000000000000000000000000000000000000..84f845b04a4e4a64a075221b5da9751a19fabb18
--- /dev/null
+++ b/src/remedi/uncertainty.py
@@ -0,0 +1,126 @@
+"""Coherent target draws, joint loss covariance, and observed-comparison calibration."""
+from dataclasses import dataclass
+import numpy as np
+from .models import design_features, _predict_head
+from .chemistry import encode, fingerprints, max_tanimoto
+
+@dataclass
+class JointLoss:
+    mean: np.ndarray
+    factor: np.ndarray
+    diagonal: np.ndarray
+    scenarios: np.ndarray
+    predictions: np.ndarray
+
+def target_draws(data, indices, count, rng, sampling=True):
+    """Draw target-response and control means, each with shape (count, features).
+
+    Sample blocks with replacement. Reuse control noise within each control
+    pool so the response difference and context retain their shared variation.
+    """
+    indices = np.asarray(indices, dtype=int)
+    if not len(indices): raise ValueError("A query needs target observations")
+    draws, controls = [], []
+    for _ in range(count):
+        chosen = rng.choice(indices, len(indices), replace=True) if sampling else indices
+        control_noise = {}
+        treated_noise = {}
+        rr, cc = [], []
+        for i in chosen:
+            cid = data.obs.iloc[i].control_id
+            if cid not in control_noise:
+                control_noise[cid] = rng.normal(size=data.response.shape[1])*np.sqrt(data.control_var[i]) if sampling else np.zeros(data.response.shape[1])
+            cn = control_noise[cid]
+            if i not in treated_noise:
+                treated_noise[i] = rng.normal(size=data.response.shape[1])*np.sqrt(data.treated_var[i]) if sampling else 0
+            tn = treated_noise[i]
+            rr.append(data.response[i]+tn-cn)
+            cc.append(data.control[i]+cn)
+        draws.append(np.mean(rr, axis=0)); controls.append(np.mean(cc, axis=0))
+    return np.asarray(draws), np.asarray(controls)
+
+def joint_losses(model, actions, targets, controls, seed=0, embeddings=None, floor=True,
+                  ablation="none"):
+    """Estimate joint endpoint losses for a common set of candidate actions.
+
+    targets and controls are aligned (S, D) draws. Return loss scenarios
+    (S, N), their mean (N,), covariance factor (N, S), and residual floor (N,).
+    """
+    rng = np.random.default_rng(seed)
+    n, s = len(actions), len(targets)
+    if s < 2: raise ValueError("At least two scenarios are required")
+    if np.asarray(controls).shape != np.asarray(targets).shape:
+        raise ValueError("Target and control draws must align")
+    # A scenario uses one head AND one query realization across every action.
+    losses, predictions = [], []
+    chemical = encode(actions.smiles, model["encoder"], embeddings or model["embeddings"])
+    for i in range(s):
+        x = design_features(chemical, actions.dose_um,
+                            np.repeat(controls[i:i+1], n, axis=0), model["kind"])
+        f = _predict_head(model["heads"][int(rng.integers(len(model["heads"])))], x)
+        losses.append(((f-targets[i])**2*model["weight"]).sum(1))
+        predictions.append(f)
+    losses, predictions = np.asarray(losses), np.asarray(predictions)
+    if ablation == "permuted":
+        losses = np.stack([rng.permutation(losses[:, a]) for a in range(n)], axis=1)
+    mean = losses.mean(0)
+    # This factor reproduces the sample covariance without a dense N-by-N matrix.
+    factor = (losses-mean).T/np.sqrt(s-1)
+    diagonal = np.zeros(n)
+    if floor:
+        residual = model["residuals"]
+        # Empirical response residuals induce query-specific scalar loss errors.
+        f = predictions.mean(0); t = targets.mean(0)
+        distance = 1-max_tanimoto(fingerprints(actions.smiles), model["train_fingerprints"])
+        bins = np.asarray([.25, .5, .75])
+        residual_bins = np.digitize(model["residual_distance"], bins)
+        for a in range(n):
+            stratum = np.flatnonzero(residual_bins == np.digitize(distance[a], bins))
+            # Blend sparse chemical-distance strata with the pooled residual sample.
+            shrinkage = len(stratum)/(len(stratum)+20.)
+            indices = rng.choice(len(residual), size=min(128, len(residual)), replace=True)
+            use_stratum = rng.random(len(indices)) < shrinkage
+            if len(stratum): indices[use_stratum] = rng.choice(stratum, int(use_stratum.sum()))
+            sampled = residual[indices]
+            perturbed = (((f[a]+sampled-t)**2-(f[a]-t)**2)*model["weight"]).sum(1)
+            diagonal[a] = max(float(np.var(perturbed, ddof=1)), 1e-10)
+    if ablation == "diagonal":
+        diagonal += np.sum(factor**2, axis=1)
+        factor = np.zeros((n, 0))
+    return JointLoss(mean, factor, diagonal, losses, predictions.mean(0))
+
+def cross_block_losses(data, action_indices, target_indices, weight):
+    """Unbiased cross-product losses only for two distinct supplied blocks."""
+    target_indices = list(target_indices)
+    if len(target_indices) < 2: raise ValueError("Two target blocks are needed")
+    u, v = target_indices[:2]
+    if data.obs.iloc[u].block == data.obs.iloc[v].block:
+        raise ValueError("Target blocks are not distinct")
+    left_controls = {data.obs.iloc[u].control_id}; right_controls = {data.obs.iloc[v].control_id}
+    values = []
+    for indices in action_indices:
+        if len(indices) < 2: raise ValueError("Two blocks per candidate are needed")
+        a, b = indices[:2]
+        if data.obs.iloc[a].block == data.obs.iloc[b].block:
+            raise ValueError("Candidate blocks are not distinct")
+        left_controls.add(data.obs.iloc[a].control_id); right_controls.add(data.obs.iloc[b].control_id)
+        values.append(np.sum((data.response[a]-data.response[u])*(data.response[b]-data.response[v])*weight))
+    if left_controls & right_controls:
+        raise ValueError("Cross-block comparisons reuse a control pool across blocks")
+    return np.asarray(values)
+
+def panel_score(joint, observed_losses, epsilon=1e-8):
+    error = joint.mean-np.asarray(observed_losses)
+    cov_diag = np.sum(joint.factor**2, axis=1)+joint.diagonal
+    variance = cov_diag[:, None]+cov_diag[None, :]-2*joint.factor@joint.factor.T
+    np.fill_diagonal(variance, 0)
+    standardized = np.abs(error[:, None]-error[None, :])/(np.sqrt(np.maximum(variance, 0))+epsilon)
+    np.fill_diagonal(standardized, 0)
+    return float(standardized.max())
+
+def calibration_radius(scores, coverage=.9):
+    scores = np.asarray(scores, float)
+    if not 0 < coverage < 1 or not len(scores) or not np.isfinite(scores).all():
+        raise ValueError("Calibration needs finite grouped scores and coverage in (0,1)")
+    # Use one score per molecule group to estimate the held-out panel quantile.
+    return float(np.quantile(scores, coverage, method="higher"))
diff --git a/tests/test_data_and_metrics.py b/tests/test_data_and_metrics.py
new file mode 100644
index 0000000000000000000000000000000000000000..8eb5ee635b51bd0cadbc22b765dc54b662872361
--- /dev/null
+++ b/tests/test_data_and_metrics.py
@@ -0,0 +1,53 @@
+import numpy as np
+import pandas as pd
+import pytest
+from scipy import sparse
+from remedi.data import normalize_counts, standardize_obs
+from remedi.splits import make_splits
+from remedi.metrics import intervention_metrics
+from remedi.demo import synthetic_data
+from remedi.io import ResponseData
+from remedi.uncertainty import cross_block_losses, calibration_radius
+
+def test_sparse_and_dense_count_normalization_agree():
+    x=np.array([[0,2,3],[0,0,0],[1,0,2]],float)
+    np.testing.assert_allclose(normalize_counts(x),normalize_counts(sparse.csr_matrix(x)).toarray())
+
+def test_missing_dose_units_fail_closed():
+    obs=pd.DataFrame({"drug":["x"],"dose":[1],"context":["c"],"block":["r"],"plate":["p"]})
+    mapping={"drug":"drug","dose":"dose","context":"context","block":"block","control_group":"plate"}
+    with pytest.raises(ValueError,match="dose_unit"):
+        standardize_obs(obs,mapping)
+
+def test_response_roundtrip_and_split_leakage(tmp_path):
+    data=synthetic_data(molecules=24)
+    data.save(tmp_path)
+    restored=ResponseData.load(tmp_path)
+    np.testing.assert_allclose(data.response,restored.response)
+    restored.obs.loc[0,"split"]="new_split"
+    with pytest.raises(ValueError,match="multiple splits"):restored.validate()
+
+def test_molecule_groups_never_cross_splits():
+    smiles=["C"*k+"c1ccccc1" for k in range(1,20)]
+    result=make_splits(pd.DataFrame({"smiles":smiles+smiles}),seed=2)
+    assert len(result)==19
+    assert result.groupby("molecule_id").split.nunique().max()==1
+    assert set(result.split)=={"train","tune","calibration","test"}
+
+def test_mixture_of_bad_endpoints_is_not_scored_as_a_good_average():
+    # Responses -1 and +1 both have squared loss 1 from target 0.
+    result=intervention_metrics([.5,.5],[1.,1.],tolerance=.1)
+    assert result["endpoint_loss"]==1.
+    assert result["success"]==0.
+
+def test_cross_block_shared_control_is_rejected():
+    data=synthetic_data(molecules=24)
+    group=data.obs.groupby(["molecule_id","dose_um","context"])
+    indices=[g.index.to_list() for _,g in group][:2]
+    data.obs.loc[indices[0][1],"control_id"]=data.obs.loc[indices[0][0],"control_id"]
+    with pytest.raises(ValueError,match="reuse a control"):
+        cross_block_losses(data,[indices[1]],indices[0],np.ones(data.response.shape[1]))
+
+def test_calibration_is_an_explicit_empirical_quantile():
+    assert calibration_radius([1.,2.,3.,4.],.8)==4.
+    with pytest.raises(ValueError): calibration_radius([], .9)
diff --git a/tests/test_optimization.py b/tests/test_optimization.py
new file mode 100644
index 0000000000000000000000000000000000000000..2fd562a5341af31ebfbba4b001a93580251cc715
--- /dev/null
+++ b/tests/test_optimization.py
@@ -0,0 +1,48 @@
+import numpy as np
+from remedi.optimization import solve_regret, regret_components, softmin, finite_minimax
+
+def test_factorized_comparator_scan_matches_dense_covariance():
+    rng = np.random.default_rng(1)
+    b = rng.normal(size=(7,3)); d = rng.uniform(size=7)
+    p = rng.dirichlet(np.ones(7)); mean = rng.normal(size=7)
+    c = b@b.T+np.diag(d)
+    expected = []
+    for e in np.eye(7):
+        z = p-e
+        expected.append(z@mean+1.2*np.sqrt(z@c@z))
+    np.testing.assert_allclose(regret_components(p,mean,b,d,1.2), expected, atol=1e-10)
+
+def test_zero_radius_has_analytic_gibbs_solution():
+    mean = np.array([1.,2.,.5]); prior = np.array([.2,.3,.5])
+    result = solve_regret(mean,np.zeros((3,1)),np.ones(3),beta=0,tau=.2,prior=prior)
+    np.testing.assert_allclose(result.probabilities,softmin(mean,.2,prior))
+
+def test_common_loss_offsets_and_common_error_factors_cancel():
+    mean = np.array([.2,.3,.8]); b = np.array([[.1],[.2],[.4]]); d = np.full(3,.01)
+    a = solve_regret(mean,b,d,tau=.1)
+    c = solve_regret(mean+123,np.column_stack([b,np.full(3,7.)]),d,tau=.1)
+    np.testing.assert_allclose(a.probabilities,c.probabilities,atol=2e-4)
+
+def test_constraint_generation_matches_full_program():
+    rng = np.random.default_rng(8)
+    mean=rng.normal(size=10); b=rng.normal(size=(10,3))*.2; d=np.full(10,.02)
+    full=solve_regret(mean,b,d,tau=.2,constraint_generation=False)
+    active=solve_regret(mean,b,d,tau=.2)
+    np.testing.assert_allclose(full.probabilities,active.probabilities,atol=3e-4)
+    assert active.constraint_violation <= 1e-6
+
+def test_conditional_bound_contains_sampled_loss_vectors():
+    rng=np.random.default_rng(3); mean=rng.normal(size=5); b=rng.normal(size=(5,2)); d=np.full(5,.05)
+    result=solve_regret(mean,b,d,beta=1.5,tau=.1)
+    factor=np.column_stack([b,np.diag(np.sqrt(d))])
+    for _ in range(30):
+        u=rng.normal(size=7);u=u/np.linalg.norm(u)*rng.uniform(0,1.5)
+        loss=mean+factor@u
+        assert result.probabilities@loss-loss.min() <= result.worst_regret+1e-6
+
+def test_finite_scenario_minimax_is_offset_invariant():
+    losses=np.array([[0.,3.,2.],[3.,0.,2.]])
+    a=finite_minimax(losses,tau=.1)
+    b=finite_minimax(losses+np.array([[100.],[-20.]]),tau=.1)
+    np.testing.assert_allclose(a,b,atol=1e-5)
+
diff --git a/tests/test_preparation.py b/tests/test_preparation.py
new file mode 100644
index 0000000000000000000000000000000000000000..b66768dc6fa45d2741a7986f88118e0cbe674ea4
--- /dev/null
+++ b/tests/test_preparation.py
@@ -0,0 +1,39 @@
+import numpy as np
+import pandas as pd
+import anndata as ad
+import pytest
+from remedi.data import prepare_h5ad
+
+
+def fixture(path, order=None):
+    rng=np.random.default_rng(3)
+    obs=[]
+    for block in ['r1','r2']:
+        for drug,smiles in [('DMSO','CS(C)=O'),('a','CCO'),('b','CCN')]:
+            for cell in range(4):
+                obs.append(dict(drug=drug,smiles=smiles,dose=100.,context='c',block=block))
+    x=ad.AnnData(rng.poisson(3,size=(len(obs),10)).astype(float),obs=pd.DataFrame(obs,index=[str(i) for i in range(len(obs))]),var=pd.DataFrame(index=[f'g{i}' for i in range(10)]))
+    if order is not None:x=x[:,order].copy()
+    x.write_h5ad(path)
+
+
+def test_external_gene_reordering_is_invariant(tmp_path):
+    first=tmp_path/'a.h5ad';second=tmp_path/'b.h5ad'
+    fixture(first);fixture(second,np.arange(10)[::-1])
+    mapping=dict(drug='drug',dose='dose',context='context',block='block',control_group='block',dose_unit='nM',smiles='smiles',control_names=['DMSO'])
+    splits=pd.DataFrame({'smiles':['CCO','CCN'],'split':['train','test']})
+    a=prepare_h5ad(first,tmp_path/'fit',mapping,splits,min_cells=2,dimensions=3)
+    b=prepare_h5ad(second,tmp_path/'transfer',mapping,splits,min_cells=2,dimensions=3,feature_model=tmp_path/'fit/cell_feature_model.joblib')
+    np.testing.assert_allclose(a.response,b.response)
+    assert a.metadata['feature_space_id']==b.metadata['feature_space_id']
+    np.testing.assert_allclose(a.obs.dose_um,.1)
+
+
+def test_external_missing_genes_rejected(tmp_path):
+    first=tmp_path/'a.h5ad';second=tmp_path/'b.h5ad'
+    fixture(first);fixture(second,np.arange(9))
+    mapping=dict(drug='drug',dose='dose',context='context',block='block',control_group='block',dose_unit='nM',smiles='smiles',control_names=['DMSO'])
+    splits=pd.DataFrame({'smiles':['CCO','CCN'],'split':['train','test']})
+    prepare_h5ad(first,tmp_path/'fit',mapping,splits,min_cells=2,dimensions=3)
+    with pytest.raises(ValueError,match='misses fitted genes'):
+        prepare_h5ad(second,tmp_path/'transfer',mapping,splits,min_cells=2,feature_model=tmp_path/'fit/cell_feature_model.joblib')