improved growth model fitting

CLAUDE.md

@@ -131,6 +131,69 @@ The hierarchical Bayesian inference engine. Key files:
 
 See **YAML Standards** below for the full conventions. The tfmodel config (`tfs_configure_config.yaml`) is generated by `tfs-configure-model` and drives `tfs-fit-model`. It contains `data`, `components`, `priors_file`, and `guesses_file` sections. Do not hand-edit it.
 
+## Simulation Internals (`simulate/`)
+
+### Overview
+
+The simulation pipeline generates ground-truth phenotypes for a synthetic TF screen experiment. The top-level entry point is `library_prediction` (`simulate/library_prediction.py`), which returns five dataframes:
+
+| Return value | Content |
+|---|---|
+| `library_df` | One row per genotype in the library |
+| `phenotype_df` | Long-form growth predictions (one row per genotype × condition) |
+| `genotype_theta_df` | Long-form theta predictions (one row per genotype × titrant_conc) |
+| `parameters_df` | One row per genotype; per-genotype Hill/theta params + dk_geno + activity |
+| `binding_theta_df` | Theta at binding concentrations for calibration genotypes (`None` if not configured) |
+
+The core calculation happens in `thermo_to_growth` (`simulate/thermo_to_growth.py`):
+
+1. **Prior-predictive theta sampling** — `sample_theta_prior` draws a `theta_gc` matrix of shape `(G, C)` where G = number of library genotypes and C = number of unique titrant concentrations. The genotype order here matches the library (`sim_data`) order.
+2. **theta_gc_override injection** — specific rows of `theta_gc` can be replaced before any further computation (see below).
+3. **Growth rate calculation** — theta is mapped to growth rates via the configured `growth_params`.
+4. **parameters_df assembly** — per-genotype Hill parameters extracted from the `theta_param` pytree, then patched with `theta_params_override` (see below).
+
+### Binding data and calibration genotypes
+
+The optional `binding_data` YAML block configures calibration genotypes for which measured binding curves are available. It has two sub-paths that can coexist:
+
+**Path 1: Simulated genotypes** (`binding_data.genotypes` list)
+- `wt` gets its natural unperturbed prior-predictive reference curve.
+- Non-wt entries are drawn via a stratified (greedy maximin) algorithm across binding concentrations, ensuring diverse coverage.
+- The selected theta values at *growth* concentrations are injected via `theta_gc_override` so that the simulated growth data is consistent with the binding data.
+
+**Path 2: Measured Hill parameters** (`binding_data.genotype_params_file`)
+- Reads a CSV with columns `genotype, theta_low, theta_high, log_hill_K, hill_n`.
+- Only supported for Hill-based theta components (`hill_geno`, `hill_mut`).
+- **`theta_low` and `theta_high` are clamped to `[1e-4, 1-1e-4]` at read time** with a `UserWarning`. This is critical: a value of e.g. `1.000004` (a common float-rounding artefact) maps to `logit ≈ +16`, making all per-mutation deltas ~13–15 σ under the `HalfNormal(1)` prior on delta scales and preventing inference from recovering reasonable theta values.
+- For `hill_mut`: the function `build_theta_gc_override_hill_mut` assembles theta for **all library genotypes** (not just the measured ones) by additively combining per-mutation logit-space deltas. The WT reference is taken from `SimPriors` defaults. Multi-mutant genotypes not directly measured in the CSV are assembled from single-mutant deltas; directly-measured multi-mutants use their CSV values directly.
+- Measured genotypes override any earlier simulated-path values in `theta_gc_override`.
+
+### theta_gc_override and theta_params_override
+
+These two dicts are the mechanism by which binding data is "pinned" into the growth simulation.
+
+**`theta_gc_override`** (`dict[str, np.ndarray]`):
+- Keys are genotype strings; values are 1-D arrays of theta at the growth titrant concentrations (sorted-unique order from `sample_df`).
+- Applied in `thermo_to_growth` *after* prior-predictive sampling but *before* noise and all downstream computations.
+- Overwrites the corresponding row of `theta_gc` in-place using `geno_to_sim_idx` (which maps genotype → its index in the original library list).
+
+**`theta_params_override`** (`dict[str, dict[str, float]]`):
+- Keys are genotype strings; values are dicts with Hill parameter keys (`theta_low`, `theta_high`, `log_hill_K`, `hill_n`).
+- Applied in `thermo_to_growth` *after* `parameters_df` is assembled from the `theta_param` pytree, overwriting the relevant columns.
+- Purpose: the `theta_param` from `sample_theta_prior` reflects the prior-predictive draw, not the override. Without this patch, `tfs_sim_parameters.csv` would show the pre-override values, making the saved parameters inconsistent with what was actually simulated.
+- Only columns already present in `parameters_df` are updated; unknown keys are silently skipped.
+- Genotype keys not found in `parameters_df` are silently skipped.
+
+### Key files
+
+| File | Role |
+|---|---|
+| `simulate/library_prediction.py` | Top-level orchestrator; assembles override dicts and calls `thermo_to_growth` |
+| `simulate/thermo_to_growth.py` | Prior-predictive sampling, theta injection, growth calculation, parameters_df assembly |
+| `simulate/binding_params.py` | CSV reading (with theta clipping), per-component override builders, binding theta assemblers |
+| `simulate/sample_theta.py` | `sample_theta_prior` (prior-predictive) and `sample_theta_stratified` (greedy maximin) |
+| `simulate/sim_data_class.py` | `SimData` container and `build_sim_data` factory |
+
 ## YAML Standards
 
 All YAML files in this codebase follow these conventions. Apply them when creating or modifying any YAML file.

docs/source/analysis.rst

@@ -53,8 +53,15 @@ conditions present in both the growth and binding data.
 
 After convergence the script updates the production ``{out_prefix}_priors.csv``
 and ``{out_prefix}_guesses.csv`` files in place (a ``.bak`` backup is written
-first), giving the full production fit a warm start. Diagnostic PDFs (one per
-genotype) and a ``{out_prefix}_calib_stats.json`` file are also written.
+first), giving the full production fit a warm start.
+
+Diagnostic artefacts from the calibration MAP run are also written (default
+prefix ``tfs_prefit``):
+
+* ``tfs_prefit_params.npz`` — MAP point estimates in constrained space.
+* ``tfs_prefit_checkpoint.pkl`` — optimizer checkpoint (can be passed to
+  ``--checkpoint_file`` to resume an interrupted calibration run).
+* ``tfs_prefit_losses.txt`` — per-epoch loss history.
 
 .. code-block:: bash
 
@@ -139,9 +146,32 @@ writes one CSV per parameter group.
 
 Outputs (default ``--out_prefix tfs_params``):
 
-* ``tfs_params_activity.csv`` — per-genotype TF activity *A* with posterior quantiles
-* ``tfs_params_theta.csv`` — inferred occupancy parameters (Hill *Kd*, *n*, etc.)
-* Additional CSVs depending on the components selected during configuration.
+* ``tfs_params_log_hill_K.csv`` — per-genotype log₁₀(*Kd*) with posterior quantiles.
+* ``tfs_params_theta_low.csv``, ``tfs_params_theta_high.csv`` — per-genotype
+  lower and upper occupancy plateaux.
+* ``tfs_params_dk_geno.csv`` — per-genotype pleiotropic growth effect.
+* Additional CSVs depending on the components selected during configuration
+  (e.g. ``tfs_params_d_log_hill_K.csv`` for per-mutation effects when using
+  ``hill_mut``; ``tfs_params_epi_*.csv`` for epistasis terms when
+  ``--epistasis`` is set).
+
+**Quantile columns**
+
+Each CSV contains metadata columns (typically ``genotype`` and
+``titrant_name``) followed by 17 quantile columns named ``q{level}``:
+
+.. code-block:: text
+
+    q0.001  q0.005  q0.01  q0.025  q0.05  q0.1  q0.159  q0.25
+    q0.5    q0.75   q0.841 q0.9    q0.95  q0.975 q0.99  q0.995  q0.999
+
+``q0.5`` is the posterior median (the recommended point estimate).  ``q0.025``
+and ``q0.975`` span the 95% credible interval; ``q0.159`` and ``q0.841``
+span the ±1 σ interval under a normal approximation.
+
+When ``tfs-extract-params`` is given a MAP checkpoint (``.pkl``) instead of
+a posterior file, only a single ``q0.5`` column is written (the MAP point
+estimate, not a true posterior quantile).
 
 Step 6: Predict Growth (``tfs-predict-growth``)
 ------------------------------------------------
@@ -206,6 +236,22 @@ Output (default ``--out_prefix tfs_cat_response``):
 * ``tfs_cat_response.csv`` — one row per (genotype, titrant_name) with
   ``best_model``, AIC weights, and fitted parameters for every model.
 
+Step 9: Summarise Fit (``tfs-summarize-fit``)
+----------------------------------------------
+
+Collects the outputs of the completed run, computes prediction quality
+statistics, and writes diagnostic plots and tables to a ``summary/``
+subdirectory.  It can be run after Steps 5–7 are complete and does not
+require the full posterior file to be retained.
+
+.. code-block:: bash
+
+    tfs-summarize-fit out/
+
+For a full guide to every output file — including how to read the theta
+correlation plots, calibration curves, and parameter-recovery scatter
+plots — see :doc:`summarize-fit`.
+
 ---
 
 .. _model-components:

docs/source/index.rst

@@ -8,9 +8,11 @@ tfscreen
    :maxdepth: 2
    :caption: Contents:
 
+   quickstart
    simulation
    process-raw
    analysis
+   summarize-fit
    grid
    ligandmpnn-features
 
@@ -30,6 +32,12 @@ The thermodynamic observable is operator occupancy (*θ*), inferred from *E. col
 
 The model in ``tfscreen.tfmodel`` is designed to jointly capture high-throughput growth data and low-throughput binding data, learning the quantitative relationship between operator occupancy and growth rate in the process.
 
+Quick Start
+-----------
+New to ``tfscreen``?  The :doc:`quickstart` page walks through a complete
+simulate-and-analyze run using the bundled example in
+``examples/simulate-and-analyze/``.
+
 Simulation
 ----------
 `tfscreen` allows you to simulate high-throughput screens starting from thermodynamic models of transcription factor binding and activity. See the :doc:`simulation` page for more details.