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

space group: use merge reduced chi^2 (not R-meas) for the consistency check

Browse Source 
The R-meas-ratio over-symmetry guard (commit 7ade6d9) missed a near-perfect
pseudo-2-fold: InsH2 (XDS R3) has a false 2-fold at operator CC 0.85 whose merged
intensities still correlate well, so R-meas inflated only ~1.3x and R32 was accepted
(ISa collapsing to 3.9). The reduced chi^2 (within-orbit scatter / sigma^2) is the
right signal - false equivalents disagree by many sigma even when they correlate well.

Measured chi^2 ratio (candidate / best subgroup) across the test set: true point groups
<= 1.63 (a cubic merge of a noisy crystal, InsI3, is the worst real case), false >= 2.26
(InsH2) and 6.0 (InsH3); threshold 2.0 separates them. InsH2 now resolves to R3 (ISa 6.1
from the false R32's 3.9); InsH3/InsI2/InsI3/cytC/lysoC/MyoB/EP0210 unchanged and correct.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-06-30 21:54:51 +02:00
co-authored by Claude Opus 4.8
parent fea2790a41
commit 636ed1e180
2 changed files with 45 additions and 37 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/SearchSpaceGroup.cpp
+36 -31
View File
@@ -173,7 +173,7 @@ SearchSpaceGroupResult SearchSpaceGroup(
// Flatten the reflections and mark which ones each stage may use. The correlation stage drops
// weak reflections; the absence stage must keep them - that is where the screw-axis signal is.
std::vector<int> H(n), K(n), L(n);
std::vector<double> I(n), IoverSigma(n);
std::vector<double> I(n), Sigma(n), IoverSigma(n);
std::vector<HKLKey> key(n);
std::vector<char> pass_absence(n, 0), pass_cc(n, 0);
@@ -181,6 +181,7 @@ SearchSpaceGroupResult SearchSpaceGroup(
const auto& r = merged[i];
H[i] = r.h; K[i] = r.k; L[i] = r.l;
I[i] = r.I;
Sigma[i] = std::isfinite(r.sigma) && r.sigma > 0 ? r.sigma : 0.0;
key[i] = Canonicalize(r.h, r.k, r.l, opt.merge_friedel);
const bool finite = std::isfinite(r.I) && std::isfinite(r.sigma) && r.sigma > 0 &&
@@ -250,14 +251,18 @@ SearchSpaceGroupResult SearchSpaceGroup(
holohedry = HolohedryRotationSet(opt.lattice_system.value());
const auto point_groups = EnumeratePointGroups(holohedry);
// R-meas of the intensities merged under a point group's rotations - how well its symmetry
// equivalents agree. A real point group keeps this low; a false operator forces non-equivalent
// reflections together and inflates it. Computed over the present (pass_cc) reflections.
auto rmeas_under = [&](const std::vector<gemmi::Op>& rotations) -> double {
std::unordered_map<HKLKey, std::pair<int, double>, HKLKeyHash> grp; // orbit rep -> (n, sum I)
// Reduced chi^2 of the intensities merged under a point group's rotations - how well its symmetry
// equivalents agree RELATIVE TO THEIR ERRORS. A real point group gives ~1; a false operator forces
// non-equivalent reflections together, so they disagree by many sigma and chi^2 blows up. This is
// more sensitive than R-meas to a strong pseudo-symmetry (where the intensities still correlate well
// - high operator CC - but not within their errors). Inverse-variance weighted mean per orbit, over
// the present (pass_cc) reflections.
auto chi2_under = [&](const std::vector<gemmi::Op>& rotations) -> double {
struct Acc { double sw = 0.0, swI = 0.0; int n = 0; };
std::unordered_map<HKLKey, Acc, HKLKeyHash> grp;
std::vector<HKLKey> rep(n);
for (size_t i = 0; i < n; ++i) {
if (!pass_cc[i])
if (!pass_cc[i] || !(Sigma[i] > 0.0))
continue;
HKLKey best = key[i];
for (const auto& op : rotations) {
@@ -268,30 +273,30 @@ SearchSpaceGroupResult SearchSpaceGroup(
}
rep[i] = best;
auto& g = grp[best];
g.first += 1; g.second += I[i];
const double w = 1.0 / (Sigma[i] * Sigma[i]);
g.sw += w; g.swI += w * I[i]; g.n += 1;
}
std::unordered_map<HKLKey, double, HKLKeyHash> absdev;
double chi2 = 0.0;
long dof = 0;
for (size_t i = 0; i < n; ++i) {
if (!pass_cc[i])
if (!pass_cc[i] || !(Sigma[i] > 0.0))
continue;
const auto& g = grp[rep[i]];
if (g.first >= 2)
absdev[rep[i]] += std::fabs(I[i] - g.second / g.first);
}
double num = 0, den = 0;
for (const auto& [k, g] : grp) {
if (g.first < 2)
if (g.n < 2)
continue;
num += std::sqrt(static_cast<double>(g.first) / (g.first - 1)) * absdev[k];
den += g.second;
const double mean = g.swI / g.sw, dev = I[i] - mean;
chi2 += dev * dev / (Sigma[i] * Sigma[i]);
}
return den > 0 ? num / den : std::numeric_limits<double>::quiet_NaN();
for (const auto& [k, g] : grp)
if (g.n >= 2)
dof += g.n - 1;
return dof > 0 ? chi2 / static_cast<double>(dof) : std::numeric_limits<double>::quiet_NaN();
};
// Operator-CC-confirmed candidates, each with its merge R-meas; rmeas_ref = the most consistent.
struct PGCand { const PointGroupInfo* pg; int order; double min_cc; double rmeas; };
// Operator-CC-confirmed candidates, each with its merge chi^2; chi2_ref = the most consistent.
struct PGCand { const PointGroupInfo* pg; int order; double min_cc; double chi2; };
std::vector<PGCand> pg_cands;
double rmeas_ref = std::numeric_limits<double>::infinity();
double chi2_ref = std::numeric_limits<double>::infinity();
for (const auto& pg : point_groups) {
bool all_present = true;
double min_cc = pg.rotations.empty() ? 1.0 : std::numeric_limits<double>::infinity();
@@ -302,22 +307,22 @@ SearchSpaceGroupResult SearchSpaceGroup(
}
if (!all_present)
continue;
const double rm = pg.rotations.empty() ? std::numeric_limits<double>::quiet_NaN()
: rmeas_under(pg.rotations);
pg_cands.push_back({&pg, static_cast<int>(pg.rotations.size()) + 1, min_cc, rm});
if (!pg.rotations.empty() && std::isfinite(rm))
rmeas_ref = std::min(rmeas_ref, rm);
const double ch = pg.rotations.empty() ? std::numeric_limits<double>::quiet_NaN()
: chi2_under(pg.rotations);
pg_cands.push_back({&pg, static_cast<int>(pg.rotations.size()) + 1, min_cc, ch});
if (!pg.rotations.empty() && std::isfinite(ch))
chi2_ref = std::min(chi2_ref, ch);
}
// Choose the largest point group that is both operator-confirmed AND self-consistent (its merge
// R-meas is not inflated past max_rmeas_ratio x the most-consistent candidate; ties -> higher min
// CC). Identity (no operators) is always consistent, so it stays the P1 fallback.
// chi^2 is not inflated past max_merge_chi2_ratio x the most-consistent candidate; ties -> higher
// min CC). Identity (no operators) is always consistent, so it stays the P1 fallback.
const PointGroupInfo* best_pg = nullptr;
int best_pg_order = 0;
double best_pg_min_cc = -2.0;
for (const auto& c : pg_cands) {
const bool consistent = c.pg->rotations.empty() || !std::isfinite(c.rmeas) ||
!std::isfinite(rmeas_ref) || c.rmeas <= rmeas_ref * opt.max_rmeas_ratio;
const bool consistent = c.pg->rotations.empty() || !std::isfinite(c.chi2) ||
!std::isfinite(chi2_ref) || c.chi2 <= chi2_ref * opt.max_merge_chi2_ratio;
if (!consistent)
continue;
if (c.order > best_pg_order || (c.order == best_pg_order && c.min_cc > best_pg_min_cc)) {
image_analysis/scale_merge/SearchSpaceGroup.h
+9 -6
View File
@@ -71,12 +71,15 @@ struct SearchSpaceGroupOptions {
double min_operator_cc = 0.5;
int min_pairs_per_operator = 20;
// Per-operator CC alone cannot tell a real weak operator from a false moderate one (a noisy
// crystal's genuine 2-fold can score below a pseudo-symmetric crystal's false one). So a point
// group is also required to be SELF-CONSISTENT: merging the intensities under it must not inflate
// R-meas beyond this factor times the most-consistent (lowest-R-meas) candidate. A false operator
// forces non-equivalent reflections together and blows R-meas up; a real one leaves it ~flat.
double max_rmeas_ratio = 1.5;
// Per-operator CC alone cannot tell a real weak operator from a false strong one (a noisy crystal's
// genuine 2-fold can score below a pseudo-symmetric crystal's near-perfect false one). So a point
// group is also required to be SELF-CONSISTENT: merging the intensities under it must not inflate the
// reduced chi^2 (within-orbit scatter / sigma^2) beyond this factor times the most-consistent
// candidate. A false operator forces non-equivalent reflections together so they disagree by many
// sigma and chi^2 blows up; a real one leaves it ~flat even when the operator CC is only moderate.
// 2.0 separates the test set: a true point group's chi^2 stays within ~1.6x the best subgroup
// (the cubic merge of a noisy crystal is the worst real case), while a false operator gives >=2.3x.
double max_merge_chi2_ratio = 2.0;
// --- Stage B: space group (screw axes / centering) ---
bool determine_space_group = true; // false: stop at the symmorphic representative