From cf7ac18eb8bed0a897d08d2a97027ad24156bb33 Mon Sep 17 00:00:00 2001 From: Filip Leonarski Date: Thu, 29 Oct 2026 06:37:29 +0100 Subject: [PATCH] RESULT-CHANGING: sparse-twin regression: a correct protein cell is no longer replaced by half of itself The full battery of rc181 39d908427 lost six protein sets to sparse-twin (08ae0fa05): 7bgt, 9pbb, 6i3j, 8xtf, 4p8r and insu_I_x06da_low_isa took a cell of half the primitive volume. Two causes, two fixes. The sub-lattice (halving) test compared the standing cell and its sub-lattice each with its own per-frame orientation polish, so they were not nested: a sub-lattice could settle on another orientation and hold spots the standing cell did not (6i3j: 4550 spots on the sub-lattice alone), and a correct cell lost to its half. It now counts on the standing cell's own fit the spots whose index falls on the nodes the sub-lattice drops, against what a wrong spindle angle puts on those nodes (half the standing cell's chance count), Poisson at SPOT_BUDGET_SIGNIFICANCE_Z (SubLatticeSuffices(dropped, chance, sub), DroppedNodeSpots). A correct cell's dropped nodes hold half its spots (thousands against tens); a supercell's hold what chance puts there (L-aspartic acid: 3-9 against 2.5-7). The twin-composite test read per-node densities that hold only for a prime index: where the two forms share a factor with n, a single crystal fills the four classes in proportions that are not their shares of the nodes, and an index-6 pair holding 1629/1572/1567/1540 spots read as two domains. The class on neither must now hold no more than chance puts there - the standing cell's chance share of its spots on the class's share of the nodes - besides the densities. And the domain must index a majority of the frames or as many as the standing cell (7bgt adopted a domain indexing 6/60 against 60/60; the beats-chance bar alone let it through). The genuine index-3 [111] composite of a twinned insulin set is still adopted, as at 297ef9dc7. Tests: SubLatticeSuffices rewritten for the nested count; TwinCompositeAt gains a case where the class on neither holds strays beyond chance. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37 --- rugnux/Rugnux.cpp | 10 ++-- rugnux/Rugnux.h | 24 ++++----- rugnux/RugnuxFirstPass.cpp | 84 ++++++++++++++++++------------- rugnux/RugnuxFirstPass.h | 4 +- rugnux/RugnuxFirstPassRescues.cpp | 43 +++++++++------- tests/RugnuxTest.cpp | 46 ++++++++++------- 6 files changed, 118 insertions(+), 93 deletions(-) diff --git a/rugnux/Rugnux.cpp b/rugnux/Rugnux.cpp index 3e91360df..ba7cb9add 100644 --- a/rugnux/Rugnux.cpp +++ b/rugnux/Rugnux.cpp @@ -137,13 +137,9 @@ 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 SubLatticeSuffices(int64_t dropped, double chance, const ValidationSpotEvidence &sub) { + return ValidationEvidenceBeatsChance(sub) + && static_cast(dropped) - chance <= SPOT_BUDGET_SIGNIFICANCE_Z * std::sqrt(std::max(1.0, chance)); } bool ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e) { diff --git a/rugnux/Rugnux.h b/rugnux/Rugnux.h index 2fdd64b38..a68da9d2d 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -377,8 +377,9 @@ std::optional LongAxisChoice(const ValidationSpotEvidence &standing, // 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. +// first, with the spots on each class and the per-node densities the verdict compares. `chance_share` +// is the share of the standing cell's spots a wrong spindle angle puts on it. 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; @@ -393,20 +394,19 @@ struct TwinCompositeEvidence { // 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; + double neither_expected = 0.0; // what chance puts on the class on neither }; std::optional TwinCompositeAt(const CrystalLattice &prim, const std::vector> &hkl, int n, - const std::vector &intensity = {}); + double chance_share, 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 -// 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); +// Whether a sub-lattice of the standing cell - the same orientation, half its nodes - explains the +// pooled validation spots as well as the standing cell does: it beats chance, and of the spots on the +// standing cell (as its own fit indexes them) those on the nodes the sub-lattice drops, `dropped`, are +// not more than what a wrong spindle angle puts on those nodes, `chance`, by SPOT_BUDGET_SIGNIFICANCE_Z +// Poisson standard errors. Where they are not, the dropped nodes hold nothing of their own, and the +// smaller cell is the hypothesis the spots support. +bool SubLatticeSuffices(int64_t dropped, double chance, 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 diff --git a/rugnux/RugnuxFirstPass.cpp b/rugnux/RugnuxFirstPass.cpp index 4c7179982..73016d10b 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -376,7 +376,7 @@ namespace { std::optional TwinCompositeAt(const CrystalLattice &prim, const std::vector> &hkl, int n, - const std::vector &intensity) { + double chance_share, 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. @@ -485,7 +485,16 @@ std::optional TwinCompositeAt(const CrystalLattice &prim, // 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; + // Two domains leave the class on neither with what chance puts there and nothing more: the spots a + // wrong spindle angle puts on the standing cell, as a share of those it indexes (`chance_share`), on + // the class's share of the nodes. The per-node densities alone do not decide it: for an n with a + // factor the two forms share, a single crystal whose lattice is a sub-lattice of the standing cell + // fills the four classes in proportions that are not their shares of the nodes (measured: an + // index-6 pair holding 1629/1572/1567/1540 spots read as a composite). + best->neither_expected = chance_share * static_cast(hkl.size()) * (1.0 - 2.0 / n + best->p_both); + const bool neither_by_chance = best->neither - best->neither_expected + <= SPOT_BUDGET_SIGNIFICANCE_Z * std::sqrt(std::max(1.0, best->neither_expected)); + best->composite = best->only2 >= 20 && best->other_density > 4.0 * best->neither_density && neither_by_chance; return best; } @@ -910,37 +919,33 @@ 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) { +// Of the frames' spots on the standing lattice `r`, as its own per-frame fit indexes them, how many sit +// on the nodes a sub-lattice drops: those whose reflection q has q . t half-integral, t the sub-lattice's +// half-translation (a vector of the same frame as r's lattice). Counted on the standing cell's fit +// alone, so the two lattices are compared on one indexing of the spots - a sub-lattice given its own +// per-frame fit is free to settle on another orientation and hold spots the standing cell does not. +int64_t Rugnux::FirstPassRun::DroppedNodeSpots(IndexAndRefine &idx, const RotationIndexerResult &r, + const Coord &t, 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}; + const float c[3] = {r.lattice.Astar() * t, r.lattice.Bstar() * t, r.lattice.Cstar() * t}; + idx.ForceRotationIndexerResult(r); + std::atomic dropped = 0; + 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) + if ((index_ice_rings || !sp.ice_ring) && (sp.h != 0 || sp.k != 0 || sp.l != 0)) { + const float x = c[0] * sp.h + c[1] * sp.k + c[2] * sp.l; + if (std::abs(x - std::round(x)) > 0.25f) + ++dropped; + } + }); + return dropped.load(); } // How deep into an image's intensity-ordered spot list this lattice is still being seen - the @@ -1784,16 +1789,19 @@ auto Rugnux::FirstPassRun::TwinCompositeHypothesis(FirstPass best) -> FirstPass for (const size_t i : kept) if (i < cloud_intensity.size()) intensity.push_back(cloud_intensity[i]); + const PooledEvidence at_standing = MeasurePooled(*indexer, *best.result); + const double chance_share = static_cast(at_standing.by_chance) + / static_cast(std::max(1, at_standing.on_lattice)); for (int n = 2; n <= TWIN_COMPOSITE_MAX_INDEX; n++) { - const auto tc = TwinCompositeAt(prim, hkl, n, intensity); + const auto tc = TwinCompositeAt(prim, hkl, n, chance_share, intensity); if (!tc) continue; logger.Info("Twin-composite hypothesis, index {}: of {} validation spots on the standing cell " "({:.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, + "together and {} on neither ({:.1f} by chance) - per node {:.1f} on the emptiest of the other " + "classes against {:.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->both, tc->neither, tc->neither_expected, tc->other_density, tc->neither_density, tc->composite ? "two domains on their coincidence lattice" : "one lattice"); if (!tc->composite) continue; @@ -1816,7 +1824,11 @@ auto Rugnux::FirstPassRun::TwinCompositeHypothesis(FirstPass best) -> FirstPass // 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)); + // And the domain must index a majority of the frames or as many as the standing cell, the bar the + // sub-lattice hypothesis asks. + const int majority = static_cast(validation.size()) / 2; + const bool adopt = (score > majority || score >= best.score) + && ValidationEvidenceBeatsChance(MeasurePooled(*indexer, cand)); const auto uc = cand.lattice.GetUnitCell(); logger.Info("Twin-composite hypothesis: the stronger domain's cell, {}-centred {} {:.2f} {:.2f} " "{:.2f} {:.1f} {:.1f} {:.1f}, indexes {}/{} validation frames - {} (the {:.0f}% of its " diff --git a/rugnux/RugnuxFirstPass.h b/rugnux/RugnuxFirstPass.h index 57ddf3be0..e0d81849d 100644 --- a/rugnux/RugnuxFirstPass.h +++ b/rugnux/RugnuxFirstPass.h @@ -133,8 +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 DroppedNodeSpots(IndexAndRefine &idx, const RotationIndexerResult &r, const Coord &t, + 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); diff --git a/rugnux/RugnuxFirstPassRescues.cpp b/rugnux/RugnuxFirstPassRescues.cpp index a9e7fbd43..99fb8793a 100644 --- a/rugnux/RugnuxFirstPassRescues.cpp +++ b/rugnux/RugnuxFirstPassRescues.cpp @@ -970,9 +970,11 @@ auto Rugnux::FirstPassRun::ShortAxisHypothesis(FirstPass best) -> FirstPass { // 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 +// of their own: of the pooled validation spots as the standing cell's own fit indexes them, those +// on the dropped nodes, against what a wrong spindle angle puts there (SubLatticeSuffices). One +// fit, so the comparison is nested: given a fit of its own per frame, a sub-lattice can settle on +// another orientation and hold spots the standing cell does not, which made a correct protein +// cell lose to its own half (measured: 4550 spots on a sub-lattice alone). 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 @@ -996,6 +998,7 @@ auto Rugnux::FirstPassRun::SubLatticeHypothesis(FirstPass best) -> FirstPass { .ToPrimitive(best.result->search_result.centering); const Coord A = prim.Vec0(), B = prim.Vec1(), C = prim.Vec2(); std::optional sub_best; + Coord sub_t; // the best sub-lattice's half-translation std::string counts; for (int u = 0; u < 2; u++) for (int v = 0; v < 2; v++) @@ -1016,10 +1019,12 @@ auto Rugnux::FirstPassRun::SubLatticeHypothesis(FirstPass best) -> FirstPass { cand.extra_lattices.clear(); const int score = CountIndexed(*indexer, cand); counts += fmt::format(" ({},{},{}) {}", u, v, w, score); - if (!sub_best || score > sub_best->score) + if (!sub_best || score > sub_best->score) { sub_best = FirstPass{std::move(cand), score, fmt::format("{} (sub-lattice {}{}{})", best.name, u, v, w), std::abs(sub.CalcVolume())}; + sub_t = t; + } } bool adopt = false; std::string verdict = " - none indexes a majority of the frames or as many as the standing cell"; @@ -1040,26 +1045,28 @@ auto Rugnux::FirstPassRun::SubLatticeHypothesis(FirstPass best) -> FirstPass { 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") + int64_t dropped = DroppedNodeSpots(*indexer, standing, sub_t, validation); + // The nodes the sub-lattice drops are half of the standing cell's, so they take half of what + // a wrong spindle angle puts on it. + const auto chance = [](const PooledEvidence &e) { return static_cast(e.by_chance) / 2.0; }; + if (NearTie("sub-lattice hypothesis", fmt::format("spots on the dropped nodes {}", dropped), + fmt::format("chance {:.1f}", chance(at_standing)), static_cast(dropped) - chance(at_standing), + std::sqrt(std::max(1.0, chance(at_standing))), "Poisson on pooled validation spots") && DecideLate()) { - const bool early = SubLatticeSuffices(only_standing, only_sub, at_standing, at_sub); + const bool early = SubLatticeSuffices(dropped, chance(at_standing), at_sub); const std::vector all = AllValidationFrames(); at_standing = MeasurePooled(*indexer, standing, all); at_sub = MeasurePooled(*indexer, sub, all); - std::tie(only_standing, only_sub) = PairedOnLattice(*indexer, standing, sub, all); + dropped = DroppedNodeSpots(*indexer, standing, sub_t, all); LogDecidedLate("sub-lattice hypothesis", "the pooled spots counted again on all the frames", early, - SubLatticeSuffices(only_standing, only_sub, at_standing, at_sub)); + SubLatticeSuffices(dropped, chance(at_standing), at_sub)); } - 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"); + adopt = SubLatticeSuffices(dropped, chance(at_standing), at_sub); + verdict = fmt::format(" - the best puts {}/{} pooled spots on itself ({} by chance); of the standing " + "cell's {}/{} ({}), {} sit on the nodes it drops against {:.1f} by chance, {}", + at_sub.on_lattice, at_sub.spots, at_sub.by_chance, at_standing.on_lattice, + at_standing.spots, at_standing.by_chance, dropped, chance(at_standing), + 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):{}{}", diff --git a/tests/RugnuxTest.cpp b/tests/RugnuxTest.cpp index 18ca97d90..66362e363 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -816,27 +816,23 @@ TEST_CASE("LongAxisChoice_PooledSpotsNotFrames", "[LongAxis]") { 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. +// The sub-lattice decision: the standing cell's spots on the nodes the sub-lattice drops, against what +// a wrong spindle angle puts on them. 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)); + // The dropped nodes hold what chance puts there, or a few more within its noise: halve. + CHECK(SubLatticeSuffices(3, 2.5, sub)); + CHECK(SubLatticeSuffices(7, 2.5, sub)); // 4.5 over, 3.29 sqrt(2.5) = 5.2 + // A class of its own (a superstructure, a pseudo-translation, or half of a correct protein cell's + // reflections): keep the standing cell. + CHECK_FALSE(SubLatticeSuffices(12, 2.5, sub)); + CHECK_FALSE(SubLatticeSuffices(4000, 70.0, 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)); + CHECK_FALSE(SubLatticeSuffices(0, 2.5, chance)); } // Two domains of a monoclinic cell twinned by a twofold about c*. With c cos(beta) = -a/4 the twin @@ -877,7 +873,7 @@ TEST_CASE("TwinCompositeAt", "[process]") { // The coincidence lattice, index 2. const CrystalLattice doubled = CrystalLattice(a, b, c * 2.0f).NiggliReduce(); - const auto at2 = TwinCompositeAt(doubled, indices(doubled, cloud), 2); + const auto at2 = TwinCompositeAt(doubled, indices(doubled, cloud), 2, 0.01); REQUIRE(at2.has_value()); CHECK(at2->composite); CHECK(at2->neither == 0); @@ -895,20 +891,34 @@ TEST_CASE("TwinCompositeAt", "[process]") { 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); + const auto weighted = TwinCompositeAt(doubled, indices(doubled, cloud), 2, 0.01, 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. 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); + const auto at10 = TwinCompositeAt(skewed, indices(skewed, cloud), 10, 0.01); REQUIRE(at10.has_value()); CHECK(at10->composite); CHECK(at10->neither == 0); CHECK(same_cell(at10->domain, reduced)); + // The same two domains with the class on neither holding far more than the chance share says it + // can (a share of 1% of the spots by chance): not two domains. + auto crowded = indices(doubled, cloud); + const size_t n_domains = crowded.size(); + // Spots on nodes of every class, as strays fall: a quarter of the nodes are on neither. + std::mt19937 rng(7); + std::uniform_int_distribution index(-6, 6); + for (int i = 0; i < 400; i++) + crowded.push_back({index(rng), index(rng), index(rng)}); + const auto with_neither = TwinCompositeAt(doubled, crowded, 2, 0.01); + REQUIRE(with_neither.has_value()); + CHECK(with_neither->neither_expected == Catch::Approx(0.01 * (n_domains + 400) * 0.25)); + CHECK_FALSE(with_neither->composite); + // One domain alone on the doubled cell: one lattice, not a composite. - const auto alone = TwinCompositeAt(doubled, indices(doubled, single), 2); + const auto alone = TwinCompositeAt(doubled, indices(doubled, single), 2, 0.01); CHECK(!(alone.has_value() && alone->composite)); }