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)); +}