diff --git a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md index 1d8bd3450..61b4bf3a1 100644 --- a/docs/CPU_DATA_ANALYSIS_INTEGRATION.md +++ b/docs/CPU_DATA_ANALYSIS_INTEGRATION.md @@ -382,9 +382,10 @@ Each batch's $B$ is fitted on **resolution-shell means**, not on single observat ### 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` = `f`, 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: +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` = `f`, a text-HKL column) for model validation (§14) and for downstream refinement. The flag is a pure function of the reflection's orbit under the **lattice holohedry** — the point group of the cell's metric, found as twin laws are (Le Page two-folds within 3° obliquity, taking the lattice of the cell's own basis vectors), which contains the merging group — together with its Friedel mate. Where the cell does not carry the merging group (a space group forced on a metric without it), the Friedel-merged (Laue) ASU index of the merging group is used instead. That gives four 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); +- **twin-law mates share one flag too**, since a twin law is a lattice symmetry the crystal lacks. Keyed on the merging group instead, nearly every free reflection's twin mate lands in the working set, and in twin refinement the free reflection's calculated intensity then carries the working-set fit. The same set is what phenix.refine generates by default (`use_lattice_symmetry`). It also makes the free set independent of the space group a file is merged in: the merged file, the P1 cross-check and a re-merge in any subgroup carry one free set (nested, where the small-data floor below lifts the fraction of one of them more than another's); - 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 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. diff --git a/image_analysis/scale_merge/Merge.cpp b/image_analysis/scale_merge/Merge.cpp index ce7b87f2e..2945ce2f7 100644 --- a/image_analysis/scale_merge/Merge.cpp +++ b/image_analysis/scale_merge/Merge.cpp @@ -389,7 +389,7 @@ std::vector MergeOnTheFly::ExportReflections() { out.emplace_back(mr); } - AssignRfreeFlags(out, space_group, scaling_settings.GetRfreeFraction()); + AssignRfreeFlags(out, space_group, scaling_settings.GetRfreeFraction(), 500, reference_cell); ApplyFrenchWilson(out, space_group); return out; } diff --git a/image_analysis/scale_merge/RfreeFlags.cpp b/image_analysis/scale_merge/RfreeFlags.cpp index 6b7b17cd2..cfc099d6e 100644 --- a/image_analysis/scale_merge/RfreeFlags.cpp +++ b/image_analysis/scale_merge/RfreeFlags.cpp @@ -8,6 +8,7 @@ #include #include "HKLKey.h" +#include "gemmi/twin.hpp" namespace { // splitmix64 bit-mix of a key -> uniform double in [0, 1). Same key -> same value, so all @@ -20,10 +21,36 @@ namespace { z = z ^ (z >> 31); return static_cast(z >> 11) * (1.0 / 9007199254740992.0); } + + + // The rotations of the lattice's own point group - its holohedry - found from the metric as for + // twin laws (Le Page two-folds, with phenix.xtriage's default obliquity of 3 degrees, so a + // pseudo-merohedral lattice counts as its higher metric; grouping too much only clusters the + // free set, grouping too little splits twin mates). The centring is taken as P: the lattice of + // the cell's own basis vectors, the same whatever group a merge is in, so the merged file, the P1 + // cross-check and any re-merge get one holohedry. On an R lattice in its hexagonal cell that is + // 6/mmm rather than -3m, which groups each present reflection with exactly its -3m mates, the + // extra ones being centring-absent. Empty when the cell does not carry the merging group itself + // (a space group forced on a metric that does not have it): there is no holohedry to key on. + std::vector LatticeHolohedry(const UnitCell &cell, const gemmi::SpaceGroup &space_group) { + std::vector rots; + for (const auto &op : gemmi::find_lattice_symmetry(static_cast(cell), 'P', 3.0).sym_ops) + rots.push_back(gemmi::Op{op.rot, {0, 0, 0}, 'x'}); + for (const auto &op : space_group.operations().sym_ops) + if (std::none_of(rots.begin(), rots.end(), [&](const gemmi::Op &r) { return r.rot == op.rot; })) + return {}; + return rots; + } + + uint64_t PackIndex(const gemmi::Op::Miller &h) { + constexpr int64_t bias = 1 << 20; + return (static_cast(h[0] + bias) << 42) | (static_cast(h[1] + bias) << 21) + | static_cast(h[2] + bias); + } } void AssignRfreeFlags(std::vector &merged, const gemmi::SpaceGroup &space_group, - double rfree_fraction, int min_free_reflections) { + double rfree_fraction, int min_free_reflections, const std::optional &cell) { for (auto &r : merged) r.rfree_flag = false; if (rfree_fraction <= 0.0 || merged.empty()) @@ -36,10 +63,34 @@ void AssignRfreeFlags(std::vector &merged, const gemmi::SpaceG // 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. + // + // The key is the orbit of the reflection under the LATTICE HOLOHEDRY, not under the merging group. + // A twin law is a lattice symmetry that the crystal lacks, so keying on the merging group puts a + // free reflection's twin mate in the working set almost every time (measured: 97-98% of the free + // reflections that have a mate), and a twin-refined R-free then reads the working set through + // I_calc. Keyed on the holohedry, twin mates share a flag, and the free set no longer depends on + // the space group a file is merged in - the merged MTZ, the P1 cross-check and a re-merge in any + // subgroup carry one free set (where the small-data floor below lifts the fraction differently in + // two of them, the smaller set is contained in the larger: each orbit has one draw). As + // phenix.refine does by default (use_lattice_symmetry). const HKLKeyGenerator laue_key(/*merge_friedel=*/true, space_group); + const std::vector holohedry = cell ? LatticeHolohedry(*cell, space_group) : std::vector{}; + auto key_of = [&](const MergedReflection &r) -> uint64_t { + if (holohedry.empty()) + return laue_key(r).pack(); + // The largest index of the orbit, Friedel mates included, names it. + gemmi::Op::Miller best{{r.h, r.k, r.l}}; + for (const auto &op : holohedry) { + const gemmi::Op::Miller h = op.apply_to_hkl({{r.h, r.k, r.l}}); + best = std::max({best, h, gemmi::Op::Miller{{-h[0], -h[1], -h[2]}}}); + } + return PackIndex(best); + }; - // Count the distinct test-eligible reflections (distinct Laue-ASU keys; mates collapse to one) so - // the fraction can be floored to a usable test-set size on small data. + // Count the distinct test-eligible reflections (distinct Laue-ASU keys of the merging group; mates + // collapse to one) so the fraction can be floored to a usable test-set size on small data. Counted + // in the merging group, not in holohedral orbits: R-free is a sum over the reflections of this file, + // and counting orbits would lift the fraction on every dataset below its holohedry. std::unordered_set distinct; distinct.reserve(merged.size()); for (const auto &r : merged) @@ -54,7 +105,7 @@ void AssignRfreeFlags(std::vector &merged, const gemmi::SpaceG const double eff_fraction = std::max(rfree_fraction, floor_fraction); for (auto &r : merged) - r.rfree_flag = UniformFromKey(laue_key(r).pack()) < eff_fraction; + r.rfree_flag = UniformFromKey(key_of(r)) < eff_fraction; } size_t ApplyReferenceFreeFlags(std::vector &merged, const gemmi::SpaceGroup &space_group, diff --git a/image_analysis/scale_merge/RfreeFlags.h b/image_analysis/scale_merge/RfreeFlags.h index fd0402577..6b00f6f7d 100644 --- a/image_analysis/scale_merge/RfreeFlags.h +++ b/image_analysis/scale_merge/RfreeFlags.h @@ -5,15 +5,21 @@ #include #include +#include #include #include "../../common/Reflection.h" +#include "../../common/UnitCell.h" #include "gemmi/symmetry.hpp" -// 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: +// Assign R-free (test-set) flags to merged reflections. With the cell given, the flag is a pure +// function of the reflection's orbit under the LATTICE HOLOHEDRY (the metric symmetry of the cell, +// which contains the merging group); without one, or when the cell does not carry the merging group, +// of its Friedel-merged (Laue) ASU key in space_group. That gives four 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; +// - with the cell, twin-law mates share one too, and every file of one crystal - whatever space +// group it is merged in - carries the same free set; // - 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: @@ -28,7 +34,8 @@ // matters more than cross-dataset identity (and a shared reference/ApplyReferenceFreeFlags is the way // to keep exact identity there). Pass 0 to disable the floor. void AssignRfreeFlags(std::vector &merged, const gemmi::SpaceGroup &space_group, - double rfree_fraction, int min_free_reflections = 500); + double rfree_fraction, int min_free_reflections = 500, + const std::optional &cell = std::nullopt); // 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 diff --git a/image_analysis/scale_merge/RotationScaleMerge.cpp b/image_analysis/scale_merge/RotationScaleMerge.cpp index 39db29683..e651313f3 100644 --- a/image_analysis/scale_merge/RotationScaleMerge.cpp +++ b/image_analysis/scale_merge/RotationScaleMerge.cpp @@ -4077,7 +4077,7 @@ RotationScaleMerge::Result RotationScaleMerge::MergeAndStats(int n_groups, bool // is not written needs neither. French-Wilson is deferred until after the anomalous accumulator // below has attached I(+)/I(-), so the two hands get their amplitudes in one pass. if (full_stats) - AssignRfreeFlags(result.merged, x.GetSpaceGroupOrP1(), rfree_fraction); + AssignRfreeFlags(result.merged, x.GetSpaceGroupOrP1(), rfree_fraction, 500, reference_cell); if (reject_count > 0) logger.Info("Merge outlier rejection: dropped {} observations", reject_count); diff --git a/tests/RfreeFlagsTest.cpp b/tests/RfreeFlagsTest.cpp index bb78a8f8b..b267034e2 100644 --- a/tests/RfreeFlagsTest.cpp +++ b/tests/RfreeFlagsTest.cpp @@ -182,3 +182,35 @@ TEST_CASE("R-free fraction of zero flags nothing", "[rfree]") { AssignRfreeFlags(v, SG(1), 0.0); for (const auto& r : v) CHECK(!r.rfree_flag); } + +TEST_CASE("R-free flags keyed on the lattice holohedry keep twin mates together", "[rfree]") { + // A P4 crystal on a tetragonal cell: the lattice is 4/mmm, and (h,k,l) -> (k,h,-l) is a twin law. + // Keyed on the holohedry, a reflection and its twin mate share a flag, and the P1 merge of the same + // crystal flags every reflection exactly as the P4 merge does. Synthetic cell. + const std::optional cell = UnitCell{47.0f, 47.0f, 63.0f, 90.0f, 90.0f, 90.0f}; + std::vector p4, p1; + for (int h = -10; h <= 10; ++h) + for (int k = -10; k <= 10; ++k) + for (int l = -10; l <= 10; ++l) { + if (h == 0 && k == 0 && l == 0) continue; + const float d = 60.0f / (1 + h * h + k * k + l * l); + p1.push_back(Refl(h, k, l, d)); + p4.push_back(Refl(h, k, l, d)); + } + AssignRfreeFlags(p4, SG(75), 0.05, /*min_free=*/0, cell); + AssignRfreeFlags(p1, SG(1), 0.05, /*min_free=*/0, cell); + + std::map, bool> flag; + for (const auto& r : p4) flag[{r.h, r.k, r.l}] = r.rfree_flag; + int twin_pairs = 0; + for (const auto& r : p4) { + CHECK(flag[{r.k, r.h, -r.l}] == r.rfree_flag); // twin mate + CHECK(flag[{-r.k, r.h, r.l}] == r.rfree_flag); // P4 mate + ++twin_pairs; + } + CHECK(twin_pairs > 0); + for (const auto& r : p1) + CHECK(flag[{r.h, r.k, r.l}] == r.rfree_flag); // one free set whatever the merging group + CHECK(FreeFraction(p4) > 0.03); + CHECK(FreeFraction(p4) < 0.08); +}