From f9ffc3d83733ae0ae07b4855ec2e04fb04545a5b Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Mon, 5 Oct 2026 10:13:50 +0200 Subject: [PATCH] SearchSpaceGroup: ask a promotion's twin-immune zone on a twinned merge, read at the other twin laws' fraction 6iu9 (deposited P3_1, merohedrally twinned at 0.4 or more by its 321 law: H 0.097, CC 0.79) went from P3_1 to P3_121 when the penalised per-frame scale smoother (495f6a97b) landed; the ingest partiality fix alone (7d6c201e9) keeps P3_1. Bisected on the two binaries: smoother-only 6iu9 P3_121, partiality-only P3_1. The smoother only removed the luck: the P3_1 answer rested on the Lorentz-filtered arm's added-operator R contrast reading 0.71 against a bound of 0.72 (the all-observation arm already passed 32 at 0.82); better scaling moved it to 0.87, no gate refused, and the twin-immune zone - which reads -401 nats, acentric - is only consulted where a gate fired. A near-perfect twin passes every agreement gate by construction (H ratio 0.97 here); only the zone can refuse it. It was confined to refusals because a genuine trigonal crystal read -46 nats with no gate firing; on the battery 5uth (genuine P3_121, twinned by a 622 law, other-law CC 0.14) reads -206. The reason: the centric density was untwinned, but a twin by a law OTHER than the zone's operators reaches the zone as it reaches every reflection, so a genuine zone reads 0.81 at a = 0.2. The operators' own twin law cannot reach their zone, so only the centric side needs it. - TwinningAnalysis: the zone's centric density is the twinned one at a caller-given fraction (weighted sum of two chi^2_1, via exp(-y) I0(y)); the control's calibration expectation is taken at the same fraction (-KL(acentric || centric_a), -0.130 at a = 0 as before). At a = 0 nothing changes. - SearchSpaceGroup: twin_fraction_outside = the fraction implied by the strongest CC of a lattice rotation outside the group, relative to the group's own mean CC, through rho = 2a(1-a)/((1-a)^2+a^2); exported for the adopted group and used by Rugnux's TwinZoneVerdict and the zone report. - New Stage A test: a candidate that does not hold every lattice rotation, on a P1 merge whose <|L|> is in the partial-twin band [0.375, 0.44), is refused when the zone over one of its index-2 subgroups reads acentric by 20 nats (the TwinZoneVerdict bound), unless zones_ambiguous. The band's lower end is the L-test gate's: below it something else compresses the zones too (a pseudo-cubic small-molecule set, cuhf2, read its 422 zone at 0.64 beside a control at 0.58 and was turned to P222 without it). Zone evidence (calibrated, at the other-law fraction): 6iu9 32 -230 / -143 (refused, P3_1); 5uth +209 / +33 (P3_121 kept); 8xtg +324 / +65; 6vww P6 +20 / +8.5; 7k1l P6 +8 / +7.5. Battery (open+inhouse, targeted, against 20261004-2350 all2-full; only beyond-noise change listed): - 6iu9: fail -> pass, P3_121 -> P3_1, R_meas 17.4% -> 14.4%, ISa 4.6 -> 5.4, R-free 0.318 -> 0.327. - unchanged: 6iu5, 6iu6, 6iu8 (P3_1), 5uth, 5j23, insu_H_x06da_twin/notwin, 6vww, 7k1l, 8xte, 8xtg, 9i80, 4bwl, 2wnq, 8c3e, 2wnn, 2xfw, 5ebi, 6p8j, 6z9g, 6rlr, 6toc, 3mc4; 321/312 controls 5lzl, 6w4h, 9gqg, 6pxb, 9z72; all 14 small-molecule sets (cuhf2 checked again after the band). - private subset (8 sets): no change beyond noise. Tests: [twinning] (new: zones read at the other twin laws' fraction), SearchSpaceGroup*. Co-Authored-By: Claude Opus 5.5 (1M context) Claude-Session: https://claude.ai/code/session_01K5K8jvPPbmCrbqnWkddTuB --- docs/CHANGELOG.md | 1 + docs/CPU_DATA_ANALYSIS_DECISIONS.md | 4 +- .../scale_merge/SearchSpaceGroup.cpp | 107 +++++++++++++- image_analysis/scale_merge/SearchSpaceGroup.h | 4 + .../scale_merge/TwinningAnalysis.cpp | 135 +++++++++++++----- image_analysis/scale_merge/TwinningAnalysis.h | 18 ++- rugnux/Rugnux.cpp | 19 ++- tests/TwinningAnalysisTest.cpp | 45 +++++- 8 files changed, 275 insertions(+), 58 deletions(-) diff --git a/docs/CHANGELOG.md b/docs/CHANGELOG.md index 26eb487a0..d2bbfc018 100644 --- a/docs/CHANGELOG.md +++ b/docs/CHANGELOG.md @@ -11,6 +11,7 @@ * Rugnux: French-Wilson amplitudes (`F`/`SIGF`) use an anisotropic Wilson prior from the fitted anisotropy tensor, as ctruncate does; intensities are unchanged. * Rugnux claims a screw axis from a short axial row of a few weak reflections, so such crystals (e.g. P2_1 with a ~30 A unique axis) are no longer written without the screw. * Rugnux keeps the screw axes it found when a higher point group is adopted after the twin-law check, instead of writing the group without screws (e.g. P 4 2 2 for P4_2 2_1 2). +* Rugnux keeps a nearly perfect merohedral twin in its true point group (e.g. P3_1 rather than P3_121) when the reflections the twin law cannot touch read acentric. * Rugnux drops a rotation reflection whose spot held a saturated (overloaded) pixel, as XDS does, and reports the count as `OBSERVATIONS_REJECTED_OVERLOAD=`. * Rugnux fits a rotation sweep's per-frame scale on the partial reflections, even with few reflections per frame, when scaling from the fulls alone leaves no measurable error model. * Rugnux names glide planes in groups without a centre of symmetry (e.g. I-42d, Fdd2), reports a group with the same absences (C2/c and Cc) as an alternative, and states how the centre was decided (`SPACE_GROUP_CENTRE`, `CENTRE_STATISTICS_*`). diff --git a/docs/CPU_DATA_ANALYSIS_DECISIONS.md b/docs/CPU_DATA_ANALYSIS_DECISIONS.md index 418fdcad3..4ff7caf17 100644 --- a/docs/CPU_DATA_ANALYSIS_DECISIONS.md +++ b/docs/CPU_DATA_ANALYSIS_DECISIONS.md @@ -48,13 +48,15 @@ Both search merges also drop the **frames whose fitted per-frame scale came out **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. That merge no longer holds the adopted group's screw-absent axial reflections, so the higher group also **keeps the screws the adopted group decided** — among the higher candidates the merge cannot separate, the one that predicts every such axial absence is taken — rather than the lowest-numbered one, which claims none. The cost is two extra merges, and only on a run that records a refusal, which is a few in a hundred. +**On a merge that reads twinned, the search asks the zone itself.** A twin of an index-2 subgroup by exactly the operators a promotion adds defeats every agreement gate once its fraction is high — the added operators then agree like real ones — so where the P1 merge's $\langle|L|\rangle$ lies in the partial-twin band (0.375 to 0.44) and the candidate does not hold every rotation of its lattice (so a twin law outside it is possible), the candidate is refused when the twin-immune zone of one of its index-2 subgroups reads acentric by 20 nats. Below 0.375 something other than a twin compresses the intensities, the zones with them, and they are not read. A $P3_1$ crystal twinned at $\alpha \ge 0.4$ by its 321 law, whose 32 promotion every agreement gate passed, reads −230 nats there. + **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. ### 13.2 Twinning check, and translational pseudo-symmetry 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. 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. +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). 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 `_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/SearchSpaceGroup.cpp b/image_analysis/scale_merge/SearchSpaceGroup.cpp index e9c0ef21d..b08ee302a 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.cpp +++ b/image_analysis/scale_merge/SearchSpaceGroup.cpp @@ -1546,6 +1546,50 @@ SearchSpaceGroupResult SearchSpaceGroup( if (opt.cell.has_value()) lattice_rotations = gemmi::find_lattice_symmetry(*opt.cell, opt.lattice_centring, LATTICE_TWIN_OBLIQUITY_DEG); + // The twin fraction the lattice operations OUTSIDE a group imply by their own correlation - what a + // twin law other than the group's operators compresses every reflection by, the zone of a genuine + // operator included (AnalyzeTwinImmuneZones' twin_fraction). The strongest of those operators' CCs, + // as a fraction of the mean CC of the group's own operators: the same noise holds a genuine + // operator's CC short of 1 and a twin law's short of what its fraction puts there. Inverted through + // rho = 2a(1-a) / ((1-a)^2 + a^2), the correlation a twin of fraction a puts between I(h) and I(Th). + // 0 where there is no lattice to ask, or no operator outside the group with pairs enough. + const auto twin_fraction_outside = [&](const auto &pg) { + if (!lattice_rotations || pg.rotations.empty()) + return 0.0; + double cc_own = 0.0; + int n_own = 0; + for (const auto &rot : pg.rotations) { + const auto &os = operator_score(rot); + if (rot.rot != gemmi::Op::identity().rot && os.n_pairs >= opt.min_pairs_per_operator) { + cc_own += os.cc; + ++n_own; + } + } + double cc_outside = 0.0; + for (const auto &op : lattice_rotations->sym_ops) { + if (op.rot == gemmi::Op::identity().rot + || std::binary_search(pg.rotation_set.begin(), pg.rotation_set.end(), RotKey(op))) + continue; + const auto &os = operator_score(op); + if (os.n_pairs >= opt.min_pairs_per_operator) + cc_outside = std::max(cc_outside, os.cc); + } + if (n_own == 0 || !(cc_own > 0.0)) + return 0.0; + const double rho = std::min(1.0, cc_outside / (cc_own / n_own)); + return 0.5 * (1.0 - std::sqrt((1.0 - rho) / (1.0 + rho))); + }; + // Whether the merge reads twinned, for the twin-immune zone gate: its <|L|> in the band a partial + // twin puts it in, from the perfect twin's 0.375 up to AnalyzeTwinning's bar of 0.44, the band the + // L-test gate below reads too. Measured once, on first use. + std::optional twinned; + const auto merge_reads_twinned = [&] { + if (!twinned) { + const auto t = AnalyzeTwinning(merged, nullptr, 20, nullptr, opt.lattice_centring); + twinned = t.l_test_pairs > 0 && t.mean_abs_l >= 0.375 && t.mean_abs_l < 0.44; + } + return *twinned; + }; for (auto& c : pg_cands) { // A genuine symmetry operator merges equivalent reflections, so it barely changes the reduced // chi^2 relative to the best subgroup - across the whole rotation-test battery every correct @@ -1761,15 +1805,16 @@ SearchSpaceGroupResult SearchSpaceGroup( // - A perfect twin is out of reach by construction: there the unmerged data already read 0.375 // and averaging cannot move them, which is the theorem that no intensity statistic separates // a perfect twin from the higher group. - bool l_refused = false; - LTestUnderMerge l_test; - if (consistent && !c.pg->rotations.empty() && c.pg->representative && lattice_rotations + const bool holds_lattice = lattice_rotations && std::all_of(lattice_rotations->sym_ops.begin(), lattice_rotations->sym_ops.end(), [&](const gemmi::Op &op) { return op.rot == gemmi::Op::identity().rot || std::binary_search(c.pg->rotation_set.begin(), c.pg->rotation_set.end(), RotKey(op)); - })) { + }); + bool l_refused = false; + LTestUnderMerge l_test; + if (consistent && !c.pg->rotations.empty() && c.pg->representative && holds_lattice) { l_test = AnalyzeLTestUnderMerge(merged, c.pg->rotations, c.pg->representative->centring_type()); constexpr double PERFECT_TWIN_L = 0.375; constexpr double TWIN_LIKE_L = 0.44; // AnalyzeTwinning's own bar for a twin-like <|L|> @@ -1782,17 +1827,61 @@ SearchSpaceGroupResult SearchSpaceGroup( if (l_refused) consistent = false; + // Twin-immune zone test, for a group that is NOT the whole symmetry of its lattice, on a merge + // that reads twinned: the other half of the twin argument. Where a lattice operation lies + // outside the candidate, a twin of one of its index-2 subgroups by exactly the operators the + // candidate adds over it is possible, and at a high twin fraction it defeats every gate above - + // the added operators then agree like real ones (H, R, chi^2 all read "real"), which is the + // theorem that no statistic of the operators' agreement separates a near-perfect twin from the + // higher group. What it cannot fake are the reflections centric in the candidate but not in the + // subgroup: they are their own twin mates, so the twin leaves them acentric and only a real + // operator makes them centric (AnalyzeTwinImmuneZones, read at the twin fraction the lattice's + // OTHER operations show, so that a twin by some other law does not make a genuine zone read + // acentric - twin_fraction_outside). Refused when a zone reads acentric by 20 nats - the bound a + // refusal's zone verdict is read on in Rugnux (TwinZoneVerdict) - unless the control is + // compressed beyond any twin, where the zones cannot separate the two (zones_ambiguous). + // - Measured: a P3_1 crystal twinned at 0.4 or more by its 321 law, whose 32 promotion every + // gate above passed (contrast 0.87, H ratio 0.97), reads -230 nats; a genuine P3_121 + // crystal twinned by a 622 law, read at its other law's fraction 0.07, +209 (-206 read + // untwinned, which is why the fraction is read). + // - Only on a merge whose L-test reads twinned (merge_reads_twinned, on the P1 merge before any + // promotion): without twinning a twin law cannot pass the gates above, and a crystal that is + // not twinned is not asked a twin's question. Below the perfect twin's 0.375 something other + // than a twin compresses the intensities - overlapping spots, or the pseudo-cubic small- + // molecule metric whose zones read 0.64 beside a control at 0.58, below acentric itself - and + // it compresses the zones with them, so they are not read there either. + constexpr double ZONE_DECISIVE_NATS = 20.0; + bool zone_refused = false; + std::string zone_operators; + double zone_evidence = 0.0; + if (consistent && !c.pg->rotations.empty() && c.pg->representative && opt.cell.has_value() + && lattice_rotations && !holds_lattice && merge_reads_twinned()) { + const auto zones = AnalyzeTwinImmuneZones(merged, *opt.cell, *c.pg->representative, nullptr, + opt.nthreads, twin_fraction_outside(*c.pg)); + if (!zones.zones_ambiguous) + for (const auto &z : zones.zones) + if (z.n > 0 && z.calibrated_evidence_nats <= -ZONE_DECISIVE_NATS + && z.calibrated_evidence_nats < zone_evidence) { + zone_refused = true; + zone_operators = z.operators; + zone_evidence = z.calibrated_evidence_nats; + } + } + if (zone_refused) + consistent = false; + if (!consistent) { // Same three sentences the highest refusal gets below, but kept per candidate for the // ledger: short, because the long form is written once for the group the user is told about. c.why = h_refused ? "H ratio" : r_refused ? "added-operator R" - : spread_refused ? "operator spread" : l_refused ? "L-test under merge" : "merge chi^2"; + : spread_refused ? "operator spread" : l_refused ? "L-test under merge" + : zone_refused ? "twin-immune zone" : "merge chi^2"; // Record the highest-order refusal so the caller can say WHY it is processing lower. if (c.order > refused_order && c.pg->representative) { refused_order = c.order; refused_pg_hm = c.pg->representative->point_group_hm(); refused_pg_rep = c.pg->representative; - refused_twin_gate = h_refused || r_refused || spread_refused || l_refused; + refused_twin_gate = h_refused || r_refused || spread_refused || l_refused || zone_refused; refused_fraction = h_refused ? std::max(0.0, 0.5 - h_added) : std::numeric_limits::quiet_NaN(); if (h_refused) @@ -1838,6 +1927,11 @@ SearchSpaceGroupResult SearchSpaceGroup( + FormatDouble(l_test.mean_abs_l_merged, 3) + " averaged over its orbits, " "on the same " + std::to_string(l_test.pairs) + " pairs - so an operator " "of it averages unequal intensities"; + else if (zone_refused) + refused_why = "the merge reads twinned, and the reflections its operators " + zone_operators + + " leave as their own twin mates read acentric (" + FormatDouble(zone_evidence, 1) + + " nats against centric, bound " + FormatDouble(ZONE_DECISIVE_NATS, 0) + + ") - so they are a twin law of the subgroup, not a symmetry"; else if (std::isfinite(c.chi2) && std::isfinite(chi2_ref)) refused_why = "merge chi^2 is " + FormatDouble(c.chi2 / chi2_ref, 2) + "x the subgroup's (bound " + FormatDouble(opt.max_merge_chi2_ratio, 2) + ")"; @@ -1892,6 +1986,7 @@ SearchSpaceGroupResult SearchSpaceGroup( // ratio of the group it forced rather than of the one Stage A would have taken. result.global_best_operator_r = std::isfinite(global_best_r) ? global_best_r : std::numeric_limits::quiet_NaN(); + result.twin_fraction_outside = twin_fraction_outside(*best_pg); for (const auto& c : pg_cands) if (c.pg == best_pg) { result.generated_point_group_adopted = c.closure; diff --git a/image_analysis/scale_merge/SearchSpaceGroup.h b/image_analysis/scale_merge/SearchSpaceGroup.h index 58f087b9b..c781ed1d1 100644 --- a/image_analysis/scale_merge/SearchSpaceGroup.h +++ b/image_analysis/scale_merge/SearchSpaceGroup.h @@ -804,6 +804,10 @@ struct SearchSpaceGroupResult { double r_over_best = std::numeric_limits::quiet_NaN(); double r_contrast = std::numeric_limits::quiet_NaN(); double global_best_operator_r = std::numeric_limits::quiet_NaN(); + // The twin fraction the lattice operations outside the ADOPTED point group imply by their own + // correlation with the intensities - what the centric side of its twin-immune zones is read at + // (AnalyzeTwinImmuneZones' twin_fraction). 0 where no cell was given. + double twin_fraction_outside = 0.0; // The intensity-weighted R of UNRELATED reflections on this merge - shell-matched pairs of // reflections no symmetry relates. The far end of the contrast scale above: where a false diff --git a/image_analysis/scale_merge/TwinningAnalysis.cpp b/image_analysis/scale_merge/TwinningAnalysis.cpp index 1788e6955..4f806f0a4 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.cpp +++ b/image_analysis/scale_merge/TwinningAnalysis.cpp @@ -616,37 +616,89 @@ 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; + // exp(-y) I0(y), I0 the modified Bessel function of order zero: its power series, and above y = 50 + // (where the series would want a hundred terms) the first terms of the asymptotic expansion. + double I0Scaled(double y) { + if (y >= 50.0) + return (1.0 + 1.0 / (8.0 * y) + 9.0 / (128.0 * y * y)) / std::sqrt(2.0 * std::numbers::pi * y); + double term = 1.0, sum = 1.0; + for (int k = 1; term > 1e-17 * sum; ++k) { + term *= y * y / (4.0 * k * k); + sum += term; + } + return std::exp(-y) * sum; + } - // 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 - // zero, reads every noisy weak reflection as decisive; flooring E^2 at s instead turns the centric - // excess of weak reflections into moderate values and read a genuine centric zone of a 1.2 A - // crystal as acentric. Integrated over u = sqrt(E^2_true), where neither density is singular, - // across +-6 s of the measurement (trapezoid). - double CentricOverAcentric(double e2, double s) { + // The Wilson densities of E^2 = u^2 as densities in u (times dE^2/du = 2u), where neither is + // singular at zero: the untwinned acentric (exponential) one, and the centric one twinned at a + // fraction a - every observed intensity (1 - a) I1 + a I2 of two independent centric reflections. + // That is the density of a weighted sum of two chi^2 variables of one degree of freedom, c1 = 1 - a + // and c2 = a: exp(-x/(2 c1)) exp(-y) I0(y) / (2 sqrt(c1 c2)), y = x (c1 - c2) / (4 c1 c2); at a = 0 + // the plain chi^2 density, at a = 0.5 the exponential. + double AcentricDensity(double u) { + return 2.0 * u * std::exp(-u * u); + } + + double CentricDensity(double u, double a) { + if (a == 0.0) + return std::sqrt(2.0 / std::numbers::pi) * std::exp(-0.5 * u * u); + const double c1 = 1.0 - a, c2 = a; + return u / std::sqrt(c1 * c2) * std::exp(-u * u / (2.0 * c1)) + * I0Scaled(u * u * (c1 - c2) / (4.0 * c1 * c2)); + } + + // A density of E^2 convolved with the measurement error of one reflection, E^2 measured with error + // s, up to a factor common to every density. Without the convolution the centric density, infinite + // at zero, reads every noisy weak reflection as decisive; flooring E^2 at s instead turns the + // centric excess of weak reflections into moderate values and read a genuine centric zone of a + // 1.2 A crystal as acentric. Integrated over u = sqrt(E^2_true) across +-6 s of the measurement + // (trapezoid). + template + double Convolved(double e2, double s, Density density) { constexpr int n = 200; const double u_lo = std::sqrt(std::max(0.0, e2 - 6.0 * s)); const double u_hi = std::sqrt(std::max(0.0, e2 + 6.0 * s)); if (!(u_hi > u_lo)) return 0.0; const double du = (u_hi - u_lo) / n; - double p_centric = 0.0, p_acentric = 0.0; + double p = 0.0; for (int i = 0; i <= n; ++i) { const double u = u_lo + i * du; - const double w = (i == 0 || i == n ? 0.5 : 1.0) * std::exp(-0.5 * std::pow((e2 - u * u) / s, 2)); - p_centric += w * std::sqrt(2.0 / std::numbers::pi) * std::exp(-0.5 * u * u); - p_acentric += w * 2.0 * u * std::exp(-u * u); + p += (i == 0 || i == n ? 0.5 : 1.0) * std::exp(-0.5 * std::pow((e2 - u * u) / s, 2)) * density(u); } + return p; + } + + // What CentricOverAcentric averages over an error-free acentric population, with the centric side + // twinned at a: -KL(acentric || centric at a), by the midpoint rule in u. At a = 0, 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 + // population reads +0.216); -0.038 at a = 0.1 and -0.012 at 0.2, the two densities closing in. + // Errors move it towards zero, so an error-free baseline takes the larger excess off a noisy + // control - the conservative side. + double AcentricEvidencePerReflection(double a) { + constexpr int n = 4000; + constexpr double u_max = 7.0; // exp(-49): nothing of the acentric density is left + const double du = u_max / n; + double sum = 0.0; + for (int i = 0; i < n; ++i) { + const double u = (i + 0.5) * du; + const double p_acentric = AcentricDensity(u); + sum += p_acentric * std::log(CentricDensity(u, a) / p_acentric) * du; + } + return sum; + } + + // ln p_centric(E^2) - ln p_acentric(E^2) for a measured E^2 with error s: the centric density at + // twin fraction a, the acentric one untwinned (see AnalyzeTwinImmuneZones). + double CentricOverAcentric(double e2, double s, double a) { + const double p_centric = Convolved(e2, s, [a](double u) { return CentricDensity(u, a); }); + const double p_acentric = Convolved(e2, s, [](double u) { return AcentricDensity(u); }); return p_centric > 0.0 && p_acentric > 0.0 ? std::log(p_centric / p_acentric) : 0.0; } + // `twin_fraction` is the compression the centric hypothesis is read at (CentricOverAcentric). TwinImmuneZone ReadZone(std::vector refl, const std::array* tncs, - size_t nthreads) { + size_t nthreads, double twin_fraction) { TwinImmuneZone z; // The evidence of each kept reflection is a pure function of its (E^2, error) and costs a few // hundred exponentials, so the pairs are collected here and evaluated together below; the sum @@ -698,7 +750,7 @@ namespace { ParallelChunks(static_cast(kept_e2_s.size()), ThreadsForWork(kept_e2_s.size(), nthreads, 1024), [&](int lo, int hi) { for (int i = lo; i < hi; ++i) - evidence[i] = CentricOverAcentric(kept_e2_s[i].first, kept_e2_s[i].second); + evidence[i] = CentricOverAcentric(kept_e2_s[i].first, kept_e2_s[i].second, twin_fraction); }); for (const double e : evidence) z.evidence_nats += e; @@ -715,7 +767,8 @@ namespace { 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); + 0.0, result.control.evidence_nats / result.control.n + - AcentricEvidencePerReflection(result.twin_fraction)); result.control.calibrated_evidence_nats = result.control.evidence_nats - result.control.n * result.control_excess_per_reflection; for (auto& z : result.zones) @@ -844,14 +897,17 @@ namespace { return strong; } - // The reflections acentric in the group: the control every zone is read beside. + // The reflections acentric in the group: the control every zone is read beside, with the evidence + // read as the zones' is (at the same twin fraction), which is what the calibration measures the + // normalisation's excess on. TwinImmuneZone ReadControl(const std::vector& strong, const gemmi::GroupOps& gops, - const std::array* tncs_vector, size_t nthreads) { + const std::array* tncs_vector, double twin_fraction, + size_t nthreads) { std::vector control; for (const auto& r : strong) if (!gops.is_reflection_centric(gemmi::Op::Miller{{r.h, r.k, r.l}})) control.push_back(&r); - TwinImmuneZone z = ReadZone(control, tncs_vector, nthreads); + TwinImmuneZone z = ReadZone(control, tncs_vector, nthreads, twin_fraction); z.operators = "acentric"; return z; } @@ -861,7 +917,8 @@ namespace { TwinImmuneZone ReadSubgroupZone(const std::vector& strong, const gemmi::GroupOps& gops, const std::vector& rotations, const std::vector& h_ops, - const std::array* tncs_vector, size_t nthreads) { + const std::array* tncs_vector, double twin_fraction, + size_t nthreads) { auto contains = [&](const std::vector& set, const gemmi::Op& op) { return std::any_of(set.begin(), set.end(), [&](const gemmi::Op& o) { return o.rot == op.rot; }); }; @@ -875,7 +932,7 @@ namespace { if (gops.is_reflection_centric(h) && !centric_in_h) zone.push_back(&r); } - TwinImmuneZone z = ReadZone(zone, tncs_vector, nthreads); + TwinImmuneZone z = ReadZone(zone, tncs_vector, nthreads, twin_fraction); for (const auto& op : rotations) if (!contains(h_ops, op)) z.operators += (z.operators.empty() ? "" : " ") + op.as_hkl().triplet(); @@ -886,8 +943,10 @@ namespace { TwinImmuneZoneResult AnalyzeTwinImmuneZone(const std::vector& p1_merged, const gemmi::UnitCell& cell, const gemmi::SpaceGroup& group, const gemmi::SpaceGroup& subgroup, - const std::array* tncs_vector) { + const std::array* tncs_vector, size_t nthreads, + double twin_fraction) { TwinImmuneZoneResult result; + result.twin_fraction = twin_fraction; result.tncs_normalised = tncs_vector != nullptr; const gemmi::GroupOps gops = group.operations(); const std::vector rotations = ProperRotations(gops); @@ -895,16 +954,19 @@ TwinImmuneZoneResult AnalyzeTwinImmuneZone(const std::vector& const auto strong = StrongReflections(p1_merged, cell, gops, result); if (strong.empty()) return result; - result.control = ReadControl(strong, gops, tncs_vector, 1); - result.zones.push_back(ReadSubgroupZone(strong, gops, rotations, h_ops, tncs_vector, 1)); + result.control = ReadControl(strong, gops, tncs_vector, twin_fraction, nthreads); + result.zones.push_back(ReadSubgroupZone(strong, gops, rotations, h_ops, tncs_vector, twin_fraction, + nthreads)); Calibrate(result); return result; } TwinImmuneZoneResult AnalyzeTwinImmuneZones(const std::vector& p1_merged, const gemmi::UnitCell& cell, const gemmi::SpaceGroup& group, - const std::array* tncs_vector, size_t nthreads) { + const std::array* tncs_vector, size_t nthreads, + double twin_fraction) { TwinImmuneZoneResult result; + result.twin_fraction = twin_fraction; result.tncs_normalised = tncs_vector != nullptr; const gemmi::GroupOps gops = group.operations(); @@ -955,12 +1017,13 @@ TwinImmuneZoneResult AnalyzeTwinImmuneZones(const std::vector& const auto strong = StrongReflections(p1_merged, cell, gops, result); if (strong.empty()) return result; - result.control = ReadControl(strong, gops, tncs_vector, nthreads); + result.control = ReadControl(strong, gops, tncs_vector, twin_fraction, nthreads); // Per hypothesis H, the zone: centric in the group, acentric in H - the reflections some added // operator sends to -h, i.e. the self-mates of that coset read as a twin law. for (const auto& h_ops : subgroups) - result.zones.push_back(ReadSubgroupZone(strong, gops, rotations, h_ops, tncs_vector, nthreads)); + result.zones.push_back(ReadSubgroupZone(strong, gops, rotations, h_ops, tncs_vector, twin_fraction, + nthreads)); Calibrate(result); return result; } @@ -977,10 +1040,14 @@ std::string TwinImmuneZonesToText(const TwinImmuneZoneResult& result) { << " 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). 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"; + << " (> 0: the added operators are real), the centric side twinned at a fraction of " + << std::setprecision(2) << result.twin_fraction << " - what the\n" + << " lattice's other twin laws show by their own correlation, and reach a genuine zone with. The\n" + << " acentric control reads " << std::setprecision(3) << AcentricEvidencePerReflection(result.twin_fraction) + << " nats per reflection when the normalisation is right; what it reads above\n" + << " that is the normalisation's (" << std::showpos << result.control_excess_per_reflection << std::noshowpos + << " here), every zone carries the same per reflection, and the calibrated\n" + << " 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 diff --git a/image_analysis/scale_merge/TwinningAnalysis.h b/image_analysis/scale_merge/TwinningAnalysis.h index 54766ffa1..6ba931957 100644 --- a/image_analysis/scale_merge/TwinningAnalysis.h +++ b/image_analysis/scale_merge/TwinningAnalysis.h @@ -139,7 +139,8 @@ struct TwinImmuneZone { // Log-likelihood ratio, centric over acentric Wilson (each convolved with the reflection's error), // summed over the zone: > 0 favours the operators being real symmetry, < 0 a twin law or // 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. + // on error-free untwinned data; less either way as the errors grow. The centric density is the one + // twinned at the fraction the caller gives (TwinImmuneZoneResult::twin_fraction). 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. @@ -159,6 +160,8 @@ struct TwinImmuneZoneResult { // 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; + // The twin fraction the centric side of every zone was read at (AnalyzeTwinImmuneZones). + double twin_fraction = 0.0; // The control reads more compressed than a PERFECT twin's acentric population (<|E^2-1|> 0.541, by // more than three standard errors). Twinning cannot do that at any fraction, so something else // averages each reflection with unrelated ones - overlapping spots of a long cell, neighbour @@ -174,16 +177,25 @@ struct TwinImmuneZoneResult { bool tncs_normalised = false; // normalised within the pseudo-translation's two phase classes }; +// `twin_fraction` is what the centric hypothesis is read at. A twin law OTHER than the operators a zone +// is read for averages the zone's reflections with unrelated ones exactly as it does every other +// reflection, so where the crystal is twinned by such a law a genuine centric zone does not read the +// untwinned centric 0.968: at a fraction a it reads 0.88 (a = 0.1), 0.81 (0.2), 0.74 (0.5), and against +// the untwinned centric density it reads acentric - a genuine 321 crystal twinned by a 622 law read its +// genuine 2-folds at -206 nats. The operators' OWN twin law cannot reach their zone, which is the whole +// argument, so the acentric side stays untwinned whatever the fraction. The caller passes the fraction +// the lattice's other twin laws imply (TwinFractionOutsideGroup); 0 reads the untwinned densities. TwinImmuneZoneResult AnalyzeTwinImmuneZones(const std::vector& p1_merged, const gemmi::UnitCell& cell, const gemmi::SpaceGroup& group, const std::array* tncs_vector = nullptr, - size_t nthreads = 1); + size_t nthreads = 1, double twin_fraction = 0.0); // 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::UnitCell& cell, const gemmi::SpaceGroup& group, const gemmi::SpaceGroup& subgroup, - const std::array* tncs_vector = nullptr); + const std::array* tncs_vector = nullptr, + size_t nthreads = 1, double twin_fraction = 0.0); std::string TwinImmuneZonesToText(const TwinImmuneZoneResult& result); diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index c368b1c5c..21131003f 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -2127,9 +2127,11 @@ namespace { // The margin is the screw-axis convention (SearchSpaceGroupOptions::min_screw_absence_evidence): // the same log-likelihood-ratio units, and the same bound of 20 nats. Read over three batteries // the zones of genuine promotions sit at +60 nats and above and the zones behind the known - // 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. + // twins at -290 nats and below, with genuine trigonal crystals at -46 and -206 nats whose + // promotion no gate refused - twinned by a law other than the 2-folds, which reaches their zone. + // Read at the twin fraction those other laws show (twin_fraction_outside) the -206 one reads + // +209, and the search itself now asks the zone of a promotion on a merge that reads twinned + // (SearchSpaceGroup's twin-immune zone test); here it is read where a gate refused. // // 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 @@ -2148,7 +2150,8 @@ namespace { if (higher.point_group_order != 2 * lower.point_group_order) return 0; const auto z = AnalyzeTwinImmuneZone(p1_merged, cell, *higher.point_group_representative, - *lower.point_group_representative, tncs_vector); + *lower.point_group_representative, tncs_vector, 1, + higher.twin_fraction_outside); Logger logger("Rugnux"); if (z.zones.empty() || z.zones[0].n == 0) { logger.Info("Twin-immune zone {} over {}: no zone reflections in the shells with signal " @@ -2161,9 +2164,10 @@ namespace { : zone.calibrated_evidence_nats >= DECISIVE_NATS ? +1 : 0; const std::string text = fmt::format( "{} over {}: operators {}; n {} <|E^2-1|> {:.3f} +- {:.3f} (control {:.3f}, {:+.3f} nats per " - "reflection above acentric), evidence {:+.1f} nats, calibrated {:+.1f} ({:.2f}-{:.2f} A{}) -> {}", + "reflection above acentric), centric read twinned at {:.2f}, 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, z.control_excess_per_reflection, + zone.standard_error, z.control.mean_abs_e2_minus_1, z.control_excess_per_reflection, z.twin_fraction, 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 && z.zones_ambiguous @@ -9544,7 +9548,8 @@ ProcessResult Rugnux::RunPipeline(RugnuxObserver *observer, bool write_output, b config_.nthreads); result.twin_immune_zones = AnalyzeTwinImmuneZones( p1.merged, gemmi::UnitCell(*result.consensus_cell), *determined, - TranslationalNCSVectorForLTest(p1_tncs), config_.nthreads); + TranslationalNCSVectorForLTest(p1_tncs), config_.nthreads, + result.space_group_search ? result.space_group_search->twin_fraction_outside : 0.0); 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 3e64e6a24..ed51ded0b 100644 --- a/tests/TwinningAnalysisTest.cpp +++ b/tests/TwinningAnalysisTest.cpp @@ -227,14 +227,17 @@ namespace { // 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) { + // `twin_super` names the group whose extra operator is the twin law. + std::vector WilsonMerge(const char *true_group, const char *twin_super, + const gemmi::UnitCell &cell, 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::SpaceGroup &super = gemmi::get_spacegroup_by_name(twin_super); const gemmi::Op twin = jfjoch_test::TwinLaw(sub, super); const gemmi::GroupOps gops = sub.operations(); const gemmi::ReciprocalAsu rasu(&sub); - const gemmi::UnitCell cell(47, 47, 63, 90, 90, 90); + const int hmax = static_cast(cell.a / 2.5) + 1, kmax = static_cast(cell.b / 2.5) + 1, + lmax = static_cast(cell.c / 2.5) + 1; auto true_intensity = [&](const gemmi::Op::Miller &hkl) { const auto asu = rasu.to_asu(hkl, gops).first; const double u1 = jfjoch_test::detail::UniformFromHkl(asu); @@ -246,9 +249,9 @@ namespace { * std::exp(-anisotropy_b * s_c * s_c / 2.0); }; std::vector out; - for (int h = -19; h <= 19; ++h) - for (int k = -19; k <= 19; ++k) - for (int l = -26; l <= 26; ++l) { + for (int h = -hmax; h <= hmax; ++h) + for (int k = -kmax; k <= kmax; ++k) + for (int l = -lmax; l <= lmax; ++l) { if (std::make_tuple(h, k, l) <= std::make_tuple(-h, -k, -l)) continue; const gemmi::Op::Miller hkl{{h, k, l}}; @@ -271,6 +274,12 @@ namespace { } return out; } + + std::vector WilsonTetragonal(const char *true_group, double twin_fraction, + double anisotropy_b = 0.0, double scale_jitter = 0.0) { + return WilsonMerge(true_group, "P 4 2 2", gemmi::UnitCell(47, 47, 63, 90, 90, 90), twin_fraction, + anisotropy_b, scale_jitter); + } } // Reflections centric in the adopted group but acentric in a subgroup are their own twin mates under @@ -345,6 +354,28 @@ TEST_CASE("Twin-immune zones are normalised for anisotropy and calibrated by the } } +// A twin by a law OTHER than the operators a zone is read for averages the zone's reflections with +// unrelated ones like every other reflection, so a genuine 321 crystal twinned by a 622 law reads its +// 2-folds' zone acentric against the untwinned centric density. Read at the twin fraction, it reads +// centric again - while a crystal twinned by the zone's own operators, which cannot reach the zone, +// stays acentric. +TEST_CASE("Twin-immune zones are read at the twin fraction of the other twin laws", "[twinning]") { + const gemmi::SpaceGroup &p321 = gemmi::get_spacegroup_by_name("P 3 2 1"); + const gemmi::SpaceGroup &p3 = gemmi::get_spacegroup_by_name("P 3"); + const gemmi::UnitCell hexagonal(51, 51, 71, 90, 90, 120); + + const auto genuine = WilsonMerge("P 3 2 1", "P 6 2 2", hexagonal, 0.1); + const auto untwinned_read = AnalyzeTwinImmuneZone(genuine, hexagonal, p321, p3); + REQUIRE(untwinned_read.zones.size() == 1); + CHECK(untwinned_read.zones[0].calibrated_evidence_nats < -20.0); + const auto twinned_read = AnalyzeTwinImmuneZone(genuine, hexagonal, p321, p3, nullptr, 1, 0.1); + CHECK(twinned_read.zones[0].calibrated_evidence_nats > 20.0); + + // Twinned by the 2-folds themselves, with no other law to read the zone at. + const auto twin = WilsonMerge("P 3", "P 3 2 1", hexagonal, 0.1); + CHECK(AnalyzeTwinImmuneZone(twin, hexagonal, p321, p3).zones[0].calibrated_evidence_nats < -20.0); +} + // Averaging each intensity over its orbit under a group leaves the L-test where it was when the group // is the crystal's, and narrows it towards the perfect-twin 0.375 when an operator of the group is // not - read on the same pairs. A perfect twin of the subgroup already reads 0.375 unmerged, so there