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

RESULT-CHANGING: rugnux first pass: a supercell of a sparse or twinned pattern comes down to the true cell

Browse Source 
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 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37
(cherry picked from commit f8544981ab647cb7c350f9d73c81d499e3d330b0)
This commit is contained in:
leonarski_f committed 2026-10-29 00:43:30 +01:00
1 parent 297ef9dc78
commit 08ae0fa05c
6 files changed
+357 -120
Download Patch File Download Diff File Expand all files Collapse all files

No files matched your search

rugnux/Rugnux.cpp
+9
View File
@@ -137,6 +137,15 @@ 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 ValidationEvidenceBeatsChance(const ValidationSpotEvidence &e) {
return ValidationEvidenceBeatsChance(e, std::erfc(SPOT_BUDGET_SIGNIFICANCE_Z / std::sqrt(2.0)));
}
rugnux/Rugnux.h
+26
View File
@@ -367,6 +367,32 @@ struct LongAxisCandidate {
std::optional<size_t> LongAxisChoice(const ValidationSpotEvidence &standing,
const std::vector<LongAxisCandidate> &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<int, 3> 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<TwinCompositeEvidence> TwinCompositeAt(const CrystalLattice &prim,
const std::vector<std::array<int, 3>> &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
rugnux/RugnuxFirstPass.cpp
+181 -79
View File
@@ -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<float>(m[i][0]) + e[1] * static_cast<float>(m[i][1])
+ e[2] * static_cast<float>(m[i][2])) / static_cast<float>(n);
return {v[0], v[1], v[2]};
}
}
std::optional<TwinCompositeEvidence> TwinCompositeAt(const CrystalLattice &prim,
const std::vector<std::array<int, 3>> &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<std::array<int, 3>> 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<int, 3> 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<int, 3>{mod(k * u), mod(k * v), mod(k * w)};
if (least)
forms.push_back(f);
}
const auto on = [&](const std::array<int, 3> &f, const std::array<int, 3> &m) {
return mod(f[0] * m[0] + f[1] * m[1] + f[2] * m[2]) == 0;
};
const auto sub_of = [&](const std::array<int, 3> &f) {
AddedClassEvidence ev;
ev.u = f[0]; ev.v = f[1]; ev.w = f[2];
return SubLattice(prim, ev, n).NiggliReduce();
};
std::vector<UnitCell> 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<bool, 3>{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<TwinCompositeEvidence> 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<double>(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<int64_t, int64_t> Rugnux::FirstPassRun::PairedOnLattice(IndexAndRefine &idx, const RotationIndexerResult &a,
const RotationIndexerResult &b,
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};
}
// 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<std::array<int, 3>> 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<int, 3> &f, const std::array<int, 3> &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<double>(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<int>(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<int>(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;
}
}
rugnux/RugnuxFirstPass.h
+4
View File
@@ -133,6 +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 MeasureBudget(IndexAndRefine &idx, const RotationIndexerResult &r);
std::unique_ptr<RotationIndexer> FeedScheme(IndexerThreadPool &pool, const std::vector<int> &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;
rugnux/RugnuxFirstPassRescues.cpp
+55 -41
View File
@@ -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<int>(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<double>(only_standing - only_sub),
std::sqrt(static_cast<double>(std::max<int64_t>(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<int> 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 {
tests/RugnuxTest.cpp
+82
View File
@@ -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<float>(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<Coord> 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<Coord> &qs) {
std::vector<std::array<int, 3>> hkl;
for (const auto &q : qs)
hkl.push_back({static_cast<int>(std::lround(s.Vec0() * q)), static_cast<int>(std::lround(s.Vec1() * q)),
static_cast<int>(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));
}