From 3171b071e6889ca2f52048d2fc6a476b679ff635 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 27 Jul 2026 21:05:16 +0200 Subject: [PATCH] Space-group search: judge a promotion against its parent, not against the error model The point-group decision moved with the AMOUNT of data at fixed physics: a partially twinned trigonal crystal was promoted into the twin's holohedry whenever the search happened to see a larger first-pass merge, and kept its true subgroup when it saw a smaller one. Simulation over 6 noise draws with only the merge multiplicity varying: the twin is promoted 0/6 at multiplicity 2 and 6/6 at 18, while the genuine control is promoted 6/6 throughout. The cause is that every existing gate is a ratio to the merge error model - b_parent grows toward the true systematic scatter as sigma shrinks with 1/sqrt(N), while b_cand is already saturated by the twin's disagreement, so the ratio slides down through a fixed veto. The parent statistic moves with data amount and the candidate statistic does not. Gate promotions on the operator disagreement H = <|I1-I2|/(I1+I2)> instead, as the ratio of the operators a promotion ADDS to the parent group's own operators on the same reflections. There is no sigma in it, so it cannot drift with the error model, and the parent normalisation cancels data quality. Measured over 27 runs, 5 promotion types and 450-1800 images: genuine symmetry 0.862-1.219, merohedral twins 1.270-2.084. On the synthetic grid it is flat across a 9x change in multiplicity - genuine pinned at 1.00, twins 3-12x the bound - which is precisely the property the old gates lacked. chi^2 and the systematic-b stay as secondary vetoes; they protect against non-crystallographic pseudo-symmetry, which is where correlation-based scoring is weak. Pick the parent carefully: 422 has two maximal subgroups of order 4, and on a tetragonal crystal twinned by 2[100] the rival (222) is CC-confirmed too and CONTAINS the twin laws, so normalising against it hides the twin among the promotion's own real operators (ratio 8.19 against the true parent, 0.78 against the rival). Where several parents tie, judge on the most damning. Also: - Report a refused promotion instead of silently processing lower. Merging a twin in the twin's holohedry averages non-equivalent reflections into each other and cannot be undone from the output; keeping the subgroup costs only redundancy. The refusal names the group and the number that caused it. - Stop the twinning report from arguing in a circle. It ran after adoption and conditioned on the adopted group, so a promotion into a holohedral Laue class made it print "no merohedral twin law exists" - the test was conditioned on the decision it should audit. Twinning is now also measured on the subgroup merge before adoption, and the post-adoption text says when its own conclusion is not authoritative. - Compare PRIMITIVE cell volumes in the first-pass scheme tie-break. A centred setting's cell is an exact integer multiple of its primitive one (a rhombohedral lattice in hexagonal axes is exactly 3x), so the integer-supercell test fired on a pure setting difference and demoted a good scheme to a threefold-smaller merge - which is what let the twin see the small merge to begin with. Rotation battery, 33 crystals: point-group agreement 30/33 -> 29/33, one crystal moved. That crystal (P422 -> P222) is the one with the known unresolved integration defect where reflections near the rotation-axis plane are wildly mis-integrated; its symmetry mates genuinely disagree, and its lower-symmetry merge is measurably better (ISa 2.72 -> 3.63, high-shell CC 75.4 -> 86.0). The threshold was not moved to accommodate it: 1.25 sits inside the measured gap and widening it would admit real twins. Separately the tie-break improved one crystal's CC1/2 from 77.7 to 84.0. Tests: a synthetic twin-fraction x multiplicity grid, which is what the search had never had - the existing tests are noise-free and exercise only Stage B absences. A NOTE ON WHAT WAS TRIED AND REJECTED, so it is not rebuilt: the obvious "physics-anchored" statistic is the disattenuated cross-validated correlation rho = corr(I_half0(h), I_half1(Rh)) / corr(I_half0, I_half1), which is 1 for real symmetry at any data quality and 2a(1-a)/((1-a)^2+a^2) for a twin. It passes the synthetic grid perfectly and FAILS ON REAL DATA IN BOTH DIRECTIONS - five false refusals of genuine symmetry on the battery, and it waves through a twin (rho 0.998) that H refuses. The reason is that cc_half correlates the two halves of the SAME reflection and so measures only random error, while cc_cross compares DIFFERENT reflections carrying different systematic error; dividing by cc_half removes the noise and leaves a systematic floor that varies by crystal AND by operator. Genuine rho measures 0.9987 on strong data and 0.73 on weak. A synthetic generator validates a statistic's arithmetic, never its premise, and this premise - that the only departure from exact symmetry is noise - is false for every real crystal. Any per-operator agreement statistic needs a same-crystal reference; an absolute threshold on one cannot be made to work by tuning. Co-Authored-By: Claude Opus 5 (1M context) --- .../scale_merge/SearchSpaceGroup.cpp | 104 ++++- image_analysis/scale_merge/SearchSpaceGroup.h | 39 ++ .../scale_merge/TwinningAnalysis.cpp | 6 + image_analysis/scale_merge/TwinningAnalysis.h | 5 + rugnux/Rugnux.cpp | 36 +- tests/CMakeLists.txt | 2 + tests/SearchSpaceGroupTwinTest.cpp | 361 ++++++++++++++++++ tests/SyntheticMergedReflections.h | 239 ++++++++++++ 8 files changed, 785 insertions(+), 7 deletions(-) create mode 100644 tests/SearchSpaceGroupTwinTest.cpp create mode 100644 tests/SyntheticMergedReflections.h diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index 4939e945..b0e64d2e 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -11,6 +11,7 @@ #include #include #include +#include #include #include #include @@ -71,6 +72,12 @@ namespace { return false; } + std::string FormatDouble(double v, int decimals) { + std::ostringstream o; + o << std::fixed << std::setprecision(decimals) << v; + return o.str(); + } + std::array RotKey(const gemmi::Op& op) { std::array out{}; for (int i = 0; i < 3; ++i) @@ -265,6 +272,17 @@ SearchSpaceGroupResult SearchSpaceGroup( s.op_triplet_hkl = op.as_hkl().triplet('h'); s.n_pairs = static_cast(x.size()); s.cc = PearsonCC(x, y); + // Sigma-free disagreement over the same pairs (see SpaceGroupOptions::max_operator_h_ratio). + double h_sum = 0.0; + int h_n = 0; + for (size_t p = 0; p < x.size(); ++p) { + const double denom = x[p] + y[p]; + if (denom > 0.0) { + h_sum += std::fabs(x[p] - y[p]) / denom; + ++h_n; + } + } + s.h_stat = h_n > 0 ? h_sum / h_n : 0.0; s.present = s.n_pairs >= opt.min_pairs_per_operator && std::isfinite(s.cc) && s.cc >= opt.min_operator_cc; return s; @@ -432,6 +450,8 @@ SearchSpaceGroupResult SearchSpaceGroup( // Operator-CC-confirmed candidates, each with its merge chi^2 and systematic-error b; chi2_ref = // the most consistent. struct PGCand { const PointGroupInfo* pg; int order; double min_class_cc; double chi2; double b_extra; }; + int refused_order = 0; + std::string refused_pg_hm, refused_why; std::vector pg_cands; double chi2_ref = std::numeric_limits::infinity(); for (const auto& pg : point_groups) { @@ -472,16 +492,59 @@ SearchSpaceGroupResult SearchSpaceGroup( // (imperfectly scaled data) and VETOES a twin whose chi^2 now looks self-consistent but whose b // balloons - the chi^2 ratio alone no longer separates them. double parent_b = -1.0; + const PointGroupInfo *parent_pg = nullptr; + // Every confirmed subgroup of the largest order below this candidate. There can be more than one + // - 422 has both 4 and 222 - and on a twinned crystal the rival is not a harmless alternative: a + // P4 crystal twinned by 2[100] has its two twin 2-folds confirmed, so 222 is CC-confirmed too and + // CONTAINS the twin laws. Normalising the H test against it hides the twin among the promotion's + // own real operators (measured on the synthetic grid: ratio 8.19 against the true parent 4, 0.78 + // against the rival 222). Which one is the true parent is exactly what is unknown here, so the + // promotion must answer to all of them. + std::vector parents; if (!c.pg->rotations.empty()) { int parent_order = 0; for (const auto& s : pg_cands) - if (s.order < c.order && s.order > parent_order + if (s.order < c.order && s.order >= parent_order && std::includes(c.pg->rotation_set.begin(), c.pg->rotation_set.end(), s.pg->rotation_set.begin(), s.pg->rotation_set.end())) { - parent_order = s.order; - parent_b = s.b_extra; + if (s.order > parent_order) { + parent_order = s.order; + parent_b = s.b_extra; + parent_pg = s.pg; + parents.clear(); + } + parents.push_back(s.pg); } } + + // Sigma-free twin test: compare the disagreement H of the operators this promotion ADDS with the + // disagreement of the parent's own operators, measured on the same reflections. A real operator + // relates equal intensities and matches the parent; a twin law relates different ones and reads + // systematically higher. Skipped when either side has too few pairs to mean anything, and when + // there is no parent group to normalise against (the first step out of P1). Where several parents + // tie (see above), the promotion is judged on the most damning of them. + double h_ratio = std::numeric_limits::quiet_NaN(); + for (const auto *parent : parents) { + double h_new = 0.0, h_par = 0.0; + int n_new = 0, n_par = 0, pairs_new = 0, pairs_par = 0; + for (const auto &rot : c.pg->rotations) { + if (rot.rot == gemmi::Op::identity().rot) + continue; + const auto &os = operator_score(rot); + if (os.n_pairs < opt.min_pairs_per_operator) + continue; + const bool in_parent = std::binary_search(parent->rotation_set.begin(), + parent->rotation_set.end(), RotKey(rot)); + if (in_parent) { h_par += os.h_stat; ++n_par; pairs_par += os.n_pairs; } + else { h_new += os.h_stat; ++n_new; pairs_new += os.n_pairs; } + } + if (n_new > 0 && n_par > 0 && pairs_new >= opt.min_pairs_for_h + && pairs_par >= opt.min_pairs_for_h && h_par > 0.0) { + const double r = (h_new / n_new) / (h_par / n_par); + if (!std::isfinite(h_ratio) || r > h_ratio) + h_ratio = r; + } + } // The chi^2 ratio is only trustworthy when the error model is calibrated. When even the best // subgroup's reduced chi^2 (chi2_ref) is far above 1 - weak, low-resolution data whose merged // sigmas are badly under-estimated - the ratio grows with point-group order for genuine high @@ -504,8 +567,29 @@ SearchSpaceGroupResult SearchSpaceGroup( && c.b_extra > std::max(parent_b, opt.min_systematic_b_for_veto) * opt.max_systematic_b_veto) consistent = false; - if (!consistent) + // The H test is a necessary condition for promotion where it can be computed: it is the only + // statistic measured to separate genuine symmetry from a merohedral twin across data amounts. + const bool h_refused = std::isfinite(h_ratio) && h_ratio > opt.max_operator_h_ratio; + if (h_refused) + consistent = false; + + if (!consistent) { + // Record the highest-order refusal so the caller can say WHY it is processing lower. + if (c.order > refused_order && c.pg->representative) { + refused_order = c.order; + refused_pg_hm = c.pg->representative->point_group_hm(); + if (h_refused) + refused_why = "operator disagreement H is " + FormatDouble(h_ratio, 2) + + "x the parent's (bound " + FormatDouble(opt.max_operator_h_ratio, 2) + + ") - the added operator relates unequal intensities, as a twin law does"; + else if (std::isfinite(c.chi2) && std::isfinite(chi2_ref)) + refused_why = "merge chi^2 is " + FormatDouble(c.chi2 / chi2_ref, 2) + + "x the subgroup's (bound " + FormatDouble(opt.max_merge_chi2_ratio, 2) + ")"; + else + refused_why = "the merge under it is not self-consistent"; + } continue; + } if (c.order > best_pg_order || (c.order == best_pg_order && c.min_class_cc > best_pg_min_cc)) { best_pg = c.pg; best_pg_order = c.order; @@ -524,6 +608,12 @@ SearchSpaceGroupResult SearchSpaceGroup( if (best_pg->representative) result.point_group_hm = best_pg->representative->point_group_hm(); + // Only report a refusal that is actually ABOVE what was adopted. + if (refused_order > best_pg_order) { + result.refused_point_group_hm = refused_pg_hm; + result.refused_reason = refused_why; + } + // --- Stage B: pick the space group within the point group --- // Without screw/centering determination, return the symmorphic representative. if (!opt.determine_space_group || best_pg->rotations.empty()) { @@ -674,6 +764,12 @@ std::string SearchSpaceGroupResultToText(const SearchSpaceGroupResult& result, size_t max_candidates_to_print) { std::ostringstream os; + if (!result.refused_point_group_hm.empty()) + os << "Higher symmetry " << result.refused_point_group_hm << " was confirmed by the operator " + "correlations but REFUSED: " << result.refused_reason << ".\n" + " Processing in the lower symmetry, which is the recoverable direction - if this is a " + "twin, merging in the higher group would average non-equivalent reflections together and " + "hide the twin law.\n"; os << "Point group: " << (result.point_group_hm.empty() ? "?" : result.point_group_hm) << " (from intensity correlations)\n"; diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index d8c35a9c..bafa48ca 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -30,6 +30,12 @@ struct SpaceGroupOperatorScore { double cc = 0.0; // correlation of I(h) with I(Rh) int n_pairs = 0; // independent reflection pairs the CC was computed from bool present = false; // operator confirmed as a real symmetry of the intensities + // Mean |I1-I2|/(I1+I2) over this operator's pairs - the disagreement the operator implies, with no + // sigma in it. Unlike chi^2 and the systematic-b, which are ratios to a merge error model that + // drifts with multiplicity, this is a property of the intensities alone; compared against an + // operator already confirmed on the same reflections it is what separates a real symmetry from a + // twin law (see max_operator_h_ratio). + double h_stat = 0.0; }; struct SpaceGroupCandidateScore { @@ -125,6 +131,30 @@ struct SearchSpaceGroupOptions { // is a non-negligible fraction of I; below the floor the increase is treated as noise, not a twin. double min_systematic_b_for_veto = 0.05; + // Promotion gate on the operator disagreement H = <|I1-I2|/(I1+I2)>, taken as the ratio of the + // operators the promotion ADDS to the operators of the parent group already confirmed on the same + // reflections. A real symmetry operator relates equal intensities, so its H matches the parent's + // (ratio ~1); a merohedral twin law relates DIFFERENT reflections mixed in proportion alpha, so its + // H is systematically larger. Measured over 27 runs spanning 5 promotion types and 450-1800 images: + // genuine symmetry 0.862-1.219, merohedral twins 1.270-2.084 - a clean gap, and unlike the chi^2 and + // systematic-b ratios (genuine 1.00-3.47 / 1.09-3.89 vs twin 1.35-3.32 / 1.77-4.76, fully + // interleaved) it does not drift with multiplicity, because there is no sigma in it and the parent + // normalisation cancels data quality. For a partial twin over the twin law is (1-2*alpha)/2, so + // the excess over the parent also estimates the twin fraction rather than being a tuned constant. + // + // The parent normalisation is what makes this work, and it is not optional. Symmetry-related + // reflections never agree exactly on real data - absorption, illumination and partiality differ + // between them - and that systematic floor varies by crystal AND by operator (a cubic 3-fold permutes + // axes, relating far-apart parts of reciprocal space, so it disagrees more than the 2-folds of its + // parent even when the symmetry is perfectly real). Measuring an added operator against the parent's + // own operators, on the same reflections, is what divides that floor out. An ABSOLUTE bound on any + // per-operator agreement statistic cannot: measured absolute values for genuine symmetry span the + // whole range from 0.99 on strong data down to 0.71 on weak, straddling every twin. + double max_operator_h_ratio = 1.25; + + // The H test needs at least this many pairs on both sides to mean anything. + int min_pairs_for_h = 200; + // Above this reduced chi^2 for the best subgroup (chi2_ref), the merged error model is treated as // badly miscalibrated (weak, low-resolution data whose sigmas are far too small): the fixed-sigma // chi^2 ratio then grows with point-group order for genuine high symmetry too and can no longer @@ -176,6 +206,15 @@ struct SearchSpaceGroupResult { std::string point_group_hm; // chosen point group, e.g. "422" std::vector operator_scores; // Stage A, all distinct operators tested std::vector candidates; // Stage B, ranked + + // A HIGHER point group whose operators the intensities confirmed (Stage A) but whose promotion the + // consistency tests refused, with the reason. Processing continues in the lower group, which is the + // safe direction: merging a twinned crystal in the twin's holohedry averages non-equivalent + // reflections into each other and is unrecoverable from the output (and makes the run report that no + // twin law exists), whereas keeping the subgroup costs only redundancy and can be promoted later. + // Empty when nothing was refused. Surfaced to the user - a silent demotion is how a twin gets missed. + std::string refused_point_group_hm; + std::string refused_reason; }; SearchSpaceGroupResult SearchSpaceGroup( diff --git a/image_analysis/scale_merge/TwinningAnalysis.cpp b/image_analysis/scale_merge/TwinningAnalysis.cpp index 92924b46..56e6fbc4 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.cpp +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -217,6 +217,12 @@ std::string TwinningAnalysisToText(const TwinningAnalysisResult& result) { os << " => Twinning suspected (estimated twin fraction ~" << result.estimated_twin_fraction << "). Statistics flag the presence of twinning, not\n" << " the twin law; confirm with a dedicated twin-law analysis.\n"; + else if (!result.merohedral_twinning_possible && result.laue_class_was_chosen_by_promotion) + os << " => Cannot rule out twinning from these numbers: the Laue class is holohedral, so no\n" + << " merohedral twin law exists WITHIN it - but this Laue class was chosen by the\n" + << " space-group search itself, and promoting into a twin's holohedry is precisely what\n" + << " a merohedral twin looks like. Judge the twinning from the subgroup statistics\n" + << " reported by the search, not from these.\n"; else if (!result.merohedral_twinning_possible) os << " => No twinning: the Laue class is holohedral, so no merohedral twin law exists\n" << " (any <|L|> below 0.5 here is a statistical artefact, not twinning).\n"; diff --git a/image_analysis/scale_merge/TwinningAnalysis.h b/image_analysis/scale_merge/TwinningAnalysis.h index 16714749..957b8efa 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.h +++ b/image_analysis/scale_merge/TwinningAnalysis.h @@ -37,6 +37,11 @@ struct TwinningAnalysisResult { // False when the Laue class is holohedral (4/mmm, 6/mmm, m-3m, rhombohedral -3m): no merohedral // twin law can exist, so a low <|L|> / second moment there is a statistical artefact, not twinning. bool merohedral_twinning_possible = true; + + // Set when this analysis ran on a merge whose space group the pipeline CHOSE by promoting past a + // subgroup. The "holohedral Laue class, so no twin law exists" conclusion is then circular - the + // promotion is exactly what a twin would have caused - so the report must not state it as a fact. + bool laue_class_was_chosen_by_promotion = false; }; // The space group (when known) is used to drop centric reflections, which follow different diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index f51c5e14..71a9a68b 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -561,7 +561,14 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b if (!result.has_value()) continue; const int score = count_indexed(idx, *result); - const double vol = std::abs(result->lattice.CalcVolume()); + // Compare PRIMITIVE volumes: two schemes can find the same lattice in different + // settings, and a centred setting's cell is an exact integer multiple of its primitive + // one - a rhombohedral lattice in hexagonal axes is exactly 3x its primitive + // rhombohedral cell. Comparing the centred volumes makes the supercell test below fire + // on that pair and "demote" a perfectly good setting to a threefold-smaller merge, which + // is enough to change the space group the search then picks. + const double vol = std::abs( + result->lattice.ToPrimitive(result->search_result.centering).CalcVolume()); const std::string &name = schemes[i].first; logger.Info("First-pass scheme '{}': indexes {}/{} validation frames", name, score, static_cast(validation.size())); @@ -575,8 +582,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // whichever ran first. When the two schemes tie on frames but their cell volumes are // related by an integer factor >=2, the larger cell is that spurious supercell and the // smaller is the true reduced cell: take it, regardless of scheme order. A near-integer - // ratio separates a real axis multiplication from a centering coincidence (a rhombohedral - // H cell vs its C2 sub-cell is 1.5x, non-integer, and is correctly left alone). + // ratio separates a real axis multiplication from a centering coincidence. Volumes are + // primitive (see above), so a pure setting difference is a ratio of 1 and never fires. bool integer_subcell = false; if (bp.result.has_value() && !clearly_more && vol > 1.0 && bp.vol > 1.0) { const bool tied = static_cast(score) >= bp.score * 0.9f - 0.5f; @@ -1128,6 +1135,10 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // test): its conventional setting does not match the indexer's primitive frame, so the two-pass must // NOT reuse it - the second pass re-searches (and re-runs the same test) instead. bool sg_reindexed = false; + // True when the point group below was chosen BY THE SEARCH (rather than given by the user) and is + // above triclinic - i.e. a promotion was made, which is what makes a later "the Laue class is + // holohedral so there is no twin law" conclusion circular. + bool promoted_point_group = false; if (!experiment_.GetGemmiSpaceGroup().has_value()) { SearchSpaceGroupOptions sg_opts; sg_opts.merge_friedel = experiment_.GetScalingSettings().GetMergeFriedel(); @@ -1139,6 +1150,22 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b sg_opts.lattice_system = end_msg.rotation_lattice_type->crystal_system; auto sg_search = SearchSpaceGroup(sm.merged, sg_opts); + // Twinning evidence measured BEFORE any promotion, on the P1/subgroup merge the search was + // given. Reported alongside the post-adoption analysis: once a higher Laue class has been + // adopted, its own twinning test can only say "no twin law exists within this class", which + // is circular - the promotion is what a twin would have caused. These numbers are not. + const auto pre_promotion_twinning = AnalyzeTwinning(sm.merged, nullptr); + if (pre_promotion_twinning.twinning_suspected) + logger.Warning("Twinning indicators BEFORE the space-group decision (subgroup merge): " + "<|L|> = {:.3f}, second moment = {:.3f} (untwinned 0.500 / 2.00) - if the " + "search promoted the point group, treat the promotion with suspicion", + pre_promotion_twinning.mean_abs_l, pre_promotion_twinning.second_moment); + promoted_point_group = !sg_search.point_group_hm.empty() && sg_search.point_group_hm != "1"; + if (!sg_search.refused_point_group_hm.empty()) + logger.Warning("Higher symmetry {} was confirmed by the operator correlations but refused: " + "{}. Processing in the lower symmetry (the recoverable direction).", + sg_search.refused_point_group_hm, sg_search.refused_reason); + // Miller-index reindex under a change of basis: (h,k,l)_conv = reindex * (h,k,l)_prim. const auto reindex_hkl = [](auto &r, const gemmi::Mat33 &m) { const double h = r.h, k = r.k, l = r.l; @@ -1338,6 +1365,9 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b const gemmi::SpaceGroup *twin_sg = twin_sg_number ? gemmi::find_spacegroup_by_number(twin_sg_number.value()) : nullptr; result.twinning = AnalyzeTwinning(sm.merged, twin_sg); + // Mark the conclusion as non-authoritative when the Laue class was reached by a promotion the + // search itself made, so the text cannot claim "no twin law exists" on its own say-so. + result.twinning.laue_class_was_chosen_by_promotion = promoted_point_group; stats_text << TwinningAnalysisToText(result.twinning) << "\n"; // Indexing-ambiguity (alternative-indexing) advisory. When the lattice metric symmetry exceeds diff --git a/tests/CMakeLists.txt b/tests/CMakeLists.txt index 9145a49b..12fc070d 100644 --- a/tests/CMakeLists.txt +++ b/tests/CMakeLists.txt @@ -75,6 +75,8 @@ ADD_EXECUTABLE(jfjoch_test HKLKeyTest.cpp TCPImagePusherTest.cpp SearchSpaceGroupTest.cpp + SearchSpaceGroupTwinTest.cpp + SyntheticMergedReflections.h XDSPluginTest.cpp MergeScaleTest.cpp RfreeFlagsTest.cpp diff --git a/tests/SearchSpaceGroupTwinTest.cpp b/tests/SearchSpaceGroupTwinTest.cpp new file mode 100644 index 00000000..43f7402b --- /dev/null +++ b/tests/SearchSpaceGroupTwinTest.cpp @@ -0,0 +1,361 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#include + +#include "../image_analysis/scale_merge/SearchSpaceGroup.h" +#include "SyntheticMergedReflections.h" +#include "gemmi/symmetry.hpp" + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +// Point-group decision on MEROHEDRALLY TWINNED data. +// +// SearchSpaceGroupTest.cpp exercises Stage B (systematic absences) on noise-free, exactly-symmetric +// intensities. The decision that actually goes wrong on real crystals is Stage A: whether the extra +// operator of the metric holohedry is a real symmetry or a twin law. That decision runs through the +// merge chi^2 gate, its systematic-b rescue and the systematic-b veto - none of which the noise-free +// set can reach, because it has no errors for a chi^2 to be reduced by. +// +// Three crystals are modelled per lattice, all with the SAME metric symmetry (a merohedral twin has +// the supergroup's metric, so the lattice cannot arbitrate - only the intensities can): +// * genuine supergroup - structure factors invariant under the supergroup; must be promoted; +// * untwinned subgroup - structure factors invariant under the subgroup only; must NOT be; +// * twinned subgroup - the same crystal at twin fraction alpha; must NOT be promoted for any +// 0 <= alpha < 0.5, because promoting averages the two twin domains into +// one intensity and the twin is then unrecoverable downstream. +// At alpha = 0.5 the twin is physically indistinguishable from real symmetry, so only "terminates +// and returns one of the two" is asserted. +// +// Each is measured through two merge-quality regimes and five merge multiplicities. Neither changes +// any physics - they change only how well the same crystal was measured and how honest its sigmas +// are - so no decision above may move with them. That is the property the harness exists to pin +// down; every case below asserts it outright. + +namespace { + using jfjoch_test::SyntheticMergeParams; + + struct TwinCrystal { + std::string name; + std::string sub; // the crystal's true space group when twinned + std::string super; // supergroup of index 2; its extra operator is the twin law + gemmi::CrystalSystem system; // metric (lattice) symmetry, as rugnux passes it from indexing + // The OTHER maximal subgroup of the supergroup of the same order as `sub`, when one exists. + // 422 has two - 4 and 222 - and only one of them is the crystal. Which one a parent-normalised + // statistic divides by decides the promotion, so the harness reports both; empty when the + // supergroup has only one maximal subgroup of that order (32 over 3) and the choice cannot arise. + std::string rival_parent; + }; + + const std::vector crystals = { + {"trigonal 3 -> 32 (R3 / R32, twin law k,h,-l)", "R 3 :H", "R 32 :H", + gemmi::CrystalSystem::Trigonal, ""}, + {"tetragonal 4 -> 422 (P4 / P422, twin law h,-k,-l)", "P 4", "P 4 2 2", + gemmi::CrystalSystem::Tetragonal, "P 2 2 2"}, + }; + + // How honest the merged sigmas are. Both regimes are ways a fitted error model misses in + // practice, and neither is a property of the crystal's symmetry. + struct MergeQuality { + std::string name; + double sigma_miscalibration; + double error_model_b; + std::optional true_systematic_b; + }; + + const std::vector merge_quality = { + // The usual case: merged sigmas come out ~1.7x too small across the board, with the error + // model's b matching the systematic scatter that is actually there. + {"sigmas 1.7x too small", 1.7, 0.05, std::nullopt}, + // The other way a fitted error model misses: the statistical sigmas come out somewhat too + // LARGE while b - the asymptotic I/sigma ceiling, ISa = 1/b - is fitted 3x too optimistic, so + // the systematic scatter present is 3x what the merged sigmas' floor admits. + {"ISa 3x too optimistic", 0.6, 0.02, 0.06}, + }; + + const std::vector multiplicities = {2, 3, 6, 9, 18}; + const std::vector twin_fractions = {0.0, 0.05, 0.10, 0.20, 0.35, 0.50}; + + struct Decision { + std::string point_group; + std::string space_group; + size_t n_merged = 0; + std::vector operators; + std::string report; + }; + + // true_group is the symmetry the structure factors have: the subgroup for the twin series, the + // supergroup for the genuine-high-symmetry control (where the twin law is a real symmetry + // operator, so the twin fraction has no effect). + Decision Decide(const TwinCrystal& c, const MergeQuality& q, + const std::string& true_group, double alpha, int multiplicity) { + SyntheticMergeParams p; + p.true_space_group = true_group; + p.twin_supergroup = c.super; + p.twin_fraction = alpha; + p.multiplicity = multiplicity; + p.sigma_miscalibration = q.sigma_miscalibration; + p.error_model_b = q.error_model_b; + p.true_systematic_b = q.true_systematic_b; + // One fixed seed for the whole harness: every case draws the same unit-normal stream, so two + // cases differ only in the knob under test and every test is reproducible. + p.seed = 20260727; + + const auto merged = jfjoch_test::GenerateSyntheticMerged(p); + REQUIRE(merged.size() > 5000); // a realistic dataset, not a handful of reflections + + SearchSpaceGroupOptions opt; // as rugnux/Rugnux.cpp sets it + opt.merge_friedel = true; + opt.lattice_system = c.system; + + const auto result = SearchSpaceGroup(merged, opt); + + Decision d; + d.point_group = result.point_group_hm; + d.space_group = result.best_space_group.has_value() + ? result.best_space_group->short_name() : "none"; + d.n_merged = merged.size(); + d.operators = result.operator_scores; + d.report = SearchSpaceGroupResultToText(result); + return d; + } + + std::string PointGroupOf(const std::string& space_group_name) { + return gemmi::get_spacegroup_by_name(space_group_name).point_group_hm(); + } + + std::string ShortNameOf(const std::string& space_group_name) { + return gemmi::get_spacegroup_by_name(space_group_name).short_name(); + } + + // The hkl triplets SearchSpaceGroup labels a space group's own rotations with. Lets a test tell the + // crystal's real symmetry operators from the twin law among result.operator_scores. + std::set OperatorTripletsOf(const std::string& space_group_name) { + std::set out; + const auto& sg = gemmi::get_spacegroup_by_name(space_group_name); + for (const auto& op : sg.operations().derive_symmorphic().sym_ops) { + if (op.rot == gemmi::Op::identity().rot) + continue; + out.insert(gemmi::Op{op.rot, {0, 0, 0}, op.notation}.as_hkl().triplet('h')); + } + return out; + } + + // The statistic the promotion is actually decided on: the mean operator disagreement + // H = <|I1-I2|/(I1+I2)> over the operators the promotion ADDS, divided by the mean over the parent + // group's own operators, measured on the same reflections. Mirrors what SearchSpaceGroup computes + // for the sub -> super step, so a test can report the margin the max_operator_h_ratio bound has. + // + // The parent normalisation is the whole design, not a detail. An ABSOLUTE per-operator bound cannot + // work: two reflections related by a real symmetry operator still disagree, because they carry + // DIFFERENT systematic error - absorption, illumination, partiality - and how much of that a + // crystal has is a property of the measurement, not of its symmetry. So a genuine operator's own + // disagreement ranges over whatever the data quality happens to be, and any fixed bound placed on + // it rejects good crystals at one end or waves twins through at the other. Dividing by the parent + // operators - already confirmed, measured on the same reflections, carrying the same systematic + // floor - cancels the data quality and leaves only the question being asked: does the ADDED + // operator relate intensities as equal as the parent's do (real symmetry), or systematically less + // equal (a twin law mixing non-equivalent reflections)? + double HRatioOfPromotion(const std::vector& operators, + const std::string& parent_group) { + const auto parent_ops = OperatorTripletsOf(parent_group); + double h_added = 0.0, h_parent = 0.0; + int n_added = 0, n_parent = 0; + for (const auto& s : operators) { + if (s.n_pairs < 200) // SearchSpaceGroupOptions::min_pairs_for_h + continue; + if (parent_ops.count(s.op_triplet_hkl) > 0) { h_parent += s.h_stat; ++n_parent; } + else { h_added += s.h_stat; ++n_added; } + } + if (n_added == 0 || n_parent == 0 || h_parent <= 0.0) + return std::numeric_limits::quiet_NaN(); + return (h_added / n_added) / (h_parent / n_parent); + } + + std::string Describe(const TwinCrystal& c, const MergeQuality& q, double alpha, int multiplicity) { + std::ostringstream os; + os << c.name << ", " << q.name << ", twin fraction " << std::fixed << std::setprecision(2) + << alpha << ", multiplicity " << multiplicity; + return os.str(); + } +} + +// Positive control: a crystal whose structure factors really do have the higher symmetry must be +// promoted to it. Guards the twin tests below against a criterion that simply never promotes. +TEST_CASE("SearchSpaceGroup promotes a genuinely high-symmetry crystal", + "[SearchSpaceGroup][twin]") { + for (const auto& c : crystals) + for (const auto& q : merge_quality) + for (int mult : multiplicities) { + DYNAMIC_SECTION(c.name << ", " << q.name << ", genuine supergroup, multiplicity " << mult) { + const auto d = Decide(c, q, c.super, 0.0, mult); + INFO(d.report); + CHECK(d.point_group == PointGroupOf(c.super)); + CHECK(d.space_group == ShortNameOf(c.super)); + } + } +} + +// Negative control: an UNTWINNED crystal of the true subgroup (alpha = 0) must not be promoted - its +// extra metric operator relates reflections that are simply not equivalent. +TEST_CASE("SearchSpaceGroup keeps an untwinned low-symmetry crystal in its subgroup", + "[SearchSpaceGroup][twin]") { + for (const auto& c : crystals) + for (const auto& q : merge_quality) + for (int mult : multiplicities) { + DYNAMIC_SECTION(Describe(c, q, 0.0, mult)) { + const auto d = Decide(c, q, c.sub, 0.0, mult); + INFO(d.report); + CHECK(d.point_group == PointGroupOf(c.sub)); + CHECK(d.space_group == ShortNameOf(c.sub)); + } + } +} + +// A partial merohedral twin must stay in its true subgroup. Promoting it averages the two twin +// domains into one intensity, which no later step can undo: the twin fraction is not recoverable and +// the merged data are simply wrong. +TEST_CASE("SearchSpaceGroup keeps a partially twinned crystal in its true subgroup", + "[SearchSpaceGroup][twin]") { + for (const auto& c : crystals) + for (const auto& q : merge_quality) + for (double alpha : {0.05, 0.10, 0.20, 0.35}) + for (int mult : multiplicities) { + DYNAMIC_SECTION(Describe(c, q, alpha, mult)) { + const auto d = Decide(c, q, c.sub, alpha, mult); + INFO(d.report); + CHECK(d.point_group == PointGroupOf(c.sub)); + CHECK(d.space_group == ShortNameOf(c.sub)); + } + } +} + +// A PERFECT (alpha = 0.5) merohedral twin produces intensities that are exactly invariant under the +// twin law: I_obs(h) = I_obs(twin h) for every reflection. No intensity statistic can tell it from a +// crystal that genuinely has the higher symmetry - the information is not in the data (it takes a +// different measurement, e.g. the |E| distribution's second moment, to even suspect it). So the only +// thing asserted here is that the search terminates and returns one of the two. +TEST_CASE("SearchSpaceGroup on a perfect merohedral twin returns one of the two symmetries", + "[SearchSpaceGroup][twin]") { + for (const auto& c : crystals) + for (const auto& q : merge_quality) + for (int mult : multiplicities) { + DYNAMIC_SECTION(Describe(c, q, 0.5, mult)) { + const auto d = Decide(c, q, c.sub, 0.5, mult); + INFO(d.report); + REQUIRE(d.space_group != "none"); + CHECK((d.point_group == PointGroupOf(c.sub) || + d.point_group == PointGroupOf(c.super))); + } + } +} + +// THE property this harness exists for. Multiplicity changes only the sigmas - the random part of a +// merged sigma averages down as 1/sqrt(n) while the systematic floor b*|I| does not - so it changes +// how well the SAME crystal is measured, never what its symmetry is. A symmetry decision that moves +// when the same crystal is merged 2x instead of 18x is a defect of the criterion, not a property of +// the data. +// +// This is what a criterion thresholded on merge chi^2 or on a systematic-b RATIO cannot deliver: both +// are ratios to an error model that multiplicity and the sigma calibration move, so the tetragonal +// 4 -> 422 twin at alpha 0.20 / 0.35 used to flip - promoted at multiplicity 2 and kept at 18 with +// under-calibrated sigmas, and the other way round with an over-optimistic ISa. The operator +// disagreement ratio H_added/H_parent holds instead because there is no sigma in it at all: it +// compares intensities with intensities, and normalising against the parent group's own operators on +// the same reflections divides out both the data quality and the systematic floor that multiplicity +// and the error model move. An absolute bound on a single operator's H would not survive this - see +// HRatioOfPromotion. +TEST_CASE("SearchSpaceGroup point-group decision does not depend on merge multiplicity", + "[SearchSpaceGroup][twin]") { + for (const auto& c : crystals) + for (const auto& q : merge_quality) + for (double alpha : twin_fractions) { + DYNAMIC_SECTION(Describe(c, q, alpha, 2) + " vs multiplicity 18") { + const auto low = Decide(c, q, c.sub, alpha, 2); + const auto high = Decide(c, q, c.sub, alpha, 18); + INFO("multiplicity 2:\n" << low.report << "\nmultiplicity 18:\n" << high.report); + CHECK(low.point_group == high.point_group); + CHECK(low.space_group == high.space_group); + } + } +} + +// Diagnostic, not run by default: ./jfjoch_test "[twin-h]" +// Prints the operator-disagreement ratio H_added/H_parent that the promotion is decided on, for the +// genuine high-symmetry crystal and for each twin fraction, across both merge-quality regimes and +// every multiplicity - i.e. how much margin the max_operator_h_ratio bound actually has, and whether +// either side of it drifts with data quality or data amount. +TEST_CASE("SearchSpaceGroup operator H ratio margins", "[.][twin-h]") { + SearchSpaceGroupOptions defaults; + std::cout << "H_added / H_parent for the sub -> super promotion; bound " + << defaults.max_operator_h_ratio << " (above = refused as a twin)\n"; + for (const auto& c : crystals) + for (const auto& q : merge_quality) + // Both normalisations where the supergroup has two maximal subgroups of the same order: + // against the crystal's true parent, and against its rival. + for (const auto& parent : c.rival_parent.empty() + ? std::vector{c.sub} + : std::vector{c.sub, c.rival_parent}) { + std::cout << "\n" << c.name << "\n merge quality: " << q.name + << " normalised against " << ShortNameOf(parent) + << (parent == c.sub ? " (the crystal's own parent)" : " (the RIVAL parent)") + << "\n"; + std::cout << " " << std::setw(22) << std::left << "true symmetry / alpha" << std::right; + for (int mult : multiplicities) + std::cout << std::setw(12) << ("mult " + std::to_string(mult)); + std::cout << "\n " << std::setw(22) << std::left << "genuine supergroup" << std::right; + for (int mult : multiplicities) + std::cout << std::setw(12) << std::fixed << std::setprecision(3) + << HRatioOfPromotion(Decide(c, q, c.super, 0.0, mult).operators, parent); + std::cout << "\n"; + for (double alpha : twin_fractions) { + std::ostringstream label; + label << "subgroup, alpha " << std::fixed << std::setprecision(2) << alpha; + std::cout << " " << std::setw(22) << std::left << label.str() << std::right; + for (int mult : multiplicities) + std::cout << std::setw(12) << std::fixed << std::setprecision(3) + << HRatioOfPromotion(Decide(c, q, c.sub, alpha, mult).operators, parent); + std::cout << "\n"; + } + } + SUCCEED(); +} + +// Diagnostic, not run by default (hidden by the [.] tag): +// ./jfjoch_test "[twin-table]" +// prints the decision for every (true symmetry, twin fraction, multiplicity) combination in both +// merge-quality regimes - the table a redesign of the point-group criterion should be judged against. +TEST_CASE("SearchSpaceGroup twin decision table", "[.][twin-table]") { + for (const auto& c : crystals) + for (const auto& q : merge_quality) { + const auto reference = Decide(c, q, c.super, 0.0, 6); + std::cout << "\n" << c.name << "\n merge quality: " << q.name + << " (true subgroup " << ShortNameOf(c.sub) + << ", supergroup " << ShortNameOf(c.super) << ", " + << reference.n_merged << " merged reflections)\n"; + std::cout << " " << std::setw(22) << std::left << "true symmetry / alpha" << std::right; + for (int mult : multiplicities) + std::cout << std::setw(12) << ("mult " + std::to_string(mult)); + std::cout << "\n " << std::setw(22) << std::left << "genuine supergroup" << std::right; + for (int mult : multiplicities) + std::cout << std::setw(12) << Decide(c, q, c.super, 0.0, mult).space_group; + std::cout << "\n"; + for (double alpha : twin_fractions) { + std::ostringstream label; + label << "subgroup, alpha " << std::fixed << std::setprecision(2) << alpha; + std::cout << " " << std::setw(22) << std::left << label.str() << std::right; + for (int mult : multiplicities) + std::cout << std::setw(12) << Decide(c, q, c.sub, alpha, mult).space_group; + std::cout << "\n"; + } + } + SUCCEED(); +} diff --git a/tests/SyntheticMergedReflections.h b/tests/SyntheticMergedReflections.h new file mode 100644 index 00000000..cc34e58d --- /dev/null +++ b/tests/SyntheticMergedReflections.h @@ -0,0 +1,239 @@ +// SPDX-FileCopyrightText: 2026 Filip Leonarski, Paul Scherrer Institute +// SPDX-License-Identifier: GPL-3.0-only + +#pragma once + +#include +#include +#include +#include +#include +#include +#include +#include +#include + +#include "../common/Reflection.h" +#include "gemmi/symmetry.hpp" +#include "gemmi/unitcell.hpp" + +// Synthetic merged intensities for the space-group / point-group search tests. +// +// The set produced here is what SearchSpaceGroup is fed in production: a P1 merge (one entry per +// Friedel-canonical hkl) with intensities, sigmas and half-set intensities. Unlike the noise-free +// set in SearchSpaceGroupTest.cpp it models the three things the point-group decision actually +// depends on: +// +// * a merohedral TWIN - the crystal's true point group is a subgroup of index 2 of the metric +// (lattice) point group, and the twin law is the operator that separates them: +// I_obs(h) = (1-alpha) I_true(h) + alpha I_true(twin h) +// I_true is a function of the TRUE (sub)group's asu, so the subgroup symmetry is exact and only +// the extra supergroup operator is broken - by (1-2 alpha) (I_true(h) - I_true(twin h)). At +// alpha = 0.5 the two are identical and the twin is indistinguishable from real symmetry. +// Setting the true group to the SUPERgroup instead gives the untwinned high-symmetry control. +// +// * the MERGE MULTIPLICITY, modelled the way the real merge behaves (Merge.h, +// SigmaWithSystematicFloor): the random part of the merged sigma averages down as +// 1/sqrt(multiplicity) while the systematic part (b*I - absorption, partiality, beam flicker; +// correlated across a reflection's repeats) does not, so the merged sigma is +// max(sigma_statistical, b*|I|). Multiplicity therefore changes the sigmas but NOT the physics, +// and no symmetry decision may depend on it. +// +// * an ERROR-MODEL MISCALIBRATION - real merged sigmas come out under-estimated (~1.7x), which is +// what pushes the merge's reduced chi^2 to ~3 and switches SearchSpaceGroup between its +// chi^2-ratio and systematic-b regimes. +// +// Intensities follow a Wilson (exponential) distribution with a resolution fall-off, so they span a +// realistic dynamic range; the "structure factor" is a hash of the asu index, not a draw from the +// RNG stream, so the same crystal is reproduced bit-for-bit whatever the multiplicity or the twin +// fraction and two runs differ only in the knob under test. + +namespace jfjoch_test { + + struct SyntheticMergeParams { + // Symmetry the structure factors actually have: the generated intensities are exactly + // invariant under it, whatever the twin fraction. + std::string true_space_group = "R 3 :H"; + + // Supergroup of index 2 over true_space_group; its extra operator is the twin law. Set it + // equal to true_space_group to model a crystal that GENUINELY has the higher symmetry - + // there is then no extra operator, twinning by a real symmetry operator is a no-op, and the + // twin fraction has no effect. + std::string twin_supergroup = "R 32 :H"; + + // Merohedral twin fraction alpha in [0, 0.5]. 0 = untwinned, 0.5 = perfect twin (whose + // intensities are exactly invariant under the twin law, hence indistinguishable from a + // crystal that really has the supergroup symmetry). + double twin_fraction = 0.0; + + // Number of observations merged into each reflection; at least 2, so both half-sets exist. + int multiplicity = 6; + + // Error model. error_model_b is the b the merge FITTED and floors its merged sigmas with + // (ISa = 1/b); true_systematic_b is the intensity-proportional systematic scatter actually + // present in the data - unset means the two agree, i.e. a perfectly fitted error model. + // sigma_miscalibration is how many times too SMALL the merged sigmas come out overall + // (> 1 = under-estimated, the usual case; < 1 = over-estimated). + double error_model_b = 0.05; + std::optional true_systematic_b; + double sigma_miscalibration = 1.7; + + // Intensity distribution: mean intensity at infinite resolution, Wilson B fall-off, and the + // background variance that keeps sigma finite for near-zero (systematically absent) + // reflections. + double mean_intensity = 8000.0; + double wilson_b_A2 = 25.0; + double background_variance = 400.0; + + double d_min_A = 2.5; + + uint32_t seed = 20260727; + }; + + namespace detail { + inline uint64_t Mix64(uint64_t x) { + x ^= x >> 33; x *= 0xff51afd7ed558ccdULL; + x ^= x >> 33; x *= 0xc4ceb9fe1a85ec53ULL; + x ^= x >> 33; return x; + } + + // Uniform in (0, 1), a pure function of the Miller index. + inline double UniformFromHkl(const gemmi::Op::Miller& hkl) { + const uint64_t x = Mix64(static_cast(hkl[0] + 512) * 0x9e3779b97f4a7c15ULL ^ + Mix64(static_cast(hkl[1] + 512)) ^ + (Mix64(static_cast(hkl[2] + 512)) << 1)); + return (static_cast(x >> 11) + 0.5) * (1.0 / 9007199254740992.0); + } + + inline std::vector> RotationSetOf(const gemmi::SpaceGroup& sg) { + std::vector> out; + for (const auto& op : sg.operations().derive_symmorphic().sym_ops) { + std::array rot{}; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + rot[i * 3 + j] = op.rot[i][j]; + out.push_back(rot); + } + std::sort(out.begin(), out.end()); + return out; + } + } + + // The twin law: a rotation of the supergroup that is not in the subgroup. Any of them gives the + // same twinned intensities (an index-2 subgroup is normal, so the coset members differ by a + // subgroup operator, which leaves I_true unchanged), so the first one found is used. + inline gemmi::Op TwinLaw(const gemmi::SpaceGroup& sub, const gemmi::SpaceGroup& super) { + const auto sub_rots = detail::RotationSetOf(sub); + for (const auto& op : super.operations().derive_symmorphic().sym_ops) { + std::array rot{}; + for (int i = 0; i < 3; ++i) + for (int j = 0; j < 3; ++j) + rot[i * 3 + j] = op.rot[i][j]; + if (!std::binary_search(sub_rots.begin(), sub_rots.end(), rot)) + return gemmi::Op{op.rot, {0, 0, 0}, op.notation}; + } + return gemmi::Op::identity(); + } + + // A cell consistent with the space group's crystal system. Synthetic throughout - the tests must + // not carry the cell of any real sample. + inline gemmi::UnitCell SyntheticCellFor(const gemmi::SpaceGroup& sg) { + switch (sg.crystal_system()) { + case gemmi::CrystalSystem::Triclinic: return {33, 37, 41, 85, 95, 105}; + case gemmi::CrystalSystem::Monoclinic: return {37, 43, 51, 90, 101, 90}; + case gemmi::CrystalSystem::Orthorhombic: return {37, 43, 51, 90, 90, 90}; + case gemmi::CrystalSystem::Tetragonal: return {47, 47, 63, 90, 90, 90}; + case gemmi::CrystalSystem::Trigonal: + case gemmi::CrystalSystem::Hexagonal: return {51, 51, 71, 90, 90, 120}; + case gemmi::CrystalSystem::Cubic: return {57, 57, 57, 90, 90, 90}; + } + return {50, 50, 50, 90, 90, 90}; + } + + inline std::vector GenerateSyntheticMerged(const SyntheticMergeParams& p) { + const gemmi::SpaceGroup& sub = gemmi::get_spacegroup_by_name(p.true_space_group); + const gemmi::SpaceGroup& super = gemmi::get_spacegroup_by_name(p.twin_supergroup); + const gemmi::Op twin = TwinLaw(sub, super); + const gemmi::UnitCell cell = SyntheticCellFor(sub); + const gemmi::GroupOps gops = sub.operations(); + const gemmi::ReciprocalAsu rasu(&sub); + + // True (untwinned) intensity: Wilson-distributed |F|^2 of the subgroup asu, with a + // resolution fall-off. Systematically absent reflections (here: the lattice centering) carry + // no intensity - they are what Stage B confirms the centering from. + auto true_intensity = [&](const gemmi::Op::Miller& hkl) -> double { + if (gops.is_systematically_absent(hkl)) + return 0.0; + const auto asu = rasu.to_asu_sign(hkl, gops).first; + const double e_squared = -std::log(detail::UniformFromHkl(asu)); // mean 1, exponential + const double d = cell.calculate_d(hkl); + return p.mean_intensity * e_squared * std::exp(-p.wilson_b_A2 / (2.0 * d * d)); + }; + + // Half-set split of the multiplicity (n0 >= n1); both halves see the same systematic error. + const int n_obs = std::max(2, p.multiplicity); + const int n_half[2] = {(n_obs + 1) / 2, n_obs / 2}; + + std::mt19937 rng(p.seed); + std::normal_distribution gauss(0.0, 1.0); + + const int hmax = static_cast(std::ceil(cell.a / p.d_min_A)) + 1; + const int kmax = static_cast(std::ceil(cell.b / p.d_min_A)) + 1; + const int lmax = static_cast(std::ceil(cell.c / p.d_min_A)) + 1; + + std::vector merged; + + for (int h = -hmax; h <= hmax; ++h) + for (int k = -kmax; k <= kmax; ++k) + for (int l = -lmax; l <= lmax; ++l) { + // One entry per Friedel pair, matching the Friedel-merged P1 set the search gets. + if (std::make_tuple(h, k, l) <= std::make_tuple(-h, -k, -l)) + continue; + const gemmi::Op::Miller hkl{{h, k, l}}; + const double d = cell.calculate_d(hkl); + if (!(d >= p.d_min_A)) + continue; + + const auto twinned = twin.apply_to_hkl(hkl); + const double i_obs = (1.0 - p.twin_fraction) * true_intensity(hkl) + + p.twin_fraction * true_intensity(twinned); + + // Statistical error of one observation, and of the merge of n of them. + const double sigma_one = std::sqrt(i_obs + p.background_variance); + // Systematic error: a property of the reflection, identical in every observation + // of it, so it survives the merge - this is what the b*|I| sigma floor models. + const double systematic = + p.true_systematic_b.value_or(p.error_model_b) * i_obs * gauss(rng); + + MergedReflection r; + r.h = h; + r.k = k; + r.l = l; + r.d = static_cast(d); + + double sum_n_i = 0.0; + for (int half = 0; half < 2; ++half) { + const double sigma_stat_half = + sigma_one / std::sqrt(static_cast(n_half[half])); + const double i_half = i_obs + systematic + sigma_stat_half * gauss(rng); + r.I_half[half] = static_cast(i_half); + r.sigma_half[half] = static_cast( + std::max(sigma_stat_half, p.error_model_b * std::abs(i_half)) / + p.sigma_miscalibration); + sum_n_i += n_half[half] * i_half; + } + + const double i_merged = sum_n_i / n_obs; + const double sigma_stat = sigma_one / std::sqrt(static_cast(n_obs)); + r.I = static_cast(i_merged); + // Merge.h SigmaWithSystematicFloor, then thrown off by the error-model + // miscalibration. + r.sigma = static_cast( + std::max(sigma_stat, p.error_model_b * std::abs(i_merged)) / p.sigma_miscalibration); + + merged.push_back(r); + } + + return merged; + } +}