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

Merge branch 'rc186-tncs-x' into rc186-merge

Browse Source
This commit is contained in:
leonarski_f committed 2026-11-12 00:27:01 +01:00
2 parents 16ab6d75bb f9a94e2c22
commit 8b6bee76dd
8 files changed
+421 -140
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

docs/CPU_DATA_ANALYSIS_DECISIONS.md
+1 -1
View File
@@ -68,7 +68,7 @@ A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$ — 0.500
Beside them, on rotation data with `--finalist-ledger`, the run reports **twin-immune zone evidence** for the operators the adopted point group adds over each of its index-2 subgroups, read on a full-resolution merge of the observations in P1 that the option makes for the purpose. Under a twin law $T$, $I_\mathrm{obs}(h)=(1-\alpha)I(h)+\alpha I(Th)$; a reflection whose twin mate is itself up to the subgroup and Friedel is untouched at every $\alpha$, and those are exactly the reflections centric in the group but acentric in the subgroup. They read centric ($\langle|E^2-1|\rangle = 0.968$) if the added operators are real and acentric (0.736) if they are a twin law or a pseudo-symmetry — the one intensity statistic that still separates the two at $\alpha = 0.5$, where every operator statistic reads "real". Each zone is normalised against its own mean in resolution bins (and within the two phase classes of a detected pseudo-translation), read only in shells with $\langle I/\sigma\rangle \ge 5$ because noise inflates every class towards centric, and reported with $n$, a standard error and the centric-over-acentric Wilson log-likelihood ratio in nats (both densities convolved with each reflection's measurement error), beside the acentric control. It is read absolutely, never as the difference to the control: a perfect twin's control (0.541) plus noise inflation would otherwise look like true symmetry. An absolute reading is only as good as the normalisation, though, and the acentric control certifies it: an acentric population reads $-0.130$ nats per reflection when the normalisation is right, and whatever the control reads above that is the normalisation's — anisotropy, a pseudo-translation, a pseudo-centring, noise all inflate every class towards centric alike — so the zone, normalised the same way, carries the same per reflection and the **calibrated** evidence has it taken off (a twinned control reads *below* the expectation, so nothing is taken off a twin). The centric side is read **twinned at the fraction the lattice's other operations show** by their own correlation — the strongest CC of a lattice rotation outside the group, relative to the mean CC of the group's own operators, inverted through $\rho = 2\alpha(1-\alpha)/((1-\alpha)^2+\alpha^2)$ — because a twin by another law reaches a genuine zone exactly as it reaches the control: a genuine 321 crystal twinned by a 622 law read its 2-folds at −206 nats against the untwinned centric density and +209 at that law's fraction of 0.07. The operators' own twin law cannot reach their zone, so the acentric side stays untwinned, and the control's expectation is taken at the same fraction. Measured, a 6/m crystal with a 67 Ų anisotropy read its control at $+0.10$ nats per reflection and the zone of a refused 622 the same, so the zone's $+124$ nats were the normalisation's; calibrated it reads $-115$, and genuine promotions keep $+200$ and above. The calibrated zone is what a refused promotion is decided on (above). It is an in-house method: published practice (Yeates' $H$-test, xtriage) excludes these reflections as uninformative about the twin fraction, and no published test was found that uses them positively; the one precedent here is a trigonal crystal whose twin-immune zone read centric and whose higher group was then confirmed by refinement.
**Translational pseudo-symmetry** — two copies of the contents of the asymmetric unit related by a pure translation that is not a lattice vector — is looked for on every merging run, because it is the classic predictor of a failed molecular replacement, and because it raises the second moment where twinning lowers it, so each can mask the other's test. The detection is the native-Patterson route of Read, Adams & McCoy (see the [references](CPU_DATA_ANALYSIS.md#references)): the largest off-origin peak of a Patterson computed from the merged intensities, as a fraction of the origin peak, with the peak vector then refined against the intensity modulation it should produce — the ratio of the strongest to the weakest bin mean of $\langle E^2\rangle$ over the phase $\mathrm{frac}(\mathbf{h}\cdot\mathbf{u})$. Both halves are scored against a null computed for the crystal at hand rather than against a fixed bound — both against the same intensities **permuted within resolution shells** — the peak against that map, the modulation against the same measurement made on the shuffled intensities and started where the measurement starts — because the noise floor of the peak statistic spans an order of magnitude across data sets, so no fixed percentage means the same thing twice. The shuffle leaves the reflection count, the $E^2$ distribution and the phase-bin populations untouched and removes only the correlation between a reflection's intensity and $\mathbf{h}\cdot\mathbf{u}$, which is the thing the modulation claims to see. Re-running the greedy search from random starting vectors over the *same* intensities does not work as a control: a greedy search started anywhere walks into a real modulation's basin, so it measures the search **and** whatever modulation the crystal carries — measured over 110 merged datasets, single draws of that control reached 199× and 6162× on crystals whose own modulation reads 4.0× and 67×, and on 28 of the 110 it came out at or above the signal it is subtracted from. Requiring both halves is what keeps the false-positive rate down; either alone over-calls by about a factor of two. A translation the merged data are **exactly** invariant under is reported as an undeclared lattice translation instead — a translation the data are exactly invariant under is a lattice vector by definition, so the centring or the cell is wrong, not the packing — and the pseudo-symmetry search continues underneath it, so a real pseudo-translation sitting under an undeclared centring is still found. The finding itself is report-only (the `TNCS_*` keys of `<prefix>_report.txt`): it gates nothing and changes no reflection, no scale and no group — but the axial-absence test (§13.1) and the $L$-test here both correct themselves against the modulation it measures.
**Translational pseudo-symmetry** — two copies of the contents of the asymmetric unit related by a pure translation that is not a lattice vector — is looked for on every merging run, because it is the classic predictor of a failed molecular replacement, and because it raises the second moment where twinning lowers it, so each can mask the other's test. The detection is the native-Patterson route of Read, Adams & McCoy (see the [references](CPU_DATA_ANALYSIS.md#references)): the largest off-origin peak of a Patterson computed from the merged intensities over 20–5 Å, as a fraction of the origin peak, with the peak vector then refined against the intensity modulation it should produce — the correlation between $E^2$ and $\cos 2\pi\,\mathbf{h}\cdot\mathbf{u}$ over the acentric reflections beyond 4 Å (or the 2000 lowest-resolution ones where a small cell has fewer there). Both halves are scored against a null computed for the crystal at hand rather than against a fixed bound — both against the same intensities **permuted within resolution shells** — the peak against that map, the modulation against the same search made on the shuffled intensities and started where the search starts — because the noise floor of the peak statistic spans an order of magnitude across data sets, so no fixed percentage means the same thing twice. The shuffle leaves the reflection count and the $E^2$ distribution untouched and removes only the correlation between a reflection's intensity and $\mathbf{h}\cdot\mathbf{u}$, which is the thing the modulation claims to see; the shells are split by the reciprocal axis $\mathbf{h}$ lies closest to, so that anisotropy, which the spherical normalisation leaves in $E^2$ as a dependence on direction, is neither read as a modulation nor left out of the null, and the search stays 15 Å from every lattice point, where the phase would sort the reflections by exactly that direction. Re-running the greedy search from random starting vectors over the *same* intensities does not work as a control: a greedy search started anywhere walks into a real modulation's basin, so it measures the search **and** whatever modulation the crystal carries. A vector that modulates is a fact about the crystal; the **pseudo-translation reported** is the one molecular replacement has to be told about, and the programs that correct for it act on a Patterson peak of a fifth of the origin — read in the band, or as the Patterson height the modulation's depth gives at the vector, fitted over every reflection of the modulation population (undiluted by a twin), which on a small cell or one whose low-resolution reflections were lost to overloads is the better measure of how much of the scattering the translation relates. Below that nothing in the data tells a weak pseudo-translation from two rotated copies that look alike at 20–5 Å (over the deposited models of a 200-dataset corpus the models' own Pattersons read 18–80 % where a pseudo-translation is there and 4–15 % where it is not), and nothing downstream needs it, so it is noted and not flagged. A peak of that size the data cannot establish — a 30 Å cell has about 150 reflections in the band, where the noise floor is itself a fifth of the origin — is carried as `POSSIBLE` rather than dropped. A translation the merged data are **exactly** invariant under is reported as an undeclared lattice translation instead — a translation the data are exactly invariant under is a lattice vector by definition, so the centring or the cell is wrong, not the packing — and the pseudo-symmetry search continues underneath it, so a real pseudo-translation sitting under an undeclared centring is still found. The finding itself is report-only (the `TNCS_*` keys of `<prefix>_report.txt`): it gates nothing and changes no reflection, no scale and no group — but the axial-absence test (§13.1) and the $L$-test here both correct themselves against the modulation it measures, and the vectors are carried onto the axes the files are written on when a model or a reference moves them after the analysis.
**The $L$-test partners are chosen so a pseudo-translation cannot silence it.** The test compares each reflection with a partner a fixed index step away, and a pseudo-translation biases $\langle|L|\rangle$ upwards unless the partner shares its modulation class. The ordinary step of 2 preserves the class of a half-integer translation — the pseudo-centering the statistic's authors call it robust to — but not of one of a third, which inflates $\langle|L|\rangle$ past the bound that is read as evidence *against* twinning and silently loses the twin call on a crystal that has one. Where a pseudo-translation is detected, the partner steps are restricted to those that preserve its class; where no step does, the statistic is dropped from the twin verdict **in both directions** — it can no longer indicate a twin and can no longer be read as proof that there is none — and the second moment decides alone. `L_TEST_VS_TNCS=` in the report says which of the three happened.
docs/RUGNUX_REPORT.md
+44 -27
View File
@@ -971,40 +971,57 @@ puts on the intensities is not in the search model. Section 4 reports it beside
two interact. The algorithm is in
[CPU/GPU data analysis ▸ Twinning and translational pseudo-symmetry](CPU_DATA_ANALYSIS_DECISIONS.md#132-twinning-check-and-translational-pseudo-symmetry).
**`TNCS_DETECTED=`** is `TRUE`, `FALSE`, `INCONCLUSIVE` or `NOT_MEASURED`. The last two are not
`FALSE`: `NOT_MEASURED` means the merge has too few reflections in 20–5 Å to compute the statistic at
all, `INCONCLUSIVE` means the Patterson was measured but too few acentric reflections remain to test
whether the vector it names modulates the intensities. Neither is a statement that the crystal has no
pseudo-symmetry.
**`TNCS_DETECTED=`** is `TRUE`, `POSSIBLE`, `FALSE`, `INCONCLUSIVE` or `NOT_MEASURED`. The last two are
not `FALSE`: `NOT_MEASURED` means the merge has too few reflections in 20–5 Å to compute the statistic
at all (or the cell is too small for a translation between molecules), `INCONCLUSIVE` means the
Patterson was measured but too few acentric reflections remain to test whether the vector it names
modulates the intensities. Neither is a statement that the crystal has no pseudo-symmetry.
A warning under the `PSEUDO_TRANSLATION` flag is given wherever `TNCS_DETECTED= TRUE`. Below a Patterson
peak of 20 % of the origin (the height phenix.xtriage flags) the detection is significant but its
modulation is small; the warning and the summary call it weak and ask to check whether molecular
replacement needs it.
A warning under the `PSEUDO_TRANSLATION` flag is given wherever `TNCS_DETECTED= TRUE` or `POSSIBLE`.
`TRUE` requires **both** of two tests, because either alone over-calls by about a factor of two:
`TRUE` is a vector that **modulates** the intensities and a Patterson peak at it of at least **20 %**
of the origin, the height at which molecular-replacement programs correct for a pseudo-translation
(Phaser's default, and the one phenix.xtriage flags). Below that nothing in the data tells a weak
pseudo-translation from two rotated copies that look alike at 20–5 Å, and nothing downstream needs
it; such a vector is noted (`TNCS_MODULATES= TRUE`) and not flagged.
- **`TNCS_PATTERSON_PEAK_PCT=`** — the largest off-origin peak of the native Patterson, as a
percentage of the origin peak, counting only peaks farther than 15 Å from any origin-equivalent
- **`TNCS_PATTERSON_PEAK_PCT=`** — the largest off-origin peak of the native Patterson over 20–5 Å, as
a percentage of the origin peak, counting only peaks farther than 15 Å from any origin-equivalent
lattice point. **`TNCS_PATTERSON_PEAK_Z=`** scores it against **`TNCS_PATTERSON_NULL_PCT=`**, the
same map recomputed with the intensities permuted within resolution shells. The null is per dataset
and not a table: the noise floor of this statistic runs from about 1% on a large merge to about 18%
and not a table: the noise floor of this statistic runs from about 1 % on a large merge to about 18 %
on a small one, so no fixed percentage separates the two populations.
- **`TNCS_MODULATION=`** — the ratio of the strongest to the weakest bin mean of ⟨E²⟩ over the phase
frac(**h**·**u**), with **u** the refined peak vector (**`TNCS_VECTOR=`**, fractional, and
**`TNCS_VECTOR_LENGTH=`** in Å). **`TNCS_MODULATION_NULL=`** is the same search started from random
vectors, so the contrast the search itself can manufacture is measured rather than assumed.
- **`TNCS_MODULATION_AVERAGED=`** / **`TNCS_MODULATION_AVERAGED_NULL=`** — written only where the
Patterson peak is significant but the vector's own modulation is not: the same test with each
reflection's class the mean of cos 2π**h**·R**u** over the lattice's rotations R (its symmetry and
every twin law the metric admits), against the same permutation null. A symmetry operator or twin
law that carries **u** to a vector other than ±**u** averages the vector's own classes - a perfect
twin does so by construction - and this is the test that still sees the modulation then.
**`TNCS_CLASSES=`** (`VECTOR` or `LATTICE_AVERAGED`) says which of the two a detection stands on; the
L-test partner steps then keep the class of every image R**u**.
- **`TNCS_VECTOR=`** (fractional; **`TNCS_VECTOR_LENGTH=`** in Å) is the peak vector refined against
the intensities, and **`TNCS_MODULATION_DEPTH_PCT=`** the depth of the modulation along it as a
Patterson height - the amplitude *a* of ⟨E²⟩ = 1 + *a* cos 2π**h**·**u** fitted over every acentric
reflection of the modulation population, times ⟨cos² 2π**h**·**u**⟩, which is what the Patterson of
an untwinned crystal reads at **u** (*a* at a half-integer vector, *a*/2 at a general one); under a
twin it is the height the crystal has rather than the diluted one the twinned Patterson shows. On a small cell, or one whose low-resolution reflections were lost to overloads, the band is
thin and this is the better measure of how much of the scattering the translation relates. `TRUE`
needs either of the two to reach 20 %.
- **`TNCS_MODULATES=`** — whether the vector modulates the intensities: **`TNCS_CORRELATION=`**, the
correlation between E² and cos 2π**h**·**u** over the acentric reflections beyond 4 Å (or the 2000
lowest-resolution ones where fewer lie there), must stand five spreads above
**`TNCS_CORRELATION_NULL=`**, the same search run on the intensities shuffled within their
resolution shells (median and spread over 20 draws). E² is normalised, and shuffled, within each
shell split by the reciprocal axis **h** lies closest to, so that anisotropy is not read as a
modulation along a vector near a lattice point. **`TNCS_MODULATION=`** is the profile's depth, the
ratio of the strongest to the weakest phase-class mean of ⟨E²⟩, and **`TNCS_MODULATION_NULL=`** what
the shuffled searches reach; descriptive, not the decision.
- **`TNCS_CORRELATION_AVERAGED=`** / **`TNCS_MODULATION_AVERAGED=`** (and their `_NULL`) — written
only where the Patterson peak is significant but the vector's own modulation is not: the same test
with each reflection's class the mean of cos 2π**h**·R**u** over the lattice's rotations R (its
symmetry and every twin law the metric admits), against the same permutation null. A symmetry
operator or twin law that carries **u** to a vector other than ±**u** averages the vector's own
classes - a perfect twin does so by construction - and this is the test that still sees the
modulation then. **`TNCS_CLASSES=`** (`VECTOR` or `LATTICE_AVERAGED`) says which of the two a
detection stands on; the L-test partner steps then keep the class of every image R**u**.
The vector is good to about 0.05 fractional. It is a starting point for a program that refines it,
not a refined result: the refinement maximises the modulation, not the accuracy of the vector.
`POSSIBLE` is a peak of 20 % or more that these data cannot establish: too few reflections in the band
for the permutation floor to lie far below it (a 30 Å cell has about 150 there) and no significant
modulation. The hypothesis is reported with what it rests on rather than dropped.
The vector is refined on the whole data, in order of **h**, and is good to about 0.01 fractional.
**`TNCS_PSEUDO_CENTRED=`** and **`TNCS_SUBLATTICE=`** are separate claims and are deliberately not
merged. A vector that is a rational translation `1/q` of the cell means the crystal is pseudo-centred;
image_analysis/scale_merge/TranslationalNCS.cpp
+209 -53
View File
@@ -42,6 +42,10 @@ namespace {
constexpr double kModulationDMin = 4.0; // coarse enough that a small error in u does not
// scramble the phase h.u
// Extended to higher resolution where fewer than this many acentric reflections lie beyond
// 4 A, as the Patterson band is: a small cell's |h| at 2.5 A is a protein cell's at 4 A, so
// the phase is no more sensitive to an error in u there.
constexpr int kMinModulationBand = 2000;
constexpr int kModulationShells = 10;
constexpr int kPhaseBins = 12;
constexpr int kMinModulationReflections = 500;
@@ -53,13 +57,25 @@ namespace {
// per trial, so the extra draws cost a few ms.
constexpr int kModulationNullDraws = 20;
// Gate constants. Read off a corpus of 128 merged datasets, on which the pair fires on 7.0% of
// them - a pathology rate, and close to the ~6% of PDB entries Read, Adams & McCoy report with a
// large off-origin Patterson peak. Either criterion ALONE fires on roughly twice as many: a 40%
// peak can carry no modulation at all, and on a 70,000-reflection dataset a 4% peak reaches
// z = 18 while modulating nothing.
// Significance. The peak against its permutation null, and the modulation against its own:
// either alone over-calls - a 40% peak can carry no modulation at all (an undeclared centring's
// shoulder, an anomalous scatterer's self-vector), and on a 70,000-reflection dataset a 4% peak
// reaches z = 18 while modulating nothing.
constexpr double kPeakZThreshold = 5.0;
constexpr double kModulationRatioThreshold = 2.5;
constexpr double kModulationZThreshold = 5.0;
// Size. A significant modulation is a fact about the crystal; a PSEUDO-TRANSLATION is the one
// molecular replacement has to be told about, and the programs that correct for it act on a
// Patterson peak of a fifth of the origin or more (Phaser's default; Read, Adams & McCoy (2013)
// Acta Cryst. D69, 176-183). Below that the two are not told apart by anything in the data: a
// 10% peak is what two copies of a quarter of the cell related by a near-translation give, and
// also what two rotated copies of a globular molecule give at 20-5 A, where any blob looks like
// any other. Over the deposited models of a 200-dataset corpus the two populations divide here:
// the models' own Patterson reads 18-80% where a pseudo-translation is there and 4-15% where
// it is not, and the calls made between 7 and 12% of the origin were all of the second kind.
constexpr double kDeclaredPeakPct = 20.0;
// A peak of the declared size standing this far above the permutation floor, with the
// modulation unable to confirm it, is reported as POSSIBLE rather than dismissed.
constexpr double kSuspectedPeakZ = 2.0;
// A modulation this deep means the suppressed class carries under 5% of the other, i.e. it is
// very nearly systematically absent - which is what a LATTICE translation of a smaller cell looks
// like measured in the larger one. Over 137 datasets exactly one crystal reaches it (at 89x, the
@@ -329,6 +345,16 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
}
};
build_near_mask();
if (std::count(near.begin(), near.end(), 0) == 0) {
std::ostringstream os;
os << "every point of the cell lies within " << kMinPeakDistanceA
<< " A of a lattice point - a small-molecule cell, with no room for a translation between molecules";
result.refusal = os.str();
fftwf_destroy_plan(plan);
fftwf_free(coeff);
fftwf_free(map);
return result;
}
auto transform = [&](const std::vector<double> &values) {
std::memset(coeff, 0, 2 * n_complex * sizeof(float));
@@ -451,9 +477,14 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
std::vector<int> mod_idx;
for (size_t i = 0; i < merged.size(); ++i) {
const auto &r = merged[i];
if (std::isfinite(r.I) && std::isfinite(r.d) && r.d >= kModulationDMin && acentric(r))
if (std::isfinite(r.I) && std::isfinite(r.d) && acentric(r))
mod_idx.push_back(static_cast<int>(i));
}
std::sort(mod_idx.begin(), mod_idx.end(),
[&](int a, int b) { return merged[a].d != merged[b].d ? merged[a].d > merged[b].d : a < b; });
const int n_beyond = static_cast<int>(std::count_if(mod_idx.begin(), mod_idx.end(),
[&](int i) { return merged[i].d >= kModulationDMin; }));
mod_idx.resize(std::max(n_beyond, std::min<int>(kMinModulationBand, static_cast<int>(mod_idx.size()))));
if (!mod_idx.empty()) {
// Shells linear in 1/d^2, and the same defensive normalisation the twinning statistics
// use: noise-only shells skipped, Wilson outliers rejected, the shell mean re-fitted once
@@ -465,11 +496,21 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
max_s = std::max(max_s, s);
}
if (max_s > min_s) {
std::vector<std::vector<int>> bins(kModulationShells);
// Shells linear in 1/d^2, each split three ways by the reciprocal axis h lies closest
// to. E^2 is normalised, and later shuffled, within these bins rather than within
// whole shells: under anisotropy a shell's E^2 depends on the direction of h, and a
// vector a little off a lattice point sorts the reflections by exactly that (the
// phase 2 pi h.(L + delta) is 2 pi h.delta), so the correlation read a modulation
// where there was a direction. A pseudo-translation's classes interleave within
// every such bin and lose nothing to it.
std::vector<std::vector<int>> bins(3 * kModulationShells);
for (int i : mod_idx) {
const double t = (1.0 / (merged[i].d * merged[i].d) - min_s) / (max_s - min_s);
bins[std::clamp(static_cast<int>(t * kModulationShells), 0, kModulationShells - 1)]
.push_back(i);
const int shell = std::clamp(static_cast<int>(t * kModulationShells), 0, kModulationShells - 1);
const double along[3] = {std::fabs(merged[i].h) * cell.ar, std::fabs(merged[i].k) * cell.br,
std::fabs(merged[i].l) * cell.cr};
const int sector = static_cast<int>(std::max_element(along, along + 3) - along);
bins[3 * shell + sector].push_back(i);
}
for (const auto &b : bins) {
if (b.empty())
@@ -511,6 +552,11 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
}
}
result.n_modulation_reflections = static_cast<int>(mod_e2.size());
// A peak of the declared size that the data cannot establish: the band holds too few reflections
// for the permutation floor to be far below it, and the modulation (below) is not significant
// either, or cannot be measured at all. The hypothesis is carried rather than dropped - the user
// is told what was seen and what it rests on.
result.suspected = result.peak_percent >= kDeclaredPeakPct && result.peak_z >= kSuspectedPeakZ;
if (result.n_modulation_reflections < kMinModulationReflections)
return result;
@@ -563,20 +609,39 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
// Coarse-to-fine local search for the vector that modulates the data most strongly, started from
// the Patterson peak. The Patterson grid is coarse and a real pseudo-translation is rarely at an
// exactly rational position, so the peak position alone understates the modulation.
//
// The objective is the correlation between E^2 and cos(2 pi h.u), not the max/min ratio of the
// phase bins the profile is reported with. A ratio of bin means is owned by its thinnest bin,
// and a vector a little off a special position fills the thin bins with the highest orders along
// one axis - so that under anisotropy, which the spherical shell normalisation leaves in E^2
// as a dependence on the direction of h, the ratio reads a modulation that is not there: 6.5x
// on one crystal whose amplitude along that vector is 0.06, and 4x on another at its own
// centring vector. The correlation reads every reflection with the same weight and neither.
//
// The search stays where the Patterson search looked: at least kMinPeakDistanceA from every
// origin-equivalent point. At a lattice vector every reflection is in the same class, and from
// just beside one the classes sort the reflections by the order of h along an axis - so a
// greedy search can walk onto a centring vector and read the anisotropy there as a modulation
// (measured: a C-centred crystal called on its own centring vector, with 99% of the origin at it).
const std::array<double, 3> step0{1.0 / nx, 1.0 / ny, 1.0 / nz};
auto between_molecules = [&](const std::array<double, 3> &v) {
return DistanceToOrigin(cell, Vec3{v[0], v[1], v[2]}, equivalents) >= kMinPeakDistanceA;
};
auto refine_over = [&](const std::vector<double> &e2v, std::array<double, 3> u,
const std::vector<Rot3> &rots) {
std::array<double, 3> s = step0;
double best = contrast_over(e2v, u, nullptr, rots);
double best = CosineCorrelation(mod_hkl, e2v, u, e2v.size(), rots);
for (int round = 0; round < 3; ++round) {
for (bool improved = true; improved;) {
improved = false;
for (int i = 0; i < 3; ++i)
// No step of 0: it is u itself, which cannot improve on its own contrast.
// No step of 0: it is u itself, which cannot improve on its own score.
for (double gstep : {-1.0, -0.5, 0.5, 1.0}) {
std::array<double, 3> v = u;
v[i] += gstep * s[i];
const double c = contrast_over(e2v, v, nullptr, rots);
if (!between_molecules(v))
continue;
const double c = CosineCorrelation(mod_hkl, e2v, v, e2v.size(), rots);
if (std::isfinite(c) && c > best + 1e-6) {
best = c;
u = v;
@@ -618,42 +683,66 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
// The search on the data is made beside its nulls rather than ahead of them: none of them reads
// another's answer, and on a large merge the one on the data alone took as long as all the nulls
// on the workers together.
// Returns the search on the data - its vector, correlation and the profile's contrast there - and
// the null's median correlation and spread, with the median contrast at the nulls' own vectors.
// The correlation is significant where it stands kModulationZThreshold spreads above what the
// search makes out of nothing.
struct SearchNull {
std::array<double, 3> u;
double correlation, contrast, null_median, null_sd, null_contrast;
};
auto search_and_null = [&](const std::vector<Rot3> &rots) {
std::vector<std::pair<std::array<double, 3>, double>> searches(kModulationNullDraws + 1);
ParallelFor(kModulationNullDraws + 1, nthreads, [&](int i) {
searches[i] = refine_over(i == 0 ? mod_e2 : permuted[i - 1], result.peak_frac, rots);
});
std::vector<double> nulls(kModulationNullDraws);
for (int draw = 0; draw < kModulationNullDraws; ++draw)
nulls[draw] = std::isfinite(searches[draw + 1].second) ? searches[draw + 1].second : 1.0;
std::vector<double> nulls(kModulationNullDraws), contrasts(kModulationNullDraws);
for (int draw = 0; draw < kModulationNullDraws; ++draw) {
nulls[draw] = searches[draw + 1].second;
const double c = contrast_over(permuted[draw], searches[draw + 1].first, nullptr, rots);
contrasts[draw] = std::isfinite(c) ? c : 1.0;
}
std::sort(nulls.begin(), nulls.end());
std::sort(contrasts.begin(), contrasts.end());
const size_t m = nulls.size();
return std::make_pair(searches[0], m % 2 ? nulls[m / 2] : 0.5 * (nulls[m / 2 - 1] + nulls[m / 2]));
const double median = m % 2 ? nulls[m / 2] : 0.5 * (nulls[m / 2 - 1] + nulls[m / 2]);
double var = 0.0;
for (double v : nulls)
var += (v - median) * (v - median);
return SearchNull{searches[0].first, searches[0].second,
contrast_over(mod_e2, searches[0].first, nullptr, rots), median,
std::sqrt(var / m), m % 2 ? contrasts[m / 2] : 0.5 * (contrasts[m / 2 - 1] + contrasts[m / 2])};
};
const auto [vector_search, vector_null] = search_and_null(no_rotations);
const auto [u_refined, best_contrast] = vector_search;
if (!std::isfinite(best_contrast))
return result;
auto significant = [](const SearchNull &s) {
return s.null_sd > 0.0 && (s.correlation - s.null_median) / s.null_sd >= kModulationZThreshold;
};
const SearchNull vector_search = search_and_null(no_rotations);
const std::array<double, 3> u_refined = vector_search.u;
for (int i = 0; i < 3; ++i)
result.vector_frac[i] = u_refined[i] - std::floor(u_refined[i]);
result.modulation_measured = true;
result.modulation = best_contrast;
result.correlation = vector_search.correlation;
result.correlation_null = vector_search.null_median;
result.correlation_null_sd = vector_search.null_sd;
result.modulation = vector_search.contrast;
result.modulation_null = vector_search.null_contrast;
contrast(mod_e2, u_refined, &result.phase_bin_mean_e2);
result.vector_length_A = DistanceToOrigin(cell, Vec3{result.vector_frac[0], result.vector_frac[1],
result.vector_frac[2]},
std::vector<Vec3>{{0.0, 0.0, 0.0}});
result.modulation_null = vector_null;
result.detected = result.peak_z >= kPeakZThreshold && result.modulation_null > 0.0 &&
result.modulation / result.modulation_null >= kModulationRatioThreshold;
result.modulates = significant(vector_search);
// The second hypothesis: classes averaged over the lattice's rotations (`modulation_averaged` in
// the header says why). The modulation model per copy pair is Read, Adams & McCoy's (2013) Acta
// Cryst. D69, 176-183; a symmetry operator or a twin law adds the pair its rotation relates, and
// the measured intensity sums them. Asked only where it can change the verdict.
// the measured intensity sums them. Asked wherever the peak is significant and the lattice has
// rotations, and taken where its classes follow the intensities better than the vector's own:
// under a perfect twin the vector's classes still correlate, only less, and the averaged ones
// are the physical ones.
std::vector<Rot3> lattice_rotations;
std::array<double, 3> u_averaged = u_refined;
if (!result.detected && result.peak_z >= kPeakZThreshold) {
if (result.peak_z >= kPeakZThreshold) {
const gemmi::GroupOps lattice = gemmi::find_lattice_symmetry(
cell, space_group ? space_group->centring_type() : 'P', kLatticeObliquityDeg);
for (const auto &op : lattice.sym_ops) {
@@ -668,20 +757,22 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
if (lattice_rotations.size() > 1) {
result.averaged_measured = true;
result.averaged_rotations = static_cast<int>(lattice_rotations.size());
const auto [averaged_search, averaged_null] = search_and_null(lattice_rotations);
std::tie(u_averaged, result.modulation_averaged) = averaged_search;
result.modulation_averaged_null = averaged_null;
result.detected_on_averaged_classes =
std::isfinite(result.modulation_averaged) && result.modulation_averaged_null > 0.0
&& result.modulation_averaged / result.modulation_averaged_null >= kModulationRatioThreshold;
result.detected = result.detected_on_averaged_classes;
const SearchNull averaged_search = search_and_null(lattice_rotations);
u_averaged = averaged_search.u;
result.correlation_averaged = averaged_search.correlation;
result.correlation_averaged_null = averaged_search.null_median;
result.correlation_averaged_null_sd = averaged_search.null_sd;
result.modulation_averaged = averaged_search.contrast;
result.modulation_averaged_null = averaged_search.null_contrast;
result.detected_on_averaged_classes = significant(averaged_search)
&& averaged_search.correlation > vector_search.correlation;
result.modulates = result.modulates || result.detected_on_averaged_classes;
}
}
const std::vector<Rot3> &class_rotations = result.detected_on_averaged_classes ? lattice_rotations
: no_rotations;
const std::array<double, 3> u_start = result.detected_on_averaged_classes ? u_averaged : u_refined;
result.near_extinct_class = result.detected && result.modulation >= kNearExtinctModulation;
// ------------------------------------------------ the vector the phase CLASSES are assigned with
//
@@ -703,9 +794,10 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
// resolution range whose objective has a local maximum every 1/|h| in each direction. Refining on
// the whole data in one go instead walks into the nearest of those.
//
// Nothing above is re-read: `modulation`, `modulation_null` and `detected` stay the gate's, so
// which crystals are called stays exactly as it was. Only where the vector points changes.
if (result.detected) {
// The gate's correlation and its null are not re-read: whether the vector modulates stays
// decided where its control was drawn. The profile and the contrast are re-read off the vector
// refined here, which is the one reported.
if (result.modulates) {
std::vector<std::array<int, 3>> cls_hkl;
std::vector<double> cls_e2;
{
@@ -774,7 +866,10 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
std::array<double, NEIGHBOURS> f{};
for (int j = 0; j < NEIGHBOURS;) {
ParallelFor(NEIGHBOURS - j, nthreads, [&](int i) {
f[j + i] = CosineCorrelation(cls_hkl, cls_e2, neighbour(j + i), n, class_rotations);
const auto v = neighbour(j + i);
f[j + i] = between_molecules(v)
? CosineCorrelation(cls_hkl, cls_e2, v, n, class_rotations)
: -std::numeric_limits<double>::infinity();
});
int k = j;
while (k < NEIGHBOURS && !(f[k] > best + 1e-9))
@@ -803,7 +898,7 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
u = u_start;
for (int i = 0; i < 3; ++i)
result.vector_frac[i] = u[i] - std::floor(u[i]);
contrast(mod_e2, u, &result.phase_bin_mean_e2);
result.modulation = contrast(mod_e2, u, &result.phase_bin_mean_e2);
result.vector_length_A = DistanceToOrigin(cell, Vec3{result.vector_frac[0], result.vector_frac[1],
result.vector_frac[2]},
std::vector<Vec3>{{0.0, 0.0, 0.0}});
@@ -811,6 +906,33 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
result.class_vectors = Images(result.vector_frac, class_rotations);
// The modulation's depth as a Patterson height: the amplitude a of <E^2> = 1 + a c(h), fitted by
// least squares over every reflection of the modulation population (c the class: cos 2 pi h.u,
// or its mean over the lattice's rotations), times <cos^2 2 pi h.u> - which is what the Patterson
// of an untwinned crystal reads at u: a at a half-integer vector, a/2 at a general one. Under a
// twin it is the height the crystal has rather than the diluted one the twinned Patterson shows.
// Where the band is thin (a small cell, low-resolution reflections lost to overloads) the band's
// own peak is a poor measure of how much of the scattering the translation relates, and a
// pseudo-centring modulating every shell to atomic resolution reads a fraction of its height there.
{
double num = 0.0, den = 0.0, cos2 = 0.0;
for (size_t i = 0; i < mod_e2.size(); ++i) {
const double c = TncsClass(mod_hkl[i][0], mod_hkl[i][1], mod_hkl[i][2], result.class_vectors);
num += (mod_e2[i] - 1.0) * c;
den += c * c;
const double cu = std::cos(2.0 * std::numbers::pi
* (mod_hkl[i][0] * result.vector_frac[0] + mod_hkl[i][1] * result.vector_frac[1]
+ mod_hkl[i][2] * result.vector_frac[2]));
cos2 += cu * cu;
}
result.modulation_depth_percent = den > 0.0 && !mod_e2.empty()
? 100.0 * num / den * cos2 / static_cast<double>(mod_e2.size()) : 0.0;
}
result.detected = result.modulates
&& std::max(result.peak_percent, result.modulation_depth_percent) >= kDeclaredPeakPct;
result.suspected = result.suspected && !result.modulates;
result.near_extinct_class = result.detected && result.modulation >= kNearExtinctModulation;
// Is the vector a rational translation of the reported cell? If it is, the cell contains a
// sub-lattice: the pseudo-translation is a candidate for a REAL lattice translation the indexing
// missed, and the reported cell is a supercell candidate rather than a settled result.
@@ -833,6 +955,27 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
return result;
}
void ChangeBasis(TranslationalNCSResult &result, const gemmi::Op &cob) {
const auto rotate = [&](std::array<double, 3> &v) {
std::array<double, 3> out{};
for (int i = 0; i < 3; ++i) {
out[i] = (cob.rot[i][0] * v[0] + cob.rot[i][1] * v[1] + cob.rot[i][2] * v[2]) / gemmi::Op::DEN;
out[i] -= std::floor(out[i]);
}
v = out;
};
rotate(result.peak_frac);
rotate(result.vector_frac);
for (auto &v : result.class_vectors)
rotate(v);
for (auto &v : result.undeclared_lattice_translations)
rotate(v);
const std::array<int, 3> p = result.commensurate_numerator;
for (int i = 0; i < 3; ++i)
result.commensurate_numerator[i] =
(cob.rot[i][0] * p[0] + cob.rot[i][1] * p[1] + cob.rot[i][2] * p[2]) / gemmi::Op::DEN;
}
std::string TranslationalNCSToText(const TranslationalNCSResult &result) {
std::ostringstream os;
os << "Translational pseudo-symmetry\n";
@@ -852,7 +995,7 @@ std::string TranslationalNCSToText(const TranslationalNCSResult &result) {
<< " pseudo-translation so strong it is nearly one. Excluded, and the search continued past it.\n";
}
os << " Native Patterson off-origin peak " << std::setprecision(1) << result.peak_percent
<< "% of the origin";
<< "% of the origin (" << result.n_reflections << " reflections, 20-5 A)";
if (result.vector_length_A > 0.0)
os << " at (" << std::setprecision(3) << result.vector_frac[0] << ", " << result.vector_frac[1]
<< ", " << result.vector_frac[2] << "), " << std::setprecision(1) << result.vector_length_A
@@ -861,27 +1004,40 @@ std::string TranslationalNCSToText(const TranslationalNCSResult &result) {
<< result.null_mean << " +- " << std::setprecision(2) << result.null_sd << "%]\n";
if (result.modulation_measured)
os << " Intensity modulation: " << std::setprecision(1) << result.modulation
<< "x between the strongest and weakest phase class [random-vector control "
<< result.modulation_null << "x]\n"
<< " [read on the search's own vector over d >= 4 A, which is what the control is drawn\n"
<< " against; the vector reported above is that one re-refined on the whole data, so it\n"
<< " points more accurately and this number is not re-read off it]\n";
<< "x between the strongest and weakest phase class [the same search on shuffled\n"
<< " intensities reaches " << result.modulation_null << "x]; correlation of E^2 with cos(2 pi h.u) "
<< std::setprecision(3) << result.correlation << " [control " << result.correlation_null << " +- "
<< result.correlation_null_sd << "]\n";
if (result.averaged_measured)
os << " Averaged over the lattice's " << result.averaged_rotations << " rotations (its symmetry and twin laws): "
<< std::setprecision(1) << result.modulation_averaged << "x [control " << result.modulation_averaged_null
<< "x]\n"
<< "x], correlation " << std::setprecision(3) << result.correlation_averaged << " [control "
<< result.correlation_averaged_null << " +- " << result.correlation_averaged_null_sd << "]\n"
<< " [a symmetry operator or a twin law R adds the copy pair related by R u, so where R u is not\n"
<< " +-u the vector's own classes are averaged with each other and its contrast understates it]\n";
if (!result.modulation_measured) {
os << " The modulation test could not be run - too few acentric reflections beyond 4 A - so\n"
<< " the peak above stands alone.\n"
<< " => Inconclusive. A Patterson peak on its own is not enough to call a pseudo-symmetry\n"
<< " (peaks of this size occur without one), and this dataset cannot supply the second\n"
<< " half of the evidence.\n";
os << " The modulation test could not be run - too few acentric reflections - so the peak\n"
<< " above stands alone.\n"
<< (result.suspected
? " => POSSIBLE translational pseudo-symmetry: a peak of the size molecular replacement\n"
" corrects for, which these data cannot establish.\n"
: " => Inconclusive. A Patterson peak on its own is not enough to call a pseudo-symmetry\n"
" (peaks of this size occur without one), and this dataset cannot supply the second\n"
" half of the evidence.\n");
return os.str();
}
if (!result.detected) {
os << " => No translational pseudo-symmetry indicated.\n";
if (result.suspected)
os << " => POSSIBLE translational pseudo-symmetry: a peak of the size molecular replacement\n"
<< " corrects for, which these data cannot establish - too few reflections in the band\n"
<< " for the noise floor to be far below it, and the modulation not significant.\n";
else if (result.modulates)
os << " => The vector modulates the intensities, but the peak is below the fifth of the\n"
<< " origin at which molecular replacement corrects for a pseudo-translation: two\n"
<< " copies of a small part of the cell, or two rotated copies alike at low\n"
<< " resolution. Noted; nothing is done with it.\n";
else
os << " => No translational pseudo-symmetry indicated.\n";
return os.str();
}
if (result.detected_on_averaged_classes)
image_analysis/scale_merge/TranslationalNCS.h
+54 -14
View File
@@ -19,28 +19,36 @@
// reliably breaks molecular replacement, because the modulation it puts on the intensities is not in
// the search model, and it is invisible in every statistic this program reported until now.
//
// TWO numbers, and BOTH must fire. Measured over a corpus of merged datasets, either one alone
// over-calls by roughly a factor of two.
// TWO measurements, and a size.
//
// - PATTERSON PEAK. The largest off-origin peak of the native Patterson computed from the merged
// intensities (shell-normalised E^2 coefficients over 20-5 A, one inverse FFT), as a percentage
// of the origin peak, counting only peaks farther than 15 A from ANY origin-equivalent lattice
// point so the origin's own shoulder and the centring vectors cannot be mistaken for a
// pseudo-translation. Answers: is there a peak above this dataset's own noise floor?
// pseudo-translation. Answers: where to look, and how much of the scattering the translation
// relates.
//
// - MODULATION DEPTH. With that peak vector refined against the intensities, the acentric
// reflections binned by the phase frac(h.u): a real pseudo-translation makes <E^2> follow a
// cosine with a deep minimum near phase 1/2, an artefact leaves it flat. Reported as the ratio
// of the strongest to the weakest bin mean. Answers: does that vector actually modulate the
// intensities - i.e. is it a translation between real copies, or a feature of the map?
// - MODULATION. With that peak vector refined against the intensities, the correlation between E^2
// and cos(2 pi h.u) over the acentric reflections: a real pseudo-translation makes <E^2> follow
// a cosine of the phase h.u, an artefact leaves it flat. Scored against the same search on
// intensities shuffled within their resolution shells. Answers: does that vector actually
// modulate the intensities - i.e. is it a translation between real copies, or a feature of the
// map? Also read as the ratio of the strongest to the weakest phase class, for the report.
//
// - SIZE. A vector that modulates is a fact about the crystal; the PSEUDO-TRANSLATION reported is
// the one molecular replacement has to be told about, which the programs correct for from a
// Patterson peak of a fifth of the origin (in the band, or as the modulation's depth over all the data).
// Below that the data do not tell a weak pseudo-translation from two rotated copies that look
// alike at 20-5 A, and nothing downstream needs it. A peak of that size the data are too few
// to establish is carried as POSSIBLE rather than dropped.
//
// NULLS ARE PER DATASET, not a table. The Patterson peak is scored against the same map recomputed
// with the intensities permuted WITHIN resolution shells (same count, same completeness, same
// resolution range, same intensity multiset; the inter-reflection correlation that makes a Patterson
// peak destroyed), which is the only way to compare datasets whose noise floor runs from ~1% to ~18%
// depending on how many reflections they have. The modulation is scored against the same refinement
// started from random vectors, so the contrast the search itself can manufacture is measured rather
// than assumed.
// on the same shuffle, so the correlation the search itself can manufacture is measured rather than
// assumed.
//
// REPORTING ONLY. Nothing here removes a reflection, changes a merge, or gates a candidate. Its one
// effect on another verdict is to hand the twinning analysis the vector it needs in order not to
@@ -88,8 +96,19 @@ struct TranslationalNCSResult {
int n_modulation_reflections = 0;
std::array<double, 3> vector_frac{{0.0, 0.0, 0.0}}; // peak vector refined against the intensities
double vector_length_A = 0.0; // shortest lattice-equivalent length of that vector
double modulation = 0.0; // max/min of <E^2> over the phase bins
double modulation_null = 0.0; // same search from random vectors (median)
// The statistic the vector is refined on and judged by: the correlation between E^2 and
// cos(2 pi h.u), against the same search on intensities shuffled within their resolution shells
// (median and spread over the draws).
double correlation = 0.0;
double correlation_null = 0.0;
double correlation_null_sd = 0.0;
// The profile's depth, for the report: max/min of <E^2> over the phase bins at the vector, and
// what the null searches reach at theirs.
double modulation = 0.0;
double modulation_null = 0.0;
double modulation_depth_percent = 0.0; // the modulation's depth over the whole modulation
// population as a Patterson height at the vector, % of
// the origin - undiluted by a twin
std::vector<double> phase_bin_mean_e2; // the profile itself, for the report
// The same test with the classes averaged over the rotations of the LATTICE (the merge's own
@@ -104,9 +123,12 @@ struct TranslationalNCSResult {
// exactly as the vector's is.
bool averaged_measured = false;
int averaged_rotations = 0; // lattice rotations averaged over
double correlation_averaged = 0.0;
double correlation_averaged_null = 0.0;
double correlation_averaged_null_sd = 0.0;
double modulation_averaged = 0.0; // max/min of <E^2> over the averaged classes
double modulation_averaged_null = 0.0;
bool detected_on_averaged_classes = false; // `detected` stands on the averaged test alone
bool detected_on_averaged_classes = false; // `modulates` stands on the averaged test alone
// The pseudo-translation's images R u over the rotations its classes are defined by: u alone
// when the call stands on the vector, every lattice rotation's image when it stands on the
@@ -115,7 +137,18 @@ struct TranslationalNCSResult {
std::vector<std::array<double, 3>> class_vectors;
// ---- verdict -------------------------------------------------------------------------------
bool detected = false; // both criteria fire
// The vector modulates the intensities: the contrast along the refined vector is significant
// against its permutation null. A fact about the crystal, whatever its size.
bool modulates = false;
// ... and its size - the band's Patterson peak, or the modulation's depth over all the data - is
// at least a fifth of the origin, the level at which molecular replacement corrects for a
// pseudo-translation (kDeclaredPeakPct in the .cpp). This is the reported pathology; a weaker
// modulation is noted and nothing is done with it.
bool detected = false;
// Not `modulates`, but a band peak of the declared size standing clear of the permutation floor:
// the data are too few to establish it (a small cell has ~150 reflections in the band), and the
// hypothesis is reported as POSSIBLE rather than dropped.
bool suspected = false;
// The refined vector is a rational translation p/q of the reported cell: the crystal is
// pseudo-centred. This is NOT by itself a statement that the cell is wrong - the suppressed class
@@ -155,3 +188,10 @@ TranslationalNCSResult AnalyzeTranslationalNCS(const std::vector<MergedReflectio
size_t nthreads = 1);
std::string TranslationalNCSToText(const TranslationalNCSResult &result);
// Carry a result across a change of basis `cob` (acting on fractional coordinates, x' = cob x, as
// CrystalSetting.h states one): every vector is rotated onto the new axes. The analysis is made on
// the axes the group was determined on, and the files can be written on others - a model's, a
// reference's, the standard setting - chosen after it; a vector left on the old axes names the
// wrong classes there.
void ChangeBasis(TranslationalNCSResult &result, const gemmi::Op &cob);
rugnux/ResultReport.cpp
+59 -35
View File
@@ -93,10 +93,6 @@ namespace {
// wording names a very strong pseudo-translation as the other reading, and refinement decides.
constexpr double LATTICE_TRANSLATION_WORDING_PCT = 90.0;
// Every detection is warned about; below this Patterson peak the wording says the modulation is small
// and may not matter to molecular replacement (phenix.xtriage flags peaks above 20% of the origin).
// On the rc173 battery the detected peaks are 8.1% on one set and 24.8-63.2% on the other 21.
constexpr double PSEUDO_TRANSLATION_WORDING_PCT = 20.0;
// The physical range of <|L|>: 0.5 untwinned, 0.375 a perfect twin, whatever the twin law. A
// reading outside it is not a twin fraction - tNCS, anisotropy, overlapping reflections of a long
@@ -2367,7 +2363,7 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
const auto &t = result.tncs;
if (!t.measurable) {
Add(s, KeyEnum("TNCS_DETECTED", "NOT_MEASURED",
{"TRUE", "FALSE", "INCONCLUSIVE", "NOT_MEASURED"}));
{"TRUE", "POSSIBLE", "FALSE", "INCONCLUSIVE", "NOT_MEASURED"}));
Add(s, Blank());
Add(s, Prose(fmt::format(
" Translational pseudo-symmetry: NOT MEASURED{}{}. A refusal is not a\n"
@@ -2376,9 +2372,10 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
} else {
Add(s, KeyEnum("TNCS_DETECTED",
t.detected ? "TRUE"
: t.suspected ? "POSSIBLE"
: t.modulation_measured ? "FALSE"
: "INCONCLUSIVE",
{"TRUE", "FALSE", "INCONCLUSIVE", "NOT_MEASURED"}));
{"TRUE", "POSSIBLE", "FALSE", "INCONCLUSIVE", "NOT_MEASURED"}));
// A different finding from a pseudo-symmetry, and reported as one: the data are
// EXACTLY invariant under these, so the centring or the cell is wrong rather than the
// molecules being packed unusually.
@@ -2405,12 +2402,21 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
fmt::format("{:.3f} {:.3f} {:.3f}", t.vector_frac[0],
t.vector_frac[1], t.vector_frac[2])));
Add(s, KeyReal("TNCS_VECTOR_LENGTH", t.vector_length_A, "{:.1f}"));
Add(s, KeyReal("TNCS_MODULATION_DEPTH_PCT", t.modulation_depth_percent, "{:.1f}"));
Add(s, KeyReal("TNCS_CORRELATION", t.correlation, "{:.3f}"));
Add(s, KeyText("TNCS_CORRELATION_NULL",
fmt::format("{:.3f} +- {:.3f}", t.correlation_null, t.correlation_null_sd)));
Add(s, KeyReal("TNCS_MODULATION", t.modulation, "{:.2f}"));
Add(s, KeyReal("TNCS_MODULATION_NULL", t.modulation_null, "{:.2f}"));
Add(s, KeyBool("TNCS_MODULATES", t.modulates));
}
// The second test, on the classes averaged over the lattice's rotations, run only where
// the vector's own test could not call a significant peak.
if (t.averaged_measured) {
Add(s, KeyReal("TNCS_CORRELATION_AVERAGED", t.correlation_averaged, "{:.3f}"));
Add(s, KeyText("TNCS_CORRELATION_AVERAGED_NULL",
fmt::format("{:.3f} +- {:.3f}", t.correlation_averaged_null,
t.correlation_averaged_null_sd)));
Add(s, KeyReal("TNCS_MODULATION_AVERAGED", t.modulation_averaged, "{:.2f}"));
Add(s, KeyReal("TNCS_MODULATION_AVERAGED_NULL", t.modulation_averaged_null, "{:.2f}"));
}
@@ -2432,17 +2438,35 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
if (t.detected)
Add(s, Prose(fmt::format(
" Translational pseudo-symmetry: INDICATED. The native Patterson has an off-origin peak\n"
" at {:.1f}% of the origin at ({:.3f}, {:.3f}, {:.3f}), {:.1f} A, and the intensities\n"
" modulate {:.1f}x along that vector{}. The molecular-replacement program needs the vector\n"
" declared: the modulation is not in the search model, and without it the search can\n"
" fail on data that are otherwise good.",
" at {:.1f}% of the origin at ({:.3f}, {:.3f}, {:.3f}), {:.1f} A (a modulation {:.1f}% deep\n"
" over all the data), and the intensities modulate {:.1f}x along it{}. The\n"
" molecular-replacement program needs the vector declared: the modulation is not in the\n"
" search model, and without it the search can fail on data that are otherwise good.",
t.peak_percent, t.vector_frac[0], t.vector_frac[1], t.vector_frac[2],
t.vector_length_A,
t.vector_length_A, t.modulation_depth_percent,
t.detected_on_averaged_classes ? t.modulation_averaged : t.modulation,
t.detected_on_averaged_classes
? " and its images under the lattice's rotations together (a twin law or a\n"
" symmetry operator that moves the vector averages its own classes)"
: "")));
else if (t.suspected)
Add(s, Prose(fmt::format(
" Translational pseudo-symmetry: POSSIBLE. The native Patterson has an off-origin peak at\n"
" {:.1f}% of the origin at ({:.3f}, {:.3f}, {:.3f}), the size molecular replacement corrects\n"
" for, but these data cannot establish it: only {} reflections lie in the 20-5 A band,\n"
" whose noise floor is {:.1f} +- {:.2f}% (z = {:.1f}), and the modulation along the vector\n"
" is not significant. Declare it to molecular replacement if the search fails without it.",
t.peak_percent, t.peak_frac[0], t.peak_frac[1], t.peak_frac[2], t.n_reflections,
t.null_mean, t.null_sd, t.peak_z)));
else if (t.modulates)
Add(s, Prose(fmt::format(
" Translational pseudo-symmetry: no indication. The vector ({:.3f}, {:.3f}, {:.3f}) modulates\n"
" the intensities, but the Patterson at it is {:.1f}% of the origin in the 20-5 A band and\n"
" the modulation {:.1f}% deep over all the data - below the fifth at which molecular\n"
" replacement corrects for a pseudo-translation, and the size two rotated copies alike at\n"
" low resolution give.",
t.vector_frac[0], t.vector_frac[1], t.vector_frac[2], t.peak_percent,
t.modulation_depth_percent)));
else if (t.modulation_measured)
Add(s, Prose(fmt::format(
" Translational pseudo-symmetry: no indication. The largest off-origin Patterson peak is\n"
@@ -2483,32 +2507,30 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
" The pseudo-symmetry numbers above were measured underneath it.",
v[0], v[1], v[2], undeclared_pct)));
}
if (t.detected && t.peak_percent < PSEUDO_TRANSLATION_WORDING_PCT)
Add(s, ProseDefaultOnly(fmt::format(
"\n The peak is significant but weak: below the {:.0f}% of the origin at which a\n"
" pseudo-translation is conventionally declared to molecular replacement.",
PSEUDO_TRANSLATION_WORDING_PCT)));
Add(s, Prose("\n" + TranslationalNCSToText(t), true));
if (t.modulation_measured && !t.detected)
if (t.modulation_measured && !t.detected && !t.suspected)
Add(s, Prose(
"\n TNCS_DETECTED= FALSE means the two tests disagreed with each other or neither fired, not\n"
" that a Patterson peak was absent: a peak is called only when it stands above THIS\n"
" dataset's own permutation noise floor (TNCS_PATTERSON_NULL_PCT, which runs from ~1% on a\n"
" large merge to ~18% on a small one) AND the vector it names is shown to modulate the\n"
" intensities. Either test alone over-calls by about a factor of two on a corpus of merged\n"
" datasets.", true));
"\n TNCS_DETECTED= FALSE means no pseudo-translation of the size molecular replacement\n"
" corrects for: a vector is called only when it is shown to modulate the intensities\n"
" (TNCS_MODULATES, the correlation of E^2 with cos 2 pi h.u against the same search on\n"
" shuffled intensities) AND its size - the Patterson peak in the 20-5 A band, or the\n"
" modulation's depth over all the data - reaches a fifth of the origin.", true));
if (t.detected)
Warn(doc, PathologyCode::PSEUDO_TRANSLATION, fmt::format(
"Translational pseudo-symmetry: Patterson off-origin peak {:.1f}% of the origin at "
"({:.3f}, {:.3f}, {:.3f}), {:.1f} A - {}",
t.peak_percent, t.vector_frac[0], t.vector_frac[1], t.vector_frac[2],
t.vector_length_A,
t.peak_percent >= PSEUDO_TRANSLATION_WORDING_PCT
? "the molecular-replacement program needs the vector declared, or it will be "
"searching against a modulation it does not model"
: "significant but weak; check whether molecular replacement needs it declared"));
"({:.3f}, {:.3f}, {:.3f}), {:.1f} A - the molecular-replacement program needs the vector "
"declared, or it will be searching against a modulation it does not model",
std::max(t.peak_percent, t.modulation_depth_percent), t.vector_frac[0], t.vector_frac[1],
t.vector_frac[2], t.vector_length_A));
else if (t.suspected)
Warn(doc, PathologyCode::PSEUDO_TRANSLATION, fmt::format(
"Possible translational pseudo-symmetry: Patterson off-origin peak {:.1f}% of the origin "
"at ({:.3f}, {:.3f}, {:.3f}), which {} reflections in the 20-5 A band cannot establish "
"(noise floor {:.1f} +- {:.2f}%) - declare it to molecular replacement if the search "
"fails without it",
t.peak_percent, t.peak_frac[0], t.peak_frac[1], t.peak_frac[2], t.n_reflections,
t.null_mean, t.null_sd));
if (undeclared)
Warn(doc, PathologyCode::LATTICE_TRANSLATION, fmt::format(
"Check the centring and the cell: the Patterson at ({:.3f}, {:.3f}, {:.3f}) is {:.0f}% of "
@@ -3731,10 +3753,12 @@ ReportDocument BuildReportDocument(const std::string &output_prefix,
}
if (t.detected)
ps += (ps.empty() ? "" : "; ") + fmt::format(
"{} (Patterson peak {:.1f}% of origin, {:.1f} A)",
t.peak_percent >= PSEUDO_TRANSLATION_WORDING_PCT
? "INDICATED" : "INDICATED, weak - check whether molecular replacement needs it",
t.peak_percent, t.vector_length_A);
"INDICATED (Patterson peak {:.1f}% of origin, {:.1f} A)",
std::max(t.peak_percent, t.modulation_depth_percent), t.vector_length_A);
else if (t.suspected)
ps += (ps.empty() ? "" : "; ") + fmt::format(
"POSSIBLE - Patterson peak {:.1f}% of origin that {} reflections cannot establish",
t.peak_percent, t.n_reflections);
checked_row("Pseudo-symmetry", {PathologyCode::PSEUDO_TRANSLATION, PathologyCode::LATTICE_TRANSLATION},
ps.empty() ? std::string("no indication") : ps);
} else {
rugnux/RugnuxScaleMerge.cpp
+3
View File
@@ -4036,6 +4036,9 @@ bool Rugnux::ScaleMergeAndSymmetry(PipelineLocals &p) {
sm.merged = std::move(*relabelled);
reseat_integrated(cob);
relabel_group(cob, *moved);
// The pseudo-translation was measured on the axes the group was determined on; its
// vectors go onto the written ones with everything else.
ChangeBasis(result.tncs, cob);
// The free set is a function of the written index (RfreeFlags.h), so that every dataset of
// a crystal form written on the same axes carries the same one: drawn again on these axes.
if (relabelled == nullptr)
tests/ResultReportTest.cpp
+20 -5
View File
@@ -1062,15 +1062,30 @@ TEST_CASE("ResultReport_TranslationalNCS", "[Diagnostics]") {
CHECK(strong.find("\nPATHOLOGY_FLAGS= LATTICE_TRANSLATION\n") != std::string::npos);
CHECK(strong.find("or a very strong pseudo-translation") != std::string::npos);
// Significant but weak: warned about, and worded as a check.
// A vector that modulates the intensities below the size molecular replacement corrects for:
// noted, not flagged.
result.tncs.undeclared_lattice_translations.clear();
result.tncs.undeclared_lattice_translation_pct.clear();
result.tncs.detected = true;
result.tncs.detected = false;
result.tncs.modulates = true;
result.tncs.peak_percent = 8.1;
result.tncs.modulation_depth_percent = 6.0;
const auto weak = RenderResultReport("prefix", "in.h5", x, result);
CHECK(weak.find("\nTNCS_DETECTED= TRUE\n") != std::string::npos);
CHECK(weak.find("\nPATHOLOGY_FLAGS= PSEUDO_TRANSLATION\n") != std::string::npos);
CHECK(weak.find("significant but weak; check whether molecular replacement needs it declared") != std::string::npos);
CHECK(weak.find("\nTNCS_DETECTED= FALSE\n") != std::string::npos);
CHECK(weak.find("\nTNCS_MODULATES= TRUE\n") != std::string::npos);
CHECK(weak.find("PSEUDO_TRANSLATION") == std::string::npos);
CHECK(weak.find("below the fifth at which molecular") != std::string::npos);
// A peak of that size the data cannot establish: carried as POSSIBLE, and warned about as one.
result.tncs.modulates = false;
result.tncs.suspected = true;
result.tncs.peak_percent = 34.6;
result.tncs.peak_z = 3.4;
result.tncs.n_reflections = 238;
const auto possible = RenderResultReport("prefix", "in.h5", x, result);
CHECK(possible.find("\nTNCS_DETECTED= POSSIBLE\n") != std::string::npos);
CHECK(possible.find("\nPATHOLOGY_FLAGS= PSEUDO_TRANSLATION\n") != std::string::npos);
CHECK(possible.find("238 reflections in the 20-5 A band cannot establish") != std::string::npos);
}
TEST_CASE("ResultReport_GeometryNotConverged", "[Diagnostics]") {
tests/TranslationalNCSTest.cpp
+31 -5
View File
@@ -83,7 +83,9 @@ TEST_CASE("A pseudo-translation is detected and its vector recovered", "[tncs]")
CHECK(r.peak_percent > 20.0);
CHECK(r.undeclared_lattice_translations.empty());
CHECK(r.peak_z > 5.0);
CHECK(r.modulation > 2.5 * r.modulation_null);
CHECK(r.modulates);
CHECK(r.correlation > r.correlation_null + 5.0 * r.correlation_null_sd);
CHECK(r.modulation_depth_percent > 20.0);
// The vector is refined against the intensities of the whole dataset, so it is good to about a
// thousandth of a cell edge - which is what the phase classes below need, and a good deal better
// than a molecular-replacement program needs as a starting point.
@@ -107,6 +109,8 @@ TEST_CASE("A crystal without a pseudo-translation stays silent", "[tncs]") {
REQUIRE(r.measurable);
CHECK_FALSE(r.detected);
CHECK_FALSE(r.suspected);
CHECK_FALSE(r.modulates);
}
// The case that would make the detector useless if it were got wrong. A centred lattice has a real,
@@ -125,10 +129,11 @@ TEST_CASE("A genuinely centred lattice is not called", "[tncs]") {
REQUIRE(sg != nullptr);
const auto merged = Synthetic(c.cell, c.sg, nullptr, 0.0);
const auto r = AnalyzeTranslationalNCS(merged, c.cell, sg);
INFO(c.sg << ": peak " << r.peak_percent << "% z " << r.peak_z << " modulation "
<< r.modulation << " vs null " << r.modulation_null);
INFO(c.sg << ": peak " << r.peak_percent << "% z " << r.peak_z << " correlation "
<< r.correlation << " vs null " << r.correlation_null << " +- " << r.correlation_null_sd);
REQUIRE(r.measurable);
CHECK_FALSE(r.detected);
CHECK_FALSE(r.suspected);
}
}
@@ -159,6 +164,27 @@ TEST_CASE("A small-molecule cell reports that it could not be measured", "[tncs]
CHECK(TranslationalNCSToText(r).find("not a statement") != std::string::npos);
}
// A small cell has ~150 reflections in the band, where the noise floor is a fifth of the origin
// itself; with a pseudo-translation that modulates only the lowest resolution (two copies alike as
// blobs, not as atoms) nothing reaches significance, and the peak is carried as POSSIBLE.
TEST_CASE("A small cell's peak the data cannot establish is carried as POSSIBLE", "[tncs]") {
const gemmi::UnitCell cell = Cell(31.0, 31.0, 81.0);
const gemmi::SpaceGroup *sg = gemmi::find_spacegroup_by_name("P 42 2 2");
const std::array<double, 3> u{0.5, 0.5, 0.5};
auto merged = Synthetic(cell, "P 42 2 2", nullptr, 0.0, 1.7);
for (auto &m : merged)
if (m.d >= 5.0)
m.I = static_cast<float>(m.I * (1.0 + 0.5 * std::cos(2.0 * std::numbers::pi
* (m.h * u[0] + m.k * u[1] + m.l * u[2]))));
const auto r = AnalyzeTranslationalNCS(merged, cell, sg);
INFO("peak " << r.peak_percent << "% z " << r.peak_z << " n " << r.n_reflections << " correlation "
<< r.correlation << " vs " << r.correlation_null << " +- " << r.correlation_null_sd);
REQUIRE(r.measurable);
CHECK(r.peak_percent > 20.0);
CHECK(r.detected != r.suspected); // one of the two, whichever the modulation reaches
CHECK(TranslationalNCSToText(r).find("pseudo-symmetry") != std::string::npos);
}
// P1 - no symmetry to expand with, and no centric reflections. The commonest case for a merge the
// space-group search declined to promote, so it must work rather than divide by zero.
TEST_CASE("P1 data are handled with no space group at all", "[tncs]") {
@@ -294,8 +320,8 @@ TEST_CASE("A pseudo-translation under a perfect twin is called on the lattice-av
for (auto &m : merged)
m.I = static_cast<float>(m.I * (1.0 + kRealisticDepth * TncsClass(m.h, m.k, m.l, images)));
const auto r = AnalyzeTranslationalNCS(merged, cell, sg);
INFO("peak " << r.peak_percent << "% z " << r.peak_z << " vector " << r.modulation << "/"
<< r.modulation_null << " averaged " << r.modulation_averaged << "/" << r.modulation_averaged_null);
INFO("peak " << r.peak_percent << "% z " << r.peak_z << " vector " << r.correlation << "/"
<< r.correlation_null << " averaged " << r.correlation_averaged << "/" << r.correlation_averaged_null);
REQUIRE(r.detected);
CHECK(r.detected_on_averaged_classes);
CHECK(r.averaged_rotations == 12);