diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 32b5f809..1c3d1a88 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -3,11 +3,11 @@ ### 1.0.0-rc.159 This is an UNSTABLE release. It includes many experimental features, as well as many AI generated fixes. We recommend using rc.152 for production use. -* rugnux: Add `--model model.pdb` - score the merged data against an atomic model and compute initial maps. It reports R-work/R-free (fitting an optimized flat bulk solvent and a per-resolution-shell scale to the observed amplitudes) and writes 2Fo-Fc / Fo-Fc electron-density maps (CCP4) plus a map-coefficient MTZ. The structure itself is not refined; the model is only re-fractionalised into the data cell. +* rugnux: Add `--model model.pdb` - score the merged data against an atomic model and compute initial maps. It reports R-work/R-free (scaling the model to the observed amplitudes with an overall scale, an anisotropic B and a flat bulk solvent - the standard few-parameter model, so a batch of maps stays directly comparable) and writes 2Fo-Fc / Fo-Fc electron-density maps (CCP4) plus a map-coefficient MTZ. The structure itself is not refined; the model is only re-fractionalised into the data cell. * rugnux: The merged reflection output now carries French-Wilson amplitudes (|F| and its sigma) next to the intensities - MTZ `F`/`SIGF`, mmCIF `_refln.F_meas_au`, and the text HKL - computed with the correct centric/acentric Wilson prior and epsilon multiplicity, so a downstream program (e.g. phenix.refine) can refine against amplitudes. The intensity columns are unchanged. -* rugnux: R-free test-set flags are now assigned deterministically and consistently across symmetry - a Bijvoet pair I(+)/I(-) is never split between the work and free sets, the free set is stratified to about 5% per resolution shell, and the assignment is reproducible run-to-run. +* rugnux: R-free test-set flags are now assigned deterministically and consistently across symmetry - a Bijvoet pair I(+)/I(-) is never split between the work and free sets, and the assignment is a reproducible per-hkl hash that depends only on the reflection index, so every dataset of one crystal form gets the same ~5% free set (what a multi-dataset campaign such as PanDDA needs). When a reference MTZ carries a `FreeR_flag` column its test set is imported instead, letting a whole campaign inherit one shared free set. * rugnux: A reference MTZ (`--reference-mtz`) can now fix the space group and cell for rotation data too (previously rejected), without being used to scale - the rotation merge stays self-consistent. When the crystal has an indexing (merohedral) ambiguity - a lattice symmetry higher than its Laue symmetry, e.g. P3/P4/P6/C2 - the reference also resolves it: each candidate reindexing (identity plus the twin-law cosets of the metric symmetry) is scored by its intensity correlation against the reference and the data are re-merged in the best-correlating one. This is a metric-preserving relabelling of hkl (the cell is unchanged) and a no-op for a holohedral crystal such as lysozyme. -* rugnux: `--model` validation now aligns the data to the model before scoring - the observed reflections are reindexed into the model's enantiomorph when the two differ only by hand (indistinguishable from merged intensities), and, for a crystal with a merohedral ambiguity, each candidate reindexing is fit to the model and the one giving the lowest R-free is kept. +* rugnux: `--model` validation now aligns the data to the model before scoring - the observed reflections are reindexed into the model's enantiomorph when the two differ only by hand (indistinguishable from merged intensities). A merohedral indexing ambiguity is resolved against the reference MTZ when one is given (so a whole campaign shares one indexing convention); only with a model and no reference does validation fall back to fitting each candidate reindexing and keeping the lowest R-free. * rugnux: De-novo symmetry - recover a genuine high-symmetry group whose data are imperfectly scaled. Such a merge's within-orbit chi² lands just past the self-consistency bound (each real symmetry step adds a little systematic scatter), right where a merohedral twin also lands, so the chi² ratio alone cannot separate them. The candidate is now rescued when the extra intensity-proportional systematic error it invokes stays small relative to the confirmed subgroup - a genuine symmetry step gains multiplicity without inflating the merge error model's b, whereas a twin forces non-equivalent reflections together and b balloons. Fixes cubic insulin (I23 instead of I222) with no change to any other crystal in the test battery, including the twins that must stay in their lower symmetry. * Docs: Document the French-Wilson amplitude estimation, R-free flagging, reference-based space-group/ambiguity resolution, and model-based validation/maps in CPU_DATA_ANALYSIS.md. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 076a14d3..b58c94cd 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -572,7 +572,9 @@ A fraction of the unique reflections (`rfree_fraction`, default 0.05) is flagged - 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. +- the hash depends only on the reflection index, **not** on this dataset's resolution range or which reflections it happens to contain, so a uniform draw takes ~`rfree_fraction` of the distinct reflections free and — crucially — **every dataset of one crystal form gets the same free set**. That cross-dataset consistency is what a multi-dataset campaign (ensemble refinement, PanDDA) requires; a per-shell stratification tied to each dataset's own $d_\mathrm{min}$ would break it. + +When a reference MTZ (`--reference-mtz`) carries a `FreeR_flag` column, its test set is **imported** instead: every merged reflection whose Laue-ASU index matches the reference takes the reference's flag (reflections absent from the reference keep the hash flag). This lets a whole fragment-screening campaign inherit one shared free set from the apo/reference dataset. The CCP4/refmac convention (test set = flag 0, including the historical 0–19 form) is assumed, with the complement taken automatically if flag 0 would be the majority (a phenix-style file where 1 marks free). ### 10.8 French–Wilson amplitudes @@ -668,15 +670,13 @@ The model electron density is sampled on a grid (IT92 X-ray form factors, with a ### 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: +A flat bulk-solvent mask around the model is transformed to $F_\mathrm{mask}$, and the model is scaled to the observed amplitudes by an overall least-squares fit of a scale $k$, an anisotropic $B$, and the flat-solvent parameters $k_\mathrm{sol}, B_\mathrm{sol}$: -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. +$ - $ - 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. +This is the standard, few-parameter scaling model used by refinement programs. A dataset-specific free-form per-resolution-shell rescale would lower this dataset's R a little, but it reshapes each map's radial amplitude profile differently, so a batch of maps would no longer be directly comparable — for a fragment-screening / PanDDA campaign, comparable maps across datasets matter more than the last bit of per-dataset R, so it is deliberately not applied. ### 14.3 R-work and R-free @@ -697,4 +697,4 @@ Two maps are formed with the model phases $\varphi_\mathrm{model}$: a $2F_o-F_c$ The model fixes a definite hand and indexing, but the merged data need not share them, so before comparison the observed reflections are brought into the model's frame. - **Enantiomorph / screw.** When the data space group is the enantiomorph of the model's (e.g. data $P4_12_12$, model $P4_32_12$; or $P3_1/P3_2$), the two are **indistinguishable from merged intensities** — $|F_\mathrm{calc}|$ is invariant under the change of hand, so R-free cannot choose between them and probing would be meaningless. The hand is therefore taken from the model: the observed reflections are reindexed by the change-of-hand operator into the model's enantiomorph. Only the map phases (the density's hand) depend on this choice. -- **Indexing (merohedral) ambiguity.** When the crystal has a merohedral ambiguity (§10.9), the observed intensities *do* differ between indexings, so here R-free is informative: the scaled model is fit to each reindexing of the data (identity plus the twin-law cosets) and the one giving the **lowest R-free** is kept for the reported statistics and the maps. This is a no-op for a holohedral crystal (no twin laws), e.g. lysozyme. +- **Indexing (merohedral) ambiguity.** When the crystal has a merohedral ambiguity (§10.9), the observed intensities *do* differ between indexings, and the right one is chosen against the best available reference. **If a reference MTZ was supplied, the data were already reindexed to agree with it** (§10.9 — by the reference-intensity correlation, at the merge stage for rotation data or per image in stills scaling), and model validation keeps that authoritative choice. **Only with a model and no reference** does validation resolve the ambiguity itself, as a fallback: the scaled model is fit to each reindexing of the data (identity plus the twin-law cosets) and the one giving the **lowest R-free** is kept. This matters for a multi-dataset campaign — a single shared reference fixes one indexing convention for every dataset, whereas an independent per-dataset lowest-R-free choice could send borderline datasets to different conventions. A no-op either way for a holohedral crystal (no twin laws), e.g. lysozyme. diff --git a/image_analysis/LoadFCalcFromMtz.cpp b/image_analysis/LoadFCalcFromMtz.cpp index a618538b..00d36b5a 100644 --- a/image_analysis/LoadFCalcFromMtz.cpp +++ b/image_analysis/LoadFCalcFromMtz.cpp @@ -69,6 +69,13 @@ ReferenceMtzData LoadReferenceMtz(const std::string& path, out.used_column = col->label; out.used_column_type = col->type; + // Optional cross-validation flags: a campaign shares ONE test set, so if the reference carries a + // FreeR flag we import it and every dataset can inherit the same free reflections. + const gemmi::Mtz::Column* free_col = nullptr; + for (const char* label : {"FreeR_flag", "FREE", "FREER", "RFREE", "R-free-flags", "FreeRflag"}) + if ((free_col = mtz.column_with_label(label, nullptr, 'I')) != nullptr) + break; + // Cell and space group the reference was recorded in, for display and the consistency check. const gemmi::UnitCell& gc = mtz.get_cell(); const bool cell_valid = gc.a > 0 && gc.b > 0 && gc.c > 0; @@ -83,6 +90,9 @@ ReferenceMtzData LoadReferenceMtz(const std::string& path, } out.reflections.reserve(static_cast(mtz.nreflections)); + std::vector raw_free; // raw FreeR value, aligned with out.reflections + if (free_col != nullptr) + raw_free.reserve(static_cast(mtz.nreflections)); const std::size_t stride = mtz.columns.size(); double d_min = std::numeric_limits::max(); @@ -110,6 +120,10 @@ ReferenceMtzData LoadReferenceMtz(const std::string& path, } out.reflections.emplace_back(r); + if (free_col != nullptr) { + const float f = (*free_col)[static_cast(i)]; + raw_free.push_back(std::isnan(f) ? -1 : static_cast(std::lround(f))); + } } if (d_max > 0.0) { @@ -117,6 +131,25 @@ ReferenceMtzData LoadReferenceMtz(const std::string& path, out.d_max = d_max; } + // Set the free flag from the raw column. The CCP4/refmac convention is that the test set is + // flag 0 (this also handles the historical 0-19 thin-shell format, where 0 is ~5%). If flag 0 + // would instead be the majority (a phenix-style file where 1 marks free), take the complement. + if (free_col != nullptr && !raw_free.empty()) { + std::size_t zeros = 0; + for (int f : raw_free) + if (f == 0) + ++zeros; + const bool free_is_zero = zeros * 2 <= raw_free.size(); + for (std::size_t j = 0; j < out.reflections.size(); ++j) { + const bool is_free = (raw_free[j] == 0) == free_is_zero; + out.reflections[j].rfree_flag = is_free; + if (is_free) + ++out.n_free; + } + out.has_free_flags = true; + out.free_column = free_col->label; + } + return out; } diff --git a/image_analysis/LoadFCalcFromMtz.h b/image_analysis/LoadFCalcFromMtz.h index 284867e2..d44c4f94 100644 --- a/image_analysis/LoadFCalcFromMtz.h +++ b/image_analysis/LoadFCalcFromMtz.h @@ -34,6 +34,12 @@ struct ReferenceMtzData { std::vector candidate_columns; double d_min = 0.0; double d_max = 0.0; + + // Cross-validation (free-R) flags, if the MTZ carried a FreeR_flag column. When present, the + // per-reflection rfree_flag above is set from it and a whole campaign can inherit one test set. + bool has_free_flags = false; + std::string free_column; // the FreeR column label used, for display + int n_free = 0; // number of reflections flagged free (test set) }; // Load reference reflections from an MTZ. With no column requested the smart default is used: diff --git a/image_analysis/scale_merge/FrenchWilson.cpp b/image_analysis/scale_merge/FrenchWilson.cpp index 0f3a5d8c..3e5a1cde 100644 --- a/image_analysis/scale_merge/FrenchWilson.cpp +++ b/image_analysis/scale_merge/FrenchWilson.cpp @@ -123,7 +123,7 @@ void ApplyFrenchWilson(std::vector &merged, int32_t space_grou const auto s = shells.GetShell(r.d); const double sigma_wilson = epsilon(r) * (s ? shell_mean[*s] : global_mean); const bool centric = gops.is_reflection_centric({{r.h, r.k, r.l}}); - const Posterior post = integrate_posterior(r.I, r.sigma, sigma_wilson, !centric, + const Posterior post = integrate_posterior(r.I, r.sigma, sigma_wilson, centric, opts.integration_points); r.F = static_cast(post.mean_F); r.sigmaF = static_cast(std::sqrt(std::max(0.0, post.mean_I - post.mean_F * post.mean_F))); diff --git a/image_analysis/scale_merge/RfreeFlags.cpp b/image_analysis/scale_merge/RfreeFlags.cpp index e039113f..875d1913 100644 --- a/image_analysis/scale_merge/RfreeFlags.cpp +++ b/image_analysis/scale_merge/RfreeFlags.cpp @@ -3,13 +3,9 @@ #include "RfreeFlags.h" -#include -#include -#include -#include +#include #include "HKLKey.h" -#include "../../common/ResolutionShells.h" namespace { // splitmix64 bit-mix of a key -> uniform double in [0, 1). Same key -> same value, so all @@ -31,49 +27,35 @@ void AssignRfreeFlags(std::vector &merged, int32_t space_group if (rfree_fraction <= 0.0 || merged.empty()) return; - // Friedel-merged (Laue) ASU key: I(+)/I(-) and every symmetry mate collapse to one key, so a - // reflection's free/work status is a pure function of this key -> mates never split. + // The flag is a pure function of the Friedel-merged (Laue) ASU key: symmetry- and Friedel- + // equivalent reflections collapse to one key and so share a flag (a Bijvoet pair I(+)/I(-) is + // never split across the work and free sets), and the draw depends only on the reflection index + // - not on this dataset's resolution range or which reflections it happens to contain. So every + // dataset of one crystal form gets the SAME free set, which is what a multi-dataset campaign + // (ensemble refinement, PanDDA) needs. A uniform hash draws ~rfree_fraction of the distinct + // reflections free; a stratified per-shell draw would be tied to the dataset and break that. const HKLKeyGenerator laue_key(/*merge_friedel=*/true, space_group_number); - - float d_min = std::numeric_limits::max(), d_max = 0.0f; - for (const auto &r : merged) - if (std::isfinite(r.d) && r.d > 0.0f) { - d_min = std::min(d_min, r.d); - d_max = std::max(d_max, r.d); - } - - // No usable resolution: fall back to a global deterministic hash (still mate-consistent). - if (!(d_min < d_max && d_min > 0.0f)) { - for (auto &r : merged) - r.rfree_flag = UniformFromKey(laue_key(r).pack()) < rfree_fraction; - return; - } - - // Stratify by resolution: within each shell take exactly round(fraction * n) of the distinct - // Laue reflections as the free set, chosen by hash rank. This guarantees ~fraction of the data - // is free in every shell, deterministically and without splitting any Bijvoet/symmetry set. - constexpr int n_shells = 20; - ResolutionShells shells(d_min * 0.999f, d_max * 1.001f, n_shells); - std::vector>> shell_keys(n_shells); // (hash, key) - std::unordered_set seen; - seen.reserve(merged.size()); - for (const auto &r : merged) { - const uint64_t key = laue_key(r).pack(); - if (!seen.insert(key).second) // count each Laue set once (Bijvoet mates share a key) - continue; - const auto shell = shells.GetShell(r.d); - if (shell) - shell_keys[*shell].push_back({UniformFromKey(key), key}); - } - - std::unordered_set free_keys; - free_keys.reserve(static_cast(rfree_fraction * seen.size()) + 1); - for (auto &keys : shell_keys) { - std::sort(keys.begin(), keys.end()); - const auto n_free = static_cast(std::llround(rfree_fraction * keys.size())); - for (size_t i = 0; i < n_free; ++i) - free_keys.insert(keys[i].second); - } for (auto &r : merged) - r.rfree_flag = free_keys.count(laue_key(r).pack()) > 0; + r.rfree_flag = UniformFromKey(laue_key(r).pack()) < rfree_fraction; +} + +size_t ApplyReferenceFreeFlags(std::vector &merged, int32_t space_group_number, + const std::vector &reference) { + // Reference free/work partition keyed by the Friedel-merged (Laue) ASU index, so it transfers + // regardless of which Bijvoet mate / symmetry equivalent each dataset happens to have measured. + const HKLKeyGenerator laue_key(/*merge_friedel=*/true, space_group_number); + std::unordered_map ref_flag; + ref_flag.reserve(reference.size()); + for (const auto &r : reference) + ref_flag[laue_key(r).pack()] = r.rfree_flag; + + size_t matched = 0; + for (auto &r : merged) { + const auto it = ref_flag.find(laue_key(r).pack()); + if (it != ref_flag.end()) { // reflections absent from the reference keep their hash flag + r.rfree_flag = it->second; + ++matched; + } + } + return matched; } diff --git a/image_analysis/scale_merge/RfreeFlags.h b/image_analysis/scale_merge/RfreeFlags.h index 6a51ed43..4645830d 100644 --- a/image_analysis/scale_merge/RfreeFlags.h +++ b/image_analysis/scale_merge/RfreeFlags.h @@ -3,16 +3,27 @@ #pragma once +#include #include #include #include "../../common/Reflection.h" -// Assign R-free (test-set) flags to merged reflections. Two guarantees: -// - All symmetry- and Friedel-equivalent reflections share one flag (the flag is a pure function -// of the Friedel-merged Laue-ASU key), so a Bijvoet pair I(+)/I(-) is never split across sets, -// and the assignment is reproducible run-to-run and independent of merge/call order. -// - The free set is stratified across resolution: ~rfree_fraction of the distinct reflections in -// every resolution shell is flagged free. +// Assign R-free (test-set) flags to merged reflections. The flag is a pure function of the +// Friedel-merged (Laue) ASU key, which gives three properties: +// - all symmetry- and Friedel-equivalent reflections share one flag, so a Bijvoet pair I(+)/I(-) +// is never split across the work and free sets; +// - the assignment is reproducible run-to-run and independent of merge/call order; +// - it depends only on the reflection index, not on the dataset's resolution range or contents, so +// every dataset of one crystal form gets the SAME free set (needed for multi-dataset campaigns: +// ensemble refinement, PanDDA). A uniform hash draws ~rfree_fraction of the distinct reflections. void AssignRfreeFlags(std::vector &merged, int32_t space_group_number, double rfree_fraction); + +// Override the R-free flags on `merged` from a reference set (e.g. a reference MTZ's FreeR_flag +// column), so an entire campaign shares one test set. A merged reflection whose Friedel-merged +// (Laue) ASU key matches a reference reflection takes the reference's flag; reflections with no +// match keep the flag they already carry (the deterministic per-hkl hash from AssignRfreeFlags). +// Returns the number of merged reflections matched to the reference (for logging). +size_t ApplyReferenceFreeFlags(std::vector &merged, int32_t space_group_number, + const std::vector &reference); diff --git a/rugnux/ModelValidation.cpp b/rugnux/ModelValidation.cpp index 935a7930..06a83847 100644 --- a/rugnux/ModelValidation.cpp +++ b/rugnux/ModelValidation.cpp @@ -6,7 +6,6 @@ #include #include #include -#include #include #include @@ -52,7 +51,8 @@ ModelValidationResult ValidateAgainstModel(const std::vector & const std::string &model_path, const std::string &output_prefix, Logger &logger, - std::optional data_space_group_number) { + std::optional data_space_group_number, + bool probe_indexing_ambiguity) { ModelValidationResult result; // --- read the atomic model --- @@ -154,7 +154,7 @@ ModelValidationResult ValidateAgainstModel(const std::vector & gemmi::AsuData> fmodel, map2fofc, mapfofc; std::unordered_map> obs_by_hkl; // hkl -> (Fobs, is_free) double r_work = 1, r_free = 1, k_sol = 0, b_sol = 0, k_overall = 0; - int n_w = 0, n_f = 0, n_shells = 0; + int n_w = 0, n_f = 0; }; auto fit_model = [&](const std::vector &obs_in) -> Fit { Fit out; @@ -188,50 +188,12 @@ ModelValidationResult ValidateAgainstModel(const std::vector & out.b_sol = scaling.b_sol; out.k_overall = scaling.k_overall; - // --- per-resolution-shell scaling on top of the overall scaling --- - // Fit a smooth scale K(1/d^2) least-squares per shell (sum(Fo*Fc)/sum(Fc^2)) and apply it to - // every reflection. This mops up the residual radial Fobs/Fmodel mismatch that a single overall - // B leaves behind. The shells are fit on the work set only, so R-free stays untouched by the fit. - { - struct ShellPoint { double s2, fo, fc; }; - std::vector pts; - pts.reserve(out.fmodel.v.size()); - for (const auto &hv : out.fmodel.v) { - auto it = out.obs_by_hkl.find(hkl_key(hv.hkl)); - if (it == out.obs_by_hkl.end() || it->second.second) continue; // skip missing + free set - if (it->second.first <= 0) continue; - pts.push_back({ucell.calculate_1_d2(hv.hkl), it->second.first, std::abs(hv.value)}); - } - std::sort(pts.begin(), pts.end(), [](const ShellPoint &a, const ShellPoint &b) { return a.s2 < b.s2; }); - int nb = std::max(6, std::min(40, static_cast(pts.size() / 300))); // ~300 refl/shell - std::vector shell_s2, shell_k; - for (int b = 0; b < nb; ++b) { - size_t i0 = pts.size() * b / nb, i1 = pts.size() * (b + 1) / nb; - double sum_fofc = 0, sum_fc2 = 0, sum_s2 = 0; - int n = 0; - for (size_t i = i0; i < i1; ++i) { - sum_fofc += pts[i].fo * pts[i].fc; - sum_fc2 += pts[i].fc * pts[i].fc; - sum_s2 += pts[i].s2; - ++n; - } - if (sum_fc2 > 0 && n > 0) { shell_s2.push_back(sum_s2 / n); shell_k.push_back(sum_fofc / sum_fc2); } - } - auto scale_at = [&](double s2) -> double { - if (shell_s2.empty()) return 1.0; - if (s2 <= shell_s2.front()) return shell_k.front(); - if (s2 >= shell_s2.back()) return shell_k.back(); - for (size_t i = 1; i < shell_s2.size(); ++i) - if (s2 <= shell_s2[i]) { - double t = (s2 - shell_s2[i - 1]) / (shell_s2[i] - shell_s2[i - 1]); - return shell_k[i - 1] * (1 - t) + shell_k[i] * t; - } - return shell_k.back(); - }; - for (auto &hv : out.fmodel.v) - hv.value *= static_cast(scale_at(ucell.calculate_1_d2(hv.hkl))); - out.n_shells = static_cast(shell_k.size()); - } + // The model is scaled to the data with an overall scale, an anisotropic B and a flat bulk + // solvent only - the standard, few-parameter model that refinement programs use. A dataset- + // specific free-form per-resolution-shell rescale would lower this dataset's R a little, but + // it reshapes each map's radial amplitude profile differently, so a batch of maps would no + // longer be directly comparable. For a fragment-screening / PanDDA campaign, comparable maps + // across datasets matter more than the last bit of per-dataset R, so it is deliberately omitted. // --- R-work / R-free and map coefficients (2Fo-Fc and Fo-Fc, model phases) --- out.map2fofc.unit_cell_ = ucell; out.map2fofc.spacegroup_ = sg; @@ -254,28 +216,33 @@ ModelValidationResult ValidateAgainstModel(const std::vector & return out; }; - // --- probe indexing (merohedral) ambiguities: keep the reindexing with the lowest R-free --- - // No-op for a holohedral crystal (no twin laws), e.g. lysozyme, where the only candidate is the - // identity. The enantiomorph/screw ambiguity is NOT probed here: |Fcalc| is the same for both - // hands, so R-free cannot distinguish them - that is resolved from the model hand above. - const auto reindex_ops = ReindexAmbiguityOperators(cell, sg->number); + // --- indexing (merohedral) ambiguity --- + // When a reference MTZ was supplied, the data were already reindexed to agree with the reference + // intensities (at the merge stage for rotation data, per image in stills scaling), and that + // choice is authoritative - we keep it. Only with a model and NO reference do we resolve the + // ambiguity here, as a fallback, by fitting each candidate reindexing and keeping the lowest + // R-free. A no-op either way for a holohedral crystal (no twin laws), e.g. lysozyme. The + // enantiomorph/screw ambiguity is never probed by R-free: |Fcalc| is the same for both hands, so + // it cannot distinguish them - that is taken from the model hand above. Fit best = fit_model(obs); - bool did_reindex = false; - for (const auto &op : reindex_ops) { - Fit cand = fit_model(ReindexReflections(obs, op)); - if (cand.r_free < best.r_free) { best = std::move(cand); did_reindex = true; } + if (probe_indexing_ambiguity) { + const auto reindex_ops = ReindexAmbiguityOperators(cell, sg->number); + bool did_reindex = false; + for (const auto &op : reindex_ops) { + Fit cand = fit_model(ReindexReflections(obs, op)); + if (cand.r_free < best.r_free) { best = std::move(cand); did_reindex = true; } + } + if (!reindex_ops.empty()) + logger.Info("Model validation: no reference - probed {} indexing solution(s) against the model; " + "{} (R-free {:.4f})", reindex_ops.size() + 1, + did_reindex ? "reindexed to the lower-R-free solution" : "kept the current indexing", + best.r_free); } - if (!reindex_ops.empty()) - logger.Info("Model validation: probed {} indexing solution(s); {} (R-free {:.4f})", - reindex_ops.size() + 1, - did_reindex ? "reindexed to the lower-R-free solution" : "kept the current indexing", - best.r_free); gemmi::AsuData> &fmodel = best.fmodel; gemmi::AsuData> &map2fofc = best.map2fofc; gemmi::AsuData> &mapfofc = best.mapfofc; std::unordered_map> &obs_by_hkl = best.obs_by_hkl; - int n_shells = best.n_shells; result.r_work = best.r_work; result.r_free = best.r_free; @@ -345,8 +312,9 @@ ModelValidationResult ValidateAgainstModel(const std::vector & result.ok = true; result.maps_prefix = output_prefix; - logger.Info("Model validation: R-work={:.4f} ({} refl) R-free={:.4f} ({} refl) [overall + {}-shell scaling]", - result.r_work, result.n_work, result.r_free, result.n_free, n_shells); + logger.Info("Model validation: R-work={:.4f} ({} refl) R-free={:.4f} ({} refl) " + "[overall + anisotropic B + bulk solvent]", + result.r_work, result.n_work, result.r_free, result.n_free); logger.Info("Model validation: bulk solvent k_sol={:.3f} b_sol={:.1f}, k_overall={:.3f}", result.k_sol, result.b_sol, result.k_overall); logger.Info("Model validation: mean 2Fo-Fc density at atom centres = {:.2f} sigma", result.mean_atom_density_sigma); diff --git a/rugnux/ModelValidation.h b/rugnux/ModelValidation.h index 379753f5..34da9bf8 100644 --- a/rugnux/ModelValidation.h +++ b/rugnux/ModelValidation.h @@ -37,9 +37,16 @@ struct ModelValidationResult { // data_space_group_number is the space group the data was merged in. If it is the enantiomorph of // the model's space group (e.g. data P4(1)2(1)2, model P4(3)2(1)2 - indistinguishable from merged // intensities), the observed reflections are reindexed into the model's hand before comparison. +// +// probe_indexing_ambiguity controls how a merohedral (alternative-indexing) ambiguity is resolved. +// Pass false when a reference (MTZ) already fixed the indexing - the data are then kept in that +// authoritative indexing. Pass true (model only, no reference) to resolve it here as a fallback, by +// fitting each candidate reindexing to the model and keeping the lowest R-free. Holohedral crystals +// (no twin laws) are unaffected either way. ModelValidationResult ValidateAgainstModel(const std::vector &merged, const UnitCell &cell, const std::string &model_path, const std::string &output_prefix, Logger &logger, - std::optional data_space_group_number = std::nullopt); + std::optional data_space_group_number = std::nullopt, + bool probe_indexing_ambiguity = true); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index ab46fdac..5ceb6a80 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -28,6 +28,7 @@ #include "../image_analysis/image_preprocessing/ImagePreprocessorCPU.h" #include "../image_analysis/image_preprocessing/ImagePreprocessorBuffer.h" #include "../image_analysis/scale_merge/Merge.h" +#include "../image_analysis/scale_merge/RfreeFlags.h" #include "../image_analysis/scale_merge/RotationScaleMerge.h" #include "../image_analysis/scale_merge/ResolutionCutoff.h" #include "../image_analysis/scale_merge/ReindexAmbiguity.h" @@ -893,6 +894,16 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { result.merge_statistics = sm.statistics; result.has_reference = !config_.reference_data.empty(); + // Inherit the campaign's shared R-free test set from the reference MTZ (overriding the + // per-hkl hash that the merge assigned), so every dataset flags the same free reflections. + if (config_.reference_has_free_flags && !config_.reference_data.empty() && !sm.merged.empty()) { + const auto sg = experiment_.GetSpaceGroupNumber().value_or(1); + const size_t matched = ApplyReferenceFreeFlags(sm.merged, static_cast(sg), + config_.reference_data); + logger.Info("R-free flags: inherited the reference test set ({} of {} merged reflections matched)", + matched, sm.merged.size()); + } + if (result.consensus_cell && write_files && config_.write_merged) { phase("Writing reflections"); WriteReflections(sm.merged, *result.consensus_cell, experiment_, sm.statistics, @@ -908,9 +919,12 @@ ProcessResult Rugnux::Run(RugnuxObserver *observer) { if (result.consensus_cell && write_files && !config_.model_path.empty()) { phase("Validating against model"); const auto data_sg = experiment_.GetSpaceGroupNumber(); + // With a reference MTZ the merohedral indexing was already resolved against it (rotation + // merge / stills scaling), so trust that; only probe indexing by R-free when model-only. ValidateAgainstModel(sm.merged, *result.consensus_cell, config_.model_path, config_.output_prefix, logger, - data_sg ? std::optional(static_cast(*data_sg)) : std::nullopt); + data_sg ? std::optional(static_cast(*data_sg)) : std::nullopt, + /*probe_indexing_ambiguity=*/config_.reference_data.empty()); } } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 1a317883..da918b96 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -67,6 +67,10 @@ struct ProcessConfig { int64_t scaling_iter = 3; std::vector reference_data; + // When the reference MTZ carried a FreeR_flag column, reference_data reflections carry its + // free/work partition and the merged reflections inherit it (one shared test set per campaign). + bool reference_has_free_flags = false; + // Diagnostic: if set, the -P rot3d combine writes the unmerged fulls here (for comparison vs XDS). std::string observation_dump_path; diff --git a/rugnux/rugnux_cli.cpp b/rugnux/rugnux_cli.cpp index ea166605..30402b28 100644 --- a/rugnux/rugnux_cli.cpp +++ b/rugnux/rugnux_cli.cpp @@ -23,6 +23,7 @@ #include "../image_analysis/UpdateReflectionResolution.h" #include "../image_analysis/WriteReflections.h" #include "../image_analysis/scale_merge/Merge.h" +#include "../image_analysis/scale_merge/RfreeFlags.h" #include "../image_analysis/scale_merge/ScaleOnTheFly.h" #include "../image_analysis/scale_merge/RotationScaleMerge.h" #include "../image_analysis/scale_merge/ResolutionCutoff.h" @@ -882,16 +883,22 @@ int main(int argc, char **argv) { logger.Info("Loaded dataset from {}", input_file); std::vector reference_data; + bool reference_has_free_flags = false; if (!ref_mtz.empty()) { try { const auto reference = LoadReferenceMtz( ref_mtz, ref_column.empty() ? std::nullopt : std::optional(ref_column)); reference_data = reference.reflections; + reference_has_free_flags = reference.has_free_flags; logger.Info("Loaded {} reference reflections from {} (column {}{}{})", reference_data.size(), ref_mtz, reference.used_column, reference.squared ? ", squared to intensity" : "", reference.default_column ? ", auto-selected" : ", user-specified"); + if (reference.has_free_flags) + logger.Info("Reference carries R-free flags (column {}): {} of {} free; the merged " + "reflections will inherit this test set", + reference.free_column, reference.n_free, reference_data.size()); if (reference.d_max > 0.0) logger.Info("Reference resolution range {:.2f} - {:.2f} A", reference.d_max, reference.d_min); if (reference.cell.has_value()) @@ -1030,6 +1037,16 @@ int main(int argc, char **argv) { std::chrono::duration(std::chrono::steady_clock::now() - scale_start).count(), merged_reflections.size()); + // Inherit the campaign's shared R-free test set from the reference MTZ (overriding the + // per-hkl hash the merge assigned), so every dataset flags the same free reflections. + if (reference_has_free_flags && !reference_data.empty() && !merged_reflections.empty()) { + const auto sg = experiment.GetSpaceGroupNumber().value_or(1); + const size_t matched = ApplyReferenceFreeFlags(merged_reflections, static_cast(sg), + reference_data); + logger.Info("R-free flags: inherited the reference test set ({} of {} merged reflections matched)", + matched, merged_reflections.size()); + } + std::cout << merged_statistics; // Space-group determination lives in the full rugnux pipeline; --scale only consumes a space @@ -1053,9 +1070,12 @@ int main(int argc, char **argv) { if (!output_prefix.empty() && !model_pdb.empty()) { const auto data_sg = experiment.GetSpaceGroupNumber(); + // With a reference MTZ the merohedral indexing was already resolved (stills per-image + // scaling); only probe indexing by R-free when model-only, with no reference. ValidateAgainstModel(merged_reflections, *experiment.GetUnitCell(), model_pdb, output_prefix, logger, - data_sg ? std::optional(static_cast(*data_sg)) : std::nullopt); + data_sg ? std::optional(static_cast(*data_sg)) : std::nullopt, + /*probe_indexing_ambiguity=*/reference_data.empty()); } return 0; } @@ -1300,6 +1320,7 @@ int main(int argc, char **argv) { config.run_scaling = run_scaling; config.scaling_iter = scaling_iter; config.reference_data = reference_data; + config.reference_has_free_flags = reference_has_free_flags; config.observation_dump_path = dump_observations; config.model_path = model_pdb; // When merging, the merged reflections (.mtz/.cif) are the wanted output; skip the large diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 5b1d2565..b10a228d 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -77,6 +77,7 @@ ADD_EXECUTABLE(jfjoch_test RfreeFlagsTest.cpp FrenchWilsonTest.cpp ReindexAmbiguityTest.cpp + LoadReferenceMtzFreeFlagsTest.cpp UnitCellTest.cpp CCTest.cpp MultiLatticeSearchTest.cpp diff --git a/tests/FrenchWilsonTest.cpp b/tests/FrenchWilsonTest.cpp index 397655a5..dab69ea9 100644 --- a/tests/FrenchWilsonTest.cpp +++ b/tests/FrenchWilsonTest.cpp @@ -64,6 +64,29 @@ TEST_CASE("French-Wilson: amplitudes are always finite and non-negative", "[fren } } +TEST_CASE("French-Wilson: centric weak reflection gets a smaller amplitude than acentric", "[french_wilson]") { + // At the same resolution (same Wilson mean Sigma) and the same near-zero intensity, the centric + // prior puts more weight near |F|=0, so the posterior <|F|> is smaller than for an acentric + // reflection (0.80*sqrt(Sigma) vs 0.89*sqrt(Sigma) at I=0). This pins the centric/acentric prior + // the right way round: if the two priors were swapped the inequality below would flip. + std::vector v; + // Uniform strong background across resolution, so every shell has ~the same Wilson mean. + for (int h = 1; h <= 10; ++h) + for (int k = 1; k <= 10; ++k) + for (int l = 1; l <= 6; ++l) + v.push_back(Refl(h, k, l, 20.0f / (0.5f + 0.1f * (h + k + l)), 1000.0f, 30.0f)); + // Two weak (I=0) probes at the SAME resolution: (3,1,0) is centric in P4 (l=0 zone), (3,1,4) is + // acentric. d is set directly, so both share a shell (hence Sigma) regardless of the cell. + v.push_back(Refl(3, 1, 0, 5.0f, 0.0f, 10.0f)); + v.push_back(Refl(3, 1, 4, 5.0f, 0.0f, 10.0f)); + ApplyFrenchWilson(v, 75); // P4 + const auto& centric = v[v.size() - 2]; + const auto& acentric = v.back(); + CHECK(centric.F > 0.0f); + CHECK(acentric.F > 0.0f); + CHECK(centric.F < acentric.F); +} + TEST_CASE("French-Wilson: unusable sigma falls back to sqrt(max(I,0))", "[french_wilson]") { auto v = Background(); v.push_back(Refl(4, 0, 0, 12.0f, 144.0f, NAN)); // no sigma diff --git a/tests/LoadReferenceMtzFreeFlagsTest.cpp b/tests/LoadReferenceMtzFreeFlagsTest.cpp new file mode 100644 index 00000000..c2c5b811 --- /dev/null +++ b/tests/LoadReferenceMtzFreeFlagsTest.cpp @@ -0,0 +1,96 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include +#include +#include +#include + +#include + +#include "../image_analysis/LoadFCalcFromMtz.h" + +namespace { + // Write a small P4(3)2(1)2 MTZ with an IMEAN column and a FreeR_flag column whose free reflections + // carry `free_value` (and the rest its complement). ~1 in 10 reflections is flagged free. + std::string WriteMtzWithFreeR(const std::string& path, int free_value) { + gemmi::Mtz mtz(true); + mtz.spacegroup = gemmi::find_spacegroup_by_number(96); + mtz.set_cell_for_all(gemmi::UnitCell(78.0, 78.0, 37.0, 90.0, 90.0, 90.0)); + mtz.add_dataset("test"); + mtz.add_column("IMEAN", 'J', -1, -1, false); + mtz.add_column("FreeR_flag", 'I', -1, -1, false); + std::vector data; + int nref = 0; + for (int h = 1; h <= 10; ++h) + for (int k = 1; k <= 10; ++k) + for (int l = 1; l <= 6; ++l) { + const bool is_free = ((h * 7 + k * 13 + l * 17) % 10 == 0); + const int flag = is_free ? free_value : (free_value == 0 ? 1 : 0); + data.insert(data.end(), {static_cast(h), static_cast(k), + static_cast(l), 100.0f, static_cast(flag)}); + ++nref; + } + mtz.nreflections = nref; + mtz.data = std::move(data); + mtz.write_to_file(path); + return path; + } + + bool ExpectedFree(int h, int k, int l) { return (h * 7 + k * 13 + l * 17) % 10 == 0; } +} + +TEST_CASE("LoadReferenceMtz imports FreeR (CCP4 convention: test set = 0)", "[reference_mtz]") { + const auto path = (std::filesystem::temp_directory_path() / "rugnux_freer_ccp4.mtz").string(); + WriteMtzWithFreeR(path, /*free_value=*/0); + const auto ref = LoadReferenceMtz(path); + std::filesystem::remove(path); + + REQUIRE(ref.has_free_flags); + CHECK(ref.n_free > 0); + CHECK(ref.n_free < static_cast(ref.reflections.size()) / 2); // free is the minority + for (const auto& r : ref.reflections) + CHECK(r.rfree_flag == ExpectedFree(r.h, r.k, r.l)); +} + +TEST_CASE("LoadReferenceMtz imports FreeR (phenix convention: test set = 1) via the complement", "[reference_mtz]") { + // Here flag 0 is the majority (work), so the loader must take the complement and treat 1 as free. + const auto path = (std::filesystem::temp_directory_path() / "rugnux_freer_phenix.mtz").string(); + WriteMtzWithFreeR(path, /*free_value=*/1); + const auto ref = LoadReferenceMtz(path); + std::filesystem::remove(path); + + REQUIRE(ref.has_free_flags); + CHECK(ref.n_free > 0); + CHECK(ref.n_free < static_cast(ref.reflections.size()) / 2); + for (const auto& r : ref.reflections) + CHECK(r.rfree_flag == ExpectedFree(r.h, r.k, r.l)); +} + +TEST_CASE("LoadReferenceMtz without a FreeR column reports no free flags", "[reference_mtz]") { + const auto path = (std::filesystem::temp_directory_path() / "rugnux_no_freer.mtz").string(); + gemmi::Mtz mtz(true); + mtz.spacegroup = gemmi::find_spacegroup_by_number(96); + mtz.set_cell_for_all(gemmi::UnitCell(78.0, 78.0, 37.0, 90.0, 90.0, 90.0)); + mtz.add_dataset("test"); + mtz.add_column("IMEAN", 'J', -1, -1, false); + std::vector data; + int nref = 0; + for (int h = 1; h <= 6; ++h) + for (int k = 1; k <= 6; ++k) + for (int l = 1; l <= 4; ++l) { + data.insert(data.end(), {static_cast(h), static_cast(k), + static_cast(l), 100.0f}); + ++nref; + } + mtz.nreflections = nref; + mtz.data = std::move(data); + mtz.write_to_file(path); + const auto ref = LoadReferenceMtz(path); + std::filesystem::remove(path); + + CHECK_FALSE(ref.has_free_flags); + CHECK(ref.n_free == 0); +} diff --git a/tests/RfreeFlagsTest.cpp b/tests/RfreeFlagsTest.cpp index 601f77f7..ce589b15 100644 --- a/tests/RfreeFlagsTest.cpp +++ b/tests/RfreeFlagsTest.cpp @@ -91,9 +91,9 @@ TEST_CASE("R-free flags are shared across symmetry equivalents", "[rfree]") { CHECK(checked > 0); } -TEST_CASE("R-free flags are stratified across resolution shells", "[rfree]") { - // Reflections in three well-separated resolution bands: each band must get some free flags, - // i.e. the free set is not clumped into one shell. +TEST_CASE("R-free flags spread across resolution", "[rfree]") { + // The uniform per-hkl hash is uncorrelated with resolution, so each of three well-separated + // resolution bands still receives a share of the free set (it is not clumped into one shell). std::vector v; for (int i = 0; i < 1000; ++i) { v.push_back(Refl(1 + i, 2, 3, 8.0f)); // low res @@ -111,6 +111,47 @@ TEST_CASE("R-free flags are stratified across resolution shells", "[rfree]") { CHECK(hi > 0); } +TEST_CASE("R-free flags are identical across datasets of one crystal form", "[rfree]") { + // The key campaign property: the free set depends only on the reflection index, not on the + // dataset's resolution extent or which reflections it contains. Two datasets with different + // resolution ranges must flag every shared reflection the same way (a per-shell stratification + // tied to each dataset's own d_min/d_max would break this). + auto wide = Grid(-15, 15); + std::vector narrow; // a lower-resolution subset + for (const auto& r : wide) + if (r.d > 6.0f) + narrow.push_back(r); + AssignRfreeFlags(wide, 96, 0.05); + AssignRfreeFlags(narrow, 96, 0.05); + + std::map, bool> flag; + for (const auto& r : wide) flag[{r.h, r.k, r.l}] = r.rfree_flag; + int checked = 0; + for (const auto& r : narrow) { + auto it = flag.find({r.h, r.k, r.l}); + REQUIRE(it != flag.end()); + CHECK(it->second == r.rfree_flag); + ++checked; + } + CHECK(checked > 0); +} + +TEST_CASE("ApplyReferenceFreeFlags imports the reference test set", "[rfree]") { + // A reference with its own free set; a dataset that starts from the per-hkl hash must, after the + // import, carry exactly the reference's flags on every reflection they share. + auto reference = Grid(-12, 12); + AssignRfreeFlags(reference, 96, 0.07); + std::map, bool> ref_flag; + for (const auto& r : reference) ref_flag[{r.h, r.k, r.l}] = r.rfree_flag; + + auto data = Grid(-12, 12); + AssignRfreeFlags(data, 96, 0.30); // deliberately a different fraction/hash split + const size_t matched = ApplyReferenceFreeFlags(data, 96, reference); + CHECK(matched == data.size()); + for (const auto& r : data) + CHECK(r.rfree_flag == ref_flag[{r.h, r.k, r.l}]); +} + TEST_CASE("R-free fraction of zero flags nothing", "[rfree]") { auto v = Grid(1, 6); AssignRfreeFlags(v, 1, 0.0);