mx/Jungfraujoch
2
1
Fork 0
Code Issues Pull Requests Actions Packages Projects Releases 72 Wiki Activity

Space-group search: score operators on normalised intensities

Browse Source 
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) <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01CHMmeM1d489zvNFT7ZMN2P
This commit is contained in:
leonarski_f
2026-08-26 00:27:13 +02:00
co-authored by Claude Opus 5
parent 4e6eb18d93
commit f6bfa54624
6 changed files with 166 additions and 19 deletions
Download Patch File Download Diff File Expand all files Collapse all files
tests/SearchSpaceGroupTest.cpp
+73
View File
@@ -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<MergedReflection> 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);
}
}