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