From 40874c2c2dc4da9f15aa89d30fdc60745087fef4 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Sun, 20 Sep 2026 08:51:02 +0200 Subject: [PATCH] Space-group search: the twin-immune zone is read on a normalisation the control certifies The zone verdict that arbitrates a refused promotion read its zone absolutely, against the centric and acentric Wilson expectations, and an absolute reading is only as good as the normalisation under it. On two refused 622 promotions - a 6/m crystal and a 312 one, whose added operators disagreed at 4.7x and 5.2x the parent's H and merged to an R_meas of 0.34 under the higher group - every class read centric: a 67 A^2 anisotropy, isotropically normalised in bins of 100, spread one shell's expected intensity over a factor of 15 between its directions, and the acentric control read 1.00 against its 0.74. The zone read the same as the control, +0.10 nats per reflection, and +124 nats over 947 reflections rescued a twin law. Two things change. The anisotropy is fitted on the acentric reflections of the shells read (ln I = c + s^T Q s, by least squares) and its deviatoric part taken out of every intensity before anything is normalised; that alone brings the controls of the crystals measured to 0.73-0.77. What no normalisation removes, the control then certifies: an acentric population reads -0.130 nats per reflection when the normalisation is right, a twinned one reads below that, so whatever the control reads above it is the normalisation's - anisotropy, a pseudo-translation, a pseudo-centring, noise all inflate every class towards centric alike - and the zone, normalised the same way, carries the same per reflection; the calibrated evidence has it taken off, and that is what the verdict reads. The two rescued twin laws now read -130 and -78 nats and stay refused, with the law named; the genuine promotions measured read +137 to +580 (a pseudo-centred orthorhombic crystal whose control reads 0.99 still +334); the trigonal and hexagonal partial twins -129 to -420 as before. The report prints the control's excess and the calibrated evidence beside the raw one. The metric-lattice re-ask carried a second copy of the two-arm rule without the zone: a P31 partial twin whose 32 the main rule had refused on the zone was promoted by that ask's Lorentz-filtered arm. It now puts the same verdict to a refusal the other arm would outvote. Co-Authored-By: Claude Fable 5.1 Claude-Session: https://claude.ai/code/session_013nW6FNRP1bBJJ8pfHiByAT --- docs/CPU_DATA_ANALYSIS_DECISIONS.md | 4 +- .../scale_merge/TwinningAnalysis.cpp | 159 ++++++++++++++---- image_analysis/scale_merge/TwinningAnalysis.h | 26 ++- rugnux/ResultReport.cpp | 3 +- rugnux/Rugnux.cpp | 70 ++++++-- tests/TwinningAnalysisTest.cpp | 70 +++++++- 6 files changed, 267 insertions(+), 65 deletions(-) diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index 8718102f0..ecc059048 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -46,7 +46,7 @@ Both search merges also drop the **frames whose fitted per-frame scale came out **Centering** is accepted when the systematically-absent class is weak relative to the present one by *either* of two floor-independent tests: its mean signed $I/\sigma$ well below the present mean, *or* its rate of individually-significant reflections well below the present class's own significant rate. The second test covers weak and low-energy data, where a positive intensity floor (background and profile leakage) lifts the absent class's mean $I/\sigma$ well above zero and, when the present class is itself weak, carries the plain mean ratio past its bound; a false centering fails both tests, its absent class being as strong as the present one. When several centerings pass they are ranked by a **Beta-tail likelihood** rather than by a count of net absences: a count lets a centering that extinguishes many more reflections out-rank one that violates far fewer, and the likelihood weighs the violations against the class each belongs to instead. The acceptance bounds above are unchanged by that ranking. -**A confirmed promotion that is refused is decided by merging under it.** Every gate that refuses a *confirmed* higher point group is a **ratio** — the added operators' disagreement against the parent's, or their $R$ against the best-agreeing operator anywhere — and both references are properties of how the crystal was mounted rather than of its symmetry. A 2-fold within a degree of the spindle records its mates on the same detector pixel half a turn later, so it carries no geometry-dependent systematic at all, and by being that clean it makes every other operator look bad against it; the best-agreeing operator can also be one of the operators under test, which collapses the ratio on worse data and raises it on better. So where a promotion is confirmed and then refused, the **merge** is asked instead of the ratio. Both groups are scaled and merged over one pinned resolution range — the adopted group's own cut — and the promotion is taken only if **neither $R_\mathrm{meas}$ nor ISa gets worse**. $R_\mathrm{meas}$ is multiplicity-corrected, so folding non-equivalent reflections together has to inflate it, and the error model is refitted per merge, so ISa says whether the extra multiplicity was bought with systematic disagreement. $CC_{1/2}$ and $\langle I/\sigma\rangle$ cannot arbitrate this: both *rise* on a false promotion too. Two details the comparison turns on: the refused point group's representative is primitive and symmorphic, so the arm merged under it takes the **adopted group's centring** — merged in the bare representative, a centred crystal is handed its centring-absent class as data, which inverts the decision — and the absence stage is asked of the adopted merge, so both arms carry the same absence classes and differ by the added rotations alone. The cost is two extra merges, and only on a run that records a refusal, which is a few in a hundred. +**A confirmed promotion that is refused is decided by the twin-immune zone of its added operators, and only where that is undecided by merging under it.** The zone (§13.2) is the one reading that does not depend on the per-frame scales or on a reference operator: read on the all-observation P1 merge where the two search arms disagree, and on the adopted group's own merge at the remerge, an index-2 promotion whose calibrated zone evidence is acentric by 20 nats stays refused — on all observations, and the first added operator is recorded as the twin law — and one whose zone is centric by 20 nats is taken. Where the zone cannot decide (a higher index, no zone reflections in the shells with signal, within the margin), the merge is asked. Every gate that refuses a *confirmed* higher point group is a **ratio** — the added operators' disagreement against the parent's, or their $R$ against the best-agreeing operator anywhere — and both references are properties of how the crystal was mounted rather than of its symmetry. A 2-fold within a degree of the spindle records its mates on the same detector pixel half a turn later, so it carries no geometry-dependent systematic at all, and by being that clean it makes every other operator look bad against it; the best-agreeing operator can also be one of the operators under test, which collapses the ratio on worse data and raises it on better. So where a promotion is confirmed and then refused, the **merge** is asked instead of the ratio. Both groups are scaled and merged over one pinned resolution range — the adopted group's own cut — and the promotion is taken only if **neither $R_\mathrm{meas}$ nor ISa gets worse**. $R_\mathrm{meas}$ is multiplicity-corrected, so folding non-equivalent reflections together has to inflate it, and the error model is refitted per merge, so ISa says whether the extra multiplicity was bought with systematic disagreement. $CC_{1/2}$ and $\langle I/\sigma\rangle$ cannot arbitrate this: both *rise* on a false promotion too. Two details the comparison turns on: the refused point group's representative is primitive and symmorphic, so the arm merged under it takes the **adopted group's centring** — merged in the bare representative, a centred crystal is handed its centring-absent class as data, which inverts the decision — and the absence stage is asked of the adopted merge, so both arms carry the same absence classes and differ by the added rotations alone. The cost is two extra merges, and only on a run that records a refusal, which is a few in a hundred. **The lattice class is re-asked where the metric knows more.** The Bravais class comes from Niggli reduction and a lattice-character lookup (§6), and that lookup can land short of the truth or beside it: near the Niggli type-I/type-II boundary it is decided by the last digits of the refined cell, and a lattice that is nearly but not exactly hexagonal matches the hexagonal character even when no point group that class can host carries the two-folds the data actually have. Two recoveries run, both settled by the intensities. Every rotation the cell metric can host beyond the named class — read off Le Page's two-fold search on the lattice itself, which measures each rotation's obliquity in a primitive basis and owes nothing to the character table — is put to the intensities as a single operator, scored exactly as the search scores its own operators, on the same reflection population and the same $E^2$ normalisation. And where the metric group is larger than the adopted class's holohedry, the merge is reindexed into the metric group's conventional cell and the space-group search is run again there, with every gate live; the reindex is committed only where the search in the new setting confirms a strictly higher point group *and* the centring the new cell describes, so a pseudo-symmetric metric leaves the answer already in hand standing. @@ -54,7 +54,7 @@ Both search merges also drop the **frames whose fitted per-frame scale came out A Padilla–Yeates $L$-test ($\langle|L|\rangle$, $\langle L^2\rangle$ — 0.500 and 0.333 untwinned, 0.375 and 0.200 for a perfect twin) and the second moment $\langle I^2\rangle/\langle I\rangle^2$ (2.0 for untwinned acentric data, 1.5 for a perfect twin) are written to the merged mmCIF as a twinning diagnostic. Both are taken on intensities divided by their resolution-shell mean, over the shells whose $\langle I/\sigma\rangle$ reaches 1, with Wilson outliers rejected; the $L$-test pairs are therefore compared on one scale even where two index steps span a steep fall-off, as in phenix.xtriage and ctruncate, and the selection is by shell, never by the individual reflection's $I/\sigma$, which would cut the weak tail and bias $\langle|L|\rangle$ down. The twin fraction is quoted from the statistic that carries the verdict — the $L$-test unless the call rests on the second moment alone — and a second moment is not turned into a fraction under a detected pseudo-translation, which inflates it. A merohedral twin law exists only where the Laue class is a proper subgroup of the lattice holohedry, so in the high-symmetry holohedral classes ($4/mmm$, $6/mmm$, $m\bar{3}m$, and $\bar{3}m$ on a rhombohedral lattice) no twin is called. The $L$-test is still read there, for a different question: merging $I(h)$ with $I(Th)$ under a false operator $T$ gives $(I(h)+I(Th))/2$ whatever the twin fraction, which has the perfect-twin distribution, so $\langle|L|\rangle$ below 0.42 in a holohedral class is reported as **an adopted operator averaging unequal intensities** — the space group is too high, or a twin law was absorbed into the point group — a warning, never a change to the space group. A genuine operator leaves the untwinned 0.5. Reflections that overlap along a very long axis narrow the distribution the same way, which is one reason it stays a warning. The same numbers measured on the P1 merge the space-group search was given, before any point group was adopted, are reported beside them (`_BEFORE_SEARCH`), with that merge's own pseudo-translation declared and the reflections its lattice centring extinguishes left out. The low-symmetry holohedral classes ($\bar{1}$, $2/m$, $mmm$) also admit no strictly merohedral law, but they stay eligible for the twin call on purpose: *pseudo*-merohedral twinning through an accidentally special metric cannot be ruled out from the symmetry alone, and those are the classes it happens in. -Beside them, on rotation data, the run reports **twin-immune zone evidence** for the operators the adopted point group adds over each of its index-2 subgroups, read on the P1 cross-check merge. 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. It is **report-only** — nothing in the space-group decision reads it. 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. +Beside them, on rotation data, the run reports **twin-immune zone evidence** for the operators the adopted point group adds over each of its index-2 subgroups, read on the P1 cross-check merge. 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). 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 `_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. diff --git a/image_analysis/scale_merge/TwinningAnalysis.cpp b/image_analysis/scale_merge/TwinningAnalysis.cpp index 0ed2afb25..2900c0a79 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.cpp +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -14,6 +14,8 @@ #include #include +#include + namespace { int64_t PackHKL(int h, int k, int l) { constexpr int64_t bias = 1 << 20; // indices assumed within +/- 2^20 @@ -488,6 +490,11 @@ namespace { // A reflection this far above its bin mean is not a Wilson draw (centric P ~ 1e-5): one such // reflection read 2.33 for a whole zone. constexpr double ZONE_OUTLIER_E2 = 20.0; + // What CentricOverAcentric averages over an error-free acentric population: with x = E^2 ~ Exp(1), + // <-ln(2 pi x)/2 + x/2> = (1 + gamma_E - ln 2 pi)/2 = -0.130 nats (a centric one reads +0.216). + // Errors move it towards zero, so an error-free baseline takes the larger excess off a noisy + // control - the conservative side. + constexpr double ACENTRIC_EVIDENCE_PER_REFLECTION = -0.1303; // ln p_centric(E^2) - ln p_acentric(E^2) for a measured E^2 with error s: both Wilson densities // convolved with the measurement error. Without the convolution the centric density, infinite at @@ -563,6 +570,18 @@ namespace { } return z; } + + // The control's excess over the acentric expectation is the normalisation's (see + // TwinImmuneZoneResult::control_excess_per_reflection); every zone carries it per reflection. + void Calibrate(TwinImmuneZoneResult& result) { + if (result.control.n > 0) + result.control_excess_per_reflection = std::max( + 0.0, result.control.evidence_nats / result.control.n - ACENTRIC_EVIDENCE_PER_REFLECTION); + result.control.calibrated_evidence_nats = + result.control.evidence_nats - result.control.n * result.control_excess_per_reflection; + for (auto& z : result.zones) + z.calibrated_evidence_nats = z.evidence_nats - z.n * result.control_excess_per_reflection; + } } namespace { @@ -576,11 +595,66 @@ namespace { return rotations; } + // An isotropic normalisation leaves an anisotropic crystal's E^2 a mixture of scales - at 2.5 A a + // 67 A^2 anisotropy spreads one shell's expected intensity over a factor of 15 between its + // directions - and every class then reads centric: the acentric control of such a crystal read + // 0.99 against its 0.74. Nor can the control calibrate that away, because a zone is a plane + // through reciprocal space and carries a different mix of directions from the control's sphere + // (a synthetic P4 twin with 60 A^2 along c*: control 0.77, zone +49 nats after the calibration). + // So the anisotropy is taken out first: ln I = c + s^T Q s on the acentric reflections of each + // shell that read above 2 sigma, by least squares, and I and sigma divided by exp(s^T Q_dev s), + // Q_dev the deviatoric part - the isotropic part is the shell mean, normalised per bin later. + // The B convention is I ~ exp(-B |s|^2 / 2) with |s| = 1/d, so delta B = 2 (Q_max - Q_min). + void RemoveAnisotropy(std::vector& strong, const std::vector& shell, + const gemmi::UnitCell& cell, const gemmi::GroupOps& gops, TwinImmuneZoneResult& result) { + auto s_of = [&](const MergedReflection& r) { + return cell.frac.mat.left_multiply(gemmi::Vec3(r.h, r.k, r.l)); + }; + auto monomials = [](const gemmi::Vec3& s) { + return Eigen::Matrix{1.0, s.x * s.x, s.y * s.y, s.z * s.z, + 2.0 * s.x * s.y, 2.0 * s.x * s.z, 2.0 * s.y * s.z}; + }; + std::vector fit; + std::vector shell_sum(*std::max_element(shell.begin(), shell.end()) + 1, 0.0); + std::vector shell_n(shell_sum.size(), 0); + for (size_t i = 0; i < strong.size(); ++i) { + const auto& r = strong[i]; + if (r.I > 2.0 * r.sigma && !gops.is_reflection_centric(gemmi::Op::Miller{{r.h, r.k, r.l}})) { + fit.push_back(i); + shell_sum[shell[i]] += r.I; + shell_n[shell[i]] += 1; + } + } + if (fit.size() < ZONE_BIN) // seven parameters want more than a bin's worth of points + return; + Eigen::Matrix ata = Eigen::Matrix::Zero(); + Eigen::Matrix aty = Eigen::Matrix::Zero(); + for (const size_t i : fit) { + const auto x = monomials(s_of(strong[i])); + ata += x * x.transpose(); + aty += x * std::log(strong[i].I / (shell_sum[shell[i]] / shell_n[shell[i]])); + } + const Eigen::Matrix c = ata.ldlt().solve(aty); + Eigen::Matrix3d q; + q << c[1], c[4], c[5], c[4], c[2], c[6], c[5], c[6], c[3]; + q -= q.trace() / 3.0 * Eigen::Matrix3d::Identity(); + const Eigen::SelfAdjointEigenSolver eig(q); + result.anisotropy_delta_b_A2 = 2.0 * (eig.eigenvalues().maxCoeff() - eig.eigenvalues().minCoeff()); + for (auto& r : strong) { + const gemmi::Vec3 s = s_of(r); + const Eigen::Vector3d v(s.x, s.y, s.z); + const double f = std::exp(v.dot(q * v)); + r.I = static_cast(r.I / f); + r.sigma = static_cast(r.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). Fills the resolution range read. - std::vector StrongReflections(const std::vector& p1_merged, - const gemmi::GroupOps& gops, - TwinImmuneZoneResult& result) { + // 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, + TwinImmuneZoneResult& result) { std::vector usable; for (const auto& r : p1_merged) { const gemmi::Op::Miller h{{r.h, r.k, r.l}}; @@ -595,7 +669,7 @@ namespace { min_s = std::min(min_s, 1.0 / (r->d * r->d)); max_s = std::max(max_s, 1.0 / (r->d * r->d)); } - std::vector strong; + std::vector strong; if (!(max_s > min_s)) return strong; auto shell_of = [&](const MergedReflection* r) { @@ -608,23 +682,40 @@ namespace { isig_sum[shell_of(r)] += r->I / r->sigma; isig_n[shell_of(r)] += 1; } + std::vector shell; for (const auto* r : usable) { const int b = shell_of(r); - if (isig_n[b] > 0 && isig_sum[b] / isig_n[b] >= ZONE_MIN_SHELL_ISIG) - strong.push_back(r); + if (isig_n[b] > 0 && isig_sum[b] / isig_n[b] >= ZONE_MIN_SHELL_ISIG) { + strong.push_back(*r); + shell.push_back(b); + } } result.d_max_A = 0.0; result.d_min_A = std::numeric_limits::infinity(); - for (const auto* r : strong) { - result.d_max_A = std::max(result.d_max_A, r->d); - result.d_min_A = std::min(result.d_min_A, r->d); + for (const auto& r : strong) { + result.d_max_A = std::max(result.d_max_A, r.d); + result.d_min_A = std::min(result.d_min_A, r.d); } + if (!strong.empty()) + RemoveAnisotropy(strong, shell, cell, gops, result); return strong; } + // The reflections acentric in the group: the control every zone is read beside. + TwinImmuneZone ReadControl(const std::vector& strong, const gemmi::GroupOps& gops, + const std::array* tncs_vector) { + 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); + z.operators = "acentric"; + return z; + } + // The zone of one subgroup hypothesis H (given by its rotations): the reflections centric in the // group but not in H - the self-mates of the coset the group adds over H. - TwinImmuneZone ReadSubgroupZone(const std::vector& strong, + TwinImmuneZone ReadSubgroupZone(const std::vector& strong, const gemmi::GroupOps& gops, const std::vector& rotations, const std::vector& h_ops, const std::array* tncs_vector) { @@ -632,14 +723,14 @@ namespace { return std::any_of(set.begin(), set.end(), [&](const gemmi::Op& o) { return o.rot == op.rot; }); }; std::vector zone; - for (const auto* r : strong) { - const gemmi::Op::Miller h{{r->h, r->k, r->l}}; - const gemmi::Op::Miller minus{{-gemmi::Op::DEN * r->h, -gemmi::Op::DEN * r->k, -gemmi::Op::DEN * r->l}}; + for (const auto& r : strong) { + const gemmi::Op::Miller h{{r.h, r.k, r.l}}; + const gemmi::Op::Miller minus{{-gemmi::Op::DEN * r.h, -gemmi::Op::DEN * r.k, -gemmi::Op::DEN * r.l}}; const bool centric_in_h = std::any_of(h_ops.begin(), h_ops.end(), [&](const gemmi::Op& op) { return op.apply_to_hkl_without_division(h) == minus; }); if (gops.is_reflection_centric(h) && !centric_in_h) - zone.push_back(r); + zone.push_back(&r); } TwinImmuneZone z = ReadZone(zone, tncs_vector); for (const auto& op : rotations) @@ -650,28 +741,25 @@ namespace { } TwinImmuneZoneResult AnalyzeTwinImmuneZone(const std::vector& p1_merged, - const gemmi::SpaceGroup& group, const gemmi::SpaceGroup& subgroup, + const gemmi::UnitCell& cell, const gemmi::SpaceGroup& group, + const gemmi::SpaceGroup& subgroup, const std::array* tncs_vector) { TwinImmuneZoneResult result; result.tncs_normalised = tncs_vector != nullptr; 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, gops, result); + const auto strong = StrongReflections(p1_merged, cell, gops, result); if (strong.empty()) return result; - 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); - result.control = ReadZone(control, tncs_vector); - result.control.operators = "acentric"; + result.control = ReadControl(strong, gops, tncs_vector); result.zones.push_back(ReadSubgroupZone(strong, gops, rotations, h_ops, tncs_vector)); + Calibrate(result); return result; } TwinImmuneZoneResult AnalyzeTwinImmuneZones(const std::vector& p1_merged, - const gemmi::SpaceGroup& group, + const gemmi::UnitCell& cell, const gemmi::SpaceGroup& group, const std::array* tncs_vector) { TwinImmuneZoneResult result; result.tncs_normalised = tncs_vector != nullptr; @@ -721,20 +809,16 @@ TwinImmuneZoneResult AnalyzeTwinImmuneZones(const std::vector& if (subgroups.empty()) return result; - const auto strong = StrongReflections(p1_merged, gops, result); + const auto strong = StrongReflections(p1_merged, cell, gops, result); if (strong.empty()) return result; - 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); - result.control = ReadZone(control, tncs_vector); - result.control.operators = "acentric"; + result.control = ReadControl(strong, gops, tncs_vector); // 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)); + Calibrate(result); return result; } @@ -743,17 +827,22 @@ std::string TwinImmuneZonesToText(const TwinImmuneZoneResult& result) { os << std::fixed; os << "Twin-immune zones of the operators the adopted group adds over each index-2 subgroup (P1 merge, " << std::setprecision(2) << result.d_max_A << "-" << result.d_min_A << " A, shells with >= " << std::setprecision(0) - << ZONE_MIN_SHELL_ISIG << (result.tncs_normalised ? ", normalised per pseudo-translation class" : "") + << ZONE_MIN_SHELL_ISIG << ", a deviatoric quadratic form of delta B " << std::setprecision(1) + << result.anisotropy_delta_b_A2 << " A^2 fitted over these shells and taken out" + << (result.tncs_normalised ? ", normalised per pseudo-translation class" : "") << ")\n" << " Reflections centric in the group but not in the subgroup are their own twin mates: centric if\n" << " the added operators are real, acentric if they are a twin law or a pseudo-symmetry. <|E^2-1|>\n" << " reads 0.968 centric, 0.736 acentric; evidence is the centric/acentric log-likelihood ratio\n" - << " (> 0: the added operators are real).\n"; + << " (> 0: the added operators are real). The acentric control reads -0.130 nats per reflection\n" + << " when the normalisation is right; what it reads above that is the normalisation's (" << std::showpos + << std::setprecision(3) << result.control_excess_per_reflection << std::noshowpos + << " here), every\n zone carries the same per reflection, and the calibrated evidence has it taken off.\n"; auto row = [&](const TwinImmuneZone& z) { os << " " << std::left << std::setw(34) << z.operators << std::right << " n " << std::setw(6) << z.n << " <|E^2-1|> " << std::setprecision(3) << z.mean_abs_e2_minus_1 << " +- " << z.standard_error - << " evidence " << std::showpos << std::setprecision(1) << z.evidence_nats << std::noshowpos - << " nats\n"; + << " evidence " << std::showpos << std::setprecision(1) << z.evidence_nats << " nats, calibrated " + << z.calibrated_evidence_nats << std::noshowpos << "\n"; }; for (const auto& z : result.zones) row(z); diff --git a/image_analysis/scale_merge/TwinningAnalysis.h b/image_analysis/scale_merge/TwinningAnalysis.h index 7165e065e..9bebb3208 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.h +++ b/image_analysis/scale_merge/TwinningAnalysis.h @@ -9,6 +9,7 @@ #include "../../common/Reflection.h" #include "gemmi/symmetry.hpp" +#include "gemmi/unitcell.hpp" // Intensity-statistics test for merohedral twinning, in the spirit of phenix.xtriage. Two // indicators on the merged intensities, both normalised per resolution shell: @@ -101,7 +102,8 @@ std::string TwinningAnalysisToText(const TwinningAnalysisResult& result); std::string TwinningVerdictLine(const TwinningAnalysisResult& result); // Twin-immune zone evidence for the operators a point group adds over each of its index-2 subgroups, -// read on a P1 merge. REPORT-ONLY. +// read on a P1 merge (or on the adopted group's own): reported for every 2-fold the adopted group +// holds, and what a refused promotion is decided on (Rugnux TwinZoneVerdict). // // Under a twin law T, I_obs(h) = (1-a) I(h) + a I(Th); a reflection whose twin mate Th is h itself, up // to the true group and Friedel, is untouched by twinning at every fraction a. For an index-2 subgroup @@ -111,7 +113,8 @@ std::string TwinningVerdictLine(const TwinningAnalysisResult& result); // statistic is blind - at a = 0.5, where I(h) = I(Th) exactly. Read ABSOLUTELY, <|E^2-1|> against // 0.968 (centric) and 0.736 (acentric); never as the difference to the control, because noise, a // pseudo-translation and anisotropy all inflate every class towards centric, and a perfect twin's -// control (0.541) plus that inflation looks like true symmetry's. +// control (0.541) plus that inflation looks like true symmetry's. But an absolute reading is only as +// good as the normalisation, and the control is what certifies it: see control_excess_per_reflection. struct TwinImmuneZone { std::string operators; // the operators G adds over H, as h,k,l triplets int n = 0; @@ -122,24 +125,39 @@ struct TwinImmuneZone { // pseudo-symmetry. About +0.2 nats per reflection for a centric zone and -0.13 for an acentric one // on error-free data; less either way as the errors grow. double evidence_nats = 0.0; + // The same less n times the control's excess (TwinImmuneZoneResult::control_excess_per_reflection): + // what a verdict is read on. + double calibrated_evidence_nats = 0.0; }; struct TwinImmuneZoneResult { std::vector zones; // one per index-2 subgroup of the group TwinImmuneZone control; // the reflections acentric in the group + // The control's evidence per reflection above what an acentric population reads (-0.130 nats on + // error-free data), floored at zero. The control is acentric whatever the truth - twinned at most, + // which reads BELOW that - so anything above it is the normalisation's: anisotropy, a + // pseudo-translation, a pseudo-centring, noise, all of which inflate every class towards centric + // alike, and the zone, normalised the same way, carries the same per reflection. Measured: a 6/m + // crystal with a 67 A^2 anisotropy read its control at +0.10 nats per reflection (centric, in an + // isotropic normalisation) and the zone of a refused 622 the same, +0.10 - so the zone's +124 + // nats were the normalisation's, and the verdict rescued a twin law. Calibrated, that zone reads + // -115 nats; the zones of genuine promotions keep +140 nats (a 90 deg tetragonal sweep) to +2000. + double control_excess_per_reflection = 0.0; + 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 }; TwinImmuneZoneResult AnalyzeTwinImmuneZones(const std::vector& p1_merged, - const gemmi::SpaceGroup& group, + const gemmi::UnitCell& cell, const gemmi::SpaceGroup& group, const std::array* tncs_vector = nullptr); // The same for ONE hypothesis: the operators `group` adds over `subgroup` (one zone in the result, // beside the control). Read where a promotion from the subgroup to the group is in question. TwinImmuneZoneResult AnalyzeTwinImmuneZone(const std::vector& p1_merged, - const gemmi::SpaceGroup& group, const gemmi::SpaceGroup& subgroup, + const gemmi::UnitCell& cell, const gemmi::SpaceGroup& group, + const gemmi::SpaceGroup& subgroup, const std::array* tncs_vector = nullptr); std::string TwinImmuneZonesToText(const TwinImmuneZoneResult& result); diff --git a/rugnux/ResultReport.cpp b/rugnux/ResultReport.cpp index 8d5f33b59..08e8a067c 100644 --- a/rugnux/ResultReport.cpp +++ b/rugnux/ResultReport.cpp @@ -588,7 +588,8 @@ ReportDocument BuildReportDocument(const std::string &output_prefix, Add(s, KeyText("TWIN_ZONE_DECISION", result.twin_zone_decision)); Add(s, Prose(" A twin gate refused this promotion, and the twin-immune zone of the operators it\n" " adds decided it (reflections a twin law cannot touch: centric if the operators\n" - " are real, acentric if they are a twin law or a pseudo-symmetry; margin 20 nats):\n" + " are real, acentric if they are a twin law or a pseudo-symmetry; read less what\n" + " the acentric control shows the normalisation to add; margin 20 nats):\n" " " + result.twin_zone_decision + ".")); } if (refused) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 1e353fcb7..92d751687 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -1739,15 +1739,24 @@ namespace { // twins at -290 nats and below, with one genuine trigonal crystal at -46 nats whose promotion // no gate refused - so the bound is not what separates them, and a refusal a gate did fire is // the only place the verdict is read. - int TwinZoneVerdict(const std::vector &p1_merged, const SearchSpaceGroupResult &higher, - const SearchSpaceGroupResult &lower, ProcessResult &result, - const std::array *tncs_vector = nullptr) { + // + // Read on the CALIBRATED evidence: the zone's less what the acentric control shows the + // normalisation to add per reflection (TwinImmuneZoneResult::control_excess_per_reflection). The + // raw evidence rescued two refused 622 promotions - a 6/m crystal and a 312 one, H 4.7x and 5.2x + // the parent's, R_meas 0.09 -> 0.34 under the higher group - whose zones read centric only + // because everything did: a 67 A^2 anisotropy and a pseudo-translation had the control at + // +0.10 nats per reflection, the zone at the same, and +124 nats over 947 reflections was the + // normalisation's. Calibrated they read -115 and -86; the genuine promotions measured keep +140 + // and above. + int TwinZoneVerdict(const std::vector &p1_merged, const gemmi::UnitCell &cell, + const SearchSpaceGroupResult &higher, const SearchSpaceGroupResult &lower, + ProcessResult &result, const std::array *tncs_vector = nullptr) { constexpr double DECISIVE_NATS = 20.0; if (!higher.point_group_representative || !lower.point_group_representative) return 0; if (higher.point_group_order != 2 * lower.point_group_order) return 0; - const auto z = AnalyzeTwinImmuneZone(p1_merged, *higher.point_group_representative, + const auto z = AnalyzeTwinImmuneZone(p1_merged, cell, *higher.point_group_representative, *lower.point_group_representative, tncs_vector); if (z.zones.empty() || z.zones[0].n == 0) { result.twin_zone_decision = fmt::format("{} over {}: no zone reflections in the shells with signal " @@ -1756,13 +1765,14 @@ namespace { return 0; } const auto &zone = z.zones[0]; - const int verdict = zone.evidence_nats <= -DECISIVE_NATS ? -1 - : zone.evidence_nats >= DECISIVE_NATS ? +1 : 0; + const int verdict = zone.calibrated_evidence_nats <= -DECISIVE_NATS ? -1 + : zone.calibrated_evidence_nats >= DECISIVE_NATS ? +1 : 0; result.twin_zone_decision = fmt::format( - "{} over {}: operators {}; n {} <|E^2-1|> {:.3f} +- {:.3f} (control {:.3f}), evidence " - "{:+.1f} nats ({:.2f}-{:.2f} A{}) -> {}", + "{} over {}: operators {}; n {} <|E^2-1|> {:.3f} +- {:.3f} (control {:.3f}, {:+.3f} nats per " + "reflection above acentric), evidence {:+.1f} nats, calibrated {:+.1f} ({:.2f}-{:.2f} A{}) -> {}", higher.point_group_hm, lower.point_group_hm, zone.operators, zone.n, zone.mean_abs_e2_minus_1, - zone.standard_error, z.control.mean_abs_e2_minus_1, zone.evidence_nats, z.d_max_A, z.d_min_A, + zone.standard_error, z.control.mean_abs_e2_minus_1, z.control_excess_per_reflection, + zone.evidence_nats, zone.calibrated_evidence_nats, z.d_max_A, z.d_min_A, z.tncs_normalised ? ", normalised per pseudo-translation class" : "", verdict < 0 ? "acentric, the added operators are a twin law or a pseudo-symmetry" : verdict > 0 ? "centric, the added operators are real" : "undecided"); @@ -5692,7 +5702,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b alt.point_group_hm, sg_search.point_group_hm); overridden(sg_search, "Lorentz-filtered", "all-observation"); zone_verdict = !sg_search.refused_point_group_hm.empty() - ? TwinZoneVerdict(sm_all.merged, alt, sg_search, result) : 0; + ? TwinZoneVerdict(sm_all.merged, gemmi::UnitCell(*result.consensus_cell), alt, + sg_search, result) : 0; if (zone_verdict < 0) { logger.Warning("Space-group search: the Lorentz-filtered merge's refusal of {} " "STANDS - the twin-immune zone of its added operators reads " @@ -5711,7 +5722,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b } if (order_of(alt) < order_of(sg_search) && !alt.refused_point_group_hm.empty() - && TwinZoneVerdict(sm_all.merged, sg_search, alt, result) < 0) { + && TwinZoneVerdict(sm_all.merged, gemmi::UnitCell(*result.consensus_cell), sg_search, alt, + result) < 0) { logger.Warning("Space-group search: the all-observation merge's refusal of {} STANDS " "- the twin-immune zone of its added operators reads acentric; keeping {}", alt.refused_point_group_hm, alt.point_group_hm); @@ -5971,14 +5983,39 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // the near-tangential measurements that make a genuine operator look like a twin. // Where it finds the larger point group, the absences are still judged on all // observations, with the point group pinned. + // And the same twin-immune zone verdict on a refusal the other arm would outvote: + // a P31 partial twin whose 32 the rule above had refused on the zone was promoted + // by this ask's filtered arm - the same override, one copy later. if (!merged_filtered.empty()) { SearchSpaceGroupOptions of = o2; of.merge_isa = merged_filtered_isa; // this arm's own error model - const auto filtered = SearchSpaceGroup(reindexed(merged_filtered), of); - if (filtered.point_group_order > s2.point_group_order) { + auto filtered = SearchSpaceGroup(reindexed(merged_filtered), of); + const auto refusal_stands = [&](const SearchSpaceGroupResult &higher, + SearchSpaceGroupResult &lower) { + if (higher.point_group_order <= lower.point_group_order + || lower.refused_point_group_hm.empty() + || TwinZoneVerdict(merged_all, *o2.cell, higher, lower, result) >= 0) + return false; + logger.Warning("Space-group search (metric lattice): the refusal of {} STANDS - the " + "twin-immune zone of its added operators reads acentric; keeping {}", + lower.refused_point_group_hm, lower.point_group_hm); + lower.refused_reason += ", and its twin-immune zone reads acentric"; + return true; + }; + if (filtered.point_group_order > s2.point_group_order && !refusal_stands(filtered, s2)) { o2.fixed_point_group = filtered.point_group_representative; s2 = SearchSpaceGroup(merged_all, o2); o2.fixed_point_group.reset(); + } else if (refusal_stands(s2, filtered)) { + o2.fixed_point_group = filtered.point_group_representative; + auto lower = SearchSpaceGroup(merged_all, o2); + o2.fixed_point_group.reset(); + lower.refused_point_group_hm = filtered.refused_point_group_hm; + lower.refused_reason = filtered.refused_reason; + lower.refused_point_group_representative = filtered.refused_point_group_representative; + lower.refused_by_twin_gate = true; + lower.refused_twin_fraction = filtered.refused_twin_fraction; + s2 = std::move(lower); } } // The centring cannot be read from the absences in this setting: reindexing @@ -6494,8 +6531,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // of this merge, at full multiplicity and production scale - where the P1 search // merge of a weak crystal has no shell with signal enough to read a zone on. const int zone = hi_sg.has_value() - ? TwinZoneVerdict(sm.merged, hi_search, sg_search, result, - TranslationalNCSVectorForLTest(pre_tncs)) : 0; + ? TwinZoneVerdict(sm.merged, gemmi::UnitCell(*result.consensus_cell), hi_search, + sg_search, result, TranslationalNCSVectorForLTest(pre_tncs)) : 0; if (zone < 0) { const std::string msg = fmt::format( "Point group {} was refused ({}), and the twin-immune zone of its added " @@ -7257,7 +7294,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b // Report-only: nothing reads it back. if (determined != nullptr && determined->number > 1) { result.twin_immune_zones = AnalyzeTwinImmuneZones( - p1.merged, *determined, TranslationalNCSVectorForLTest(result.tncs)); + p1.merged, gemmi::UnitCell(*result.consensus_cell), *determined, + TranslationalNCSVectorForLTest(result.tncs)); if (!result.twin_immune_zones->zones.empty()) logger.Info("{}", TwinImmuneZonesToText(*result.twin_immune_zones)); } diff --git a/tests/TwinningAnalysisTest.cpp b/tests/TwinningAnalysisTest.cpp index 06c6b8d60..81dbfad9f 100644 --- a/tests/TwinningAnalysisTest.cpp +++ b/tests/TwinningAnalysisTest.cpp @@ -223,8 +223,12 @@ namespace { // intensity from the acentric distribution, which is all the searches it was written for need, but // the zone test reads exactly the difference. E^2 is chi^2 with one degree of freedom for a // reflection centric in the true group (a squared Box-Muller normal), exponential otherwise, keyed - // on the asymmetric-unit index so the true group's symmetry is exact. - std::vector WilsonTetragonal(const char *true_group, double twin_fraction) { + // on the asymmetric-unit index so the true group's symmetry is exact. `anisotropy_b` is an extra + // Debye-Waller B along c* alone (A^2) - the same in both twin domains, a twin law being a lattice + // operation - and `scale_jitter` a log-normal factor of that width on every observed intensity, + // a nuisance with no direction that no normalisation removes. + std::vector WilsonTetragonal(const char *true_group, double twin_fraction, + double anisotropy_b = 0.0, double scale_jitter = 0.0) { const gemmi::SpaceGroup &sub = gemmi::get_spacegroup_by_name(true_group); const gemmi::SpaceGroup &super = gemmi::get_spacegroup_by_name("P 4 2 2"); const gemmi::Op twin = jfjoch_test::TwinLaw(sub, super); @@ -237,7 +241,9 @@ namespace { const double u2 = jfjoch_test::detail::UniformFromHkl({{asu[0] + 300, asu[1] + 300, asu[2] + 300}}); const double c = std::cos(2.0 * std::numbers::pi * u2); const double e2 = gops.is_reflection_centric(hkl) ? -2.0 * std::log(u1) * c * c : -std::log(u1); - return 1000.0 * e2 * std::exp(-20.0 / (2.0 * std::pow(cell.calculate_d(hkl), 2))); + const double s_c = hkl[2] / cell.c; + return 1000.0 * e2 * std::exp(-20.0 / (2.0 * std::pow(cell.calculate_d(hkl), 2))) + * std::exp(-anisotropy_b * s_c * s_c / 2.0); }; std::vector out; for (int h = -19; h <= 19; ++h) @@ -254,8 +260,12 @@ namespace { r.k = k; r.l = l; r.d = static_cast(d); - r.I = static_cast((1.0 - twin_fraction) * true_intensity(hkl) - + twin_fraction * true_intensity(twin.apply_to_hkl(hkl))); + const double u3 = jfjoch_test::detail::UniformFromHkl({{h + 700, k + 700, l + 700}}); + const double u4 = jfjoch_test::detail::UniformFromHkl({{h + 900, k + 900, l + 900}}); + const double jitter = std::exp(scale_jitter * std::sqrt(-2.0 * std::log(u3)) + * std::cos(2.0 * std::numbers::pi * u4)); + r.I = static_cast(((1.0 - twin_fraction) * true_intensity(hkl) + + twin_fraction * true_intensity(twin.apply_to_hkl(hkl))) * jitter); r.sigma = static_cast(0.01 * r.I + 1.0); out.push_back(r); } @@ -268,8 +278,9 @@ namespace { // stay acentric under it and are centric only where the added operators are real. TEST_CASE("Twin-immune zones tell a perfect twin from real symmetry", "[twinning]") { const gemmi::SpaceGroup &p422 = gemmi::get_spacegroup_by_name("P 4 2 2"); + const gemmi::UnitCell tetragonal(47, 47, 63, 90, 90, 90); - const auto real = AnalyzeTwinImmuneZones(WilsonTetragonal("P 4 2 2", 0.0), p422); + const auto real = AnalyzeTwinImmuneZones(WilsonTetragonal("P 4 2 2", 0.0), tetragonal, p422); REQUIRE(real.zones.size() == 3); // over 4, over 222, over the diagonal 222 for (const auto &z : real.zones) { CAPTURE(z.operators); @@ -280,11 +291,56 @@ TEST_CASE("Twin-immune zones tell a perfect twin from real symmetry", "[twinning // A perfect P4 twin: whichever subgroup is taken as the parent, the operators added over it are // not all real, and the zone reads acentric. - const auto twin = AnalyzeTwinImmuneZones(WilsonTetragonal("P 4", 0.5), p422); + const auto twin = AnalyzeTwinImmuneZones(WilsonTetragonal("P 4", 0.5), tetragonal, p422); REQUIRE(twin.zones.size() == 3); for (const auto &z : twin.zones) { CAPTURE(z.operators); CHECK(z.mean_abs_e2_minus_1 == Catch::Approx(0.736).margin(0.08)); CHECK(z.evidence_nats < 0.0); + CHECK(z.calibrated_evidence_nats < 0.0); + } + // Wilson data normalise cleanly: the control reads the acentric expectation and the + // calibration takes nothing off. + CHECK(real.control_excess_per_reflection < 0.01); + CHECK(twin.control_excess_per_reflection < 0.01); +} + +// The zone is read absolutely, so it is only as good as the normalisation. A strong anisotropy +// leaves every class of an isotropically normalised merge reading centric, the acentric control +// included, and a zone - a plane - more so than the control's sphere; it is fitted and taken out +// before anything is read. What no normalisation removes, the control certifies: its excess over the +// acentric expectation is taken off every zone, so a partial twin's zone reads acentric and a real +// operator's centric, as they do on clean data. +TEST_CASE("Twin-immune zones are normalised for anisotropy and calibrated by the acentric control", + "[twinning]") { + const gemmi::SpaceGroup &p422 = gemmi::get_spacegroup_by_name("P 4 2 2"); + const gemmi::SpaceGroup &p4 = gemmi::get_spacegroup_by_name("P 4"); + const gemmi::UnitCell tetragonal(47, 47, 63, 90, 90, 90); + + SECTION("a 60 A^2 anisotropy along c* is taken out") { + const auto twin = AnalyzeTwinImmuneZone(WilsonTetragonal("P 4", 0.2, 60.0), tetragonal, p422, p4); + REQUIRE(twin.zones.size() == 1); + CHECK(twin.anisotropy_delta_b_A2 == Catch::Approx(60.0).margin(15.0)); + CHECK(twin.control.mean_abs_e2_minus_1 < 0.76); // a partial twin's control: at or below acentric + CHECK(twin.zones[0].calibrated_evidence_nats < -20.0); + + const auto real = AnalyzeTwinImmuneZone(WilsonTetragonal("P 4 2 2", 0.0, 60.0), tetragonal, p422, p4); + REQUIRE(real.zones.size() == 1); + CHECK(real.control.mean_abs_e2_minus_1 == Catch::Approx(0.736).margin(0.05)); + CHECK(real.zones[0].calibrated_evidence_nats > 20.0); + } + + SECTION("a scale jitter no normalisation removes is calibrated off") { + const auto real = AnalyzeTwinImmuneZone(WilsonTetragonal("P 4 2 2", 0.0, 0.0, 0.6), tetragonal, p422, p4); + REQUIRE(real.zones.size() == 1); + CHECK(real.control.mean_abs_e2_minus_1 > 0.80); // inflated: the normalisation is wrong + CHECK(real.control_excess_per_reflection > 0.02); + CHECK(real.zones[0].calibrated_evidence_nats > 20.0); + + // The twin's control is deflated by the twinning as much as the jitter inflates it, so + // nothing is taken off - and its zone, acentric, still reads so. + const auto twin = AnalyzeTwinImmuneZone(WilsonTetragonal("P 4", 0.2, 0.0, 0.6), tetragonal, p422, p4); + REQUIRE(twin.zones.size() == 1); + CHECK(twin.zones[0].calibrated_evidence_nats < -20.0); } }