diff --git a/image_analysis/IndexAndRefine.cpp b/image_analysis/IndexAndRefine.cpp index eb04815f4..45f1632d8 100644 --- a/image_analysis/IndexAndRefine.cpp +++ b/image_analysis/IndexAndRefine.cpp @@ -167,6 +167,20 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetryRotat static_cast(std::lround(6.0f * (tw->t * l.Bstar()))), static_cast(std::lround(6.0f * (tw->t * l.Cstar())))}; outcome.twin_overlap_n = 6 * tw->n; + // The other domain's reflection at the node of h, read in this domain's basis: the index of + // R^-1 q, q the node, R the rotation from this domain to the other. + if (tw->rotation) { + const Coord a[3] = {l.Vec0(), l.Vec1(), l.Vec2()}, as[3] = {l.Astar(), l.Bstar(), l.Cstar()}; + std::array p{}; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + p[3 * i + j] = (*tw->rotation * a[i]) * as[j]; + outcome.twin_partner = p; + outcome.twin_domain_share = 1.0f - tw->fraction; + const auto right = [](const Coord &x, const Coord &y) { return std::abs(angle_deg(x, y) - 90.0f) < 0.5f; }; + for (int i = 0; i < 3; i++) + outcome.twin_axis_kept[i] = right(a[i], a[(i + 1) % 3]) && right(a[i], a[(i + 2) % 3]); + } } } return outcome; @@ -663,9 +677,34 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg, std::vector integrated = integrate(prediction.GetReflections(), nrefl, msg.number); if (const int n = outcome.twin_overlap_n; n > 0) { const auto &w = outcome.twin_overlap_w; + const auto &p = outcome.twin_partner; + // Kept only where the coincidence is exact: a twin by reticular PSEUDO-merohedry puts the + // other domain's node beside this one's rather than on it, by an offset that grows with the + // index, so a partner that rounds to the reflection itself proves nothing in general. Along an + // axis perpendicular to the other two it does - the twin takes that axis onto itself or its + // reverse exactly - and those are the rows a screw is read on. + const auto exact = [&](const Reflection &r) { + const int i_axis = r.h != 0 ? 0 : r.k != 0 ? 1 : 2; + const bool axial = (r.h != 0) + (r.k != 0) + (r.l != 0) == 1; + if (!p || !axial || !outcome.twin_axis_kept[i_axis]) + return false; + const auto row = [&](int i) { + return std::lround((*p)[3 * i] * r.h + (*p)[3 * i + 1] * r.k + (*p)[3 * i + 2] * r.l); + }; + const long h = row(0), k = row(1), l = row(2); + return (h == r.h && k == r.k && l == r.l) || (h == -r.h && k == -r.k && l == -r.l); + }; std::erase_if(integrated, [&](const Reflection &r) { - return (r.h * w[0] + r.k * w[1] + r.l * w[2]) % n == 0; + return (r.h * w[0] + r.k * w[1] + r.l * w[2]) % n == 0 && !exact(r); }); + // Such a reflection holds the whole crystal's intensity, every other one this domain's share of + // it: scaled by that share, it is on the others' scale. + for (auto &r : integrated) + if ((r.h * w[0] + r.k * w[1] + r.l * w[2]) % n == 0) { + r.I *= outcome.twin_domain_share; + r.sigma *= outcome.twin_domain_share; + r.var_bkg *= outcome.twin_domain_share * outcome.twin_domain_share; + } } msg.integrated_reflections = integrated.size(); auto integration_end_time = std::chrono::steady_clock::now(); diff --git a/image_analysis/IndexAndRefine.h b/image_analysis/IndexAndRefine.h index 5eaae7f43..75787533e 100644 --- a/image_analysis/IndexAndRefine.h +++ b/image_analysis/IndexAndRefine.h @@ -131,6 +131,13 @@ class IndexAndRefine { // overlap's index) are left out of integration. std::array twin_overlap_w{}; int twin_overlap_n = 0; + // Where the twin rotation is known: the twin partner of reflection h, in the candidate's basis, + // is the nearest integer to twin_partner * h (row-major 3x3); twin_axis_kept[i] where axis i is + // perpendicular to the other two, so that a reflection along it whose partner is itself or its + // Friedel mate coincides with that partner exactly. + std::optional> twin_partner; + std::array twin_axis_kept{}; + float twin_domain_share = 1.0f; // this domain's share of the crystal (1 - the twin fraction) explicit IndexingOutcome(const DiffractionExperiment& experiment_ref) : experiment(experiment_ref) {} diff --git a/image_analysis/rotation_indexer/RotationIndexer.h b/image_analysis/rotation_indexer/RotationIndexer.h index e17e22315..5825f5606 100644 --- a/image_analysis/rotation_indexer/RotationIndexer.h +++ b/image_analysis/rotation_indexer/RotationIndexer.h @@ -40,10 +40,15 @@ struct RotationIndexerResult { // Twinning by reticular merohedry (rugnux's twin-composite step): a second domain whose lattice // shares the reflections h of this one with h . t divisible by n, t a vector of this lattice (real // space, in the frame of `lattice`, so it holds through any change of basis). Those reflections - // carry both domains' intensity, and are not integrated. + // carry both domains' intensity, and are not integrated - but for the axial reflections whose twin + // partner is the reflection itself or its Friedel mate, which carry this domain's intensity alone + // (IndexAndRefine) - scaled by this domain's share. `rotation` takes this domain's lattice onto the + // other's (same frame as `lattice`), where it is known. struct TwinOverlap { Coord t; int n = 0; + std::optional rotation; + float fraction = 0.0f; // the other domain's share of the crystal (the twin fraction) }; std::optional twin_overlap; }; diff --git a/rugnux/Rugnux.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..0aaca5272 100644 --- a/rugnux/Rugnux.h +++ b/rugnux/Rugnux.h @@ -367,6 +367,40 @@ 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 + // The other domain's share of the crystal (the twin fraction): of the summed intensity of the spots + // on one domain alone, the share on the other's - the two classes sample the same reflections of + // one crystal. Zero where no intensities were given. + double fraction = 0.0; + double other_density = 0.0, neither_density = 0.0; + bool composite = false; // two domains on their coincidence lattice + CrystalLattice domain; // the stronger domain's cell, Niggli-reduced + // The proper rotation taking `domain` onto the other domain's cell, where their reduced bases + // agree to one (the twin law, up to the lattice's own rotations). + std::optional rotation; +}; +std::optional TwinCompositeAt(const CrystalLattice &prim, + const std::vector> &hkl, int n, + const std::vector &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); + // 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..4c7179982 100644 --- a/rugnux/RugnuxFirstPass.cpp +++ b/rugnux/RugnuxFirstPass.cpp @@ -273,18 +273,24 @@ namespace { // // Following the reduction of an axis multiple to a centring condition, and its measurement as an // absence test, from this project's own lattice-arbiter analysis. + // `kept`, where given, receives the position in `cloud` of each index returned. std::vector> IndexedMillerIndices(const CrystalLattice &latt, - const std::vector &cloud, float tol) { + const std::vector &cloud, float tol, + std::vector *kept = nullptr) { const Coord a = latt.Vec0(), b = latt.Vec1(), c = latt.Vec2(); const float tol_sq = tol * tol; std::vector> hkl; hkl.reserve(cloud.size()); - for (const Coord &s : cloud) { + for (size_t i = 0; i < cloud.size(); i++) { + const Coord &s = cloud[i]; const float fh = a * s, fk = b * s, fl = c * s; // Coord operator* = dot = Miller index const float h = std::rint(fh), k = std::rint(fk), l = std::rint(fl); const float dh = fh - h, dk = fk - k, dl = fl - l; - if (dh * dh + dk * dk + dl * dl < tol_sq) + if (dh * dh + dk * dk + dl * dl < tol_sq) { hkl.push_back({static_cast(h), static_cast(k), static_cast(l)}); + if (kept) + kept->push_back(i); + } } return hkl; } @@ -342,10 +348,147 @@ 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 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. + 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); + if (intensity.size() == hkl.size()) { + double i1 = 0.0, i2 = 0.0; + for (size_t m = 0; m < hkl.size(); m++) { + const bool o1 = on(best->f1, hkl[m]), o2 = on(best->f2, hkl[m]); + if (o1 && !o2) + i1 += intensity[m]; + else if (o2 && !o1) + i2 += intensity[m]; + } + if (i1 + i2 > 0.0) + best->fraction = i2 / (i1 + i2); + } + // The rotation from the stronger domain's reduced basis to the other's, the other's with the pair of + // axes reversed (or none) that makes it a rotation - the one nearest to orthogonal. + { + const CrystalLattice other = sub_of(best->f2); + const RotMatrix from = RotMatrix(best->domain.Vec0(), best->domain.Vec1(), best->domain.Vec2()).invert(); + double best_err = 0.05; + for (const auto &[s0, s1, s2] : {std::array{1, 1, 1}, {1, -1, -1}, {-1, 1, -1}, {-1, -1, 1}}) { + const RotMatrix r = RotMatrix(other.Vec0() * s0, other.Vec1() * s1, other.Vec2() * s2) * from; + double err = 0.0; + for (int i = 0; i < 3; i++) + for (int j = 0; j < 3; j++) + err = std::max(err, std::abs(r.Column(i) * r.Column(j) - (i == j ? 1.0 : 0.0))); + if (err < best_err) { + best_err = err; + best->rotation = r; + } + } + } + // The share of the nodes on both sub-lattices, counted over one period of the forms: 1/n^2 + // for a prime n, more where the two forms share a factor with n. + int n_both = 0; + 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 +540,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 +910,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,89 +1747,58 @@ 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; const CrystalLattice prim = best.result->lattice .ToPrimitive(best.result->search_result.centering).NiggliReduce(); - const auto hkl = IndexedMillerIndices(prim, ValidationCloud(*best.result), index_tol); + std::vector kept; + const auto hkl = IndexedMillerIndices(prim, ValidationCloud(*best.result), index_tol, &kept); if (hkl.empty()) return best; - 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; + // The spots' intensities in the cloud's order (ValidationCloud), for the twin fraction. + std::vector cloud_intensity, intensity; + for (const int ordinal : validation) + if (const auto it = spot_cache.find(ordinal); it != spot_cache.end()) + for (const auto &sp : it->second) + cloud_intensity.push_back(sp.intensity); + for (const size_t i : kept) + if (i < cloud_intensity.size()) + intensity.push_back(cloud_intensity[i]); + for (int n = 2; n <= TWIN_COMPOSITE_MAX_INDEX; n++) { + const auto tc = TwinCompositeAt(prim, hkl, n, intensity); + 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 +1806,36 @@ 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(), - .n = n}; + .t = tc->f2[0] * prim.Vec0() + tc->f2[1] * prim.Vec1() + tc->f2[2] * prim.Vec2(), + .n = n, + .rotation = tc->rotation, + .fraction = static_cast(tc->fraction)}; const int score = CountIndexed(*indexer, cand); + // The composite is the evidence that the domain's cell is this crystal's; what is left to ask + // of the cell is that it is a lattice at all - that the validation spots it holds beat what a + // 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)", + "of integration{}; twin fraction {:.2f} from the intensity of the spots on one domain alone)", cand.search_result.centering, gemmi::crystal_system_str(cand.search_result.system), uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, score, static_cast(validation.size()), adopt ? "adopting it" : "keeping the standing cell", - 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, + tc->rotation ? fmt::format(", but for the axial ones whose twin partner is the reflection itself " + "or its Friedel mate - the twin law a rotation of {:.1f} deg", + std::acos(std::clamp((tc->rotation->Column(0).x + tc->rotation->Column(1).y + + tc->rotation->Column(2).z - 1.0f) / 2.0f, -1.0f, 1.0f)) + * 180.0f / static_cast(M_PI)) + : std::string(), tc->fraction); if (adopt) { best = FirstPass{std::move(cand), score, fmt::format("{} (twin domain)", best.name), - std::abs(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..5b842beba 100644 --- a/tests/RugnuxTest.cpp +++ b/tests/RugnuxTest.cpp @@ -752,3 +752,100 @@ 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)); + // The twin law, a twofold up to the lattice's own rotations (about c*, or about a, which the + // monoclinic twofold about b takes it to): either reverses b, so 0k0 is its own partner's Friedel + // mate and is kept. + REQUIRE(at2->rotation.has_value()); + const RotMatrix &r = *at2->rotation; + CHECK(r.Column(0).x + r.Column(1).y + r.Column(2).z == Catch::Approx(-1.0f).margin(1e-3)); + CHECK(((r * b) + b).Length() < 1e-3f * b.Length()); + // The twin fraction from the intensity of the spots on one domain alone: domain two at half the + // intensity of domain one is a third of the crystal (which domain comes first is the count's tie). + std::vector intensity; + for (size_t i = 0; i < cloud.size(); i++) + intensity.push_back(i % 2 ? 0.5f : 1.0f); + const auto weighted = TwinCompositeAt(doubled, indices(doubled, cloud), 2, intensity); + REQUIRE(weighted.has_value()); + CHECK(std::min(weighted->fraction, 1.0 - weighted->fraction) == Catch::Approx(1.0 / 3.0).margin(1e-3)); + + // A skewed index-5 supercell of it: the domains are index-10 sub-lattices, among them forms with no + // coefficient prime to 10. + 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)); +}