From f6bfa546244586f58a3026e7377f21fe28b53f6d Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 26 Aug 2026 00:27:13 +0200 Subject: [PATCH] Space-group search: score operators on normalised intensities Symmetry operators were scored by Pearson correlation on raw merged intensities. Both members of a symmetry pair sit at the same |s|, so the resolution fall-off is variance the two arms share exactly, and it inflates the correlation of true and false operators alike. The clearest demonstration is the test this commit adds: a synthetic data set with no symmetry at all - a radial fall-off times an independent per- reflection factor - is assigned point group 432 by the shipped code, with all 23 rotations confirmed at 0.632 to 0.656. The existing suite passes identically before and after, because nothing covered this. The new case fails 64 of its 100 assertions on the old scoring and passes on the new. On real data the same effect had the gate leaking: on one cubic crystal three pseudo-symmetric operators scored 0.506 to 0.517, above the 0.5 bound, so they were confirmed and 432 had to be refused further downstream by the twin-law and systematic-b guards. Scored on E-squared they read 0.283 to 0.298 and exactly the eleven genuine rotations of 23 are confirmed. The normalisation has to be over the reflections the correlation actually pairs. Reusing the existing normalised array is worse than doing nothing: it is normalised over the set the absence tests use, whose surviving fraction is itself resolution-dependent, and the coupling to the resolution cut rises from 0.086 to 0.262 against a raw baseline of 0.086. Normalised over the paired set it falls to 0.023. The twin-law H statistic keeps its own vectors on raw intensities. It shares the pair arrays with the correlation, and normalising in place moves it by up to 12% against a bound whose window is 5.5% wide. Verified rather than assumed: two instrumented binaries print the same H to twelve significant figures while the correlation differs. min_operator_cc goes 0.5 to 0.30. Normalised correlations run lower, and the observed window on rugnux's own search merges is 0.298 to 0.351; 0.35 is too high, because one crystal's weakest genuine operator reads 0.351. The headroom between a crystal's weakest true operator and its own measured false -operator floor widens on 14 of 14 crystals, median 0.430 to 0.619, and the worst operational margin goes from 0.031 to 0.051. Battery, twice, against a baseline reproducible to zero: the space group is identical on all 38 crystals, per-shell merging is 0 better and 0 worse across all 380 shells, and the merged mmCIF is byte-identical on 38 of 38 - on the 12 crystals whose integration radius now adapts as well as the 25 that do not. The arm is live rather than inert: all 313 operator correlations move while every pair count and every H value stays bit-identical, and on one cubic crystal three operators that raw intensities confirmed are rejected, with the space group unchanged. Following Padilla and Yeates (2003) Acta Cryst. D59, 1124-1130 for why a resolution-normalised statistic is the right one for a symmetry test. Co-Authored-By: Claude Opus 5 (1M context) Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P --- docs/ACKNOWLEDGEMENT.md | 5 +- docs/CHANGELOG.md | 1 + docs/CPU_DATA_ANALYSIS.md | 4 +- .../scale_merge/SearchSpaceGroup.cpp | 66 ++++++++++++++--- image_analysis/scale_merge/SearchSpaceGroup.h | 36 +++++++-- tests/SearchSpaceGroupTest.cpp | 73 +++++++++++++++++++ 6 files changed, 166 insertions(+), 19 deletions(-) diff --git a/docs/ACKNOWLEDGEMENT.md b/docs/ACKNOWLEDGEMENT.md index 4211a4e9..fde8e15b 100644 --- a/docs/ACKNOWLEDGEMENT.md +++ b/docs/ACKNOWLEDGEMENT.md @@ -59,7 +59,10 @@ G. Winter, D. G. Waterman et al., "Robust background modelling in DIALS" (2016), 1912-1921 [doi:10.1107/S1600576716013595](https://doi.org/10.1107/S1600576716013595). **[POINTLESS](https://www.ccp4.ac.uk/)** (CCP4) — the space-group search. Stage A scores each -candidate rotation operator by the correlation of I(h) with I(Rh); the screw-axis test scores a +candidate rotation operator by the correlation of I(h) with I(Rh) on **resolution-normalised** +intensities (E²), as POINTLESS does — both arms of a symmetry pair sit at the same |s|, so on raw +intensities the resolution fall-off is variance shared between them and lifts a false operator's +correlation as much as a true one's; the screw-axis test scores a predicted-absent class against the rest of its own axial row rather than against a global mean or a fixed cut, and lets confidence fall away with the number of axial reflections instead of refusing below a count. P. Evans, "Scaling and assessment of data quality" (2006), Acta Cryst. D62, 72-82 diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 2e148d10..211cbbe5 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -9,6 +9,7 @@ This is an UNSTABLE release. It includes many experimental features, as well as * rugnux: a lattice centring the data could not test - the crystal was integrated on the primitive sub-cell, so the reflections the centring extinguishes were never measured - is marked `UNTESTED` in the space-group candidate table and, where it is adopted, is warned about as coming from the lattice metric rather than from the intensities. * rugnux: the space-group search prints the twin-law disagreement H for every operator it tested and the H ratio of the point group it adopted, on every run, instead of only when that ratio refuses a promotion. * rugnux: `--mode scale` works on a `_process.h5` whose space group came from re-seating the lattice; a file written before this stops with a message naming the two cells instead of failing inside the merge. +* rugnux: de-novo space-group determination scores each symmetry operator on resolution-normalised intensities (E²) rather than on raw merged intensities, so the operator threshold means the same thing on every crystal and no longer moves with the search resolution cut. * rugnux: on a pattern too dense for the widened radius - where neighbouring reflections leave more than 1.1% of the reflections without a background ring - the second integration pass goes back to the fixed 4 px radius, and says so. * rugnux: on rotation data the integration signal radius is set from how wide the crystal's own spots are, measured in the pre-scan, instead of the fixed 4 px; `--adaptive-integration-radius=off` restores the fixed radius, and an explicit `--integration-radius` still overrides both. * rugnux: a directional diffraction limit that is the edge of the measured data rather than the crystal's own limit is marked as such - with a `<` in the report and in `ANISOTROPY_D_MIN_CENSORED`, and in the mmCIF - so `ANISOTROPY_D_MIN_SPREAD` is not read as a measurement when it is a lower bound. diff --git a/docs/CPU_DATA_ANALYSIS.md b/docs/CPU_DATA_ANALYSIS.md index 22d22f4a..ab73d1ae 100644 --- a/docs/CPU_DATA_ANALYSIS.md +++ b/docs/CPU_DATA_ANALYSIS.md @@ -935,12 +935,14 @@ A **dataset-wide** Wilson $B$ is also estimated over the merged reflections — ### 13.1 Space-group determination -When no space group is supplied, a POINTLESS-like search scores Laue-group symmetry (CC of $I(h)$ vs $I(Rh)$ plus merge self-consistency) and detects screw/centering absences from the $P1$-merged intensities. Three tests gate a promotion to higher symmetry, all aimed at the merohedral twin, whose twin law forces non-equivalent reflections together and so mimics symmetry: +When no space group is supplied, a POINTLESS-like search scores Laue-group symmetry (CC of $E^2(h)$ vs $E^2(Rh)$ — the intensities normalised by the mean of their own resolution shell — plus merge self-consistency) and detects screw/centering absences from the $P1$-merged intensities. Three tests gate a promotion to higher symmetry, all aimed at the merohedral twin, whose twin law forces non-equivalent reflections together and so mimics symmetry: 1. **Merge self-consistency** ($\chi^2$ under the candidate group, relative to the confirmed subgroup). On its own this is not sufficient: it is a ratio to an error model that moves with the *amount* of data — the parent's systematic term grows as $\sigma$ shrinks with $1/\sqrt{N}$, while a twin's is already saturated — so its verdict depends on how much data the search saw. 2. **Error-model $b$** (the intensity-proportional systematic). A genuine symmetry step gains multiplicity without inflating $b$; merging a twin law's extra operator inflates it. A $\chi^2$-passing promotion is vetoed when $b$ rises past a bound relative to the confirmed subgroup. 3. **Operator disagreement**, a sigma-free statistic $H=\mathrm{median}\,|I_1-I_2|/(I_1+I_2)$, formed as the ratio of the operators a promotion *adds* to the parent's own, measured on the same reflections. Normalising against the parent divides out the systematic floor that symmetry mates carry on real data, which varies by crystal and by operator; a median is used because a twin perturbs every pair whereas a badly-measured minority perturbs only the tail. Where a candidate has several parents of the same order, it is judged against the worst of them, since a rival subgroup can itself contain the twin laws. +The correlation is on **resolution-normalised** intensity $E^2 = I/\langle I\rangle(\text{shell})$, normalised over exactly the reflections the correlation pairs. Both members of a symmetry pair lie at the same $|s|$, so on raw $I$ the resolution fall-off is variance shared perfectly between the two arms and appears as a positive correlation for *any* pairing at all: a shell-matched random pairing — the exact null for a metrically-allowed false operator — scores a median 0.31 across the rotation battery, and on one crystal 0.53 — above the bound the correlation is tested against. That floor varies more from crystal to crystal (spread 0.46) than the whole true/false gap is wide (0.38), so an absolute bound on the raw statistic is a different test on every crystal; and it moves with the search resolution cut, which is what made that cut a symmetry-deciding parameter. Normalised, the floor has a median of 0.015, never exceeds 0.06, and barely moves with the cut. + The operator correlations are taken on reflections above an $I/\sigma$ cut, and that cut is **capped at the merge's own $I/\sigma$ quantile** rather than applied as a fixed number. On a search merge whose ISa is below 3, a fixed cut of 3 selects nothing at all, leaving every operator correlation undefined and collapsing the point group to 1. The cap keeps at least the strongest quarter and is inert — the cut stays exactly 3.0 — on a healthy merge. Several space groups may share an absence pattern exactly. Where they do, the search scores them identically and **all of them are named** in the result rather than one being reported as the answer: some are enantiomorph pairs, which merged intensities cannot distinguish in principle, and others differ only by a screw condition that the centering condition already implies, so the screw has no observable signature at all. The representative reported first is the lowest space-group number, which is a convention and not a measurement. diff --git a/image_analysis/scale_merge/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index e8c37ad3..fd328b0d 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -276,15 +276,16 @@ SearchSpaceGroupResult SearchSpaceGroup( pass_cc[i] = pass_absence[i] && IoverSigma[i] >= present_cut && (opt.min_i_over_sigma <= 0 || IoverSigma[i] >= opt.min_i_over_sigma); - // Resolution-normalised intensity E^2 = I / (shell), from equal-count resolution shells over - // the reflections the absence test uses. Lets the absence test judge "present" by intensity - // magnitude, not by a possibly under-estimated sigma (see present_e_squared). - std::vector Esq(n, 0.0); - { + // Resolution-normalised intensity E^2 = I / (shell), from equal-count resolution shells over a + // given subset of the merge. Which subset matters: E^2 is only free of the resolution fall-off on + // the population it was normalised over, so a caller has to normalise over the reflections it will + // actually use (see Ecc below). + auto shell_normalised = [&](const std::vector& subset) { + std::vector E(n, 0.0); std::vector order; order.reserve(n); for (size_t i = 0; i < n; ++i) - if (pass_absence[i]) + if (subset[i]) order.push_back(i); std::sort(order.begin(), order.end(), [&](size_t a, size_t b) { return merged[a].d > merged[b].d; }); // low res -> high res @@ -297,9 +298,14 @@ SearchSpaceGroupResult SearchSpaceGroup( sum += I[order[j]]; const double mean = (hi > lo) ? sum / static_cast(hi - lo) : 0.0; for (size_t j = lo; j < hi; ++j) - Esq[order[j]] = mean > 0.0 ? I[order[j]] / mean : 0.0; + E[order[j]] = mean > 0.0 ? I[order[j]] / mean : 0.0; } - } + return E; + }; + + // Over the reflections the absence test uses, so it can judge "present" by intensity magnitude + // rather than by a possibly under-estimated sigma (see present_e_squared). + const std::vector Esq = shell_normalised(pass_absence); // Overlap guard (Stage A / correlation only): drop the extreme resolution-normalised-E tail, which // on a two-lattice crystal is the one-sided overlap contamination that poisons the operator CC. @@ -309,6 +315,33 @@ SearchSpaceGroupResult SearchSpaceGroup( if (pass_cc[i] && Esq[i] > opt.max_e_squared_for_cc) pass_cc[i] = false; + // The operator correlation is scored on E^2, not on raw I. Both members of a symmetry pair sit at + // the same |s|, so the resolution fall-off is variance shared perfectly between the two arms of + // every pair: a Pearson CC on raw I measures the fall-off as well as the symmetry, and the fall-off + // lifts a FALSE operator's CC as much as a true one's. Measured over the rotation battery with + // shell-matched random pairing as the null for a metrically-allowed false operator, that raw-I noise + // floor has a median of 0.31 and reaches 0.53 on one crystal - above the old bound of 0.5 outright - + // and it varies more between crystals (spread 0.46) than the whole true/false gap is wide (0.38), so + // an absolute bound on it was a different test on every crystal. It also moves with the search + // resolution cut, by a median 0.09 and up to 0.23, which is what made that cut decide symmetries. + // Normalised, that floor has a median of 0.015 and a maximum of 0.06, and moves by a median 0.03. + // Following POINTLESS, which likewise scores each symmetry element on normalised intensities + // (Evans, Acta Cryst. D62, 72-82 (2006)). + // + // Normalised over pass_cc - the reflections the correlation actually pairs - and NOT over the + // pass_absence set Esq uses. pass_cc keeps only the stronger reflections and the fraction it keeps + // itself falls with resolution, so an E^2 built on pass_absence still carries a resolution trend + // inside the pass_cc subset; reusing Esq here makes the coupling WORSE than raw I (floor movement + // 0.14). Esq stays as it is: the absence tests need their own set, and the E^2 cap above helps + // DEFINE pass_cc, so normalising that over pass_cc would be circular. + // + // The price is a new dependence in place of the old one: the CC now moves with whatever defines + // pass_cc, i.e. with present_i_over_sigma and the MIN_PRESENT_FRACTION cap above. On the battery + // that is inert - the cap fires on no crystal that has an operator to score - but the cap exists + // for weak merges, so a merge on which it fires decides its point group on a differently + // normalised statistic. + const std::vector Ecc = shell_normalised(pass_cc); + std::unordered_map key_to_index; key_to_index.reserve(n * 2); for (size_t i = 0; i < n; ++i) @@ -321,7 +354,8 @@ SearchSpaceGroupResult SearchSpaceGroup( auto score_operator = [&](const gemmi::Op& op, std::vector& visited, uint32_t& epoch) -> SpaceGroupOperatorScore { ++epoch; - std::vector x, y; + std::vector x, y; // raw merged I of each pair, for the H statistic + std::vector ex, ey; // the same pairs as E^2, for the correlation for (size_t i = 0; i < n; ++i) { if (!pass_cc[i] || visited[i] == epoch) continue; @@ -337,14 +371,21 @@ SearchSpaceGroupResult SearchSpaceGroup( continue; x.push_back(I[i]); y.push_back(I[j]); + ex.push_back(Ecc[i]); + ey.push_back(Ecc[j]); visited[i] = epoch; visited[j] = epoch; } SpaceGroupOperatorScore s; s.op_triplet_hkl = op.as_hkl().triplet('h'); s.n_pairs = static_cast(x.size()); - s.cc = PearsonCC(x, y); + s.cc = PearsonCC(ex, ey); // Sigma-free disagreement over the same pairs (see SpaceGroupOptions::max_operator_h_ratio). + // On RAW I, deliberately - which is why the pairs are collected twice. The shell divisor cancels + // in |I1-I2|/(I1+I2) exactly for a pair whose members share a shell, but not for one that + // straddles a shell boundary: measured over the battery's operators, normalising moves H by a + // median 0.08% but by 3.7% at p95 and 12% at worst, and max_operator_h_ratio has been decided on + // a margin of 1.8%. H is calibrated on raw I and stays there. std::vector hv; hv.reserve(x.size()); for (size_t p = 0; p < x.size(); ++p) { @@ -375,8 +416,9 @@ SearchSpaceGroupResult SearchSpaceGroup( // Conjugate rotations (symmetry-equivalent within the point group) relate symmetry-equivalent // reflection sets, so on real data their CCs cluster; a noisy crystal can push one class member - // below min_operator_cc while the class is unmistakably present (e.g. one cubic 3-fold at 0.48 - // among siblings at 0.53-0.66). Judge each conjugacy class by its mean CC, not its weakest + // below min_operator_cc while the class is unmistakably present (this was first seen on a cubic + // crystal whose three 3-folds spread over 0.13 of CC with only the weakest below the bound). + // Judge each conjugacy class by its mean CC, not its weakest // member, so a genuine high-symmetry point group is not lost to one marginal operator. chi2_under // (below) remains the safety net against a truly false promotion. Returns {all classes present, // worst class-mean CC}. diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index 278bc897..91077afb 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -15,7 +15,8 @@ // of POINTLESS (Evans 2006): // // Stage A - point group (Laue) symmetry. Every candidate rotation operator is scored once by -// the correlation of I(h) with I(Rh). The chosen point group is the largest one all +// the correlation of E^2(h) with E^2(Rh), on resolution-normalised intensities as +// POINTLESS does. The chosen point group is the largest one all // of whose operators are confirmed (high CC). A wrong operator scores ~0, so this is // self-pruning - the unit cell is not needed. // @@ -28,7 +29,7 @@ struct SpaceGroupOperatorScore { std::string op_triplet_hkl; // reciprocal-space triplet of the rotation, e.g. "-h,-k,l" - double cc = 0.0; // correlation of I(h) with I(Rh) + double cc = 0.0; // correlation of E^2(h) with E^2(Rh) - resolution-normalised int n_pairs = 0; // independent reflection pairs the CC was computed from bool present = false; // operator confirmed as a real symmetry of the intensities // MEDIAN |I1-I2|/(I1+I2) over this operator's pairs - the disagreement the operator implies, with no @@ -91,15 +92,40 @@ struct SearchSpaceGroupOptions { // stage only. A second lattice deposits intensity on one reciprocal position but not its symmetry // mate, so an overlap-contaminated reflection is a one-sided E^2 outlier that poisons an operator's // I(h)/I(Rh) correlation (it flipped a pseudo-merohedral P2_1 case to P1 under -A: excluding the E^2>9 tail, ~0.4% - // of reflections, lifted the 2-fold CC 0.33->0.53 back over the gate). Clean Wilson-distributed data + // of reflections, lifted the 2-fold CC 0.33->0.53 back over the gate - values on the raw-I scale the + // operator CC used before it was normalised). Clean Wilson-distributed data // almost never reaches E^2=9 (P(E^2>9) ~ 0.01-0.3%), so this removes essentially nothing there and // only trims the overlap tail. 0 disables. Absences are unaffected (they need the weak tail). + // This E^2 is the one normalised over pass_absence, NOT the one the correlation is scored on: the + // cut helps DEFINE pass_cc, so normalising it over pass_cc would be circular. double max_e_squared_for_cc = 9.0; // --- Stage A: point group --- - // A rotation is accepted as a real symmetry when its I(h)/I(Rh) correlation reaches this over + // A rotation is accepted as a real symmetry when its E^2(h)/E^2(Rh) correlation reaches this over // at least min_pairs_per_operator independent pairs. - double min_operator_cc = 0.5; + // + // The correlation is on RESOLUTION-NORMALISED intensity (see SearchSpaceGroup.cpp), which is what + // sets this value. On raw I the two arms of every pair share the whole resolution fall-off, so a + // completely false operator still scores a large positive CC - measured over the rotation battery + // by shell-matched random pairing, that noise floor runs 0.05 to 0.53 across crystals, a SPREAD + // (0.46) wider than the entire true/false gap (0.38). One absolute number therefore meant a + // different test on every crystal: the old bound of 0.5 on raw I corresponds, crystal by crystal, + // to a normalised threshold anywhere from -0.03 to 0.47. On the normalised statistic the floor has + // a median of 0.015 and never exceeds 0.06, so a single number finally means the same thing on + // every crystal. + // + // 0.30 is the midpoint of the gap the battery leaves: the weakest genuine conjugacy-class mean is + // 0.35 and the strongest false candidate 0.23 (both on the pessimistic arm, which carries rugnux's + // own search merge onto the normalised scale). Scanning the constant against the battery: 0.20 lets + // a cubic over-promotion in, 0.35 loses a genuine monoclinic 2-fold, and 0.25-0.30 changes nothing. + // It is also the median of the per-crystal threshold the old 0.5 already imposed (0.309), so it is + // the behaviour-preserving choice at the median crystal and only redistributes strictness at the + // tails - away from the strong crystals, which were being handed margin out of their own noise + // floor, and towards the weak ones, which had none. + double min_operator_cc = 0.30; + // NOTE this is under-calibrated for any threshold in this range: a Pearson CC over 20 null pairs + // has sd ~0.24, so a 1.3-sigma fluctuation clears the bound. Raising it (or making the bound + // pair-count-aware) is a separate change with its own calibration. int min_pairs_per_operator = 20; // Per-operator CC alone cannot tell a real weak operator from a false strong one (a noisy crystal's diff --git a/tests/SearchSpaceGroupTest.cpp b/tests/SearchSpaceGroupTest.cpp index 3d3e7ba3..0426e964 100644 --- a/tests/SearchSpaceGroupTest.cpp +++ b/tests/SearchSpaceGroupTest.cpp @@ -294,3 +294,76 @@ TEST_CASE("SearchSpaceGroup weighs a screw's absences by evidence, not by how ma CHECK(result.best_space_group->short_name() == "P2"); } } + +// The operator correlation is on resolution-normalised E^2, not on raw I (see SearchSpaceGroup.cpp). +// Both members of a symmetry pair sit at the same |s|, so on raw intensities the resolution fall-off +// is variance shared perfectly between the two arms of every pair and reads as a correlation for ANY +// pairing at all. These two cases pin that down from both sides. +TEST_CASE("SearchSpaceGroup operator correlation reads symmetry, not the resolution fall-off", + "[SearchSpaceGroup]") { + // Intensities that are a smooth function of resolution times an INDEPENDENT per-reflection + // factor: a Wilson-like fall-off with no symmetry in it whatsoever. + auto radial_only = [](int hmax) { + std::vector merged; + for (int h = -hmax; h <= hmax; ++h) + for (int k = -hmax; k <= hmax; ++k) + for (int l = -hmax; l <= hmax; ++l) { + if ((h == 0 && k == 0 && l == 0) || std::make_tuple(-h, -k, -l) < std::make_tuple(h, k, l)) + continue; + const double d = CalcSyntheticD(h, k, l); + const double falloff = std::exp(-30.0 / (d * d)); + // Deterministic, independent of any symmetry mate: reuse the hash on the raw index. + const double jitter = SyntheticIntensityFromAsu(gemmi::Op::Miller{{h, k, l}}) / 350.0; + const double I = 1.0e5 * falloff * jitter; + merged.push_back(MergedReflection{ + .h = h, .k = k, .l = l, .I = I, .sigma = I / 20.0, .d = d}); + } + return merged; + }; + + SearchSpaceGroupOptions opt; + opt.merge_friedel = true; + + SECTION("a fall-off with no symmetry in it confirms no operator") { + const auto result = SearchSpaceGroup(radial_only(8), opt); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.operator_scores.size() > 1); + for (const auto& s : result.operator_scores) { + INFO("operator " << s.op_triplet_hkl); + CHECK(s.n_pairs >= opt.min_pairs_per_operator); + CHECK(s.cc < opt.min_operator_cc); + CHECK_FALSE(s.present); + } + CHECK(result.point_group_hm == "1"); + } + + SECTION("a real operator under the same fall-off is confirmed, and does not move with the cut") { + // Same fall-off, but the intensities now carry a genuine monoclinic 2-fold. + const gemmi::SpaceGroup& sg = gemmi::get_spacegroup_by_name("P 1 2 1"); + const gemmi::ReciprocalAsu rasu(&sg); + const gemmi::GroupOps gops = sg.operations(); + auto merged = radial_only(8); + for (auto& r : merged) { + const auto [asu, plus] = rasu.to_asu_sign(gemmi::Op::Miller{{r.h, r.k, r.l}}, gops); + const double falloff = std::exp(-30.0 / (r.d * r.d)); + r.I = 1.0e5 * falloff * SyntheticIntensityFromAsu(asu) / 350.0; + r.sigma = r.I / 20.0; + } + auto two_fold_cc = [&](double d_min) { + SearchSpaceGroupOptions o = opt; + o.d_min_limit_A = d_min; + const auto result = SearchSpaceGroup(merged, o); + INFO(SearchSpaceGroupResultToText(result)); + REQUIRE(result.point_group_hm == "2"); + double cc = -2.0; + for (const auto& s : result.operator_scores) + if (s.present) + cc = s.cc; + REQUIRE(cc > opt.min_operator_cc); + return cc; + }; + // The whole point of normalising: how much of the fall-off is inside the merge no longer + // moves the operator's score, so the search resolution cut cannot decide the symmetry. + CHECK(std::fabs(two_fold_cc(0.0) - two_fold_cc(6.0)) < 0.05); + } +}