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

rugnux: recover genuine high symmetry in de-novo space-group search

Browse Source 
The POINTLESS-style symmetry search discarded a real point group on two
hard rotation-test crystals:

- Ins_I_3 (I23): one of the eight cubic 3-fold operators scored CC 0.483,
  just under the 0.5 per-operator floor, so I23 fell to I222 and merging
  collapsed (3-fold orientation ambiguity; CC1/2 1%, R-meas 558%). Gate A
  now scores each CONJUGACY CLASS of a point group's rotations by its mean
  CC - conjugate operators are symmetry-equivalent and must stand or fall
  together - so the cubic 3-fold class (mean 0.57) is confirmed. -> I23,
  CC1/2 99.9%. Degenerate (few-pair / NaN) class members are skipped
  rather than vetoing the class.

- Benas_3/7 (F432): the chi^2 self-consistency gate rejected the true
  cubic group because the weak-data error model is badly miscalibrated
  (reduced chi^2 ~150), which inflates the chi^2 ratio with point-group
  order for genuine high symmetry too. The ratio bound is now widened by
  log10(chi2_ref) so a broken error model does not veto real symmetry,
  while a false operator whose ratio is far worse than its subgroup's is
  still caught. The bound stays a per-candidate chi^2 test - chi^2 still
  arbitrates and is never bypassed. Verified: F432 (chi2_ref 151, ratio
  3.2 < bound 4.2) accepted; a calibrated pseudo/twin R32 (Ins_H_3:
  chi2_ref 2.6, ratio 4.2 > bound 2.4) still rejected, stays R3.

SearchSpaceGroup.{cpp,h} only. Validated on the 24-crystal rotation-test
battery: exactly these two space groups change (both now match the
reference); zero regressions - every pseudo-symmetry P1 case and every
correct assignment unchanged.

Co-Authored-By: Claude Opus 4.8 (1M context) <noreply@anthropic.com>
This commit is contained in:
leonarski_f
2026-07-09 21:14:18 +02:00
co-authored by Claude Opus 4.8
parent 4f3153d6c7
commit 8c288c8697
2 changed files with 69 additions and 15 deletions
Download Patch File Download Diff File Expand all files Collapse all files
image_analysis/scale_merge/SearchSpaceGroup.cpp
+66 -15
View File
@@ -271,6 +271,54 @@ SearchSpaceGroupResult SearchSpaceGroup(
return op_cache.emplace(rk, score_operator(op)).first->second;
};
// 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
// 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}.
auto point_group_present = [&](const std::vector<gemmi::Op>& rots) -> std::pair<bool, double> {
const size_t m = rots.size();
std::vector<int> cls(m, -1);
int n_cls = 0;
for (size_t i = 0; i < m; ++i) {
if (cls[i] >= 0)
continue;
cls[i] = n_cls;
for (size_t j = i + 1; j < m; ++j)
if (cls[j] < 0)
for (const auto& p : rots)
if ((p * rots[i] * p.inverse()).rot == rots[j].rot) {
cls[j] = n_cls;
break;
}
++n_cls;
}
bool ok = true;
double worst_mean = 1.0;
for (int c = 0; c < n_cls; ++c) {
// Average over the class members that actually have enough pairs to score; a single
// low-multiplicity / degenerate (NaN) operator in an otherwise strong class is skipped,
// not allowed to veto the class. The class must still have at least one scored member.
double sum_cc = 0.0;
int n_valid = 0;
for (size_t i = 0; i < m; ++i)
if (cls[i] == c) {
const auto& s = operator_score(rots[i]);
if (s.n_pairs >= opt.min_pairs_per_operator && std::isfinite(s.cc)) {
sum_cc += s.cc;
++n_valid;
}
}
const double mean_cc = n_valid > 0 ? sum_cc / n_valid : 0.0;
worst_mean = std::min(worst_mean, mean_cc);
if (n_valid == 0 || mean_cc < opt.min_operator_cc)
ok = false;
}
return {ok, worst_mean};
};
std::optional<RotationSet> holohedry;
if (opt.lattice_system.has_value())
holohedry = HolohedryRotationSet(opt.lattice_system.value());
@@ -319,41 +367,44 @@ SearchSpaceGroupResult SearchSpaceGroup(
};
// 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; };
struct PGCand { const PointGroupInfo* pg; int order; double min_class_cc; double chi2; };
std::vector<PGCand> pg_cands;
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();
for (const auto& op : pg.rotations) {
const auto& s = operator_score(op);
all_present = all_present && s.present;
min_cc = std::min(min_cc, s.cc);
}
if (!all_present)
const auto [present, min_class_cc] = point_group_present(pg.rotations);
if (!present)
continue;
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});
pg_cands.push_back({&pg, static_cast<int>(pg.rotations.size()) + 1, min_class_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
// 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.
// chi^2 is not inflated past the miscalibration-widened bound below; ties -> higher min class 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) {
// Widen the chi^2 ratio bound by the error-model miscalibration. When even the best candidate's
// reduced chi^2 (chi2_ref) is far above 1 (weak data / uncorrected decay), the ratio grows with
// point-group order for genuine high symmetry too, so the fixed bound wrongly rejects it. Adding
// log10(chi2_ref) lets the bound exceed max_merge_chi2_ratio by ~the order of magnitude of the
// miscalibration - enough to keep the true group (F432 weak data: chi2_ref~150, ratio 3.2 vs
// bound ~4.2) while STILL rejecting a false operator whose ratio is far worse than its subgroup's
// (a twin/pseudo R32 on calibrated data: chi2_ref~2.6, ratio 4.2 vs bound ~2.4). The bound stays
// a per-candidate chi^2 test, so chi^2 still arbitrates between candidates - it is never bypassed.
const double allowed_ratio = opt.max_merge_chi2_ratio + std::log10(std::max(1.0, chi2_ref));
const bool consistent = c.pg->rotations.empty() || !std::isfinite(c.chi2) ||
!std::isfinite(chi2_ref) || c.chi2 <= chi2_ref * opt.max_merge_chi2_ratio;
!std::isfinite(chi2_ref) || c.chi2 <= chi2_ref * allowed_ratio;
if (!consistent)
continue;
if (c.order > best_pg_order || (c.order == best_pg_order && c.min_cc > best_pg_min_cc)) {
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;
best_pg_min_cc = c.min_cc;
best_pg_min_cc = c.min_class_cc;
}
}
image_analysis/scale_merge/SearchSpaceGroup.h
+3
View File
@@ -79,6 +79,9 @@ struct SearchSpaceGroupOptions {
// 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.
// This bound is the calibrated-error baseline; SearchSpaceGroup widens it by log10(chi2_ref) so a
// badly miscalibrated error model (weak data / uncorrected decay, where the ratio grows with
// point-group order for real symmetry too) does not spuriously reject a genuine high-symmetry group.
double max_merge_chi2_ratio = 2.0;
// --- Stage B: space group (screw axes / centering) ---