Merge sparse-twin into rc181: a first-pass supercell on sparse or twinned data comes down to the true cell (paired sub-lattice test, twin composites to index 10), the 2_1 screw kept from exactly coinciding twin overlaps

Co-Authored-By: Claude Opus 5.5 <noreply@anthropic.com>
Claude-Session: https://claude.ai/code/session_01HVeDQqN9qRD57vnShzTC37

image_analysis/IndexAndRefine.cpp

@@ -167,6 +167,20 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetryRotat
                                       static_cast<int>(std::lround(6.0f * (tw->t * l.Bstar()))),
                                       static_cast<int>(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<float, 9> 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<Reflection> 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();

image_analysis/IndexAndRefine.h

@@ -131,6 +131,13 @@ class IndexAndRefine {
         // overlap's index) are left out of integration.
         std::array<int, 3> 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<std::array<float, 9>> twin_partner;
+        std::array<bool, 3> 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) {}

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<RotMatrix> rotation;
+        float fraction = 0.0f;   // the other domain's share of the crystal (the twin fraction)
     };
     std::optional<TwinOverlap> twin_overlap;
 };

rugnux/Rugnux.cpp

@@ -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

@@ -367,6 +367,40 @@ 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
+    // 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<RotMatrix> rotation;
+};
+std::optional<TwinCompositeEvidence> TwinCompositeAt(const CrystalLattice &prim,
+                                                     const std::vector<std::array<int, 3>> &hkl, int n,
+                                                     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);
+
 // 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

@@ -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<std::array<int, 3>> IndexedMillerIndices(const CrystalLattice &latt,
-                                                        const std::vector<Coord> &cloud, float tol) {
+                                                        const std::vector<Coord> &cloud, float tol,
+                                                        std::vector<size_t> *kept = nullptr) {
         const Coord a = latt.Vec0(), b = latt.Vec1(), c = latt.Vec2();
         const float tol_sq = tol * tol;
         std::vector<std::array<int, 3>> 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<int>(h), static_cast<int>(k), static_cast<int>(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<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 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.
+    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);
+    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<float, 3>{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<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 +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<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,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<size_t> 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<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;
+    // The spots' intensities in the cloud's order (ValidationCloud), for the twin fraction.
+    std::vector<float> 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<float>(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<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)",
+                    "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<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,
+                    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<float>(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;
         }
     }

rugnux/RugnuxFirstPass.h

@@ -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

@@ -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

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