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

rugnux: twin-immune zone takes a pseudo-translation's modulation out per shell, not by splitting the zone

Browse Source 
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 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQsU5a8AJUdJPgWevAK
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed 2026-11-12 07:02:36 +01:00
1 parent 8f18e4175f
commit 23adbcb6b3
3 files changed
+104 -54
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

image_analysis/scale_merge/TwinningAnalysis.cpp
+77 -53
View File
@@ -811,7 +811,7 @@ namespace {
// <I/sigma> 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<const MergedReflection*> refl, const std::vector<std::array<double, 3>>* tncs,
size_t nthreads, double twin_fraction) {
TwinImmuneZone ReadZone(std::vector<const MergedReflection*> 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<std::pair<double, double>> 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<const MergedReflection*> 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<size_t>(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<double>(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<size_t>(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<double>(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<double> 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<MergedReflection>& strong, const std::vector<int>& shell,
const gemmi::GroupOps& gops, const std::vector<std::array<double, 3>>& tncs) {
const int n_shells = *std::max_element(shell.begin(), shell.end()) + 1;
std::vector<std::array<double, 2>> sum(n_shells, {0.0, 0.0});
std::vector<std::array<int, 2>> n(n_shells, {0, 0});
std::vector<int> 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<float>(strong[i].I / f);
strong[i].sigma = static_cast<float>(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<MergedReflection> StrongReflections(const std::vector<MergedReflection>& p1_merged,
const gemmi::UnitCell& cell, const gemmi::GroupOps& gops,
const std::vector<std::array<double, 3>>* tncs_vector,
TwinImmuneZoneResult& result) {
std::vector<const MergedReflection*> 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<MergedReflection>& strong, const gemmi::GroupOps& gops,
const std::vector<std::array<double, 3>>* tncs_vector, double twin_fraction,
size_t nthreads) {
double twin_fraction, size_t nthreads) {
std::vector<const MergedReflection*> 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<MergedReflection>& strong,
const gemmi::GroupOps& gops, const std::vector<gemmi::Op>& rotations,
const std::vector<gemmi::Op>& h_ops,
const std::vector<std::array<double, 3>>* tncs_vector, double twin_fraction,
size_t nthreads) {
const std::vector<gemmi::Op>& h_ops, double twin_fraction, size_t nthreads) {
auto contains = [&](const std::vector<gemmi::Op>& 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<MergedReflection>&
const gemmi::GroupOps gops = group.operations();
const std::vector<gemmi::Op> rotations = ProperRotations(gops);
const std::vector<gemmi::Op> 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<MergedReflection>&
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;
image_analysis/scale_merge/TwinningAnalysis.h
+1 -1
View File
@@ -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
tests/TwinningAnalysisTest.cpp
+26
View File
@@ -13,6 +13,7 @@
#include <vector>
#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<std::array<double, 3>> 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<MergedReflection> 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<float>(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