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; + } +}