From 08ae0fa05c3fd2a34dce8d1f447dda4e26e8c159 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 21:08:06 +0100 Subject: [PATCH 1/2] RESULT-CHANGING: rugnux first pass: a supercell of a sparse or twinned pattern comes down to the true cell On a sparse rotation sweep (a few real spots per frame among many stray ones) the FFT, whose floor is 10 A, returns a cell that is a multiple of the true one - doubled along every axis, or the coincidence lattice of two twin domains at 5 or 10 times the volume - and the frame counts that judged the alternatives favour it, because its extra nodes catch the stray spots that put frames over the bar. Measured on the public L-aspartic acid set (APEX II, P21 5.12/6.92/7.60, a ~50:50 non-merohedral twin, twofold about c*, domains coinciding on h even): every one of its seven runs merged in an 8x or 10x P1/C2 cell or not in the true one. - Sub-lattice hypothesis: decided by a paired test on the pooled validation spots - the spots on the standing cell alone against those on its index-2 sub-lattice alone, less what the extra nodes take at a wrong spindle angle (SubLatticeSuffices, McNemar). The sub-lattice is adopted unless the standing cell's extra nodes hold spots of their own; before, it had to beat the standing cell by 3.3 sigma on an unpaired comparison, which a sub-lattice of a supercell cannot do. Asked again of each adopted sub-lattice (index 8 at most). The nominee needs a majority of the frames OR as many frames as the standing cell (it needed both). - Twin-composite hypothesis: asked at every level of that descent, before the halving (two domains take fewer nodes than any halving of their composite, and a halving can split it); indices 2..10, every index-n sub-lattice as a form mod n (SubLattice gains the general case of a form with no coefficient prime to n); only congruent pairs are compared, congruence allowing the axis reversals that turn two angles into their supplements (a monoclinic domain's reduced cell comes out with beta 100 or 80 deg); the domain is adopted when it beats chance on the pooled spots, not on a majority of frames. A domain cell that is itself a multiple is halved on, its overlap form held. laspartic runs 01-07 (CPU), base 856ad5106 -> this: all seven P2 5.11-5.14/6.91-6.93/7.58-7.60 beta 99.7-100.5 (base: 10x P1 x2, 4x C2 x3, true cell x2); SHELXL R1 against COD 1577465 3.6-7.2% on all seven (base: no cell match x5, 17.5% and 19.4% where the cell was right but the twin overlap was merged in). The screw is not determined yet: the overlap exclusion takes 0k0. Tests: SubLatticeSuffices, TwinCompositeAt (synthetic monoclinic twin at index 2 and through a skewed index-10 supercell). Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 (cherry picked from commit f8544981ab647cb7c350f9d73c81d499e3d330b0) --- rugnux/Rugnux.cpp | 9 ++ rugnux/Rugnux.h | 26 +++ rugnux/RugnuxFirstPass.cpp | 260 +++++++++++++++++++++--------- rugnux/RugnuxFirstPass.h | 4 + rugnux/RugnuxFirstPassRescues.cpp | 96 ++++++----- tests/RugnuxTest.cpp | 82 ++++++++++ 6 files changed, 357 insertions(+), 120 deletions(-) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 90fa50985..d700c24ec 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -137,6 +137,15 @@ std::optional LongAxisChoice(const ValidationSpotEvidence &standing, return best; } +// Following McNemar (1947) Psychometrika 12, 153-157 +bool SubLatticeSuffices(int64_t only_standing, int64_t only_sub, const ValidationSpotEvidence &standing, + const ValidationSpotEvidence &sub) { + const double gain = static_cast(only_standing - only_sub) + - static_cast(std::max(0, standing.by_chance - sub.by_chance)); + const double sigma = std::sqrt(static_cast(std::max(1, only_standing + only_sub))); + return ValidationEvidenceBeatsChance(sub) && gain <= SPOT_BUDGET_SIGNIFICANCE_Z * sigma; +} + bool ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e) { return ValidationEvidenceBeatsChance(e, std::erfc(SPOT_BUDGET_SIGNIFICANCE_Z / std::sqrt(2.0))); } diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index d5defa3b7..0297c3c2f 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -367,6 +367,32 @@ struct LongAxisCandidate { std::optional LongAxisChoice(const ValidationSpotEvidence &standing, const std::vector &candidates); +// The twin-composite test at one index n (see Rugnux::FirstPassRun::TwinCompositeHypothesis): of the +// index-n sub-lattices of the standing cell `prim` - the kernels of the forms f.h = 0 (mod n) - the +// pair with congruent cells whose union holds the most of the indexed spots `hkl`, the stronger +// first, with the spots on each class and the per-node densities the verdict compares. None where no +// two index-n sub-lattices have congruent cells. +struct TwinCompositeEvidence { + std::array f1{}, f2{}; + int both = 0, only1 = 0, only2 = 0, neither = 0; + double p_both = 0.0; // the share of the nodes on both sub-lattices + double other_density = 0.0, neither_density = 0.0; + bool composite = false; // two domains on their coincidence lattice + CrystalLattice domain; // the stronger domain's cell, Niggli-reduced +}; +std::optional TwinCompositeAt(const CrystalLattice &prim, + const std::vector> &hkl, int n); + +// Whether a sub-lattice of the standing cell - the same orientation, fewer nodes - explains the +// pooled validation spots as well as the standing cell does. The two are compared on the same spots, so +// the comparison is paired: `only_standing` spots sit on the standing cell and not on the sub-lattice, +// `only_sub` the other way round. The standing cell's extra nodes hold spots of their own when the +// first exceeds the second, less what the extra nodes take at a wrong spindle angle, by more than the +// paired standard error (SPOT_BUDGET_SIGNIFICANCE_Z). Where they do not, the extra nodes are evidence +// of nothing, and the smaller cell, if it beats chance, is the hypothesis the spots support. +bool SubLatticeSuffices(int64_t only_standing, int64_t only_sub, const ValidationSpotEvidence &standing, + const ValidationSpotEvidence &sub); + // The goniometer rotation scale - the factor by which the stage turned relative to the angles stored // in the file - walked to the fit's fixed point and decided on the whole sweep. The fit only sees the // frames the angles it was measured at still track, and a stage at the wrong rate is exactly what diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp index c50c0514e..a836008a4 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -342,10 +342,116 @@ namespace { v[i] = t; return {v[0], v[1], v[2]}; } - return latt; + // No coefficient prime to n (n = 10, form (2,5,0)): the lattice the three axes and the + // translation span, reduced to a basis by integer row elimination, in units of an axis / n. + int m[4][3] = {{n, 0, 0}, {0, n, 0}, {0, 0, n}, {c[0], c[1], c[2]}}; + for (int j = 0; j < 3; j++) + for (bool done = false; !done;) { + int p = j; + for (int r = j; r < 4; r++) + if (m[r][j] != 0 && (m[p][j] == 0 || std::abs(m[r][j]) < std::abs(m[p][j]))) + p = r; + std::swap(m[j], m[p]); + done = true; + for (int r = j + 1; r < 4; r++) { + const int q = m[r][j] / m[j][j]; + for (int k = 0; k < 3; k++) + m[r][k] -= q * m[j][k]; + done = done && m[r][j] == 0; + } + } + Coord v[3]; + for (int i = 0; i < 3; i++) + v[i] = (e[0] * static_cast(m[i][0]) + e[1] * static_cast(m[i][1]) + + e[2] * static_cast(m[i][2])) / static_cast(n); + return {v[0], v[1], v[2]}; } } +std::optional TwinCompositeAt(const CrystalLattice &prim, + const std::vector> &hkl, int n) { + const auto mod = [n](int x) { return (x % n + n) % n; }; + // Each index-n sub-lattice once: the forms (u,v,w) mod n whose kernel has index n (no factor + // shared by all of u, v, w and n), each the least of its multiples by the units mod n. + std::vector> forms; + for (int u = 0; u < n; u++) + for (int v = 0; v < n; v++) + for (int w = 0; w < n; w++) { + if (std::gcd(std::gcd(std::gcd(u, v), w), n) != 1) + continue; + const std::array f = {u, v, w}; + bool least = true; + for (int k = 2; k < n && least; k++) + if (std::gcd(k, n) == 1) + least = f <= std::array{mod(k * u), mod(k * v), mod(k * w)}; + if (least) + forms.push_back(f); + } + const auto on = [&](const std::array &f, const std::array &m) { + return mod(f[0] * m[0] + f[1] * m[1] + f[2] * m[2]) == 0; + }; + const auto sub_of = [&](const std::array &f) { + AddedClassEvidence ev; + ev.u = f[0]; ev.v = f[1]; ev.w = f[2]; + return SubLattice(prim, ev, n).NiggliReduce(); + }; + std::vector cells; + for (const auto &f : forms) + cells.push_back(sub_of(f).GetUnitCell()); + // Reversing two axes of a cell turns two of its angles into their supplements and leaves the + // lattice as it was, and a reduced cell with angles near 90 deg can come out either way: one + // domain's monoclinic cell with beta 100 deg, the other's with 80. + const auto congruent = [](const UnitCell &u1, const UnitCell &u2) { + if (std::abs(u1.a - u2.a) >= 0.02 * u1.a || std::abs(u1.b - u2.b) >= 0.02 * u1.b + || std::abs(u1.c - u2.c) >= 0.02 * u1.c) + return false; + const auto close = [](double x, double y, bool supplement) { + return std::abs(x - (supplement ? 180.0 - y : y)) < 2.0; + }; + for (const auto &[sa, sb, sg] : {std::array{false, false, false}, {true, true, false}, + {true, false, true}, {false, true, true}}) + if (close(u1.alpha, u2.alpha, sa) && close(u1.beta, u2.beta, sb) && close(u1.gamma, u2.gamma, sg)) + return true; + return false; + }; + // The congruent pair whose union holds the most spots, the stronger domain first. + std::optional best; + for (size_t i = 0; i < forms.size(); i++) + for (size_t j = i + 1; j < forms.size(); j++) { + if (!congruent(cells[i], cells[j])) + continue; + TwinCompositeEvidence t{.f1 = forms[i], .f2 = forms[j]}; + for (const auto &m : hkl) { + const bool o1 = on(forms[i], m), o2 = on(forms[j], m); + (o1 && o2 ? t.both : o1 ? t.only1 : o2 ? t.only2 : t.neither)++; + } + if (!best || t.neither < best->neither) + best = t; + } + if (!best) + return best; + if (best->only2 > best->only1) { + std::swap(best->f1, best->f2); + std::swap(best->only1, best->only2); + } + best->domain = sub_of(best->f1); + // The share of the nodes on both sub-lattices, counted over one period of the forms: 1/n^2 + // for a prime n, more where the two forms share a factor with n. + int n_both = 0; + for (int h = 0; h < n; h++) + for (int k = 0; k < n; k++) + for (int l = 0; l < n; l++) + n_both += on(best->f1, {h, k, l}) && on(best->f2, {h, k, l}); + best->p_both = n_both / static_cast(n * n * n); + // Spots per node: each domain alone holds 1/n - p_both of the nodes, both p_both, neither + // 1 - 2/n + p_both. One crystal reads the emptiest of the first three about as dense as + // neither; two domains read it tens of times denser. + best->other_density = std::min(best->only2 / (1.0 / n - best->p_both), best->both / best->p_both); + best->neither_density = best->neither / (1.0 - 2.0 / n + best->p_both); + best->composite = best->only2 >= 20 && best->other_density > 4.0 * best->neither_density; + return best; +} + // The first pass of a rotation run, as the ordered list of its decisions: index the first part of the // sweep on two sampling schemes, put the standing lattice through the rescues and the hypotheses that // may replace it, accept or refuse it on the pooled validation spots, and leave it forced on the @@ -397,8 +503,7 @@ bool Rugnux::FirstPassRotationIndexing(PipelineLocals &p) { run.RequireLattice(best); best = run.ShortAxisHypothesis(std::move(best)); - best = run.SubLatticeHypothesis(std::move(best)); - best = run.TwinCompositeHypothesis(std::move(best)); + best = run.SubLatticeHypothesis(std::move(best)); // and the twin-composite hypothesis if (!p.geometry_prepass && prepass_result_.has_value()) best = run.Pass1Fallback(std::move(best)); // The last resort, and ONLY on a run that has no usable lattice: the pass-1 fallback above @@ -768,6 +873,39 @@ auto Rugnux::FirstPassRun::MeasurePooled(IndexAndRefine &idx, const RotationInde return e; } +// The paired form of PooledOnLattice, for two lattices compared on the same spots: how many of the +// frames' spots `a` puts on itself and `b` does not, and the other way round. Each frame's spots are +// fitted to each lattice as the pooled count fits them. +std::pair Rugnux::FirstPassRun::PairedOnLattice(IndexAndRefine &idx, const RotationIndexerResult &a, + const RotationIndexerResult &b, + const std::vector &frames) { + PrefetchSpots(frames); + const bool index_ice_rings = rx.experiment_.GetIndexingSettings().GetIndexIceRings(); + const auto on = [&](const RotationIndexerResult &r) { + idx.ForceRotationIndexerResult(r); + std::vector> flags(frames.size()); + ParallelFor(static_cast(frames.size()), + std::min(frames.size(), rx.config_.nthreads), [&](int i) { + DataMessage m{}; + m.number = frames[i]; + m.spots = spot_cache.at(frames[i]); + idx.CountSpotsOnLattice(m, validation_settings); + // As in PooledOnLattice: the spots the count took are the ones it gave a Miller index. + for (const auto &sp : m.spots) + flags[i].push_back((index_ice_rings || !sp.ice_ring) && (sp.h != 0 || sp.k != 0 || sp.l != 0)); + }); + return flags; + }; + const auto on_a = on(a), on_b = on(b); + int64_t only_a = 0, only_b = 0; + for (size_t i = 0; i < frames.size(); i++) + for (size_t k = 0; k < std::min(on_a[i].size(), on_b[i].size()); k++) { + only_a += on_a[i][k] && !on_b[i][k]; + only_b += on_b[i][k] && !on_a[i][k]; + } + return {only_a, only_b}; +} + // How deep into an image's intensity-ordered spot list this lattice is still being seen - the // measured spot budget. Same frames, same per-image path as CountIndexed above, but scoring the // spots rather than the frames: each frame's spots are tallied by rank and the run keeps the @@ -1572,17 +1710,25 @@ auto Rugnux::FirstPassRun::IndexSchemes() -> FirstPass { // domain's own cell, and it indexes more frames and more spots than a domain's own cell does - so // neither the frame count nor the occupancy of the classes it adds tells it from a superstructure. // What does is the class of nodes on neither domain. Two domains put every reflection on the union of -// their own index-n sub-lattices K1 and K2 and leave the (1 - 1/n)^2 of the nodes on neither empty but -// for chance coincidences, while the nodes on both (the overlap) and on each alone are populated. One -// crystal cannot do that: whatever sub-lattice its reflections favour (a pseudo-translation), the -// nodes off K1 that lie on K2 and those that lie on neither are the same cosets of it and as dense in -// spots, and a crystal whose spots avoid one coset of a sub-lattice thins the overlap with the class -// on neither alike. So the class on neither must be far emptier than every other class. The two -// sub-lattices are the same lattice in two orientations, so their cells must also agree. Where all -// of this holds the stronger domain's own cell is adopted, and the other domain is left to the -// leftover-lattice census, which names it. +// their own index-n sub-lattices K1 and K2 and leave the nodes on neither empty but for chance +// coincidences, while the nodes on both (the overlap) and on each alone are populated. One crystal +// cannot do that: whatever sub-lattice its reflections favour (a pseudo-translation), the nodes off K1 +// that lie on K2 and those that lie on neither are the same cosets of it and as dense in spots, and a +// crystal whose spots avoid one coset of a sub-lattice thins the overlap with the class on neither +// alike. So the class on neither must be far emptier than every other class. The two sub-lattices are +// the same lattice in two orientations, so their cells must also agree, and only such pairs are +// compared. Where all of this holds the stronger domain's own cell is adopted, and the other domain is +// left to the leftover-lattice census, which names it. +// +// The index is not only 2 or 3. A twin by reticular pseudo-merohedry - a monoclinic crystal twinned +// about a reciprocal axis that is nearly a lattice row, the domains coinciding on alternate layers - +// has no exact coincidence lattice, and the FFT then returns whatever multiple of the cell absorbs +// both domains' spots best, measured 5 and 10 times the domain's volume. Every index-n sub-lattice of +// the standing cell is the kernel of a linear form mod n when n has no square factor, so the forms +// mod n cover those indices; n runs to TWIN_COMPOSITE_MAX_INDEX. // Following Friedel's twin index (Friedel 1926, Lecons de cristallographie) and the obverse/reverse -// reflection classes of Herbst-Irmer & Sheldrick (2002) Acta Cryst. B58, 477-481 +// reflection classes of Herbst-Irmer & Sheldrick (2002) Acta Cryst. B58, 477-481; twinning by +// reticular pseudo-merohedry after Herbst-Irmer & Sheldrick (1998) Acta Cryst. B54, 443-449 auto Rugnux::FirstPassRun::TwinCompositeHypothesis(FirstPass best) -> FirstPass { if (rx.cancelled_ || !best.result.has_value() || best.score <= 0) return best; @@ -1591,70 +1737,21 @@ auto Rugnux::FirstPassRun::TwinCompositeHypothesis(FirstPass best) -> FirstPass const auto hkl = IndexedMillerIndices(prim, ValidationCloud(*best.result), index_tol); if (hkl.empty()) return best; - for (const int n : {2, 3}) { - // Each index-n sub-lattice once: the form (u,v,w) mod n with its first non-zero entry 1. - std::vector> forms; - for (int u = 0; u < n; u++) - for (int v = 0; v < n; v++) - for (int w = 0; w < n; w++) - if ((u ? u : v ? v : w) == 1) - forms.push_back({u, v, w}); - const auto on = [n](const std::array &f, const std::array &m) { - return ((f[0] * m[0] + f[1] * m[1] + f[2] * m[2]) % n + n) % n == 0; - }; - struct Classes { int both = 0, only1 = 0, only2 = 0, neither = 0; }; - const auto classes = [&](size_t i, size_t j) { - Classes c; - for (const auto &m : hkl) { - const bool o1 = on(forms[i], m), o2 = on(forms[j], m); - (o1 && o2 ? c.both : o1 ? c.only1 : o2 ? c.only2 : c.neither)++; - } - return c; - }; - // The pair whose union holds the most spots, the stronger domain first. - size_t bi = 0, bj = 1; - Classes c = classes(0, 1); - for (size_t i = 0; i < forms.size(); i++) - for (size_t j = i + 1; j < forms.size(); j++) { - const Classes t = classes(i, j); - if (t.neither < c.neither) { - c = t; bi = i; bj = j; - } - } - if (c.only2 > c.only1) { - std::swap(bi, bj); - std::swap(c.only1, c.only2); - } - const auto sub_of = [&](size_t f) { - AddedClassEvidence ev; - ev.u = forms[f][0]; ev.v = forms[f][1]; ev.w = forms[f][2]; - return SubLattice(prim, ev, n).NiggliReduce(); - }; - const CrystalLattice sub = sub_of(bi); - const UnitCell u1 = sub.GetUnitCell(), u2 = sub_of(bj).GetUnitCell(); - const bool congruent = std::abs(u1.a - u2.a) < 0.02 * u1.a && std::abs(u1.b - u2.b) < 0.02 * u1.b - && std::abs(u1.c - u2.c) < 0.02 * u1.c && std::abs(u1.alpha - u2.alpha) < 2.0 - && std::abs(u1.beta - u2.beta) < 2.0 && std::abs(u1.gamma - u2.gamma) < 2.0; - // Spots per node: each domain alone holds 1/n - 1/n^2 of the nodes, both 1/n^2, neither - // (1 - 1/n)^2. One crystal reads the emptiest of the first three about as dense as neither (up - // to 1.5 times, the scatter of picking the emptiest pair out of 7 or 13 forms); two domains read - // it tens of times denser. - const double other_density = std::min(c.only2 / (1.0 / n - 1.0 / (n * n)), - c.both * static_cast(n * n)); - const double neither_density = c.neither / ((1.0 - 1.0 / n) * (1.0 - 1.0 / n)); - const bool composite = congruent && c.only2 >= 20 && other_density > 4.0 * neither_density; + for (int n = 2; n <= TWIN_COMPOSITE_MAX_INDEX; n++) { + const auto tc = TwinCompositeAt(prim, hkl, n); + if (!tc) + continue; logger.Info("Twin-composite hypothesis, index {}: of {} validation spots on the standing cell " - "({:.0f} A^3), the sub-lattices ({},{},{}) and ({},{},{}) hold {} and {} alone, {} " + "({:.0f} A^3), the congruent sub-lattices ({},{},{}) and ({},{},{}) hold {} and {} alone, {} " "together and {} on neither - per node {:.1f} on the emptiest of the other classes against " - "{:.1f} on neither{} => {}", n, hkl.size(), best.vol, - forms[bi][0], forms[bi][1], forms[bi][2], forms[bj][0], forms[bj][1], forms[bj][2], - c.only1, c.only2, c.both, c.neither, other_density, neither_density, - congruent ? "" : ", cells not congruent", - composite ? "two domains on their coincidence lattice" : "one lattice"); - if (!composite) + "{:.1f} on neither => {}", n, hkl.size(), best.vol, + tc->f1[0], tc->f1[1], tc->f1[2], tc->f2[0], tc->f2[1], tc->f2[2], tc->only1, tc->only2, + tc->both, tc->neither, tc->other_density, tc->neither_density, + tc->composite ? "two domains on their coincidence lattice" : "one lattice"); + if (!tc->composite) continue; RotationIndexerResult cand = *best.result; - cand.search_result = LatticeSearch(sub); + cand.search_result = LatticeSearch(tc->domain); cand.lattice = cand.search_result.conventional; cand.unconstrained.reset(); cand.extra_lattices.clear(); @@ -1662,23 +1759,28 @@ auto Rugnux::FirstPassRun::TwinCompositeHypothesis(FirstPass best) -> FirstPass // stronger domain's reflections h on the weaker's sub-lattice are those with h . t divisible // by n, t the weaker sub-lattice's form taken as a vector of the standing cell. cand.twin_overlap = RotationIndexerResult::TwinOverlap{ - .t = forms[bj][0] * prim.Vec0() + forms[bj][1] * prim.Vec1() + forms[bj][2] * prim.Vec2(), + .t = tc->f2[0] * prim.Vec0() + tc->f2[1] * prim.Vec1() + tc->f2[2] * prim.Vec2(), .n = n}; const int score = CountIndexed(*indexer, cand); + // The composite is the evidence that the domain's cell is this crystal's; what is left to ask + // of the cell is that it is a lattice at all - that the validation spots it holds beat what a + // wrong spindle angle puts on it. Not a frame count: on a sparse pattern with two domains in + // the beam a frame holds a handful of each domain's spots, and a bar on frames is a bar on + // how sparse the pattern is. + const bool adopt = score > 0 && ValidationEvidenceBeatsChance(MeasurePooled(*indexer, cand)); const auto uc = cand.lattice.GetUnitCell(); - const bool adopt = score > static_cast(validation.size()) / 2; logger.Info("Twin-composite hypothesis: the stronger domain's cell, {}-centred {} {:.2f} {:.2f} " - "{:.2f} {:.1f} {:.1f} {:.1f}, indexes {}/{} validation frames - {} (the 1/{} of its " + "{:.2f} {:.1f} {:.1f} {:.1f}, indexes {}/{} validation frames - {} (the {:.0f}% of its " "reflections the other domain overlaps, h.({:.2f},{:.2f},{:.2f}) divisible by {}, is left out " "of integration)", cand.search_result.centering, gemmi::crystal_system_str(cand.search_result.system), uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, score, static_cast(validation.size()), adopt ? "adopting it" : "keeping the standing cell", - n, cand.twin_overlap->t * cand.lattice.Astar(), cand.twin_overlap->t * cand.lattice.Bstar(), - cand.twin_overlap->t * cand.lattice.Cstar(), n); + 100.0 * tc->p_both * n, cand.twin_overlap->t * cand.lattice.Astar(), + cand.twin_overlap->t * cand.lattice.Bstar(), cand.twin_overlap->t * cand.lattice.Cstar(), n); if (adopt) { best = FirstPass{std::move(cand), score, fmt::format("{} (twin domain)", best.name), - std::abs(sub.CalcVolume())}; + std::abs(tc->domain.CalcVolume())}; break; } } diff --git a/rugnux/RugnuxFirstPass.h b/rugnux/RugnuxFirstPass.h index 001277b8b..57ddf3be0 100644 --- a/rugnux/RugnuxFirstPass.h +++ b/rugnux/RugnuxFirstPass.h @@ -133,6 +133,8 @@ struct Rugnux::FirstPassRun { const std::vector &frames); PooledEvidence MeasurePooled(IndexAndRefine &idx, const RotationIndexerResult &r); PooledEvidence MeasurePooled(IndexAndRefine &idx, const RotationIndexerResult &r, const std::vector &frames); + std::pair PairedOnLattice(IndexAndRefine &idx, const RotationIndexerResult &a, + const RotationIndexerResult &b, const std::vector &frames); int64_t MeasureBudget(IndexAndRefine &idx, const RotationIndexerResult &r); std::unique_ptr FeedScheme(IndexerThreadPool &pool, const std::vector &ordinals, const DiffractionExperiment *x = nullptr); @@ -157,6 +159,8 @@ struct Rugnux::FirstPassRun { // and for two pooled spot shares the binomial standard error ValidationEvidencePrefers uses. // Each returns whether the decision is a near tie. static constexpr double NEAR_TIE_SIGMA = 2.0; + // The largest index of a twin's coincidence lattice the twin-composite hypothesis asks about. + static constexpr int TWIN_COMPOSITE_MAX_INDEX = 10; std::string PassTag() const; // A frame count as the log gives it: "k/N", or "no lattice" for a pass that found none. std::string Frames(int score) const; diff --git a/rugnux/RugnuxFirstPassRescues.cpp b/rugnux/RugnuxFirstPassRescues.cpp index 238e6867c..a9e7fbd43 100644 --- a/rugnux/RugnuxFirstPassRescues.cpp +++ b/rugnux/RugnuxFirstPassRescues.cpp @@ -967,19 +967,31 @@ auto Rugnux::FirstPassRun::ShortAxisHypothesis(FirstPass best) -> FirstPass { // // The frame count only nominates: a coarser lattice asks less accuracy of every spot, so a // sub-lattice of the true cell can index as many frames as the cell itself, or a few more -// where the cell's own score is depressed. What decides is the comparison this pass already -// makes between a cell and its axis harmonic: the share of the validation frames' pooled -// spots each puts on itself over and above what a wrong spindle angle puts there, each -// against its own null (ValidationEvidencePrefers). A real superstructure keeps its cell on -// that count, because the class it adds holds reflections; a cell built out of two crystals -// loses it, because it has to place every spot of the one it describes twice as accurately -// with an orientation fitted to both. The nominee must also index a majority of the frames, -// the bar the false-promotion test above asks of a cell that "works", and no fewer frames -// than the standing cell: a standing cell that indexes almost nothing is not this crystal's -// lattice, its sub-lattices are not hypotheses about it, and the rescues below are what -// such a run needs. +// where the cell's own score is depressed, and on a sparse pattern a supercell can index more +// frames than the true cell because its extra nodes catch the stray spots. What decides is +// whether the standing cell's extra nodes - the half of them off the sub-lattice - hold spots +// of their own. The two lattices are counted on the same pooled validation spots, so the test +// is paired: the spots on the standing cell alone against those on the sub-lattice alone, less +// what the extra nodes take at a wrong spindle angle (SubLatticeSuffices). A real +// superstructure keeps its cell, because the class it adds holds reflections; a cell whose +// extra nodes hold nothing beyond that is explained as well by half the nodes, and the smaller +// cell is the one the spots support. Asked again of an adopted sub-lattice, so a cell doubled +// along every axis comes down to the true one a halving at a time. The nominee must also index +// a majority of the frames, the bar the false-promotion test above asks of a cell that +// "works", or as many as the standing cell does: a standing cell whose sub-lattices index almost +// nothing is not this crystal's lattice, and the rescues below are what such a run needs. auto Rugnux::FirstPassRun::SubLatticeHypothesis(FirstPass best) -> FirstPass { - if (!rx.cancelled_ && best.result.has_value() && best.score > 0) { + // Asked again of each sub-lattice it adopts: a cell doubled along two or three axes is a sub-lattice + // of index 4 or 8, reached one halving at a time. Each cell on the way is first asked whether it is + // the coincidence lattice of two twin domains: two domains of a cell take fewer nodes than any + // halving of their composite, and a halving can split the composite and lose one domain. A domain's + // cell can itself be a multiple of the true one, and is halved on as any other; the overlap it + // carries is a vector of every finer lattice in its orientation, so it holds through the halving. + for (int level = 0;; level++) { + if (best.result.has_value() && !best.result->twin_overlap.has_value()) + best = TwinCompositeHypothesis(std::move(best)); + if (level == 3 || rx.cancelled_ || !best.result.has_value() || best.score <= 0) + break; const CrystalLattice prim = best.result->lattice .ToPrimitive(best.result->search_result.centering); const Coord A = prim.Vec0(), B = prim.Vec1(), C = prim.Vec2(); @@ -1010,52 +1022,54 @@ auto Rugnux::FirstPassRun::SubLatticeHypothesis(FirstPass best) -> FirstPass { std::abs(sub.CalcVolume())}; } bool adopt = false; - std::string verdict = " - none indexes a majority of the frames and as many as the standing cell"; + std::string verdict = " - none indexes a majority of the frames or as many as the standing cell"; const int majority = static_cast(validation.size()) / 2; - bool nominated = sub_best && sub_best->score >= best.score && sub_best->score > majority; - // Carried only where the gate is what decides: a sub-lattice that indexes fewer frames than - // the standing cell is not nominated whatever the gate says. - if (sub_best && NearTieGate("sub-lattice majority gate", "best sub-lattice", sub_best->score, majority) - && sub_best->score >= best.score && DecideLate()) { + const bool as_many = sub_best && sub_best->score >= best.score; + bool nominated = sub_best && (as_many || sub_best->score > majority); + if (sub_best && !as_many + && NearTieGate("sub-lattice majority gate", "best sub-lattice", sub_best->score, majority) + && DecideLate()) { const int late = CountIndexed(*indexer, *sub_best->result, late_validation); - const int standing_late = CountIndexed(*indexer, *best.result, late_validation); const bool early = nominated; - nominated = sub_best->score + late >= best.score + standing_late - && LateOverGate(sub_best->score, late, majority); + nominated = LateOverGate(sub_best->score, late, majority); LogDecidedLate("sub-lattice majority gate", - fmt::format("best sub-lattice {}/{} of the later frames, the standing cell {}/{}", late, - late_validation.size(), standing_late, late_validation.size()), + fmt::format("best sub-lattice {}/{} of the later frames", late, late_validation.size()), early, nominated); } if (nominated) { - PooledEvidence at_standing = MeasurePooled(*indexer, *best.result); - PooledEvidence at_sub = MeasurePooled(*indexer, *sub_best->result); - if (NearTiePooled("sub-lattice hypothesis", "standing cell", at_standing, "sub-lattice", at_sub) + const RotationIndexerResult &standing = *best.result, &sub = *sub_best->result; + PooledEvidence at_standing = MeasurePooled(*indexer, standing); + PooledEvidence at_sub = MeasurePooled(*indexer, sub); + auto [only_standing, only_sub] = PairedOnLattice(*indexer, standing, sub, validation); + if (NearTie("sub-lattice hypothesis", fmt::format("standing cell alone {}", only_standing), + fmt::format("sub-lattice alone {}", only_sub), static_cast(only_standing - only_sub), + std::sqrt(static_cast(std::max(1, only_standing + only_sub))), + "paired on pooled validation spots") && DecideLate()) { - const bool early = ValidationEvidencePrefers(at_standing, at_sub); + const bool early = SubLatticeSuffices(only_standing, only_sub, at_standing, at_sub); const std::vector all = AllValidationFrames(); - at_standing = MeasurePooled(*indexer, *best.result, all); - at_sub = MeasurePooled(*indexer, *sub_best->result, all); - LogDecidedLate("sub-lattice hypothesis", "the pooled spots counted again on all the frames", - early, ValidationEvidencePrefers(at_standing, at_sub)); + at_standing = MeasurePooled(*indexer, standing, all); + at_sub = MeasurePooled(*indexer, sub, all); + std::tie(only_standing, only_sub) = PairedOnLattice(*indexer, standing, sub, all); + LogDecidedLate("sub-lattice hypothesis", "the pooled spots counted again on all the frames", early, + SubLatticeSuffices(only_standing, only_sub, at_standing, at_sub)); } - adopt = ValidationEvidencePrefers(at_standing, at_sub); - verdict = fmt::format(" - the best puts {}/{} pooled spots on itself ({} by chance) against " - "the standing cell's {}/{} ({}), {}", - at_sub.on_lattice, at_sub.spots, at_sub.by_chance, - at_standing.on_lattice, at_standing.spots, at_standing.by_chance, - adopt ? "adopting it" : "keeping the standing cell"); + adopt = SubLatticeSuffices(only_standing, only_sub, at_standing, at_sub); + verdict = fmt::format(" - the best puts {}/{} pooled spots on itself ({} by chance) against the standing " + "cell's {}/{} ({}); {} spots on the standing cell alone, {} on the sub-lattice " + "alone, {}", at_sub.on_lattice, at_sub.spots, at_sub.by_chance, + at_standing.on_lattice, at_standing.spots, at_standing.by_chance, only_standing, + only_sub, adopt ? "adopting it" : "keeping the standing cell"); } logger.Info("Sub-lattice hypothesis: the standing cell ({:.0f} A^3) indexes {} validation " "frames; its index-2 sub-lattices (half-translation u,v,w):{}{}", best.vol, Frames(best.score), counts, verdict); - if (adopt) - best = std::move(*sub_best); + if (!adopt) + break; + best = std::move(*sub_best); } return best; } - -// Pass 1's lattice refit at the rotation scale this pass runs at (see Pass1Fallback): cell and // orientation, at pass 1's detector geometry and with the detector and the axis held, on the spots of // the validation frames of the whole sweep. auto Rugnux::FirstPassRun::RefitAtScale(RotationIndexerResult pass1_lattice) -> RotationIndexerResult { diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 7114e22f3..e4d2024ac 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -752,3 +752,85 @@ TEST_CASE("LongAxisChoice_PooledSpotsNotFrames", "[LongAxis]") { CHECK_FALSE(LongAxisChoice({}, {{ev(50, 60), false}}).has_value()); CHECK_FALSE(LongAxisChoice(standing, {}).has_value()); } + +// The sub-lattice decision is paired: the spots on the standing cell alone against those on its +// sub-lattice alone, less what the standing cell's extra nodes take at a wrong spindle angle. +TEST_CASE("SubLatticeSuffices", "[process]") { + // A sub-lattice that beats chance by far. + const ValidationSpotEvidence sub{.spots = 1400, .on_lattice = 340, .by_chance = 5, .reflections = 300, + .reflections_by_chance = 30}; + ValidationSpotEvidence standing = sub; + standing.on_lattice = 300; + // The standing cell's extra nodes hold fewer spots than the sub-lattice's better fit wins: halve. + CHECK(SubLatticeSuffices(8, 52, standing, sub)); + // A few more on the standing cell, within the paired noise (gain 10, sqrt(130) = 11.4): halve. + CHECK(SubLatticeSuffices(70, 60, standing, sub)); + // A class of its own (a superstructure, a pseudo-translation): keep the standing cell. + CHECK_FALSE(SubLatticeSuffices(120, 20, standing, sub)); + // What the extra nodes take by chance is not evidence for them. + standing.by_chance = 30; + CHECK(SubLatticeSuffices(80, 30, standing, sub)); + // A sub-lattice that does not beat chance is no hypothesis. + const ValidationSpotEvidence chance{.spots = 1400, .on_lattice = 10, .by_chance = 9, .reflections = 10, + .reflections_by_chance = 45}; + CHECK_FALSE(SubLatticeSuffices(0, 0, standing, chance)); +} + +// Two domains of a monoclinic cell twinned by a twofold about c*. With c cos(beta) = -a/4 the twin +// takes c to c + a/2, so the domains coincide on the reflections with h even and their coincidence +// lattice doubles c; any supercell of that one also indexes both. A synthetic cell. +TEST_CASE("TwinCompositeAt", "[process]") { + const float beta = std::acos(-1.0f / 6.0f) * 180.0f / static_cast(M_PI); + const CrystalLattice domain(6.0f, 8.0f, 9.0f, 90.0f, beta, 90.0f); + const Coord a = domain.Vec0(), b = domain.Vec1(), c = domain.Vec2(); + const Coord u = domain.Cstar() * (1.0f / domain.Cstar().Length()); + const auto twin = [&](const Coord &q) { return u * (2.0f * (q * u)) - q; }; + std::vector cloud, single; + for (int h = -4; h <= 4; h++) + for (int k = -4; k <= 4; k++) + for (int l = -4; l <= 4; l++) { + if (h == 0 && k == 0 && l == 0) + continue; + const Coord q = domain.Astar() * h + domain.Bstar() * k + domain.Cstar() * l; + if (q.Length() > 1.0f) + continue; + cloud.push_back(q); + cloud.push_back(twin(q)); + single.push_back(q); + } + const auto indices = [](const CrystalLattice &s, const std::vector &qs) { + std::vector> hkl; + for (const auto &q : qs) + hkl.push_back({static_cast(std::lround(s.Vec0() * q)), static_cast(std::lround(s.Vec1() * q)), + static_cast(std::lround(s.Vec2() * q))}); + return hkl; + }; + const auto same_cell = [](const CrystalLattice &x, const CrystalLattice &y) { + const auto ux = x.GetUnitCell(), uy = y.GetUnitCell(); + return std::abs(ux.a - uy.a) < 0.01f && std::abs(ux.b - uy.b) < 0.01f && std::abs(ux.c - uy.c) < 0.01f + && std::abs(std::abs(std::cos(ux.beta * M_PI / 180.0)) - std::abs(std::cos(uy.beta * M_PI / 180.0))) < 1e-3; + }; + const CrystalLattice reduced = domain.NiggliReduce(); + + // The coincidence lattice, index 2. + const CrystalLattice doubled = CrystalLattice(a, b, c * 2.0f).NiggliReduce(); + const auto at2 = TwinCompositeAt(doubled, indices(doubled, cloud), 2); + REQUIRE(at2.has_value()); + CHECK(at2->composite); + CHECK(at2->neither == 0); + CHECK(at2->p_both == Catch::Approx(0.25)); + CHECK(same_cell(at2->domain, reduced)); + + // A skewed index-5 supercell of it: the domains are index-10 sub-lattices, among them forms with no + // coefficient prime to 10. + const CrystalLattice skewed = CrystalLattice(a * 2.0f + b, b * 2.0f - a, c * 2.0f).NiggliReduce(); + const auto at10 = TwinCompositeAt(skewed, indices(skewed, cloud), 10); + REQUIRE(at10.has_value()); + CHECK(at10->composite); + CHECK(at10->neither == 0); + CHECK(same_cell(at10->domain, reduced)); + + // One domain alone on the doubled cell: one lattice, not a composite. + const auto alone = TwinCompositeAt(doubled, indices(doubled, single), 2); + CHECK(!(alone.has_value() && alone->composite)); +} From f1f663a33e52c3ce9631f5d8b87729f0e5292035 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Wed, 28 Oct 2026 22:16:36 +0100 Subject: [PATCH 2/2] RESULT-CHANGING: twin-composite domain: axial reflections the twin maps onto themselves are kept, on the domain's scale The twin-composite hypothesis leaves every reflection the other domain overlaps out of integration. On a monoclinic crystal twinned about c* (domains coinciding on h even) that takes the whole 0k0 row, so the screw along b could not be read and a P21 crystal was written as P2 with the screw undetermined. Of the overlapped reflections, those along a cell axis perpendicular to the other two whose twin partner is the reflection itself or its Friedel mate coincide with it EXACTLY - the twin takes that axis onto itself or its reverse - and are kept. Elsewhere the coincidence of a twin by reticular pseudo-merohedry is only approximate (the partner's node sits beside the reflection by an offset that grows with the index): keeping every overlapped reflection whose partner rounds to itself or its Friedel mate was tried and raised SHELXL R1 on L-aspartic acid from 3.6-7.2% to 5.7-14.1%. A kept reflection holds the whole crystal's intensity and every other one its domain's share, so it is scaled by that share: 1 - the twin fraction, the other domain's share of the summed intensity of the spots on one domain alone. Unscaled, the three 0k0 reflections of run 02 alone raised R1 from 3.6 to 4.6%. The twin rotation is the rotation between the two domains' reduced bases (TwinCompositeEvidence:: rotation, carried in RotationIndexerResult::TwinOverlap); the partner's index is read through it. laspartic runs 01-07 (CPU), against f8544981a: P21 on 01/02/03/04/06 (was P2, screw undetermined); 05 and 07 stay P2 - their 0k0 row holds too few reflections or none (rugnux's own warning). SHELXL R1 vs COD 1577465: 6.24/3.92/4.50/6.48/7.20/5.13/5.31% (was 5.48/3.63/4.28/5.65/7.20/5.22/5.31%). Test: TwinCompositeAt checks the twin rotation (a twofold reversing b) and the twin fraction. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 (cherry picked from commit 15bd8c39a8726c00daee9e4ec9b9e6f0d42e0124) --- image_analysis/IndexAndRefine.cpp | 41 ++++++++++- image_analysis/IndexAndRefine.h | 7 ++ .../rotation_indexer/RotationIndexer.h | 7 +- rugnux/Rugnux.h | 10 ++- rugnux/RugnuxFirstPass.cpp | 73 ++++++++++++++++--- tests/RugnuxTest.cpp | 15 ++++ 6 files changed, 141 insertions(+), 12 deletions(-) diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index eb04815f4..45f1632d8 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -167,6 +167,20 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetryRotat static_cast(std::lround(6.0f * (tw->t * l.Bstar()))), static_cast(std::lround(6.0f * (tw->t * l.Cstar())))}; outcome.twin_overlap_n = 6 * tw->n; + // The other domain's reflection at the node of h, read in this domain's basis: the index of + // R^-1 q, q the node, R the rotation from this domain to the other. + if (tw->rotation) { + const Coord a[3] = {l.Vec0(), l.Vec1(), l.Vec2()}, as[3] = {l.Astar(), l.Bstar(), l.Cstar()}; + std::array p{}; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + p[3 * i + j] = (*tw->rotation * a[i]) * as[j]; + outcome.twin_partner = p; + outcome.twin_domain_share = 1.0f - tw->fraction; + const auto right = [](const Coord &x, const Coord &y) { return std::abs(angle_deg(x, y) - 90.0f) < 0.5f; }; + for (int i = 0; i < 3; i++) + outcome.twin_axis_kept[i] = right(a[i], a[(i + 1) % 3]) && right(a[i], a[(i + 2) % 3]); + } } } return outcome; @@ -663,9 +677,34 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, std::vector integrated = integrate(prediction.GetReflections(), nrefl, msg.number); if (const int n = outcome.twin_overlap_n; n > 0) { const auto &w = outcome.twin_overlap_w; + const auto &p = outcome.twin_partner; + // Kept only where the coincidence is exact: a twin by reticular PSEUDO-merohedry puts the + // other domain's node beside this one's rather than on it, by an offset that grows with the + // index, so a partner that rounds to the reflection itself proves nothing in general. Along an + // axis perpendicular to the other two it does - the twin takes that axis onto itself or its + // reverse exactly - and those are the rows a screw is read on. + const auto exact = [&](const Reflection &r) { + const int i_axis = r.h != 0 ? 0 : r.k != 0 ? 1 : 2; + const bool axial = (r.h != 0) + (r.k != 0) + (r.l != 0) == 1; + if (!p || !axial || !outcome.twin_axis_kept[i_axis]) + return false; + const auto row = [&](int i) { + return std::lround((*p)[3 * i] * r.h + (*p)[3 * i + 1] * r.k + (*p)[3 * i + 2] * r.l); + }; + const long h = row(0), k = row(1), l = row(2); + return (h == r.h && k == r.k && l == r.l) || (h == -r.h && k == -r.k && l == -r.l); + }; std::erase_if(integrated, [&](const Reflection &r) { - return (r.h * w[0] + r.k * w[1] + r.l * w[2]) % n == 0; + return (r.h * w[0] + r.k * w[1] + r.l * w[2]) % n == 0 && !exact(r); }); + // Such a reflection holds the whole crystal's intensity, every other one this domain's share of + // it: scaled by that share, it is on the others' scale. + for (auto &r : integrated) + if ((r.h * w[0] + r.k * w[1] + r.l * w[2]) % n == 0) { + r.I *= outcome.twin_domain_share; + r.sigma *= outcome.twin_domain_share; + r.var_bkg *= outcome.twin_domain_share * outcome.twin_domain_share; + } } msg.integrated_reflections = integrated.size(); auto integration_end_time = std::chrono::steady_clock::now(); diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 5eaae7f43..75787533e 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -131,6 +131,13 @@ class IndexAndRefine { // overlap's index) are left out of integration. std::array twin_overlap_w{}; int twin_overlap_n = 0; + // Where the twin rotation is known: the twin partner of reflection h, in the candidate's basis, + // is the nearest integer to twin_partner * h (row-major 3x3); twin_axis_kept[i] where axis i is + // perpendicular to the other two, so that a reflection along it whose partner is itself or its + // Friedel mate coincides with that partner exactly. + std::optional> twin_partner; + std::array twin_axis_kept{}; + float twin_domain_share = 1.0f; // this domain's share of the crystal (1 - the twin fraction) explicit IndexingOutcome(const DiffractionExperiment& experiment_ref) : experiment(experiment_ref) {} diff --git a/image_analysis/rotation_indexer/RotationIndexer.h b/image_analysis/rotation_indexer/RotationIndexer.h index e17e22315..5825f5606 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.h +++ b/image_analysis/rotation_indexer/RotationIndexer.h @@ -40,10 +40,15 @@ struct RotationIndexerResult { // Twinning by reticular merohedry (rugnux's twin-composite step): a second domain whose lattice // shares the reflections h of this one with h . t divisible by n, t a vector of this lattice (real // space, in the frame of `lattice`, so it holds through any change of basis). Those reflections - // carry both domains' intensity, and are not integrated. + // carry both domains' intensity, and are not integrated - but for the axial reflections whose twin + // partner is the reflection itself or its Friedel mate, which carry this domain's intensity alone + // (IndexAndRefine) - scaled by this domain's share. `rotation` takes this domain's lattice onto the + // other's (same frame as `lattice`), where it is known. struct TwinOverlap { Coord t; int n = 0; + std::optional rotation; + float fraction = 0.0f; // the other domain's share of the crystal (the twin fraction) }; std::optional twin_overlap; }; diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 0297c3c2f..0aaca5272 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -376,12 +376,20 @@ struct TwinCompositeEvidence { std::array f1{}, f2{}; int both = 0, only1 = 0, only2 = 0, neither = 0; double p_both = 0.0; // the share of the nodes on both sub-lattices + // The other domain's share of the crystal (the twin fraction): of the summed intensity of the spots + // on one domain alone, the share on the other's - the two classes sample the same reflections of + // one crystal. Zero where no intensities were given. + double fraction = 0.0; double other_density = 0.0, neither_density = 0.0; bool composite = false; // two domains on their coincidence lattice CrystalLattice domain; // the stronger domain's cell, Niggli-reduced + // The proper rotation taking `domain` onto the other domain's cell, where their reduced bases + // agree to one (the twin law, up to the lattice's own rotations). + std::optional rotation; }; std::optional TwinCompositeAt(const CrystalLattice &prim, - const std::vector> &hkl, int n); + const std::vector> &hkl, int n, + const std::vector &intensity = {}); // Whether a sub-lattice of the standing cell - the same orientation, fewer nodes - explains the // pooled validation spots as well as the standing cell does. The two are compared on the same spots, so diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp index a836008a4..4c7179982 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -273,18 +273,24 @@ namespace { // // Following the reduction of an axis multiple to a centring condition, and its measurement as an // absence test, from this project's own lattice-arbiter analysis. + // `kept`, where given, receives the position in `cloud` of each index returned. std::vector> IndexedMillerIndices(const CrystalLattice &latt, - const std::vector &cloud, float tol) { + const std::vector &cloud, float tol, + std::vector *kept = nullptr) { const Coord a = latt.Vec0(), b = latt.Vec1(), c = latt.Vec2(); const float tol_sq = tol * tol; std::vector> hkl; hkl.reserve(cloud.size()); - for (const Coord &s : cloud) { + for (size_t i = 0; i < cloud.size(); i++) { + const Coord &s = cloud[i]; const float fh = a * s, fk = b * s, fl = c * s; // Coord operator* = dot = Miller index const float h = std::rint(fh), k = std::rint(fk), l = std::rint(fl); const float dh = fh - h, dk = fk - k, dl = fl - l; - if (dh * dh + dk * dk + dl * dl < tol_sq) + if (dh * dh + dk * dk + dl * dl < tol_sq) { hkl.push_back({static_cast(h), static_cast(k), static_cast(l)}); + if (kept) + kept->push_back(i); + } } return hkl; } @@ -369,7 +375,8 @@ namespace { } std::optional TwinCompositeAt(const CrystalLattice &prim, - const std::vector> &hkl, int n) { + const std::vector> &hkl, int n, + const std::vector &intensity) { const auto mod = [n](int x) { return (x % n + n) % n; }; // Each index-n sub-lattice once: the forms (u,v,w) mod n whose kernel has index n (no factor // shared by all of u, v, w and n), each the least of its multiples by the units mod n. @@ -435,6 +442,36 @@ std::optional TwinCompositeAt(const CrystalLattice &prim, std::swap(best->only1, best->only2); } best->domain = sub_of(best->f1); + if (intensity.size() == hkl.size()) { + double i1 = 0.0, i2 = 0.0; + for (size_t m = 0; m < hkl.size(); m++) { + const bool o1 = on(best->f1, hkl[m]), o2 = on(best->f2, hkl[m]); + if (o1 && !o2) + i1 += intensity[m]; + else if (o2 && !o1) + i2 += intensity[m]; + } + if (i1 + i2 > 0.0) + best->fraction = i2 / (i1 + i2); + } + // The rotation from the stronger domain's reduced basis to the other's, the other's with the pair of + // axes reversed (or none) that makes it a rotation - the one nearest to orthogonal. + { + const CrystalLattice other = sub_of(best->f2); + const RotMatrix from = RotMatrix(best->domain.Vec0(), best->domain.Vec1(), best->domain.Vec2()).invert(); + double best_err = 0.05; + for (const auto &[s0, s1, s2] : {std::array{1, 1, 1}, {1, -1, -1}, {-1, 1, -1}, {-1, -1, 1}}) { + const RotMatrix r = RotMatrix(other.Vec0() * s0, other.Vec1() * s1, other.Vec2() * s2) * from; + double err = 0.0; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + err = std::max(err, std::abs(r.Column(i) * r.Column(j) - (i == j ? 1.0 : 0.0))); + if (err < best_err) { + best_err = err; + best->rotation = r; + } + } + } // The share of the nodes on both sub-lattices, counted over one period of the forms: 1/n^2 // for a prime n, more where the two forms share a factor with n. int n_both = 0; @@ -1734,11 +1771,21 @@ auto Rugnux::FirstPassRun::TwinCompositeHypothesis(FirstPass best) -> FirstPass return best; const CrystalLattice prim = best.result->lattice .ToPrimitive(best.result->search_result.centering).NiggliReduce(); - const auto hkl = IndexedMillerIndices(prim, ValidationCloud(*best.result), index_tol); + std::vector kept; + const auto hkl = IndexedMillerIndices(prim, ValidationCloud(*best.result), index_tol, &kept); if (hkl.empty()) return best; + // The spots' intensities in the cloud's order (ValidationCloud), for the twin fraction. + std::vector cloud_intensity, intensity; + for (const int ordinal : validation) + if (const auto it = spot_cache.find(ordinal); it != spot_cache.end()) + for (const auto &sp : it->second) + cloud_intensity.push_back(sp.intensity); + for (const size_t i : kept) + if (i < cloud_intensity.size()) + intensity.push_back(cloud_intensity[i]); for (int n = 2; n <= TWIN_COMPOSITE_MAX_INDEX; n++) { - const auto tc = TwinCompositeAt(prim, hkl, n); + const auto tc = TwinCompositeAt(prim, hkl, n, intensity); if (!tc) continue; logger.Info("Twin-composite hypothesis, index {}: of {} validation spots on the standing cell " @@ -1760,7 +1807,9 @@ auto Rugnux::FirstPassRun::TwinCompositeHypothesis(FirstPass best) -> FirstPass // by n, t the weaker sub-lattice's form taken as a vector of the standing cell. cand.twin_overlap = RotationIndexerResult::TwinOverlap{ .t = tc->f2[0] * prim.Vec0() + tc->f2[1] * prim.Vec1() + tc->f2[2] * prim.Vec2(), - .n = n}; + .n = n, + .rotation = tc->rotation, + .fraction = static_cast(tc->fraction)}; const int score = CountIndexed(*indexer, cand); // The composite is the evidence that the domain's cell is this crystal's; what is left to ask // of the cell is that it is a lattice at all - that the validation spots it holds beat what a @@ -1772,12 +1821,18 @@ auto Rugnux::FirstPassRun::TwinCompositeHypothesis(FirstPass best) -> FirstPass logger.Info("Twin-composite hypothesis: the stronger domain's cell, {}-centred {} {:.2f} {:.2f} " "{:.2f} {:.1f} {:.1f} {:.1f}, indexes {}/{} validation frames - {} (the {:.0f}% of its " "reflections the other domain overlaps, h.({:.2f},{:.2f},{:.2f}) divisible by {}, is left out " - "of integration)", + "of integration{}; twin fraction {:.2f} from the intensity of the spots on one domain alone)", cand.search_result.centering, gemmi::crystal_system_str(cand.search_result.system), uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, score, static_cast(validation.size()), adopt ? "adopting it" : "keeping the standing cell", 100.0 * tc->p_both * n, cand.twin_overlap->t * cand.lattice.Astar(), - cand.twin_overlap->t * cand.lattice.Bstar(), cand.twin_overlap->t * cand.lattice.Cstar(), n); + cand.twin_overlap->t * cand.lattice.Bstar(), cand.twin_overlap->t * cand.lattice.Cstar(), n, + tc->rotation ? fmt::format(", but for the axial ones whose twin partner is the reflection itself " + "or its Friedel mate - the twin law a rotation of {:.1f} deg", + std::acos(std::clamp((tc->rotation->Column(0).x + tc->rotation->Column(1).y + + tc->rotation->Column(2).z - 1.0f) / 2.0f, -1.0f, 1.0f)) + * 180.0f / static_cast(M_PI)) + : std::string(), tc->fraction); if (adopt) { best = FirstPass{std::move(cand), score, fmt::format("{} (twin domain)", best.name), std::abs(tc->domain.CalcVolume())}; diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index e4d2024ac..5b842beba 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -820,6 +820,21 @@ TEST_CASE("TwinCompositeAt", "[process]") { CHECK(at2->neither == 0); CHECK(at2->p_both == Catch::Approx(0.25)); CHECK(same_cell(at2->domain, reduced)); + // The twin law, a twofold up to the lattice's own rotations (about c*, or about a, which the + // monoclinic twofold about b takes it to): either reverses b, so 0k0 is its own partner's Friedel + // mate and is kept. + REQUIRE(at2->rotation.has_value()); + const RotMatrix &r = *at2->rotation; + CHECK(r.Column(0).x + r.Column(1).y + r.Column(2).z == Catch::Approx(-1.0f).margin(1e-3)); + CHECK(((r * b) + b).Length() < 1e-3f * b.Length()); + // The twin fraction from the intensity of the spots on one domain alone: domain two at half the + // intensity of domain one is a third of the crystal (which domain comes first is the count's tie). + std::vector intensity; + for (size_t i = 0; i < cloud.size(); i++) + intensity.push_back(i % 2 ? 0.5f : 1.0f); + const auto weighted = TwinCompositeAt(doubled, indices(doubled, cloud), 2, intensity); + REQUIRE(weighted.has_value()); + CHECK(std::min(weighted->fraction, 1.0 - weighted->fraction) == Catch::Approx(1.0 / 3.0).margin(1e-3)); // A skewed index-5 supercell of it: the domains are index-10 sub-lattices, among them forms with no // coefficient prime to 10.