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