From 23adbcb6b30a4e7a0e77cade8a80d1dbfb9aaa5d Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 12 Nov 2026 07:02:36 +0100 Subject: [PATCH] rugnux: twin-immune zone takes a pseudo-translation's modulation out per shell, not by splitting the zone With a pseudo-translation declared, ReadZone split the zone into the two cos(2 pi h.u) classes and normalised each in bins of 100 reflections. The zone a refused promotion is decided on is read on the adopted group's merge and can hold only a couple of hundred reflections, so each class got a single bin spanning the whole resolution range: the fall-off read as spread and the zone read centric. A twinned 321 crystal (8c3e) whose pseudo-translation is now detected read its 622-over-32 zone at +0.9 nats as one population and +21.4 split, past the 20-nat bound, and its twin law 622 was adopted (single sweep and multisweep). The modulation is now divided out of every strong reflection per resolution shell, by its class's mean over the shell's mean, both measured on the acentric reflections, and the zone and control are then read as one population, as without a pseudo-translation. Synthetic test: a 321 twin (30%, 622 law) with a 0.6 modulation at (2/3,1/3,1/2), 248 zone reflections: -11.2 nats before, -25.8 after. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK --- .../scale_merge/TwinningAnalysis.cpp | 130 +++++++++++------- image_analysis/scale_merge/TwinningAnalysis.h | 2 +- tests/TwinningAnalysisTest.cpp | 26 ++++ 3 files changed, 104 insertions(+), 54 deletions(-) diff --git a/image_analysis/scale_merge/TwinningAnalysis.cpp b/image_analysis/scale_merge/TwinningAnalysis.cpp index 9939e4aa1..fbf1239f1 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.cpp +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -811,7 +811,7 @@ namespace { // of 2.5 and ~0 at 16, measured), so only shells above this are read. At 5 the acentric // control of four reference merges reads 0.65-0.74. constexpr double ZONE_MIN_SHELL_ISIG = 5.0; - // Each class is normalised against its own mean in resolution bins of this many reflections, so a + // Each zone is normalised against its own mean in resolution bins of this many reflections, so a // zone that is a plane through reciprocal space is not judged on another direction's fall-off. constexpr size_t ZONE_BIN = 100; // A reflection this far above its bin mean is not a Wilson draw (centric P ~ 1e-5): one such @@ -898,52 +898,42 @@ namespace { } // `twin_fraction` is the compression the centric hypothesis is read at (CentricOverAcentric). - TwinImmuneZone ReadZone(std::vector refl, const std::vector>* tncs, - size_t nthreads, double twin_fraction) { + TwinImmuneZone ReadZone(std::vector refl, size_t nthreads, double twin_fraction) { TwinImmuneZone z; // The evidence of each kept reflection is a pure function of its (E^2, error) and costs a few // hundred exponentials, so the pairs are collected here and evaluated together below; the sum // is then taken over them in this same order. std::vector> kept_e2_s; - // A pseudo-translation u multiplies the intensity by 1 + f cos(2 pi h.u): normalising within - // each sign of the cosine keeps the two populations from reading as one broad distribution. - std::vector classes[2]; - for (const auto* r : refl) { - const double c = tncs == nullptr ? 1.0 : TncsClass(r->h, r->k, r->l, *tncs); - classes[c >= 0.0 ? 0 : 1].push_back(r); - } double sum = 0.0, sum2 = 0.0; - for (auto& cls : classes) { - std::sort(cls.begin(), cls.end(), [](const auto* a, const auto* b) { return a->d > b->d; }); - const size_t n_bins = std::max(1, cls.size() / ZONE_BIN); - for (size_t b = 0; b < n_bins; ++b) { - const size_t first = b * cls.size() / n_bins, last = (b + 1) * cls.size() / n_bins; - double mean = 0.0; - for (size_t i = first; i < last; ++i) - mean += cls[i]->I; - mean /= static_cast(last - first); - if (mean <= 0.0) - continue; - double kept = 0.0; - int n_kept = 0; - for (size_t i = first; i < last; ++i) - if (cls[i]->I / mean <= ZONE_OUTLIER_E2) { - kept += cls[i]->I; - ++n_kept; - } - if (n_kept == 0 || kept <= 0.0) - continue; - mean = kept / n_kept; - for (size_t i = first; i < last; ++i) { - const double e2 = cls[i]->I / mean; - if (e2 > ZONE_OUTLIER_E2) - continue; - const double v = std::fabs(e2 - 1.0); - sum += v; - sum2 += v * v; - ++z.n; - kept_e2_s.emplace_back(e2, cls[i]->sigma / mean); + std::sort(refl.begin(), refl.end(), [](const auto* a, const auto* b) { return a->d > b->d; }); + const size_t n_bins = std::max(1, refl.size() / ZONE_BIN); + for (size_t b = 0; b < n_bins; ++b) { + const size_t first = b * refl.size() / n_bins, last = (b + 1) * refl.size() / n_bins; + double mean = 0.0; + for (size_t i = first; i < last; ++i) + mean += refl[i]->I; + mean /= static_cast(last - first); + if (mean <= 0.0) + continue; + double kept = 0.0; + int n_kept = 0; + for (size_t i = first; i < last; ++i) + if (refl[i]->I / mean <= ZONE_OUTLIER_E2) { + kept += refl[i]->I; + ++n_kept; } + if (n_kept == 0 || kept <= 0.0) + continue; + mean = kept / n_kept; + for (size_t i = first; i < last; ++i) { + const double e2 = refl[i]->I / mean; + if (e2 > ZONE_OUTLIER_E2) + continue; + const double v = std::fabs(e2 - 1.0); + sum += v; + sum2 += v * v; + ++z.n; + kept_e2_s.emplace_back(e2, refl[i]->sigma / mean); } } std::vector evidence(kept_e2_s.size()); @@ -1045,11 +1035,46 @@ namespace { } } + // A pseudo-translation u multiplies the intensity by 1 + f cos(2 pi h.u), so each shell holds two + // populations, one on each sign of the cosine, of different mean. Each is divided by its own mean + // over the shell's mean, both measured on the shell's acentric reflections, before any zone is read. + // Splitting the zone itself into the two classes instead left a zone of a few hundred reflections + // one resolution bin per class, normalised across the whole fall-off, and it read centric for that: + // a twinned 321 crystal's zone of 227 reflections read +0.9 nats as one population and +21.4 + // split, and 622 - its twin law - was adopted. + void RemoveTncsModulation(std::vector& strong, const std::vector& shell, + const gemmi::GroupOps& gops, const std::vector>& tncs) { + const int n_shells = *std::max_element(shell.begin(), shell.end()) + 1; + std::vector> sum(n_shells, {0.0, 0.0}); + std::vector> n(n_shells, {0, 0}); + std::vector cls(strong.size()); + for (size_t i = 0; i < strong.size(); ++i) { + const auto& r = strong[i]; + cls[i] = TncsClass(r.h, r.k, r.l, tncs) >= 0.0 ? 0 : 1; + if (!gops.is_reflection_centric(gemmi::Op::Miller{{r.h, r.k, r.l}})) { + sum[shell[i]][cls[i]] += r.I; + n[shell[i]][cls[i]] += 1; + } + } + for (size_t i = 0; i < strong.size(); ++i) { + const int b = shell[i]; + if (n[b][0] == 0 || n[b][1] == 0) + continue; + const double shell_mean = (sum[b][0] + sum[b][1]) / (n[b][0] + n[b][1]); + const double f = sum[b][cls[i]] / n[b][cls[i]] / shell_mean; + if (f > 0.0) { + strong[i].I = static_cast(strong[i].I / f); + strong[i].sigma = static_cast(strong[i].sigma / f); + } + } + } + // The reflections a zone test can be read on: present in the group, general position, in the // shells with signal enough (ZONE_MIN_SHELL_ISIG), copied with the anisotropy taken out of I and // sigma (RemoveAnisotropy). Fills the resolution range read. std::vector StrongReflections(const std::vector& p1_merged, const gemmi::UnitCell& cell, const gemmi::GroupOps& gops, + const std::vector>* tncs_vector, TwinImmuneZoneResult& result) { std::vector usable; for (const auto& r : p1_merged) { @@ -1094,6 +1119,8 @@ namespace { } if (!strong.empty()) RemoveAnisotropy(strong, shell, cell, gops, result); + if (!strong.empty() && tncs_vector != nullptr) + RemoveTncsModulation(strong, shell, gops, *tncs_vector); return strong; } @@ -1101,13 +1128,12 @@ namespace { // read as the zones' is (at the same twin fraction), which is what the calibration measures the // normalisation's excess on. TwinImmuneZone ReadControl(const std::vector& strong, const gemmi::GroupOps& gops, - const std::vector>* tncs_vector, double twin_fraction, - size_t nthreads) { + double twin_fraction, size_t nthreads) { std::vector control; for (const auto& r : strong) if (!gops.is_reflection_centric(gemmi::Op::Miller{{r.h, r.k, r.l}})) control.push_back(&r); - TwinImmuneZone z = ReadZone(control, tncs_vector, nthreads, twin_fraction); + TwinImmuneZone z = ReadZone(control, nthreads, twin_fraction); z.operators = "acentric"; return z; } @@ -1116,9 +1142,7 @@ namespace { // group but not in H - the self-mates of the coset the group adds over H. TwinImmuneZone ReadSubgroupZone(const std::vector& strong, const gemmi::GroupOps& gops, const std::vector& rotations, - const std::vector& h_ops, - const std::vector>* tncs_vector, double twin_fraction, - size_t nthreads) { + const std::vector& h_ops, double twin_fraction, size_t nthreads) { auto contains = [&](const std::vector& set, const gemmi::Op& op) { return std::any_of(set.begin(), set.end(), [&](const gemmi::Op& o) { return o.rot == op.rot; }); }; @@ -1132,7 +1156,7 @@ namespace { if (gops.is_reflection_centric(h) && !centric_in_h) zone.push_back(&r); } - TwinImmuneZone z = ReadZone(zone, tncs_vector, nthreads, twin_fraction); + TwinImmuneZone z = ReadZone(zone, nthreads, twin_fraction); for (const auto& op : rotations) if (!contains(h_ops, op)) z.operators += (z.operators.empty() ? "" : " ") + op.as_hkl().triplet(); @@ -1151,11 +1175,11 @@ TwinImmuneZoneResult AnalyzeTwinImmuneZone(const std::vector& const gemmi::GroupOps gops = group.operations(); const std::vector rotations = ProperRotations(gops); const std::vector h_ops = ProperRotations(subgroup.operations()); - const auto strong = StrongReflections(p1_merged, cell, gops, result); + const auto strong = StrongReflections(p1_merged, cell, gops, tncs_vector, result); if (strong.empty()) return result; - result.control = ReadControl(strong, gops, tncs_vector, twin_fraction, nthreads); - result.zones.push_back(ReadSubgroupZone(strong, gops, rotations, h_ops, tncs_vector, twin_fraction, + result.control = ReadControl(strong, gops, twin_fraction, nthreads); + result.zones.push_back(ReadSubgroupZone(strong, gops, rotations, h_ops, twin_fraction, nthreads)); Calibrate(result); return result; @@ -1214,15 +1238,15 @@ TwinImmuneZoneResult AnalyzeTwinImmuneZones(const std::vector& if (subgroups.empty()) return result; - const auto strong = StrongReflections(p1_merged, cell, gops, result); + const auto strong = StrongReflections(p1_merged, cell, gops, tncs_vector, result); if (strong.empty()) return result; - result.control = ReadControl(strong, gops, tncs_vector, twin_fraction, nthreads); + result.control = ReadControl(strong, gops, twin_fraction, nthreads); // Per hypothesis H, the zone: centric in the group, acentric in H - the reflections some added // operator sends to -h, i.e. the self-mates of that coset read as a twin law. for (const auto& h_ops : subgroups) - result.zones.push_back(ReadSubgroupZone(strong, gops, rotations, h_ops, tncs_vector, twin_fraction, + result.zones.push_back(ReadSubgroupZone(strong, gops, rotations, h_ops, twin_fraction, nthreads)); Calibrate(result); return result; diff --git a/image_analysis/scale_merge/TwinningAnalysis.h b/image_analysis/scale_merge/TwinningAnalysis.h index 4630fd175..f5eb7d4e4 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.h +++ b/image_analysis/scale_merge/TwinningAnalysis.h @@ -233,7 +233,7 @@ struct TwinImmuneZoneResult { double anisotropy_delta_b_A2 = 0.0; // the anisotropy taken out of the intensities before normalising double d_max_A = 0.0; // resolution range of the shells read double d_min_A = 0.0; - bool tncs_normalised = false; // normalised within the pseudo-translation's two phase classes + bool tncs_normalised = false; // the pseudo-translation's modulation divided out per shell }; // `twin_fraction` is what the centric hypothesis is read at. A twin law OTHER than the operators a zone diff --git a/tests/TwinningAnalysisTest.cpp b/tests/TwinningAnalysisTest.cpp index f2e9b1cdc..fdfbbeeac 100644 --- a/tests/TwinningAnalysisTest.cpp +++ b/tests/TwinningAnalysisTest.cpp @@ -13,6 +13,7 @@ #include #include "../image_analysis/scale_merge/TwinningAnalysis.h" +#include "../image_analysis/scale_merge/TranslationalNCS.h" #include "SyntheticMergedReflections.h" namespace { @@ -431,6 +432,31 @@ TEST_CASE("Twin-immune zones are read at the twin fraction of the other twin law CHECK(AnalyzeTwinImmuneZone(twin, hexagonal, p321, p3).zones[0].calibrated_evidence_nats < -20.0); } +// A pseudo-translation splits every shell into two populations of different mean; the zone is read +// with that modulation divided out. A thin zone - a few hundred reflections - split into the two +// classes and normalised per class had one resolution bin per class, read the fall-off as spread, and +// a 321 crystal twinned by its 622 law read its 2-folds' zone centric. +TEST_CASE("Twin-immune zone of a twin with a pseudo-translation", "[twinning][tncs]") { + const gemmi::SpaceGroup &p622 = gemmi::get_spacegroup_by_name("P 6 2 2"); + const gemmi::SpaceGroup &p321 = gemmi::get_spacegroup_by_name("P 3 2 1"); + const gemmi::UnitCell hexagonal(58, 58, 42, 90, 90, 120); + const std::vector> tncs{{{2.0 / 3.0, 1.0 / 3.0, 0.5}}}; + + // As the verdict reads it: on the 321 merge, one reflection per unique. + const gemmi::ReciprocalAsu asu(&p321); + std::vector merge; + for (auto r : WilsonMerge("P 3 2 1", "P 6 2 2", hexagonal, 0.3)) + if (asu.is_in({{r.h, r.k, r.l}})) { + r.I = static_cast(r.I * (1.0 + 0.6 * TncsClass(r.h, r.k, r.l, tncs))); + merge.push_back(r); + } + + const auto z = AnalyzeTwinImmuneZone(merge, hexagonal, p622, p321, &tncs); + REQUIRE(z.zones.size() == 1); + CHECK(z.zones[0].n < 400); + CHECK(z.zones[0].calibrated_evidence_nats < -20.0); +} + // Averaging each intensity over its orbit under a group leaves the L-test where it was when the group // is the crystal's, and narrows it towards the perfect-twin 0.375 when an operator of the group is // not - read on the same pairs. A perfect twin of the subgroup already reads 0.375 unmerged, so there