Merge branch 'twin-overlap' into rc175: the twin overlap is left out of integration where the twin-composite step adopted a domain

Co-Authored-By: Claude Opus 5.5 (1M context) <noreply@anthropic.com>

image_analysis/IndexAndRefine.cpp

@@ -134,6 +134,15 @@ IndexAndRefine::IndexingOutcome IndexAndRefine::DetermineLatticeAndSymmetryRotat
         outcome.symmetry.centering = result->search_result.centering;
         outcome.symmetry.niggli_class = result->search_result.niggli_class;
         outcome.symmetry.crystal_system = result->search_result.system;
+        if (const auto &tw = result->twin_overlap) {
+            // t is a vector of the lattice, so in a centred cell's basis its coordinates are multiples
+            // of 1/2 or 1/3: in sixths they are integers, and the test is then modulo 6n.
+            const auto &l = result->lattice;
+            outcome.twin_overlap_w = {static_cast<int>(std::lround(6.0f * (tw->t * l.Astar()))),
+                                      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;
+        }
     }
     return outcome;
 }
@@ -587,7 +596,13 @@ void IndexAndRefine::QuickPredictAndIntegrate(DataMessage &msg,
     // caller's callback binds it to the right image (GPU-resident buffer, host buffer, or the assembled
     // FPGA image read straight on the CPU).
     auto integration_start_time = std::chrono::steady_clock::now();
-    const std::vector<Reflection> integrated = integrate(prediction.GetReflections(), nrefl, msg.number);
+    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;
+        std::erase_if(integrated, [&](const Reflection &r) {
+            return (r.h * w[0] + r.k * w[1] + r.l * w[2]) % n == 0;
+        });
+    }
     msg.integrated_reflections = integrated.size();
     auto integration_end_time = std::chrono::steady_clock::now();
     msg.integration_time_s = std::chrono::duration<float>(integration_end_time - integration_start_time).count();

image_analysis/IndexAndRefine.h

@@ -85,6 +85,11 @@ class IndexAndRefine {
             .crystal_system = gemmi::CrystalSystem::Triclinic
         };
         bool beam_center_updated = false;
+        // The rotation result's twin overlap (RotationIndexerResult::TwinOverlap) in the candidate's
+        // basis, in sixths: reflections with h*w[0] + k*w[1] + l*w[2] divisible by n (6 times the
+        // overlap's index) are left out of integration.
+        std::array<int, 3> twin_overlap_w{};
+        int twin_overlap_n = 0;
 
         explicit IndexingOutcome(const DiffractionExperiment& experiment_ref)
             : experiment(experiment_ref) {}

image_analysis/rotation_indexer/RotationIndexer.cpp

@@ -244,6 +244,7 @@ void RotationIndexer::RunIndexing() {
     std::unique_lock ul(m);
     if (!axis_)
         return;
+    twin_overlap_.reset();
 
     const KeyDigest key = Digest(InputKey());
     std::optional<IndexingOutcome> kept;
@@ -886,6 +887,7 @@ std::optional<RotationIndexerResult> RotationIndexer::GetLattice() const {
         .axis = axis_,
         .unconstrained = unconstrained_,
         .tilt_walk = tilt_walk_,
+        .twin_overlap = twin_overlap_,
     };
 }
 
@@ -903,6 +905,7 @@ void RotationIndexer::ForceResult(const RotationIndexerResult &result) {
     axis_ = result.axis;
     unconstrained_ = result.unconstrained;
     tilt_walk_ = result.tilt_walk;
+    twin_overlap_ = result.twin_overlap;
 }
 
 bool RotationIndexer::AccumulationFull() const {
@@ -912,6 +915,7 @@ bool RotationIndexer::AccumulationFull() const {
 
 void RotationIndexer::ForceLattice(const CrystalLattice &lattice) {
     indexed_lattice = lattice;
+    twin_overlap_.reset();
     const gemmi::SpaceGroup &sg = experiment.GetSpaceGroupOrP1();
     search_result_ = LatticeSearchResult{
         .niggli_class = 0, // Since Niggli class was not searched for, we don't know which one

image_analysis/rotation_indexer/RotationIndexer.h

@@ -37,6 +37,15 @@ struct RotationIndexerResult {
         bool refused = false;
     };
     std::optional<TiltWalk> tilt_walk;
+    // 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.
+    struct TwinOverlap {
+        Coord t;
+        int n = 0;
+    };
+    std::optional<TwinOverlap> twin_overlap;
 };
 
 // What a memo verification throws where a reused result differs from its recomputation - its own
@@ -92,6 +101,7 @@ private:
     std::vector<CrystalLattice> extra_lattices_;
     std::shared_ptr<RotationIndexerResult> unconstrained_;
     std::optional<RotationIndexerResult::TiltWalk> tilt_walk_;
+    std::optional<RotationIndexerResult::TwinOverlap> twin_overlap_;   // only ever forced, see ForceResult
 
     IndexerThreadPool &indexer_;
 

rugnux/Rugnux.h

@@ -404,6 +404,9 @@ struct ProcessResult {
     // number a few dozen against the merge's tens of thousands, and a bulk statistic therefore cannot
     // see them at all - which is why the two-pass quality guard counts them separately.
     int64_t search_merge_axial_reflections = 0;
+    // The index n of a twin overlap the pass left out of integration (RotationIndexerResult::
+    // TwinOverlap), 0 where there was none: its merges then hold (n - 1)/n of the lattice's nodes.
+    int twin_overlap_index = 0;
     // How much signal that merge holds: reflections at I/sigma >= 2, and how many reflections it holds
     // at all. This is what the two-pass quality guard compares the passes on - see RefinedPassIsWorse
     // for why it is these and not CC1/2.

rugnux/RugnuxFirstPass.cpp

@@ -1553,14 +1553,24 @@ auto Rugnux::FirstPassRun::TwinCompositeHypothesis(FirstPass best) -> FirstPass
         cand.lattice = cand.search_result.conventional;
         cand.unconstrained.reset();
         cand.extra_lattices.clear();
+        // The overlap, the reflections on both sub-lattices, holds both domains' intensity: the
+        // 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};
         const int score = CountIndexed(*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 - {}",
+                    "{:.2f} {:.1f} {:.1f} {:.1f}, indexes {}/{} validation frames - {} (the 1/{} of its "
+                    "reflections the other domain overlaps, h.({:.2f},{:.2f},{:.2f}) divisible by {}, is left out "
+                    "of integration)",
                     cand.search_result.centering, gemmi::crystal_system_str(cand.search_result.system),
                     uc.a, uc.b, uc.c, uc.alpha, uc.beta, uc.gamma, score,
-                    static_cast<int>(validation.size()), adopt ? "adopting it" : "keeping the standing cell");
+                    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);
         if (adopt) {
             best = FirstPass{std::move(cand), score, fmt::format("{} (twin domain)", best.name),
                              std::abs(sub.CalcVolume())};

rugnux/RugnuxImagePass.cpp

@@ -557,6 +557,7 @@ void Rugnux::ImagePass(PipelineLocals &p) {
             };
             result.rotation_lattice_found = true;
             result.rotation_lattice_type = end_msg.rotation_lattice_type;
+            result.twin_overlap_index = rot->twin_overlap ? rot->twin_overlap->n : 0;
             // What the equality this class asserts is worth on a REFINED cell. The indexer refines the
             // same candidate a second time with nothing held (its pseudo-symmetry guard) and hands that
             // cell over whenever the constrained one won; it comes back in the primitive reduced

rugnux/RugnuxPasses.cpp

@@ -122,9 +122,10 @@ namespace {
     // of reflections at a given Miller index the same count on both sides. A pass that re-indexes at
     // the refined geometry can settle on a centred setting of the same lattice the other pass held
     // primitive; the cells are then different bases and only the counts taken over all reflections
-    // survive the change.
+    // survive the change. Nor is it the same count where one pass left a twin overlap out.
     bool SameSetting(const ProcessResult &pass1, const ProcessResult &pass2) {
         if (pass1.consensus_centering != pass2.consensus_centering
+            || pass1.twin_overlap_index != pass2.twin_overlap_index
             || !pass1.consensus_cell || !pass2.consensus_cell)
             return false;
         const UnitCell &c1 = *pass1.consensus_cell, &c2 = *pass2.consensus_cell;
@@ -193,8 +194,14 @@ std::string RefinedPassIsWorse(const ProcessResult &pass1, const ProcessResult &
     // pass on the crystal's own lattice and the run went back to the supercell, while that merge
     // held more reflections at <= -2 than at >= 2. On a real superstructure the extra classes hold
     // signal, and it is counted.
-    const auto signal = [](const ProcessResult &r) {
-        return static_cast<double>(r.search_merge_strong_reflections - r.search_merge_negative_reflections);
+    // A pass that left a twin overlap of index n out of integration merged (n - 1)/n of the
+    // lattice's nodes, and its counts are taken up to the whole lattice.
+    const auto whole_lattice = [](const ProcessResult &r) {
+        return r.twin_overlap_index > 0 ? r.twin_overlap_index / (r.twin_overlap_index - 1.0) : 1.0;
+    };
+    const auto signal = [&](const ProcessResult &r) {
+        return static_cast<double>(r.search_merge_strong_reflections - r.search_merge_negative_reflections)
+               * whole_lattice(r);
     };
     // ... and the loss has to be a loss rather than a dilution. The wider disk pulls weaker
     // reflections into the merge and the error model is then refitted over them, so a pass can hold
@@ -208,8 +215,8 @@ std::string RefinedPassIsWorse(const ProcessResult &pass1, const ProcessResult &
     // the smaller lattice is scaled by the index n, which is an integer - the same lattice, whose
     // volume two refinements put a percent or two apart, is index 1 and compares as it always did.
     const auto coverage_does_not_grow = [&] {
-        double n1 = static_cast<double>(pass1.search_merge_reflections);
-        double n2 = static_cast<double>(pass2.search_merge_reflections);
+        double n1 = static_cast<double>(pass1.search_merge_reflections) * whole_lattice(pass1);
+        double n2 = static_cast<double>(pass2.search_merge_reflections) * whole_lattice(pass2);
         const double v1 = PrimitiveCellVolume(pass1), v2 = PrimitiveCellVolume(pass2);
         if (v1 > 0.0 && v2 > 0.0) {
             if (v1 >= v2)

tests/RugnuxTest.cpp

@@ -395,6 +395,15 @@ TEST_CASE("RefinedPassIsWorse", "[process]") {
     with_rows.search_merge_axial_reflections = 16;
     lost_rows.search_merge_axial_reflections = 8;
     CHECK(!RefinedPassIsWorse(with_rows, lost_rows).empty());
+
+    // A refined pass that left a twin overlap of index 3 out of integration merged two thirds of the
+    // nodes: two thirds of the signal is no loss, and its axial rows are not the header pass's rows.
+    ProcessResult twin = pass(cell, 'P', 34000, 700, 132000);
+    twin.twin_overlap_index = 3;
+    twin.search_merge_axial_reflections = 10;
+    CHECK(RefinedPassIsWorse(with_rows, twin).empty());
+    twin.search_merge_strong_reflections = 28000;
+    CHECK(!RefinedPassIsWorse(with_rows, twin).empty());
 }
 
 // The beam-centre arbitration between a first pass at the file's centre and one at the measured