docs: document French-Wilson amplitudes, R-free flags and model-based maps

Update CPU_DATA_ANALYSIS.md for the new merge/analysis capabilities:
  - §10.7 R-free test-set flags: Friedel-merged Laue-ASU key (Bijvoet-safe),
    deterministic splitmix64 selection, resolution stratification;
  - §10.8 French-Wilson amplitudes: posterior-mean |F| from I and sigma with the
    Wilson prior (centric/acentric, epsilon multiplicity, strong-reflection
    short-circuit), written as MTZ F/SIGF and mmCIF F_meas_au;
  - §14 model-based validation (rugnux --model): Fcalc, flat bulk solvent + overall
    and per-shell scaling, R-work/R-free vs a model, and 2Fo-Fc / Fo-Fc maps.
Also refresh the scope list and references, and correct the old "intensities only"
note (the output now carries amplitudes too).

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>

docs/CPU_DATA_ANALYSIS.md

@@ -14,7 +14,9 @@ This document describes the crystallographic algorithms implemented in Jungfrauj
 8. Bragg integration by either 2D box summation or profile fitting (Kabsch, reference-free),
 9. scaling and merging,
 10. merge-level error modelling and outlier rejection,
-11. auxiliary statistics (Wilson plot, ⟨I/σ(I)⟩, CC1/2, CCref).
+11. auxiliary statistics (Wilson plot, ⟨I/σ(I)⟩, CC1/2, CCref),
+12. amplitude estimation (French–Wilson) and R-free test-set flagging,
+13. optional model-based validation: R-free against a supplied model and 2Fo−Fc / Fo−Fc electron-density maps.
 
 ## References
 
@@ -25,6 +27,9 @@ The methods are inspired and reuising solutions implemented in:
 - T. A. White et al., CrystFEL method papers (spot finding, three‑ring integration, serial/still diffraction processing concepts).
 - J. Kieffer & J. P. Wright, "PyFAI: a Python library for high performance azimuthal integration on GPU", *Powder Diffraction* **28** (2013), S339-S350 (detector geometry definition, azimuthal integration)
 - H. Powell, "The Rossmann Fourier autoindexing algorithm in MOSFLM", *Acta Cryst.* **D55** (1999), 1690-1695 (FFT indexing)
+- S. French & K. Wilson, "On the treatment of negative intensity observations", *Acta Cryst.* **A34** (1978), 517-525 (Bayesian amplitude estimation from intensities).
+- A. T. Brünger, "Free R value: a novel statistical quantity for assessing the accuracy of crystal structures", *Nature* **355** (1992), 472-475 (R-free cross-validation).
+- M. Wojdyr, "GEMMI: A library for structural biology", *J. Open Source Softw.* **7** (2022), 4200 (model / structure-factor / map machinery used in §14).
 (list is not exhaustive)
 
 ## 1. Geometry, reciprocal-space mapping, and basic quantities
@@ -561,6 +566,34 @@ After scale-fulls, two **optional correction surfaces** can be fitted on the com
 
 Both surfaces are **cross-validated**: fitted on even-numbered frames and kept only if they improve the held-out odd-frame symmetry-equivalent agreement by a clear margin (and vice versa). A surface fitted to noise where its systematic is absent therefore does not generalize and is discarded — an opt-in correction never adds scatter.
 
+### 10.7 R-free test-set flags
+
+A fraction of the unique reflections (`rfree_fraction`, default 0.05) is flagged as a **free (test) set**, written to the output (MTZ `FreeR_flag`, mmCIF `_refln.status_free`, a text-HKL column) for model validation (§14) and for downstream refinement. The flag is a pure function of the reflection's **Friedel-merged (Laue) ASU index**, which gives three properties:
+
+- all symmetry- and Friedel-equivalent reflections share one flag — in particular a Bijvoet pair $I(+)/I(-)$, kept as two separate merged rows in anomalous mode, is **never split** across the work and free sets (which would bias R-free);
+- the free/work decision is a deterministic hash of that key, so the same reflection always lands in the same set — reproducible run-to-run and independent of the order in which observations were merged;
+- the set is **stratified by resolution**: within each of 20 shells (equal-volume in $1/d^2$), exactly $\mathrm{round}(\text{fraction}\cdot n)$ of the distinct reflections — those with the smallest hash — are flagged free, so ~5 % of the data is free in every shell.
+
+### 10.8 French–Wilson amplitudes
+
+The last step of the merge estimates a Bayesian structure-factor amplitude $|F|$ for each unique reflection from its intensity $I$ and error $\sigma$, so the output carries amplitudes alongside intensities (a naïve $\sqrt{\max(I,0)}$ turns every weak or negative measurement into a biased — or zero — amplitude). With the Wilson prior for the true intensity $J\ge 0$ at that resolution,
+
+$
+P_\mathrm{acentric}(J) \propto e^{-J/\Sigma},\qquad
+P_\mathrm{centric}(J) \propto J^{-1/2}\,e^{-J/2\Sigma},
+$
+
+and a Gaussian likelihood $\mathcal{N}(I;J,\sigma^2)$, the posterior mean amplitude and its uncertainty are
+
+$
+\langle |F|\rangle = \frac{\int_0^\infty \sqrt{J}\,\mathcal{N}(I;J,\sigma^2)\,P(J)\,\mathrm{d}J}{\int_0^\infty \mathcal{N}(I;J,\sigma^2)\,P(J)\,\mathrm{d}J},\qquad
+\sigma_F = \sqrt{\langle J\rangle - \langle|F|\rangle^2}.
+$
+
+The prior mean is $\Sigma = \varepsilon\,\langle I/\varepsilon\rangle_\mathrm{shell}$, where $\varepsilon$ is the reflection's epsilon (symmetry-enhancement) multiplicity and $\langle I/\varepsilon\rangle$ is the Wilson mean in its resolution shell (so reflections on symmetry elements, and each shell, are treated correctly). Strong reflections ($I>4\sigma$) short-circuit to $|F|=\sqrt{I}$, where the French–Wilson bias is negligible; a reflection with an unusable $I/\sigma$ falls back to $\sqrt{\max(I,0)}$. The integral is evaluated numerically with a log-shift for stability.
+
+Amplitudes are written as MTZ `F`/`SIGF`, mmCIF `_refln.F_meas_au`/`F_meas_sigma_au`, and appended to the text HKL, alongside the intensity columns. The **same** $|F|$ feed the model-validation step (§14), so the reflection file and the maps use one consistent set of amplitudes.
+
 ---
 
 ## 11. Mosaicity and “profile radius” monitoring
@@ -613,4 +646,40 @@ A linear regression of $\log\langle I\rangle$ vs $1/d^2$ provides an estimate of
 - **Twinning check.** A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ (taken per resolution shell with noise-only shells skipped and Wilson outliers rejected, so a single strong reflection in a collapsed-mean shell cannot skew it) are written to the merged mmCIF as a twinning diagnostic. Twinning is only flagged in Laue classes where a merohedral twin law can exist; the holohedral high-symmetry classes ($4/mmm$, $6/mmm$, $m\bar{3}m$, and $\bar{3}m$ on a rhombohedral lattice) are exempt, so a low $\langle|L|\rangle$ there is reported as a statistical artefact rather than twinning.
 - **Outlier rejection.** Merging applies an optional per-observation median-based $N\sigma$ cut (default 6σ for `rot3d`) and an optional per-crystal $\Delta\mathrm{CC}_{1/2}$ image rejection (`--reject-delta-cchalf`, CrystFEL-style, off by default). The same $N\sigma$ cut is fed back into the error model: after an initial $a,b$ fit the parameters are re-fit once on the reflections that survive rejection (dropping any whose squared deviation exceeds $N\sigma^2\,[a\,\sigma^2 + (b\,\langle I\rangle)^2]$), so the calibrated errors describe the reflections that actually enter the merge rather than the pre-rejection pool.
 - **Automatic resolution cutoff.** By default the reported/written high-resolution limit is trimmed where $\mathrm{CC}_{1/2}$ falls off (logistic, target 0.30); `--scaling-high-resolution` overrides it and `--resolution-cutoff off` disables it.
-- **Intensities only.** The merged output carries intensities (mmCIF `intensity_meas`, MTZ `IMEAN`/`SIGIMEAN`); it does not convert to amplitudes $|F|$ (no French–Wilson / truncate step) — do that downstream.
+- **Amplitudes and intensities.** The merged output carries both intensities (mmCIF `intensity_meas`, MTZ `IMEAN`/`SIGIMEAN`) and French–Wilson amplitudes (mmCIF `F_meas_au`, MTZ `F`/`SIGF`; §10.8), so a downstream program can refine against either.
+
+---
+
+## 14. Model-based validation: R-free against a model and electron-density maps
+
+Offline (`rugnux --model model.pdb`) the merged data can be scored against a supplied atomic model and **initial** electron-density maps computed — enough to confirm that a model fits the data and to inspect the density, not a substitute for refinement. **The structure itself is not refined**; the model is only re-fractionalized into the data unit cell (a rigid cell adjustment, so a deposited model with a slightly different cell still lines up), and the observed amplitudes are the French–Wilson $|F|$ from §10.8, so the R-free and the maps use exactly the same amplitudes as the written reflection file. The model, structure-factor, bulk-solvent and FFT machinery is provided by GEMMI.
+
+### 14.1 Model structure factors
+
+The model electron density is sampled on a grid (IT92 X-ray form factors, with a Refmac-compatible Gaussian blur chosen for the grid spacing) and Fourier-transformed to structure factors $F_\mathrm{calc}(hkl)$ up to the data resolution.
+
+### 14.2 Bulk solvent and scaling
+
+A flat bulk-solvent mask around the model is transformed to $F_\mathrm{mask}$, and the model is scaled to the observed amplitudes in two stages:
+
+1. an overall least-squares fit of a scale $k$, an anisotropic $B$, and the flat-solvent parameters $k_\mathrm{sol}, B_\mathrm{sol}$, giving
+
+   $
+   F_\mathrm{model} = k\,e^{-\mathbf{h}^\top \mathbf{B}\,\mathbf{h}/4}\left(F_\mathrm{calc} + k_\mathrm{sol}\,e^{-B_\mathrm{sol}\,s^2}\,F_\mathrm{mask}\right),\quad s^2 = 1/4d^2;
+   $
+
+2. a smooth **per-resolution-shell scale** $K(1/d^2)$ (least squares $\sum F_o F_c/\sum F_c^2$ per shell, fit on the work reflections only) applied on top, which absorbs the residual radial $F_o/F_\mathrm{model}$ mismatch a single overall $B$ cannot.
+
+### 14.3 R-work and R-free
+
+Crystallographic R-factors are reported over the work and free sets (the §10.7 flags):
+
+$
+R = \frac{\sum \big|\,|F_o| - |F_\mathrm{model}|\,\big|}{\sum |F_o|},
+$
+
+with R-free the same sum restricted to the free set — an unbiased measure of how well the model explains data it was not scaled against.
+
+### 14.4 Electron-density maps
+
+Two maps are formed with the model phases $\varphi_\mathrm{model}$: a $2F_o-F_c$ map, coefficients $(2|F_o|-|F_\mathrm{model}|)\,e^{i\varphi_\mathrm{model}}$, and an $F_o-F_c$ difference map, $(|F_o|-|F_\mathrm{model}|)\,e^{i\varphi_\mathrm{model}}$, each inverse-Fourier-transformed to a real-space CCP4 map (`<prefix>_2fofc.ccp4`, `<prefix>_fofc.ccp4`). A map-coefficient MTZ (`<prefix>_maps.mtz`: `FP`, `FC`, `PHIC`, `FWT`/`PHWT`, `DELFWT`/`PHDELWT`, `FREE`) is written alongside so the maps can be reopened or rebuilt in Coot / PyMOL. These are unweighted difference coefficients (no $\sigma_A$ / figure-of-merit weighting), which is why they are described as *initial* maps.