mx/Jungfraujoch
3
2 Fork 0
Code Issues Pull Requests Actions 3 Packages Projects Releases 88 Wiki Activity

RESULT-CHANGING: sparse-twin regression: a correct protein cell is no longer replaced by half of itself

Browse Source 
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 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37
This commit is contained in:
leonarski_fandClaude Opus 5.5 committed 2026-11-01 15:00:03 +01:00
1 parent 5c5ee93672
commit cf7ac18eb8
6 files changed
+118 -93
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

rugnux/Rugnux.cpp
+3 -7
View File
@@ -137,13 +137,9 @@ std::optional<size_t> 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<double>(only_standing - only_sub)
- static_cast<double>(std::max<int64_t>(0, standing.by_chance - sub.by_chance));
const double sigma = std::sqrt(static_cast<double>(std::max<int64_t>(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<double>(dropped) - chance <= SPOT_BUDGET_SIGNIFICANCE_Z * std::sqrt(std::max(1.0, chance));
}
bool ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e) {
rugnux/Rugnux.h
+12 -12
View File
@@ -377,8 +377,9 @@ std::optional<size_t> 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<int, 3> 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<RotMatrix> rotation;
double neither_expected = 0.0; // what chance puts on the class on neither
};
std::optional<TwinCompositeEvidence> TwinCompositeAt(const CrystalLattice &prim,
const std::vector<std::array<int, 3>> &hkl, int n,
const std::vector<float> &intensity = {});
double chance_share, const std::vector<float> &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
rugnux/RugnuxFirstPass.cpp
+48 -36
View File
@@ -376,7 +376,7 @@ namespace {
std::optional<TwinCompositeEvidence> TwinCompositeAt(const CrystalLattice &prim,
const std::vector<std::array<int, 3>> &hkl, int n,
const std::vector<float> &intensity) {
double chance_share, const std::vector<float> &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<TwinCompositeEvidence> 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<double>(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<int64_t, int64_t> Rugnux::FirstPassRun::PairedOnLattice(IndexAndRefine &idx, const RotationIndexerResult &a,
const RotationIndexerResult &b,
const std::vector<int> &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<int> &frames) {
PrefetchSpots(frames);
const bool index_ice_rings = rx.experiment_.GetIndexingSettings().GetIndexIceRings();
const auto on = [&](const RotationIndexerResult &r) {
idx.ForceRotationIndexerResult(r);
std::vector<std::vector<bool>> flags(frames.size());
ParallelFor(static_cast<int>(frames.size()),
std::min<size_t>(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<int64_t> dropped = 0;
ParallelFor(static_cast<int>(frames.size()),
std::min<size_t>(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<double>(at_standing.by_chance)
/ static_cast<double>(std::max<int64_t>(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<int>(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 "
rugnux/RugnuxFirstPass.h
+2 -2
View File
@@ -133,8 +133,8 @@ struct Rugnux::FirstPassRun {
const std::vector<int> &frames);
PooledEvidence MeasurePooled(IndexAndRefine &idx, const RotationIndexerResult &r);
PooledEvidence MeasurePooled(IndexAndRefine &idx, const RotationIndexerResult &r, const std::vector<int> &frames);
std::pair<int64_t, int64_t> PairedOnLattice(IndexAndRefine &idx, const RotationIndexerResult &a,
const RotationIndexerResult &b, const std::vector<int> &frames);
int64_t DroppedNodeSpots(IndexAndRefine &idx, const RotationIndexerResult &r, const Coord &t,
const std::vector<int> &frames);
int64_t MeasureBudget(IndexAndRefine &idx, const RotationIndexerResult &r);
std::unique_ptr<RotationIndexer> FeedScheme(IndexerThreadPool &pool, const std::vector<int> &ordinals,
const DiffractionExperiment *x = nullptr);
rugnux/RugnuxFirstPassRescues.cpp
+25 -18
View File
@@ -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<FirstPass> 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<double>(only_standing - only_sub),
std::sqrt(static_cast<double>(std::max<int64_t>(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<double>(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<double>(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<int> 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):{}{}",
tests/RugnuxTest.cpp
+28 -18
View File
@@ -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<float> 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<int> 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));
}